Carbon migration and metagenomic characteristics during anaerobic digestion of rice straw

Background Considerable interest has been expressed in the development of anaerobic digestion (AD) of straw to solve the environmental problems caused by the dumping and burning of straw and to generate clean energy. However, the poor biodegradability of straw and the low efficiency of energy generation achieved during its AD are problematic. Studying the parameter changes involved in the process of AD is helpful for clarifying its micro-mechanisms and providing a theoretical basis for improving its efficiency. Currently, most research into process parameters has focused on gas production, methane content, pH, and volatile fatty acid (VFA) content; limited research has focused on carbon migration and functional gene changes during the AD of straw. Results Carbon migration and changes in metagenomic characteristics during the AD of rice straw (RS) were investigated. Accumulated biogas production was 388.43 mL/g VS. Carbon in RS was consumed, and the amount of carbon decreased from 76.28 to 36.83 g (conversion rate 51.72%). The degree of hydrolysis rapidly increased during the first 5 days, and a large amount of carbon accumulated in the liquid phase before migrating into the gas phase. By the end of AD, the amount of carbon in the liquid and gas phases was 2.67 and 36.78 g, respectively. According to our metagenomic analysis, at the module level, the abundance of M00357, M00567, M00356, and M00563 (the modules related to the generation of methane) during AD were 51.23–65.43%, 13.96–26.88%, 16.44–22.98%, and 0.83–2.40%, respectively. Methyl-CoM, 5-methyl-5,6,7,8-tetrahydromethanopterin, and Acetyl-CoA were important intermediates. Conclusions Carbon was enriched in the liquid phase for the first 5 days and then gradually consumed, and most of the carbon was transferred to the gas phase by the end of AD. In this study, AD proceeded mainly via aceticlastic methanogenesis, which was indicated to be a dominant pathway in methane metabolism. Batch AD could be divided into three stages, including initiation (days 1–5), adaptation (days 6–20), and stabilization (days 21–50), according to biogas production performance, carbon migration, and metagenomic characteristics during AD.


