Thermodynamic Control on Biogeography and Functioning of Rare-Biosphere Propionate Syntrophs in Paddy Field Soils

Global biogeochemical processes are not only gauged by dominant taxa of soil microbiome but also depend on the critical functions of “rare biosphere” members. Here we evaluated the biogeographical pattern of “rare biosphere” propionate-oxidizing syntrophs in 113 paddy soil samples collected across eastern China. Results The relative abundance, functioning capacity and growth potential of propionate-oxidizing syntrophs were analyzed to provide a panoramic view of syntroph biogeographical distribution at the continental scale. The relative abundances of four syntroph genera, Syntrophobacter , Pelotomaculum , Smithella and Syntrophomonas were significantly greater at the warm low latitudes than at the cool high latitudes. Correspondingly, the functioning potential of propionate degradation was greater in the low latitude soils compared with the high latitude soils. The slow rate of propionate degradation in high latitude soils resulted in a greater fold change in increase of the relative abundance, probably due to the growth rate-yield tradeoff relationship. The mean annual temperature (MAT) is the most important factor shaping the biogeographical pattern of propionate-oxidizing syntrophs, with the next factor to be the total S content (TS) in soil. suggest that related to change, which the tension propionate oxidation is leveraged with the increase MAT. The TS effect is likely due to that some propionate syntrophs can facultatively perform sulfate respiration.


Results
The relative abundance, functioning capacity and growth potential of propionate-oxidizing syntrophs were analyzed to provide a panoramic view of syntroph biogeographical distribution at the continental scale. The relative abundances of four syntroph genera, Syntrophobacter, Pelotomaculum, Smithella and Syntrophomonas were significantly greater at the warm low latitudes than at the cool high latitudes. Correspondingly, the functioning potential of propionate degradation was greater in the low latitude soils compared with the high latitude soils. The slow rate of propionate degradation in high latitude soils resulted in a greater fold change in increase of the relative abundance, probably due to the growth rate-yield tradeoff relationship. The mean annual temperature (MAT) is the most important factor shaping the biogeographical pattern of propionate-oxidizing syntrophs, with the next factor to be the total S content (TS) in soil.

Conclusions
We suggest that the effect of MAT is related to the Gibbs free energy change, in which the endergonic tension of propionate oxidation is leveraged with the increase of MAT. The TS effect is likely due to that some propionate syntrophs can facultatively perform sulfate respiration.

