Microbial Interactions and Roles in Soil Fertility in Seasonal Freeze-Thaw Periods under Different Straw Returning Strategies

: With the increase of world food demand, the intensity of cultivated land use also increased. To improve soil nutrient concentrations and crop yield, several straw returning techniques have been developed. Studies have shown that straw returning is beneﬁcial to soil, but few studies have focused on the relationship between microbes and fertility in seasonal freeze-thaw periods. A two-year cropland experiment was set up that comprised three different straw return strategies, namely covering tillage with straw return for two years (CS), rotary tillage and straw return for two years (RS), rotary covering tillage with straw return (ﬁrst year covering and the second year rotary tillage) (CRS), and conventional tillage with no straw return (CK). Illumina Miseq high throughput sequencing of 16S rRNA was applied to assess bacteria community structure. The relationship between bacteria community structure and changes in soil fertility induced by different straw incorporating during seasonal trends was studied. Our results showed that soil bacterial communities varied signiﬁcantly during the soil seasonal freeze-thaw period in the northwest of Jilin province, China, and were inﬂuenced, to some extent, by the different straw returning procedures. Multidimensional analysis revealed that total phosphorus (TP), available nitrogen (AN), and total nitrogen (TN) were the major drivers of bacterial community structure. The co-occurrence network was divided into several modules. Notably, the major bacterial modules varied signiﬁcantly in different sampling periods and different treatments. These results suggested that speciﬁc bacterial groups could contribute to soil fertility in relation to environmental ﬂuctuations. Some bacterial groups (e.g., Pyrinomonadales , Rhizobiales , Sphingomonadales , and Xanthomonadales , in order level) were directly linked with speciﬁc environmental factors, indicating the key roles of these groups in soil fertility. In summary, the soil bacterial communities varied signiﬁcantly during the freeze-thaw period and might play important roles in the degradation of straw. Thus, the straw return could enhance soil fertility.


Introduction
The soil microbiome is one of the most complex and dynamic microbiomes on earth [1]. Environmental conditions significantly influenced microbial community structure in the soil ecosystems. It has been revealed that an optimum condition for each microbe is the best for its growth and activities [2]. Seasonal soil temperature and moisture variations

Site Description
The field experiment began in early 2017 at the experimental site of the Jilin Academy of Agricultural Sciences (45 • 69 N, 122 • 86 E), located in the northwest of Jilin province, China ( Figure 1). The test area has a temperate continental monsoon climate, with an average annual precipitation of 399.9 mm, an average annual temperature of 5.2 • C, an average active accumulated temperature above 10 • C of 2996.2 • C·day, an average frost-free period of 144 days, and an average annual sunshine duration of 2915 h. The temperature and precipitation of sampling durations were shown in Figure S1. The soil type is light chernozem, and the essential chemical properties are as following TN 1.04 g kg −1 , TP 0.31 g kg −1 , total potassium (TK) 24.79 g kg −1 , SOC 11.60 g kg −1 , and pH 6.40.
Agriculture 2021, 11, x FOR PEER REVIEW 3 of 16 and precipitation of sampling durations were shown in Figure S1. The soil type is light chernozem, and the essential chemical properties are as following TN 1.04 g kg −1 , TP 0.31 g kg −1 , total potassium (TK) 24.79 g kg −1 , SOC 11.60 g kg −1 , and pH 6.40.