Background
In recent years, anaerobic digestion (AD) technology has been recognized as an environmentally friendly method to control greenhouse gas emissions and produce clean energy in the form of biogas [1]. AD technology can address the environmental issues caused by the dumping and burning of straw and exhibits immense potential for the production of biogas with high methane content [2]. The efficient conversion of straw to methane is limited by the recalcitrant lignocellulosic structure of straw, resulting in the low susceptibility of its lignocelluloses to hydrolysis [3]. However, these problems of poor biodegradability and low efficiency of AD when using straw have been widely researched [3][4][5], whereas research into the fundamental process changes involved in AD is limited.
To elucidate the mechanism of AD, many researchers have focused on process parameters as AD proceeds. For instance, Boe et al. monitored the process parameters of biogas production, methane content, pH value, and volatile fatty acid (VFA) content in the AD of cow dung and analyzed the acetic acid and propionic acid content; the results indicated that a high organic loading rate could inhibit the AD process [6]. Lv et al. monitored biogas production, methane content, pH value, VFAs, and other process parameters in a continuous stirred tank reactor (CSTR) used for the AD of corn silage and found that ammonia inhibition occurred in the reactor [7]. Li et al. investigated biogas production, biogas components, pH values, VFAs, ammonia nitrogen concentration, and other process parameters in a CSTR system used for the AD of vegetable waste and found that the reaction was inhibited due to its high organic loading rate [8]. However, most of these studies of process parameters focused on CSTRs, and the main parameters investigated were biogas production, methane content, pH value, VFAs, and ammonia nitrogen concentration. Carbon is an important component of straw. In AD, solid carbon in straw is first hydrolyzed and transformed into a liquid state before being transformed into CH 4 and CO 2 by microorganisms. The proportion of carbon converted differs among the three phases of the AD process. The conversion law of carbon elements in the AD process can provide a basis for the further clarification of the law of AD of straw and should be investigated. However, relatively few studies have investigated carbon migration and changes in parameters during AD of straw.
As promoters of AD, microorganisms play a key role in this process. Many researchers have conducted considerable research into the microbial community involved in the AD system. Chen et al. (2017) reported that Euryarchaeota and Firmicutes formed the dominant microbial communities at pH levels of 7-8 at high temperatures, with abundances of 58-60% and 19-23%, respectively, and the abundance of Firmicutes reached 60% when the pH increased to 9 [9]. Cai et al. discovered that the addition of Fe, Co, Ni, Mo, Se, or Mn during the AD of straw affected the relative abundance of bacteria and archaea in the system [4]. Zhou et al. [10] demonstrated that Firmicutes were the dominant bacteria in a two-phase system, whereas Proteobacteria were the dominant bacteria in a one-phase system [11]. However, studies of microbial communities during AD have mainly focused on comparisons among different AD conditions, and no reports have investigated the microscopic mechanisms during batch AD processes. According to the literature, the succession of microbial communities during the AD process is very complex; therefore, metagenomics may be used as an alternative approach. Wang et al. conducted a metagenomic analysis of the AD of weed silage, and the results showed that the composition, structure, and abundance of dominant bacterial communities in an AD system changed under different environmental conditions [11]. Fontana et al. performed a metagenomic study of one-and two-phase thermophilic AD of cheese wastewater; the results showed that the microbial community structures were significantly different between the two types of AD reactors and that increased microbial diversity was observed in the two-phase AD reactors [12]. Zhang et al. found that a high concentration of ammonia nitrogen was beneficial for the accumulation of VFAs but can inhibit Acetyl-CoA and Methyl-CoM gene expression, thereby significantly inhibiting the methane metabolism pathway [13]. Hu et al. analyzed the biodegradation pathway in an optimized activated carbon-containing AD reactor using Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis and found that the abundance of the pathways involved in the conversion of propionic acid into acetic acid was increased in the activated carbonenhanced reactor and that the pathways involved in lipid and methane metabolism were also enhanced [14]. However, most metagenomic analyses exclusively focused on microbial communities during digestion; in addition, the succession of metabolic pathways during the AD process has not been revealed to date.
The laws of carbon migration and knowledge of the functional gene changes during the AD of straw will provide a new perspective for the study of the micro-mechanisms that occur in this process. Therefore, the present study focused on carbon migration and metagenomic characteristics during the AD of rice straw (RS) to elucidate this process from the perspectives of the conversion of materials and functional characteristics. The aim of the present study was to present a thorough understanding of the AD of straw and provide a theoretical reference for further research into improving the performance of AD of straw.