Background
Soil microbiome plays the vital role in regulating global biogeochemistry. It has been recently documented that albeit immense biodiversity of soil microbiome, only a few hundreds of phylotypes are prevalent across global soils and these dominant phylotypes display distinct biogeographical distribution according to their habitat preferences [1]. The factors shaping the biogeographic patterns of global soil microbiome include climatic factors, edaphic properties and biological interactions [2][3][4][5][6].
Identification of soil dominant taxa, their geographical distributions and the controlling factors provide a "most-wanted" list for cultivation and genomic scrutiny that shall help decoding biogeochemical mechanisms and pave a way toward developing global Earth system models that are soil contextdependent [1,7,8]. The dominants-hunting approaches, however, will inevitably overlook the function of "rare biosphere" specialists [9]. It has been understood that biogeochemical processes are not only gauged by dominant taxa that operate through biomass turnover but also depend on critical functions of specialist members [7]. A typical example is microbial syntrophs, which plays a critical role in anaerobic decomposition of organic matter and methanogenesis but have only low relative abundance in natural environments like paddy field soils [10,11].
Complex organic matter in anoxic environments is degraded by an anaerobic food chain comprising primary and secondary fermentors, homoacetogens and methanogens [12,13]. A series of shortchain fatty acids and alcohols (SCFAs) are transiently accumulated as intermediate products during the process [12,[14][15][16]. Syntrophic cooperation between secondary fermentors (i.e. the syntrophs) and methanogens is essential to degrade these chemicals, which otherwise cannot be broken down by either guilds alone [13,17,18]. Propionic acid is one of major intermediates accounting for up to 30% of the total CH 4 formation in paddy soils [14,19]. The process of propionate degradation is thermodynamically the most unfavorable among common SCFAs like ethanol and butyrate, making propionate degraders being the slowest in growth and sensitive to environmental conditions [17,20,21]. Identifying propionate-oxidizing key taxa and delineating their biogeographical distribution are therefore important for better understanding and predicting organic matter decomposition and C cycling in anoxic environments including natural wetlands and paddy field soils.
Paddy field soils are human-managed ecosystems and an important contributor to global CH 4 emissions [22,23]. About 90% of paddy fields are located in Asia with 20% of that in China [24]. Rice production in China is dominated by irrigated systems, which is mainly distributed in the lowland of eastern China, extending from warm sub-tropics at 18°N latitude to cool temperate regions at 50°N [25]. The long history of rice cultivation makes the region an ideal site for microbial biogeography investigation that is independent of large scale variations of vegetation and soil. In the present study, we collected 113 paddy field soils from different regions across eastern China with a latitudinal 4 distance of 3 689 kilometers. The potentials of propionate degradation in each soil were determined through anaerobic incubations in laboratory. The composition and abundance of propionate syntrophs in original soils and soil samples after anaerobic incubations were investigated using high throughput sequencing of bacterial 16S rRNA genes. We show here that temperature and the total S content in soil are the most important factors shaping the biogeographic distribution and functioning of propionate-oxidizing syntrophs in paddy field soils at the continental scale.

Soil sampling and analysis of basic soil properties
One hundred thirteen soil samples were collected from paddy fields across eastern China in the summer from July to September of 2017. At each site, five soil sub-samples were taken at a depth of 0-15 cm and mixed thoroughly. The collected paddy soils were placed in ziplock plastic bags and transported to lab at 4 °C within 36 h. Standard soil testing procedures were followed to measure soil pH, cation exchange capacity (CEC), organic matter (OM), dissolved organic carbon (DOC), total nitrogen (TN), total phosphorus (TP), total sulfur (TS), total Fe (TFe) and fourteen additional properties [26]. Climatic variables, including mean annual temperature (MAT) and mean annual precipitation (MAP) at each sampling site were obtained using site coordinates from the WorldClim database (www.worldclim.org).

Soil slurry preparation and anaerobic incubation
To determine the potentials of propionate degradation, soil samples were subjected to anaerobic incubation under laboratory conditions with propionate as a sole substrate. Soil samples were suspended in autoclaved degassed water at a soil-to-water ratio of 1:5 (1 g dry weight soil plus 5 ml autoclaved degassed water). Aliquots (45 ml) of homogenized soil slurry were transferred into 120 ml sterile serum bottles which were capped with black butyl rubber stoppers and aluminum crimps. The sodium propionate was added into the bottles to a final concentration of 10 mM. The bottles were vigorously shaken by hand to homogenize soil slurries, flushed with N 2 for 5 min and then incubated under the dark at 30 °C without shaking. Experiment was carried out in triplicate.

5
The gaseous samples (200 μL) were regularly collected from the headspace of incubation bottles using pressure-lock analytical syringe (Baton Rouge, LA). The CH 4 concentration was analyzed by Agilent 7890B gas chromatograph (Agilent Technologies, USA) equipped with flame ionization detector (FID) and thermal conductivity detector (TCD) [27]. For the analysis of fatty acids, 0.8 ml of liquid samples were collected from the bottles with sterile syringes and centrifuged for 10 min at 6870 g at 4 °C. The supernatants were filtered through 0.22 μm cellulose membrane filters and acidified by 18.4 mM H 2 SO 4 . Acetate, propionate and butyrate were analyzed by HPLC (Agilent 1260, USA) equipped with Zorbax SB-AQ C18 column and UV detector at 210 nm [28].

Soil DNA extraction and Illumina sequencing of 16 S rRNA genes
Microbial DNA samples were extracted both from original soil samples and from soil samples after anaerobic incubations. For the original soil samples, 0.5 g of fresh sample was used for DNA extraction using the FastDNA SPIN Kit (MP Biomedicals, USA) according to the manufacturer's protocol [26]. For the soil samples after anaerobic incubations, 2 ml of soil slurries were collected at the time when at least 90% of propionate was consumed in each soil incubation. The date of sampling differed among soil samples because the time for propionate consumption varied markedly. The collected soil slurries were transferred to the sterile centrifugal tube and centrifuged for 10 min at 6 870 g at 4 °C.
Total genomic DNA was extracted from the precipitated samples using the FastDNA SPIN Kit. DNA samples extracted in triplicate for each soil sample were pooled in equal-volume and then used for high-throughput sequencing on an Illumina Miseq platform [29][30][31]. Primer set Ba338F (5'-ACTCCTACGGGAGGCAGCAG-3') and 806R (5'-GGACTACHVGGGTWTCTAAT-3') was employed to amplify V3-V4 regions of the bacterial 16S rRNA genes. Sequences were split into groups according to taxonomy and assigned to operational taxonomy units (OTUs) at a 3% dissimilarity level using the UPARSE pipeline [32]. OTUs with fewer than two sequences were removed, and representative sequence of each OTU was assigned to a taxonomic lineage by RDP classifier against the SILVA database.

6
To correct for sampling effort in DNA sequencing (number of analyzed sequences per sample), the sequence number was rarefied at 27 812 high-quality sequences per sample for the original soil samples and 18 441 high-quality sequences per sample for the samples after anaerobic incubations, respectively [26]. Bacterial α-diversity (Shannon index) was calculated using QIIME (http://qiime.org/index.html). Bacterial β-diversity was estimated according to the Bray-Curtis distance between samples. Metacommunity co-occurrence network analysis at OTU or genus levels was based on Pearson's correlation coefficient between microbial taxa and visualized by Gephi 0.9.2 [33]. The circle phylogenetic tree was built with FastTree using a GTR model and visualized by iTOL after the sequence cutoff and gap trimming treatment by Silva Incremental Aligner (SINA v1.2.11) and trimAl, respectively [34,35]. The phylogenetic trees were constructed using the neighbor-joining algorithm with 1 000 Bootstrap replications in MEGA 7.0.14.
To build predictive maps of the spatial distributions of the bacterial syntroph taxa, we used a Kriging interpolation method with the "automap" package in R to estimate the relative abundance and the changes of syntrophs in paddy soils. Kriging interpolations were performed with the best-fitted variogram model through cross validation, which was performed using autoKrige.cv in the "automap" package by setting the parameter "nfold" to 113, the number of soil samples in analysis [36]. The correlation between estimated values of sampling site in cross validation and the observed values at the corresponding sites were displayed with the Pearson correlation coefficient and P value in the map. The predictive maps were also build for the functional potential of propionate degradation, which was evaluated based on the time lapse for CH 4 production from propionate oxidation and the maximum rate of methanogenesis during anaerobic incubation and for the Gibbs free-energy changes (△G 0' ) of propionate oxidation, which were calculated by setting standard conditions except temperature.
Regression and correlation analyses were used to delineate the relationship between environmental factors, syntroph abundances and functional potentials. Boruta algorithm implemented in R package "boruta" was used to identify the relevant environmental factors and estimate their importances to 7 the relative abundance of propionate syntrophs in paddy soils [37].
The processing and visualization of all data were conducted in the R environment (v3.6.0; http://www.r-project.org/) using the following packages unless otherwise indicated: ggplot2 [38], Boruta [37], corrplot [39], automap [40], pheatmap [41] and psych [42]. and up to C 10 fatty acids [43]. The inclusion of Syntrophomonas in analysis was due to the fact that Smithella, which utilize the C 6 dismutation pathway, can release butyrate as an intermediate product into environment, which is then metabolized by Syntrophomonas [44]. A total of 30 OTUs belonging to five syntroph genera were retrieved.

Biogeographical distribution of propionate syntrophs
The total abundance of all five syntroph genera together (referred to synTotal) ranged from 0.01 to 0.34% across 113 soil samples. This low abundance certificates that syntrophs in paddy soils belong to a subset of "rare biosphere" community (defined as individual species or OTUs accounting for <0.1% of the total abundance). Albeit "rare" in relative abundance, the synTotal displayed a distinct geographical distribution, showing the highest abundance in the warm low latitude soils and gradually decreasing towards the high latitude regions (Fig. 1A). The conventional geographical division separates China into the south (low latitude) and the north (high latitude) regions according to the Qinling Mountains-Huaihe River Line (latitude ≈ 32°), a critical geographical boundary for climate, landform and soil conditions [45]. Accordingly, we summarized our soil samples into the low latitude and the high latitude groups, respectively. The mean relative abundance of the low latitude group is 8 significantly greater (by 1.84-fold) than that of the high latitude group (bottom-left of Fig. 1A). In line with the synTotal, the α-diversity of syntrophs, estimated based on OTU richness, significantly decreased with the increase of latitude (Fig. 1B). We then identified the key environmental factors related to the synTotal distribution through Spearman correlation analysis. MAT and TS (total soil S) were identified as the two most important factors followed by NH 4 + -N, soil OM, MBC and TN (Fig. 1C).
Soil CEC and pH showed the negative correlations. The importance index estimated based on Boruta algorithm reiterated that MAT and TS were the two most important factors linking to biogeographical distribution of synTotal (Fig. 1C).
To compare the distribution of different syntrophs, the relative abundances of five syntroph genera Syntrophobacter showed the greatest mean relative abundance followed by Syntrophomonas and Smithella, while Pelotomaculum and Desulfotomaculum showed very low abundances (Additional file 1: Supplementary Figure 3F).

Biogeographical feature of syntroph functioning potentials
Next we evaluated the biogeographical feature of syntroph functioning potentials. For this purpose, fresh soil samples were incubated under anaerobic conditions with addition of 10 mM propionate. We observed the CH 4 production until >90% of the added propionate consumed. Three distinct patterns of methanogenesis were identified according to the time lapse required for CH 4 production from propionate oxidation (Fig. 2). The pattern I comprised 45 soil samples with the time lapse of 13 d to 27 d, representing the fast rate group of syntrophic metabolisms ( Fig. 2A). A greater proportion of soil 9 samples tends to be located at the low latitudes in this group. The pattern II, including 51 samples, represented the median rate of syntrophic metabolisms (28 d to 43 d), which did not show a distinct latitudinal tendency (Fig. 2B). The pattern III comprising 17 soil samples required 44 d to 82 d for CH 4 production from propionate oxidation and hence represented the slow rate group. The soils of this group were distributed mainly in the cool high latitude regions (Fig. 2C). The analyses of propionate consumption confirmed the separation of soil samples into three groups (Fig. 3A). Acetate was a major intermediate detected in all samples (Fig. 3B). Butyrate was detected in 28.3% of soil samples (detection limit of 0.05 mM) (Fig. 3C), with the highest concentration occurring in the pattern III soils (Fig. 3C). Taking all soil samples together, the time lapse for methanogenesis from propionate degradation displays an explicit biogeographical pattern, being significantly slower in the high latitude soils than in the low latitude soils (Fig. 4A). The linear least squares regression revealed that the time lapse for methanogenesis was significantly negatively correlated with the relative abundance of synTotal in original soils, MAT, TS, and other edaphic factors including soil OM, MBC, available Fe, Cu and Mn (Fig. 4C). The functional potential of propionate degradation can also be inferred from the maximum rate of methanogenesis (Fig. 4B), which supports the biogeographical tendency revealed by the time lapse (i.e. the shorter the time lapse, the greater the maximum rate).

Shift of microbial community during the anaerobic incubation
At the end of anaerobic incubations, the structure of microbial community in soil samples was  Figure 4C). The co-occurrence network analysis of the top 10% OTUs confirmed that the OTUs for the low and high latitude regions tended to group separately (Additional file 1: Supplementary Figure 4D). These results indicate that the bacterial community shifted significantly during anaerobic incubation, while the separation of metacommunities into the low latitude and the high latitude groups remained robust.
The relative abundance of syntrophs except Desulfotomaculum increased markedly after the incubation with the maximum abundance of synTotal reaching to 12.5%, indicating significant growth.
We estimated the increase of individual syntrophs by calculating the logarithmic (log 2 ) fold change (R) of the relative abundance over the incubation. The relative abundance of Desulfotomaculum, which was low at the beginning, did not change over the incubation, indicating no response to propionate. Other four genera showed significant increases but to different extents (Fig. 5).
Pelotomaculum showed the greatest increase, followed by Smithella and Syntrophomonas, while Syntrophobacter exhibited the least increase (Fig. 5F). Strikingly, the relative increases of syntrophs showed an opposite geographic tendency compared with their abundances in original soils (Fig. 2).
The log 2 R values for Syntrophobacter, Syntrophomonas and Smithella increased with the increase of latitude (Fig. 5A,B,D). The values for Pelotomaculum also increased with latitude, though showing the highest value at the middle latitude (Fig. 5C). These results indicate that propionate syntrophs in anaerobic incubation grew to a greater extend in soil samples from the high latitudes than those from the low latitudes.

Discussion
Given their low abundance in paddy soils in situ but playing the critical role in oxidizing intermediate products during anaerobic decomposition of organic matter, propionate-oxidizing syntrophs represent the typical example of rare but important taxa of soil microbiome. Here we show that the relative abundance of propionate syntrophs is significantly correlated with MAT (Fig. 1). Based on the time lapse for methanogenesis from propionate oxidation, we further reveal that the biogeographical pattern of syntroph functioning potentials (Figs. 2,3) is in accordance with their relative abundances in soils. It has been documented that within the physiological range, the microbial diversity, metabolic activity and population growth rates in terrestrial ecosystems increase exponentially with temperature [46]. Hence, climate warming is expected to accelerate temporal turnover and divergent succession of microbial communities in soils [47,48]. However, such a temperature impact can be universal to all microbes. Syntrophs are known to have a specific lifestyle by living at thermodynamic limit. In order to obtain a deeper insight into temperature effect, we calculated the standard Gibbs free-energy changes (△G 0' ) for the reaction of propionate oxidation [CH 3 CH 2 COO − (aq) + 2H 2 O(l)-CH 3 COO − (aq) + 3H 2 (g) + CO 2 (g)] (Fig. 6A). Temperature was assumed to be an only variable with other variables set to standard conditions and the △G 0' was calculated as described [49]. We found that the total S in soil was the second important factor. This was possibly due to the fact that many of propionate syntrophs like Syntrophobacter fumaroxidans and Pelotomaculum thermopropionicum are facultative sulfate reducers [50][51][52]. In the presence of sulfate these organisms tend to perform anaerobic sulfate respiration instead of syntrophy with methanogens for maximizing energy conservation [52][53][54][55][56]. Notably, while the relative abundances of Syntrophobacter, Smithella and Syntrophomonas in paddy soils were significantly correlated with TS, that of Pelotomaculum did not (Fig. 1C). Therefore, individual syntrophs have different responses to TS.
Besides TS, the relative abundances of syntrophs were also positively correlated with the contents of soil OM, total N, MBC, and in addition the functioning potential was correlated with the contents of trace elements Fe, Cu and Mn (Fig. 4C). Since a high soil OM content is usually associated with high TN, MBC and trace element contents and even TS, these edaphic factors are possibly co-variables, and the positive correlations possibly indicate the favorable living conditions for most microbes 12 including syntrophs.
Anaerobic incubation caused a marked shift of soil bacterial community (Additional file 1: Supplementary Fig. 4). Specifically, the relative abundances of syntrophic populations increased by up to 284-fold over the incubation. The increase in relative abundances, however, differs among five syntroph genera. Pelotomaculum showed the greatest growth, followed by Smithella and Syntrophomonas while that of Syntrophobacter and Desulfotomaculum being the least (Fig. 5). These results indicate that Pelotomaculum species are the most important syntrophs metabolizing propionate in paddy soils across a wide geographical location. Smithella in combination with Syntrophomonas are the next important group, while Syntrophobacter play only a minor role and that of Desulfotomaculum is negligible. Pelotomaculum and Smithella are known to use the methylmalonyl-CoA (MMC) pathway and the C 6 dismutation route, respectively [20,44,57,58].
Smithella are considered to tolerate a stricter thermodynamic condition than the MMC-utilizing syntrophs [59]. In supporting this idea, we found that the maximum growth of Pelotomaculum happened at the middle latitudes while that of Smithella and Syntrophomonas occurred at the further north (Fig. 5).
The increase in relative abundances, expressed as log 2 R, tends to be greater in the high latitude soils relative to the low latitude soils (Fig. 5), in contrast to the geographical distribution of relative abundances in original soils (Fig. 1). Two plausible explanations may cause this result. First, with identical quantity of propionate consumed, a greater fold change in relative abundance could be expected for soils having the lower initial abundances. Second, the anaerobic incubations were stopped when > 90% of propionate were consumed, which made the duration of incubation varying from 13 d to 82 d across different soils (Fig. 2). With identical quantity of substrate consumption, a longer duration would be expected to produce a greater biomass (yield) [60]. Such a rate-yield tradeoff relationship has been demonstrated in E. coli [61,62] and methanogens [63]. It remains, however, unknown if such a relationship also occurs for propionate syntrophs in paddy soils in situ.

Conclusions
In the present study, we demonstrate that both the relative abundance and functioning potential of

Ethics approval and consent to participate
Not applicable

Consent for publication
Not applicable

Availability of data and material
The datasets supporting the conclusions of this article are available in the NCBI Sequence Read Archive under BioProject PRJNA544819 and PRJNA601098 that are publicly accessible at https://www.ncbi.nlm.nih.gov. R codes on the statistical analyses are available at https://github.com/jinyidan/Rare-Biosphere-Propionate-Syntrophs-in-Paddy-Field-Soils.

Competing interests
The authors declare that they have no competing interests.    for the low latitude regions (South) and the high latitude regions (North), respectively. The red asterisks denote the significant differences between the "South" and "North" groups: * 0.01 < P ≤ 0.05, ** P ≤ 0.01.

Figure 6
The geographical feature of standard Gibbs free energy change for propionate degradation.