Experimental Design and Sample Collection
A control and three different treatments were designed as follows: (i) conventional tillage-no straw returning (NT, as CK); total maize straw removed from the field, then rotary cultivating two times by using rotary cultivator (YTO Group Corporation, China, working depth was about 10-15 cm); (ii) covering + straw returning for two years continuously (CS); the straws were used to cover the field; (iii) rotary tillage + straw returning for two years continuously (RS); the straw was crushed and extruded into small pieces with less than 5 cm, then rotary cultivating two times by using IGQN-200K-QY rotary cultivator with a working depth of 10-15 cm; (iv) covering + rotary tillage + straw returning (CRS, first year CS and second year RS, straw incorporation treatment was repeated every year). Each treatment plot was 33 m long and 20 m wide in triplicate. In the present experiment, 9600 kg ha −1 of above-ground parts of maize were incorporated into the field. Maize (Zea mays L. cultivar Xiangyu '998′) was mono-cultured in early May at a density of 70,000 plants ha −1 , and then harvested in October. The experimental field was managed in accordance with conventional management methods.
In order to evaluate the soil condition after two consecutive tillage seasons, soil samples were gathered randomly from every treated plot from August 2018 to May 2019. The specific sampling times were the growing season (20 August 2018), freezing period (12 December 2018 and 19 March 2019), thawing period (2 April 2019), and later thawing period (3 May 2019). We collected 20 samples from four tillage methods, including conventional tillage and straw returning cropland (CK, CS, RS, and CRS). Samples were collected from five different points per plot at a depth between 5 cm to 25 cm and then formed a mixed sample, representing one of the three replicates of each treatment. Debris of each mixed sample (plant tissues, rocks, and roots) was removed by sieving through a 2-mm mesh in the field (Ren et al. 2016) and kept in a sterilized plastic bag. Each sample was divided into two subsamples: one was for chemical analyses by air drying (n = 3); and the

Experimental Design and Sample Collection
A control and three different treatments were designed as follows: (i) conventional tillage-no straw returning (NT, as CK); total maize straw removed from the field, then rotary cultivating two times by using rotary cultivator (YTO Group Corporation, China, working depth was about 10-15 cm); (ii) covering + straw returning for two years continuously (CS); the straws were used to cover the field; (iii) rotary tillage + straw returning for two years continuously (RS); the straw was crushed and extruded into small pieces with less than 5 cm, then rotary cultivating two times by using IGQN-200K-QY rotary cultivator with a working depth of 10-15 cm; (iv) covering + rotary tillage + straw returning (CRS, first year CS and second year RS, straw incorporation treatment was repeated every year). Each treatment plot was 33 m long and 20 m wide in triplicate. In the present experiment, 9600 kg ha −1 of above-ground parts of maize were incorporated into the field. Maize (Zea mays L. cultivar Xiangyu '998 ) was mono-cultured in early May at a density of 70,000 plants ha −1 , and then harvested in October. The experimental field was managed in accordance with conventional management methods.
In order to evaluate the soil condition after two consecutive tillage seasons, soil samples were gathered randomly from every treated plot from August 2018 to May 2019. The specific sampling times were the growing season (20 August 2018), freezing period (12 December 2018 and 19 March 2019), thawing period (2 April 2019), and later thawing period (3 May 2019). We collected 20 samples from four tillage methods, including conventional tillage and straw returning cropland (CK, CS, RS, and CRS). Samples were collected from five different points per plot at a depth between 5 cm to 25 cm and then formed a mixed sample, representing one of the three replicates of each treatment. Debris of each mixed sample (plant tissues, rocks, and roots) was removed by sieving through a 2-mm mesh in the field (Ren et al. 2016) and kept in a sterilized plastic bag. Each sample was divided into two subsamples: one was for chemical analyses by air drying (n = 3); and the other was stored at −80 • C, then sent to Meiji Biomedical Technology Co., Ltd. (Shanghai, China) for Illumina MiSeq sequencing analysis.

Environmental Parameters
The soil samples TN and available nitrogen (AN) concentrations were determined using the Kjeldahl method with H 2 SO 4 + H 2 O 2 digestion and the alkaline KMnO 4 method, respectively [20]. TP and available phosphorus (AP) were analyzed by using the molyb-Agriculture 2021, 11, 779 4 of 15 date colorimetric method [21]. TK and available K (AK) were measured by the flame photometer method [22]. The soil pH was determined with a pH meter at a soil: water ratio of 1:5 (DELTA 320, Mettler Toledo instruments Co., Ltd.(Shanghai, China)). The SOC concentration was measured through heating with an oil bath based on the Wlakley-Black chromic acid wet oxidation method [12].

DNA Extraction and Sequencing
The total microbial DNA was extracted by using the E.Z.N.A. ® soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA) following the manufacture's protocol from the samples. The primer information was according to reference [23], and the PCR amlification and purification method were shown in Text S1.

Bacterial Community Analysis
The raw sequences of 16S genes were processed to calculate the Amplicon Sequence Variants (ASVs) by using DADA2 (version 1.8), on the basis of the pipeline tutorial 1.8 (https://benjjneb.github.io/dada2/tutorial_1_8.html (accessed on 14 October 2018)) in R [24]. The Silva database 138 was used to align and classify the sequences of the 16S rRNA gene [25]. Archaea, chloroplast, eukaryote, and mitochondria sequences were removed from the results after classifying the sequences. Furthermore, ASVs that contained only singletons, doubletons, and tripletons were removed from the dataset of the 16S rRNA gene. Each ASV of the 16S gene was defined as each different bacterial genotype in this study. The accompanying metadata and raw sequences are deposited in the Sequence Read Archive (SRA) of the NCBI with the project accession number SRP296995.

Correlation-Based Co-Occurrence Network and Topological Features
To reveal the relationships among different bacteria in relation to environmental fluctuations, a correlation-based co-occurrence network was established on the basis of Spearman's rank correlation coefficient (ρ), which was calculated with R software. ASVs that met the following standard for every given data set was used to construct the network: (1) detected in 40% of the samples and (2) a relative proportion of >0.01% for at least one sample. Only positive correlations that had a p-value < 0.001 were selected to construct the network. The network was visualized by using Cytoscape 3.8.2 ( U.S. National Institute of General Medical Sciences (NIGMS) [26]. Network topological features (average shortest path length (L), average node degree, clustering coefficient (C), diameter, and modularity) were calculated using the "Networkanalyzer" plugin in Cytoscape and R software (package: igraph) [27,28]. A random undirected network with equal amounts of nodes and edges as the co-occurrence network was constructed with the Erdős-Rényi model using the Network Randomizer plugin in Cytoscape. Modules were identifying with the Louvain algorithm using R software (package: igraph) [27]. The small-world coefficient (σ) was worked out to investigate the small-world property of the networks (in other words, the degree of clustering and shortness of paths between nodes) [29]. All networks were visualized with Cytoscape 3.8.2 [26].

Variation Patterns of Treatment Modules
To explore the variation patterns of major modules under different treatments, we first normalized the relative abundance of ASVs in each module (including the control) using feature scaling. To calculate the transition patterns of each module, the normalized ASVs were averaged together after normalization [30]. Then, the obtained values of each treatment module (CRS, CS, and RS) were divided by those of the control (CK) to obtain the variation patterns of treatment modules.

Statistical Analysis
The samples were ordinated by using a non-metric multidimensional scaling (NMDS) analysis with Bray-Curtis distances that are based on their dissimilarity using the "metaMDS" function in Vegan [31]. The assemblage-environment relationships were defined by fitting vectors onto the ordination space using the Vegan function "envfit". The consequence of the fitted vectors was evaluated using a permutation procedure (permutation999). The significance of soil chemical properties was analyzed with SPSS 22.0 software (IBM, Armonk, NY, USA)with two-way ANOVA.

Chemical Properties of Soil
The chemical characteristics of soil samples from the control (CK) and three treatments (CS, RS, CRS) are shown in Table 1. The results of two-way ANOVA showed that control and different treatments, sampling periods, and their interaction had significant effects on most soil chemical properties (TN, AN, TP, AP, AK, pH, and SOC), except for TK ( Table 2). At the beginning of the experimental period (August 2018, the crop growing season), the concentrations of SOC and TN in CS were higher than those in CK, reaching 15.89 g kg −1 and 1.53 g kg −1 , respectively. Therefore, the SOC concentrations in RS and CRS were 91% and 81% in CK, respectively, and the TN concentrations in RS and CRS were 89% and 82% in CK, respectively. The concentration of AN was the same in CS and CK and was lower in RS and CRS than in CK. In the case of TP and AP, the concentrations were lower in CS than in CK. The pH value of CK was lower than in the straw returning treatments. The results differed before and after freeze-thawing (from December 2018 to May 2019), the off-growing season. The concentration of SOC was 12.65 g kg −1 in CS in December 2018; it increased by 3.24 g kg −1 and reached 15.89 g kg −1 in May 2019, which was higher than the changes in CK (2.84 g kg −1 ). Meanwhile, the increases of RS and CRS were lower than in CK, which were 1.22 g kg −1 and 0.29 g kg −1 , respectively. In general, the concentration increments of nitrogen (TN and AN) and phosphorus (TP and AP) in the straw treatments (CS, RS, and CRS) were smaller than in CK. The increase in TN concentration in CK was 0.33 g kg −1 , while the TN concentrations increased by 0.12 g kg −1 , 0.04 g kg −1 , and 0.02 g kg −1 in RS, CS, and CRS, respectively. The AN concentration in CS decreased by 8.81 g kg −1 , while in CK it increased by 30.84 g kg −1 . In the case of soil phosphorus concentrations, the TP and AP concentration increased in RS and CK, while they decreased in CS and CRS. The pH values of the treatments (CS, RS, and CRS) were 6.33, 6.31, and 6.38, respectively; these pH values were all higher than that of CK (6.04).

Bacterial Community Composition
At the class level, Alphaproteobacteria dominated in both the control and the straw return treatments during the experiment, while no apparent shifting patterns were observed (Figure 2A). Actinobacteria tended to increase toward the end of the experiment, especially in the control. This group comprised 35.5% of the control but 16.7-27.4% in the treatments in May. Gammaproteobacteria was approximately 1.5-3 times more enriched in the control than in the treatments. It tended to decrease to one-fifth of the proportion in May at the end of the experiment in the control. However, this group showed little variation in CS during the experiment. Thermoleophilia increased sharply until March and decreased again in both the control and the treatments. In winter, Blastocatellia made up 11.3% of the control, and this proportion was approximately twice as high as those in CS and RS. the control than in the treatments. It tended to decrease to one-fifth of the proportion in May at the end of the experiment in the control. However, this group showed little variation in CS during the experiment. Thermoleophilia increased sharply until March and decreased again in both the control and the treatments. In winter, Blastocatellia made up 11.3% of the control, and this proportion was approximately twice as high as those in CS and RS.
At the order level, Rhizobiales always dominated in both the control (CK: 7.4-12.0%) and treatments (CS: 8.3-13.7%; CRS: 7.6-10.9%; RS: 6.9-12.8%) ( Figure 2B). Xanthomonadales in the control was up to eight times higher than in the treatments in August 2018. Vicinamibacterales, Pyrinomonadales, and Rubrobacterales increased approximately 2-3 times compared to those in the control in CS and CRS, but showed no apparent variations in RS compared to the control. The proportions of Sphingomonadales, Propionibacteriales, Micromonosporales, Frankiales, and Streptomycetales decreased to 30% of those in the control in CS and CRS in May.  At the order level, Rhizobiales always dominated in both the control (CK: 7.4-12.0%) and treatments (CS: 8.3-13.7%; CRS: 7.6-10.9%; RS: 6.9-12.8%) ( Figure 2B). Xanthomonadales in the control was up to eight times higher than in the treatments in August 2018. Vicinamibacterales, Pyrinomonadales, and Rubrobacterales increased approximately 2-3 times compared to those in the control in CS and CRS, but showed no apparent variations in RS compared to the control. The proportions of Sphingomonadales, Propionibacteriales, Micromonosporales, Frankiales, and Streptomycetales decreased to 30% of those in the control in CS and CRS in May.

Multidimensional Analysis and the Potential Drivers of Bacterial Community Composition
NMDS constructed with the Bray-Curtis distance of obtained ASVs variations in bacterial community composition between samples and seasonal trends (Figure 3). Both ANOSIM and PERMANOVA also indicated that bacterial community structure showed strong seasonality (ANOSIM, R = 0.60, p = 0.001 and PERMANOVA, pseudo-F = 3.47, p = 0.01). A multidimensional analysis was performed to identify the key environmental drivers of bacterial communities (Figure 3). After fitting the environmental parameters to the NMDS ordination, it was revealed that AN, TN, and TP were significantly associated with the NMDS axes of bacterial community composition.

Multidimensional Analysis and the Potential Drivers of Bacterial Community Composition
NMDS constructed with the Bray-Curtis distance of obtained ASVs variations in bacterial community composition between samples and seasonal trends (Figure 3). Both ANOSIM and PERMANOVA also indicated that bacterial community structure showed strong seasonality (ANOSIM, R = 0.60, p = 0.001 and PERMANOVA, pseudo-F = 3.47, p = 0.01). A multidimensional analysis was performed to identify the key environmental drivers of bacterial communities (Figure 3). After fitting the environmental parameters to the NMDS ordination, it was revealed that AN, TN, and TP were significantly associated with the NMDS axes of bacterial community composition.

Co-Occurrence Networks, the Modular Structures, and the Taxonomic Features
A co-occurrence network was constructed to investigate the co-occurrence modes within bacterial ASVs related to environmental fluctuations. After filtering (ρ ≥ 0.6 and p < 0.01), we obtained a total of 4892 correlations (edges) that were statistically significant among 649 variables (nodes) ( Figure 4A,B). The topological features of the co-occurrence networks are depicted in Table S1. Compared with the Erdös-Rényi random network, the cooccurrence networks (or real network) showed higher values of the clustering coefficient (C) (0.36), the average path length (L) (4.65), and the modularity (0.49). The small-world coefficient of the co-occurrence network was 8.97, which was much higher than 1, indicating that the bacterial network showed small-world properties with modular structures.
To explore how module composition varied under different treatments at each sampling time point, we calculated the variation patterns of modules (Var-M) under different treatments against the control ( Figure 4C). Therefore, the Var-M could 1) represent the microbial transitions of a specific module in different sampling time points against the control and 2) indicate which distinctive module dominated under which treatment against the control. Compared to the control, the module I was only enriched in March