Biogas production performance
The biogas-producing performance of the AD of RS is shown in Fig. 1. Three peaks in gas production were noted during the AD process. The first peak in daily biogas production occurred on day 3, the second peak occurred on day 15, and the third peak occurred on day 27. During the first peak, the volume of gas produced was 26.81 mL/g VS; however, the methane content was low at 33.37%. During the second peak, the volume of gas produced was 14.12 mL/g VS, and the methane content was 59.73%. During the third peak, the volume of gas produced was less than that produced in the previous two peaks at 11.95 mL/g VS with the methane content of 56.30%. These results are consistent with those reported by Hu et al., who showed that three peaks appeared as the digestion process proceeded with the production of biogas decreasing in sequence [15]. The methane content gradually increased from 33.37% on day 1 to a maximum value of 64.59% by day 14; thereafter, the levels began to decline, reaching a stable state after day 20 and remaining within the range of 49.89% to 56.51%, which is highlighted in Fig. 1. The cumulative volume of gas produced by the end of the AD was 388.43 mL/g VS, which is similar to the results reported by Dai et al. [16].
The pH level, ammonia nitrogen content, and total alkalinity are important indexes that can be used to evaluate the stability of anaerobic digestive systems and microbial metabolism. The pH, ammonia nitrogen content, and total alkalinity of the liquid discharged following AD were measured. The pH of the liquid was 7.33, and the alkalinity was 4000 mg/L. Both of these values were within the normal range (pH > 6.8, alkalinity > 2000 mg/L). However, the ammonia nitrogen content was 358 mg/L, which did not exceed the tolerance range of anaerobic microorganisms (2 g/L). The methane production process was not inhibited [17], and the system reached a stable state following AD.
The biogas-producing performance of the AD of RS determined in the present study was similar to other results reported in the literature [18,19]. These results indicate that the AD process was successful, and subsequent analysis based on this study is justified. If the analysis was performed by considering 5-day units, the AD process can be divided into three stages, which are indicated by the vertical dotted line in Fig. 1. The daily biogas production during the first unit of digestion time exhibited large fluctuations. The methane content gradually increased, and the first rapid increase in cumulative biogas production occurred. This increase may be attributed to the fact that bacteria are more active than archaea at the beginning of AD, leading to the accumulation of a large amount of acidic substances and inhibiting the progression of methane production. During the next three units of digestion time, the range in the fluctuations of daily biogas production decreased. The methane content increased to the maximum observed value in the present study and subsequently decreased to 56.1%, and the cumulative biogas production rapidly increased for the second time. This process potentially occurred due to the gradual acclimatization of the microorganisms to the system and stabilization of the synergy among the microorganisms. During the remainder of the unit digestion time, the daily biogas production decreased following the transitory increasing period and finally approached zero. The methane content tended to be stable, and the cumulative biogas production increased gradually after the third period of rapid increase. This finding could be attributed to the relatively stabilized microbial composition; however, the amounts of components that could be used for biological purposes in the RS raw material gradually decreased, and finally AD ceased.

Distribution of carbon among the three phases
Total organic carbon (TOC) reflects the amount of carbon in the liquid phase during AD because almost all the carbon in the liquid phase exists in the form of organic carbon [20]. The amount of carbon in the gas phase is indicated by the amount of carbon in the form of CH 4 and CO 2 . Carbon in the solid phase is contained in the straw; the quantity can be calculated by subtracting the quantity in the liquid and gas phases from the initial quantity of total carbon. Based on the above assumptions, the distribution of carbon in the gas, liquid, and solid phases during AD is shown in Fig. 2a. The amount of carbon in the liquid phase rapidly increased at the beginning of AD, reaching a maximum value of 18.81 g on day 5. Thereafter, the level gradually declined until it was maintained at approximately 3 g after day 30. At the end of AD, the amount of carbon in the gas phase reached 36.78 g, whereas only 2.67 g remained in the liquid phase, indicating that most of the carbon in the liquid phase that had been hydrolyzed from the solid phase was converted into the gas phase. The amount of carbon in the solid phase changed from 76.28 to 36.83 g with a conversion rate of 51.72%, which was quite similar to the VS degradation rate of 54.57%.
Carbon exists in two forms in the gas phase, CH 4 and CO 2 , and a floating diagram of carbon enriched in the target product CH 4 is shown in Fig. 2b. At the beginning of AD, carbon enriched in CH 4 was reduced by 96.55% compared with that enriched in the by-product CO 2 . As AD progressed, carbon enriched in CH 4 continually increased. By the end of AD, on day 50, 18.46 g of carbon was enriched in CH 4 ; i.e., 0.76% higher than that enriched in the by-product CO 2 .
The degree of hydrolysis (DH) and degree of gasification (DG) during AD of RS along with data fitting results obtained using the first-order dynamics model, modified Gompertz model, and Cone model are shown in Fig. 3. The DH sharply increased from 8.20% to 36.49% over the first 5 days, and then the rate of increase in the DH tended to level off. At the end of AD on day 50, the DH was 51.72%. Conversely, the DG steadily increased throughout AD; however, this increasing trend slightly slowed during the final stages of AD. By the end of day 50 of AD, the DG was 48.22%. Relevant data from these dynamic models were calculated and analyzed, and the parameters estimated by the different models are shown in Table 1.
According to the data fitting results generated by each model, the maximum predicted DH for the AD of RS was in the range of 45.67% to 53.49%, and the Cone model predictions were the closest to the actual situation. The Cone model generated the highest R 2 value (0.960) and exhibited a good correlation with the experimental data. Moreover, the Cone model exhibited the lowest values of root mean square prediction error (RMSPE) and the Akaike information criterion (AIC), demonstrating that it was the best of the three models for predicting the DH associated with the AD of RS.
The maximum predicted DG from the AD of RS was in the range of 49.54% to 71.56%. The modified Gompertz model output was the closest to the actual situation; the predicted maximum DG in the other two models was higher than the predicted maximum DH in the Cone model, which is not consistent with the actual condition in AD. The R 2 value of the modified Gompertz model was 0.9891, demonstrating that it was the best model for predicting the DG of the AD of RS.

Metagenomic analysis
A total of 48.81 GB of read data were generated, including 447985312 original sequences and 67645782112 original basic groups. After sequence quality control, 443321486 sequences and 66703997126 basic groups remained. High-quality sequences and basic groups comprised > 98% of the total. A total of 7484761 gene sequences were predicted, and 3447695 non-redundant gene catalogs were constructed. In specific, 533761 KEGG genes and 2684 genera were annotated.
The composition of the pathways on levels 1 and 2 is shown in Fig. 4a. Metabolism was the dominant pathway on level 1, and Amino Acid Metabolism was the dominant pathway on level 2, which is consistent with the results reported by Li et al. [21]. Carbohydrate Metabolism was the second most dominant pathway on level 2, which is consistent with the chemical components of straw.
The pathway with the highest abundance in Carbohydrate Metabolism was Glycolysis/Gluconeogenesis (ko00010, 1.71%-1.99%), which was consistent with the conclusions obtained by Soares et al. [22]. Moreover, ko00620, ko00520, ko00630, and ko00020 were pathways with high abundance on level 3. This finding indicates that Carbohydrate Metabolism is a complex process that involves a series of complicated steps that require catalysis by multiple enzymes and is accompanied by the generation of ATP and NADH [23].
At the Module level, four Modules are related to the generation of methane: M00357, M00567, M00356, and M00563. M00357 uses acetate as a substrate. M00567 uses H 2 or CO 2 as a substrate. M00356 uses methanol as a substrate. M00563 uses methylamine, dimethylamine, and trimethylamine as substrates. The abundance of these four modules is shown in Fig. 4b. M00357 accounted for the highest proportion with abundance values in the range of 51.23 to 65.43% followed by M00567 with an abundance of 16.96 to 26.88%, M00356 with an abundance of 16.44 to 22.98%, and M00563 with the lowest abundance of 0.83 to 2.40%. This finding indicated that the methane production process in the present study mainly involved aceticlastic methanogenesis. The abundance of archaea at the genus level is shown in Fig. 4c. Methanosaeta was the most abundant archaeal genus, which uses acetic acid as a substrate [24], followed by Methanosarcina and Methanobacterium. These findings were consistent with the results of the KEGG functional analysis and further indicated that the methanogenic process in the present study mainly involved aceticlastic methanogenesis. These results are consistent with the results reported by Chen et al. [25].
The abundance of Methyl-CoM slightly increased between days 1 and 5, and fluctuations subsequently occurred. The minimum value was reached on day 20 before Methyl-CoM began to continuously increase. This trend is strongly correlated with the trend of biogas production, except for the increase in the abundance of Methyl-CoM and the decrease in biogas production after 25 days. This finding may be due to the consumption of the substrate despite the increase in the capacity for methanation in the AD system. On the other hand, the stage at which the abundance of Methyl-CoM steadily increased was the stage at which the methane content in the biogas was becoming relatively stable. Therefore, the abundance of Methyl-CoM exhibited a strong correlation with biogas production and methane content in the AD system and may be a useful indicator to determine whether the AD system has entered the stage of stable methane production.
CO was produced during reaction R09096 with the generation of 5-methyl-5,6,7,8-tetrahydromethanopterin. CO can be converted to CO 2 via a reaction involving EC 1.2.7.4, which is highly abundant in the AD of RS, and the relative abundances of EC  10, were present at a low abundance; thus, a limited amount of CO 2 enters the hydrotropic methanogenic pathway, resulting in high CO 2 content in biogas and low efficiency of methane production.
In addition to acetic acid, small molecule organic compounds (SMOC), such as propionic acid and butyric acid, are also important intermediate products in the hydrolysis and acidification stages of AD. The level at which these compounds are present can directly affect AD efficiency. At the enzyme level, 31, 7, and 5 biological processes associated with acetic acid, propionic acid, and butyric acid metabolism, respectively, were identified in the present study, and their abundances are shown in Table 3. The highest abundance was observed for enzymes associated with acetic acid metabolism followed by enzymes associated with propionic acid metabolism. A total of 20 enzymes associated with acetic acid production were identified. The total abundance of these enzymes increased significantly between days 1 and 5, decreased slightly between days 5 and 15, and then increased again until it reached the maximum value on day 25 before decreasing to the minimum value by day 50. The first 15 days showed a trend of an initial increase followed by a decrease, indicating that the production of acetic acid continued for 1-5 days. Subsequently, the microorganisms present could not consume all the acetic acid due to its increasing accumulation, resulting in the inhibition of acetic acid production and a decrease in the total amount of acetic acid produced. During the middle stages of AD, the concentration of acetic acid-generating enzymes increased, and the accumulated acetic acid was utilized by microorganisms, leading to the more stable production of methane.