Co-Occurrence Networks, the Modular Structures, and the Taxonomic Features
A co-occurrence network was constructed to investigate the co-occurrence modes within bacterial ASVs related to environmental fluctuations. After filtering (ρ ≥ 0.6 and p < 0.01), we obtained a total of 4892 correlations (edges) that were statistically significant among 649 variables (nodes) ( Figure 4A,B). The topological features of the co-occurrence networks are depicted in Table S1. Compared with the Erdös-Rényi random network, the co-occurrence networks (or real network) showed higher values of the clustering coefficient (C) (0.36), the average path length (L) (4.65), and the modularity (0.49). The small-world coefficient of the co-occurrence network was 8.97, which was much higher than 1, indicating that the bacterial network showed small-world properties with modular structures.
To explore how module composition varied under different treatments at each sampling time point, we calculated the variation patterns of modules (Var-M) under different treatments against the control ( Figure 4C). Therefore, the Var-M could (1) represent the microbial transitions of a specific module in different sampling time points against the control and (2)     The different modules comprised different microbial compositions taxonomically ( Figure 4D). For instance, Rhizobiales dominated in modules I and II, while Gaiellales dominated in Module V. Micrococcales, Propionibacteriales, and Frankiales dominated in Module II, but Pyrinomonadales dominated in Module III. Rubrobacterales was more enriched in Module IV, but Solirubrobacterales was more enriched in Module V.

Distinctive Bacterial Associations with Environmental Parameters
The sub-network was extracted from the microbial association network to explore direct links between environmental factors and microbes ( Figure 5). The sub-network was composed of eight environmental and 68 bacterial nodes as well as 208 edges. The bacteria belong to several phyla, including Proteobacteria, Gemmatimonadota, Acidobacteriota, Verrucomicrobiota, Actinobacteriota, Chloroflexi, and others, in the subnetwork. Different colors designated the different bacterial phyla; the names of phyla Acidobacteriota, Actinobacteriota, Verrucomicrobiota were not validly published but proposed recently [32]. Notably, Rhizobiales showed correlations with multiple environmental factors (AP, TP, AK, TK, AN, and TN), while Sphingomonadales (belonged to Proteobacteria) showed correlations with the abovementioned environmental factors and SOC. Xanthomonadales (belonged to Proteobacteria) correlated with AP, AK, TN, and SOC. Interestingly, Pyrinomonadales (which belonged to Acidobacteriota) only correlated with pH but not other environmental factors.

Distinctive Bacterial Associations with Environmental Parameters
The sub-network was extracted from the microbial association network to explore direct links between environmental factors and microbes ( Figure 5). The sub-network was composed of eight environmental and 68 bacterial nodes as well as 208 edges. The bacteria belong to several phyla, including Proteobacteria, Gemmatimonadota, Acidobacteriota, Verrucomicrobiota, Actinobacteriota, Chloroflexi, and others, in the subnetwork. Different colors designated the different bacterial phyla; the names of phyla Acidobacteriota, Actinobacteriota, Verrucomicrobiota were not validly published but proposed recently [32]. Notably, Rhizobiales showed correlations with multiple environmental factors (AP, TP, AK, TK, AN, and TN), while Sphingomonadales (belonged to Proteobacteria) showed correlations with the abovementioned environmental factors and SOC. Xanthomonadales (belonged to Proteobacteria) correlated with AP, AK, TN, and SOC. Interestingly, Pyrinomonadales (which belonged to Acidobacteriota) only correlated with pH but not other environmental factors.