Three stages of batch anaerobic digestion
According to the biogas production performance, carbon migration, and metagenomic characteristics during the AD of RS, certain rules related to the changes in carbon content, the microbial community, and functions of this community during batch AD were observed. The beta-diversity and functions of the microbial community during the different phases of the AD of RS were analyzed. Principal component analysis (PCA) of the microbial community (Fig. 6a) showed that components in principle component axis low. This finding suggested significant differences in the microbial community structure at the different stages of the AD of RS. Furthermore, the PCA of KEGG functions (Fig. 6b) revealed significant differences in the  metabolic pathways at different stages of digestion, indicating that the pathways changed following succession in the microbial community, enabling adaptation to changes in the substrate and the digestive environment at different stages. The microbial community and functional composition of A5, A10, A15, and A20 demonstrated similar characteristics, and A25 and A50 were similar but clearly different from A0. Therefore, the entire AD could be divided into three stages: the initiation stage (days 1-5), adaptation stage (days [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], and stabilization stage (days 21-50). The result of the iPath pathway analysis of the three stages is shown in Fig. 6c According to the changing trends in the three stages, the unique pathways in the initiation stage included 21 pathways, such as ko00130, ko00360, ko00940, ko0059, and ko00330; the microorganisms corresponding to these pathways were inhibited or eliminated during AD and eventually disappeared. Numerous new pathways, including ko00600, ko100, ko00561, ko00564, and ko00565, appeared during the adaptation stage; however, these pathways were not adapted to an AD system and were eventually eliminated. No unique pathways were identified in the stabilization stage, indicating that the AD system had reached a stable state by this stage and that it is reasonable to divide the AD into three stages.

Conclusions
From the perspective of carbon migration, carbon in RS was consumed (76.28-36.83 g) with a conversion rate of 51.72%. In the first 5 days, DH rapidly increased, and a large amount of carbon accumulated in the liquid phase and subsequently was converted into gas. By the end of the AD of RS, the amount of carbon in the liquid and gas phases was 2.67 and 36.78 g, respectively. From the perspective of metagenomic characteristics, the Methyl-CoM, 5-methyl-5,6,7,8-tetrahydromethanopterin and Acetyl-CoA were important intermediates during the AD of RS. According to the changes of carbon migration and metagenomic characteristics, the batch AD of RS could be divided into three stages (initiation, adaptation, and stabilization).

Methods
Mesophilic AD of RS was performed, and the carbon migration and changes in metagenomic characteristics during the process were investigated to study the mechanisms involved in the AD of RS.

Feedstock and inoculum
The RS was collected from Ji County (Tianjin, China). After natural air-drying, the RS was cut into 3-to 4-cmlong sections using a straw chopper. These sections were then smashed to the size of 20-mesh using a pulverizer (YSW-180). The inoculum was obtained from a biogas station in Shunyi district (Beijing). The characteristics of the raw materials and the inoculum are shown in Table 4.

Experimental equipment
Glass bottles with a volume of 5 L were connected to a gas-collecting device by latex and glass tubes and used as batch reactors. The discharge interface was connected by a latex tube and sealed with a water stopper.

Experimental methods
Raw materials and the inoculum were added to the reactor; the organic loading rate of the pretreated RS was 50 gTS/L, and that of the inoculum was 20 gTS/L. Pure water was added to adjust the volume to the working volume of the AD reactor (4 L), and the initial pH in the reactor was adjusted to a value between 6.8 and 7.2 by adding Ca(OH) 2 . The AD reactor device was placed in a constant-temperature water bath at 35 °C ± 1 °C to allow mesophilic AD. The AD process was allowed to continue for 50 days. Three parallel experiments were established for each group of experiments.

Analytical methods
Daily biogas production was measured using the water displacement method [26]. The biogas composition was determined using gas chromatography (SP-2100, BeiFenRuiLi, Beijing, China). The American Public Health Association standard method was used to determine total solids (TS) and volatile solids (VS) [26]. The total organic carbon (TOC) content was measured by the multiple dilution method using

Degree of hydrolysis and degree of gasification
The specific quantity of carbon in the cumulative gas produced together with the amount of carbon in the liquid of the AD system and the total quantity of carbon in the RS was defined as the DH of the straw during the AD process. The specific quantity of carbon in the cumulative gas produced together with the total quantity of carbon in the RS was defined as the DG of the straw during the AD process. A brief description of DH and DG is shown in Eqs. 1 and 2: where C(g) is the amount of carbon in the cumulative gas produced (unit: g); C(l) the amount of carbon in the liquid phase (unit: g); and C(t) the amount of the carbon in the RS material (unit: g).

Dynamic analysis
A first-order reaction kinetics model [27], the modified Gompertz model [28], and the Cone model [29] are the most widely used models for the analysis of AD of complex organic matter. These models are typically used to analyze the degree and rate of organic matter degradation. DH and DG were fitted by the three models.
The formula of the first-order dynamic model is shown in Eq. 3: where P (t) is the degree of hydrolysis or of gasification of the material on day t (%); P m the maximum degree of (1) DH(%) = C(g) + C(l) C(t) × 100% (2) DG(%) = C(g) C(t) × 100% (3) P (t) − P m (1 − exp (−kt)) hydrolysis/gasification (%); k the first-order reaction constant (day −1 ); and t the AD time (days). The modified Gompertz model formula is shown in Eq. 4: where P (t) is the degree of hydrolysis/gasification of the material on day t (%); t is the digestion time (d); P m is the maximum degree of hydrolysis/gasification (%); R m represents the maximum hydrolysis/gasification rate (%); λ is the delay time (days); and e is the natural constant, 2.71828.
The Cone model formula is shown in Eq. 5: where P (t) is the degree of hydrolysis/gasification of the material on day t (%); P m is the maximum degree of hydrolysis/gasification (%); k represents the rate constant (d −1 ); n is the shape factor (dimensionless); and t is the digestion time (days). The applicability of the models was evaluated by the root mean square percentage error (RMSPE) and the Akaike information criterion (AIC) [30].
The RMSPE was calculated using Eq. 6: where Pv i represents the degree of hydrolysis/gasification (%); Mv i represents the measured degree of hydrolysis/ gasification (%); and n is the number of measurements. AIC was calculated using Eq. 7 and was used to represent the degree of information missing from the model: where RSS represents the sum of the squares of the residuals; n is the number of data points; and N is the number of model parameters.