Discussion
Straw return could improve soil organic carbon, which could increase nutrient concentrations (including total nitrogen and total phosphorus) as the by-products of microbial metabolism, and shift the soil microbial community in agro-ecosystems. In this study, we obtained the following key results: (1) soil bacterial community structure varied significantly during the frozen and thawed periods and were mostly influenced by the concentrations of nitrogenous and phosphorus compounds; (2) network analysis revealed that the dominated bacterial modules varied in relation to the different treatments as well as in different sampling time points; and (3) the bacterial groups that were directly linked to the nutrient concentrations could contribute to the enhancement of soil fertility.

Discussion
Straw return could improve soil organic carbon, which could increase nutrient concentrations (including total nitrogen and total phosphorus) as the by-products of microbial metabolism, and shift the soil microbial community in agro-ecosystems. In this study, we obtained the following key results: (1) soil bacterial community structure varied significantly during the frozen and thawed periods and were mostly influenced by the concentrations of nitrogenous and phosphorus compounds; (2) network analysis revealed that the dominated bacterial modules varied in relation to the different treatments as well as in different sampling time points; and (3) the bacterial groups that were directly linked to the nutrient concentrations could contribute to the enhancement of soil fertility.
Microbial community variations and ecological functions in soil ecosystems could be caused by the seasonal fluctuations of environmental parameters [33,34]. During the winter, the soil microbes were reported to have more capacity to degrade organic materials, such as plant detritus, than during other seasons. Still, the utilization of these materials was lower in winter [34]. In addition, the shift of microbial community during the frozen and thawed periods could be one of the key points in the summer nitrogen cycle in the soil ecosystem [33]. In our test, the bacterial communities showed strong seasonality during the experimental period (Figures 2 and 3), agreeing with previous studies [33,34]. They also varied in order level, to some extent, in different treatments and in the control (Figure 2). Rhizobiales occupied approximately 9.5 ± 1.8% of the microbial community in this study; accordingly, they dominated both the control and the treatments during the sampling period. Among Rhizobiales ASVs, 23~41% of ASVs were assigned Bradyrhizobium, Mesorhizobium, Microvirga, and Phyllobacterium. These bacteria were reported to have the capacity to fix nitrogen, and their abundance and communities were influenced by the surrounding carbon and nitrogen sources [35,36]. Therefore, Rhizobiales, one of the dominant groups in this study, could play crucial roles in nitrogen cycling in both the control and treatment systems. Sphingomonadales, by contrast, decreased during the frozen period, while they increased again in the warmer season. This group of bacteria was one of the plant growth promoters in the plant rhizosphere [37]. Sphingomonas species could promote plant growth via producing plant hormones including abscisic acid, gibberellins, indole-3-acetic acid (IAA), salicylic acid (SA), and zeatin [38][39][40]. Therefore, the increase of Sphigomonadales during the warmer season could be due to their interactions with the rhizosphere of plants and further promote plant growth.
Many previous studies have revealed the roles of N and P sources in soil microbial community structure and their enzyme activities [6][7][8][9][10]. Consistently, the multidimensional analysis showed that these two factors were the major drivers of the bacterial community structure in this study. The straw itself contained large amounts of biodegradable organic materials, including cellulose, lignin, and pentosan. These organic materials could be good carbon, nitrogen, and phosphorus sources for specific groups of bacteria. After the biodegradation of the added straw in the treatments, the concentrations of nitrogen and phosphorus increased, which affected the local microbial community structure further. In addition, the concentrations of N and P were higher in the summer and lowered in the winter. The high concentrations of N and P could have been partly caused by fertilizer application in early summer for crop cultivation. The high concentrations of N and P could be utilized and transformed into easy-to-use materials by specific groups of bacteria and further utilized by the crops to increase crop yields.
Consistent with previous reports in soil ecosystems [41][42][43][44], the co-occurrence patterns of bacterial communities in this study displayed modular structure with small-world properties, indicating that the interspecies interactions among bacteria were nonrandom but highly organized in accordance with their demands. Bacterial communities were divided into six major modules; specific modules dominated in different treatment systems as well as in the frozen and thawed periods compared with the control. For instance, Module V was more enriched in December in all three treatments, but Module VI was more enriched in treatment CS. Solirubrobacterales and Gaiellales, which belong to Actinobacteria, were two dominant groups in Module V, and Rhizobiales and Solirubrobacterales were dominated in Module VI. Solirubrobacterales and Gaiellales have been mentioned to have the capacity to degrade lignin [45]. A recent study pointed out that mainly fungi are involved in decomposing refractory substances in summer, while bacteria replace them in this role in winter in forest soils [46]. Since straw contains a large amount of lignin, which is also hard to degrade, we assumed that these bacteria could contribute to lignin degradation during the frozen season. Module II was more enriched in the warmer season than in the cold season. One of the major bacterial components of this module was Pyrinomonadales. The reason was unknown, but members belonging to this group seemed to be mesophiles or thermophiles [47,48], which could explain why they were more enriched in the warmer season. Module IV was highly enriched in the treatments with comparing the control in May and August during the crop growing period. Among many bacteria, members assigned to Actinobacteria, such as Rubrobacterales, dominated in this period. Past studies have shown that Actinobacteria can contribute to degrading organic matters, including cellulose and phytoplankton detritus [49,50]. The organic/inorganic materials released as their by-products could be used as a good energy source for crop growth. The treatment soils contained high amounts of straw that could be used as substrates and degraded by Actinobacteria, which explained the higher abundance of Actinobacteria in the spring and summer in the treatments than in the control. To sum up, the separation of the module formation could be derived from the differences in bacterial physiological and ecological characteristics in the ecosystem, which further forced bacterial species to establish interconnectivity with specific groups of bacteria and dominated in certain crop growing periods and/or frozen and thawed periods.
Environmental factors were directly linked with certain groups of bacteria in this study ( Figure 5). ASVs assigned to Pyrinomonadales were only correlated with pH, indicating that the abundance of these mesophiles or thermophiles [47,48] was highly susceptible to pH variations, although the reasons remain unclear. Rhizobiales and Sphingomonadales correlated with most of the environmental factors, including the concentration of nitrogen sources. These two groups of bacteria were reported to contribute to the nitrogen cycle and acted as plant growth promoters, discussed in detail above. Because they could convert nitrogen (gas) to other forms of nitrogenous compounds and further influence the local N concentrations, we assumed that these two groups could play key roles in ecosystem soil fertility. Xanthomonadales correlated with N, P, and K sources as well as SOC. Huang et al. (2019) reported that Xanthomonadales increased with the concentrations of bio-organic fertilizer in cucumber farmland. Their biodegradation of bio-organic fertilizers could increase N, P, and K concentrations in the soil. Therefore, the direct links of Xanthomonadales with the abovementioned environmental factors suggested their key roles in soil fertility.

Conclusions
In our study, the soil bacterial communities varied significantly during the experimental period, and they also varied, to some extent, in different treatments. Among the environmental factors, AN, TN, and TP were the key drivers of bacterial community structure. The co-occurrence network constructed with bacterial and environmental data could be divided into several modules and showed nonrandom but small-world properties. Notably, the dominated bacterial modules varied in different months and different treatments, suggesting that specific microbial groups could contribute to soil fertility in relation to environmental fluctuations. Rhizobiales, Pyrinomonadales, Sphingomonadales, and Xanthomonadales showed correlations to specific nutrient concentrations and pH, indicated that they could play important roles in the changes in soil fertility. In this experiment, we elucidated the short-term responses of microbes to different straw returning procedures. As the relationships between soil fertility, microorganisms, and environmental factors remain unclear, a long-term study on these relationships will be conducted in the future.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/agriculture11080779/s1, Text S1: The information of primers and the method of PCR amplification and purification. Figure S1: Monthly temperature and precipitation of the test area during the sampling period. Table S1: Topological features and statistics of the microbial association networks.