Differential Assembly and Shifts of the Rhizosphere Bacterial Community by a Dual Transgenic Glyphosate-Tolerant Soybean Line with and without Glyphosate Application

Modern agriculture has gained significant economic benefits worldwide with the use of genetically modified (GM) technologies. While GM crops provide convenience to humans, their biosafety has attracted increasing concern. In this study, the Illumina MiSeq was used to perform a high-throughput sequencing of the V3-V4 hypervariable regions of 16S rRNA gene (16S rDNA) amplicons to compare the rhizosphere bacterial communities of the EPSPS/GAT dual transgenic glyphosate-tolerant soybean line Z106, its recipient variety ZH10, and Z106 with glyphosate application (Z106G) during flowering, seed filling, and maturing stages under field settings. At each of the three stages, the alpha and beta diversity of rhizosphere bacterial communities revealed no significant differences between ZH10, Z106, and Z106G. However, some bacterial taxa demonstrated a greater proportional contribution, particularly the nitrogen-fixing rhizobium Ensifer fredii, in the rhizospheric soil of Z106 at the seed filling and maturing stages, when compared to ZH10 and Z106G. The present study therefore suggests that the EPSPS/GAT dual transgenic line Z106 and exogenous glyphosate application have a minimal effect on the composition of the soybean rhizosphere bacterial community but have no impact on the structure of the rhizosphere microbial community during a single planting season.


Introduction
Following the turn of the 21st century, biotech/genetically modified (GM) crop technology was widely adopted in a number of developed and developing countries, primarily in the cultivation of canola, maize, cotton, and soybean [1,2]. China, a leading country in planting GM crops, planted 2.94 million hectares of GM crops in 2018, with 2.93 million hectares of insect-resistant cotton and 9600 hectares of virus-resistant papaya [3,4]. In China, the commercialization of GM soybean is still being considered, owing to the potential for biosafety concerns [5], such as the influence of glyphosate-tolerant transgenic soybean with herbicide application on the soil ecological environment and microorganisms [6]. Glyphosate-based herbicides are the most widely used herbicides in the world, and they are widely employed on glyphosate-tolerant GM crops [7,8]. Glyphosate has a chelating property that allows it to be rapidly absorbed and translocated in plant tissues via a variety of action sites, resulting in the lack of metabolism in plants and changes in functions or enzymes in carbon metabolism and nitrogen metabolism, and it is toxic to soil organisms other than plants. It is worth considering the ecological risk assessment of glyphosate-tolerant GM crops and glyphosate application in the field [7,[9][10][11][12].
To improve the methodological approach in this study, we used high-throughput sequencing of 16S rDNA (V3-V4 hypervariable region) amplicons via the Illumina MiSeq, a technology that has been widely used for the analysis of the structure and function of the rhizosphere bacterial community [17][18][19], and analyzed the rhizosphere bacterial communities between Z106, its recipient Zhonghuang 10 (ZH10), and Z106 with glyphosate application (Z106G) at the flowering, seed filling, and maturing stages.

Plant Materials, Field Design, and Sampling Methods
In this study, the soybean (Glycine max L. Merr.) line Z106, which carries the EPSPS and GAT glyphosate-tolerant transgenes, its recipient soybean line ZH10, and Z106 treated with 900 g·hectare −1 glyphosate at the seedling stage on July 24 in 2017 (Z106G) were selected as three plant material treatments [20]. The soybeans were planted in the Shunyi Experimental Field Station of the Crop Research Institute (Chinese Academy of Agricultural Sciences) in Shunyi district, Beijing, China (N 40.23 • , E 116.56 • ) in July 2017. The soil type used is aquic cinnamon soil according to the China Soil Database (http://vdb3.soil.csdb.cn/ (accessed on 15 December 2017)). The experimental field was divided into 9 plots (15 m × 10 m per plot) using a randomized block design, with three replications for each of the three treatments. Without any other fertilizers, [NH 4 ] 2 HPO 4 300 kg·hectare −1 and KCl 150 kg·hectare −1 were applied prior to planting [13].
We collected the experimental samples at the flowering (R1-R2) stage on August 16, seed filling (R5-R6) stage on September 11, and maturing (R8) stage on October 23 in 2017, and soil before planting was collected on June 30 in 2017. At three growth stages, the surrounding soil or rhizospheric soil of each soybean treatment was sampled from three different plots with five sampling points each, in which were distributed diagonal crossover, and composite samples from five sampling points per plot were made and treated as one replicate. F, S, and M represent the soil sampled at the flowering, seed filling, and maturing stages, respectively, while BS and RS represent surrounding soil that was loosely adhered to the soybean's roots and rhizospheric soil that was collected by brushing off soil that was tightly adhered to the root surface, respectively. The BS without marked growth stages represents soil before planting (Supplementary Table S1). All 63 soil samples were selected for 21 different groups, with 3 replications in each group. For basic physicochemical analysis, parts of the soil samples were dried in the air, while the remainder were stored at −80 • C for DNA extraction [21,22]. The entire plant at the three stages, including pods, was sampled and dried in the air, after which the sample was ground into a powder by a disintegrator for plant analysis.

Basic Physicochemical Properties of Soils and Plants: Key Soil Enzyme Activities
The carbon and nitrogen contents of the dried rhizospheric soil and dried plant tissue samples at the three stages were analyzed using the Dumas Combustion Method by Vario MICRO cube elemental analyzer (Elementar Analysensysteme GmbH, Langenselbold, Germany) in the Modern Analysis Center of Nanjing University. The pH value and water content of fresh rhizospheric soil samples were tested under laboratory conditions by the potentiometric method and stoving method [13], respectively. Furthermore, the activities of key soil enzymes involved in the nitrogen cycle and carbon cycle, such as Solid-Urease (S-UE), Solid-Nitrate Reductase (S-NR), Solid-Nitrite Reductase (S-NiR), and Solid-Sucrase (S-SC), were tested for dried surrounding soil and dried rhizospheric soil using the corresponding kit developed in 2017 (Beijing Solarbio Science & Technology Co., Ltd., Beijing, China).

Culturable Nitrogen-Fixing Bacteria in Rhizospheric Soil and Nodulation Effect of Soil Nitrogen-Fixing Bacteria
To count the number of culturable nitrogen-fixing bacteria in soil, 1.0 g of fresh rhizospheric soil was diluted to the concentrations of 10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 , 10 −6 , 10 −7 , and 10 −8 in sterile water. To determine the appropriate bacterial concentration for counting, 0.1 mL different concentrations of bacterial solution were plated evenly on nitrogen-fixing rhizobium medium (Qingdao Hope Biol-Technology Co., Ltd., Qingdao, China). After 48 h incubation at 28 • C, the optimal concentration of nitrogen-fixing bacteria colony distribution was observed, then the 10 −5 diluent was selected as the most suitable concentration. Finally, colonies were enumerated and the number of culturable nitrogenfixing bacteria per gram of soil was calculated as follows: Number of culturable nitrogen-fixing bacteria per gram of soil = Colonies number × Dilution concentration/Weight of fresh soil Soybean roots and nodules were collected and stored in sealed bags at the flowering, seed filling, and maturing stages to study the nodulation effect of soil nitrogen-fixing bacteria. The nodules were carefully detached from the roots, rinsed, and dried in the air, and the number and fresh weight of nodules were then recorded. Each sampling plot had 20 individual soybean plants. The number and fresh weight of root nodules per plant were calculated as follows: Number In this study, each biological replicate of the total metagenomic DNA was extracted from approximately 0.6 g soil by using the PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA), following the methods as Lu et al. described with slight modifications [23,24]. The extracted DNA solution was stored at −80 • C for subsequent experiments.
We used the forward primer 341F (5 -ACTCCTACGGGAGGCAGCAG-3 ) and the reverse primer 806R (5 -GGACTACHVGGGTWTCTAAT-3 ) to amplify the V3-V4 region of the 16S rDNA by Polymerase Chain Reaction (PCR) [25,26]. The qualified libraries were subjected to PCR amplification, product purification, library quality determinations, quantification, and high-throughput sequencing of the qualified libraries on the Illumina MiSeq platform (Illumina, San Diego, CA, USA) using the MiSeq Reagent Kit by BGI Tech Solutions (Wuhan, China) as described previously by Lu et al. (2017). The sequencing clean data of 63 samples were submitted to the Sequence Read Archive (SRA) with the accession number PRJNA734728.  Table S2). The sequencing data statistics were provided by BGI Tech Solutions.

OTU Analysis of 16S rDNA Amplicon Sequencing Data
The clean tags were clustered into operational taxonomic units (OTUs) by using UP-ARSE (v7.0.1090) [27], and minor modification was conducted as previously described by Lu et al. [21,22]. OTU counts were normalized after species annotation and phylogenetic tree construction, according to the minimum sample number sequence. Mitochondrion/ chloroplast-related sequences (phylum Cyanobacteria, family Mitochondria) were eliminated [28,29], to avoid the interference of mitochondrial and chloroplast genes. The OTU rank abundance curve was drawn to reflect the richness and evenness of species (Supplementary Figure S1). The Pan and Core OTU analysis was used to observe the increase in total species and decrease in common shared species, respectively (Supplementary Figure S2). All of the operations described above were performed on the I-Sanger platform (http: //www.i-sanger.com (accessed on 30 March 2021)) [26].

Alpha Diversity, Beta Diversity, Taxonomy, Functional Analysis, and Statistical Analysis
Alpha and beta diversity analyses were conducted based on OTUs and species annotation results between different treatments at each of the three stages [30].
Alpha diversity analysis calculated by Mothur (v1.30.2) was applied to express the complexity of species diversity for a sample by different indices, such as sobs, Chao1, ACE, Shannon, and Simpson, then Graphpad Prism 8.4.2 was used to make the boxplot of each index.
Taxa clustering was performed based on the relative abundance of each taxon. The bar plot drawn by R software (v3.3.1) showed the dominant species at different taxonomic levels in each sample and relative abundances of these species. A ternary plot was drawn by GGTERN (v3.3.0) to display the ratio relationship of the different attributes of three variables [28,33]. The Circos plot was drawn by Circos-0.67-7 to describe the correspondence between samples and species [34].
The Cluster of Orthologous Groups of proteins (COG) function classification was performed using the software Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) (v2.2.0), which was based on the optimization and improvement of PICRUSt1 [35], according to the database EggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups, http://eggnog.embl.de/ (accessed on 30 March 2021)) and KEGG (Kyoto Encyclopedia of Genes and Genomes).
For statistical analysis, the Student's t-test for different indices of alpha diversity was conducted. Analysis of similarities (ANOSIM) and permutational MANOVA (PERMANOVA/ Adonis) were performed using the vegan package of R software (v3.3.1) based on the Bray-Curtis, WUF, and UUF distance metrics with the database Silva132/16S_bacteria [26]. We processed the above data on the I-Sanger platform. One-way ANOVA performed by Graphpad Prism 8.4.2 was used to compare basic physicochemical properties, key soil enzyme activities, culturable nitrogen-fixing bacteria, the nodulation effect, and relative abundances of rhizosphere bacterial communities and COG functional classification of ZH10, Z106, and Z106G. The Tukey post hoc test was used after one-way ANOVA.

Basic Physicochemical Properties of Rhizospheric Soil and Plants
At the flowering stage, a significant difference in rhizospheric soil pH was observed between ZH10 and Z106, but no significant difference in water content was identified between the three soybean treatments at any of the three stages, as displayed in Table 1. The rhizospheric soil carbon content differed significantly between the three treatments at the flowering stage, while the rhizospheric soil nitrogen content differed substantially at the flowering and maturing stages. The plant carbon and nitrogen content did not differ significantly across the three treatments at any of the three stages.

Key Enzyme Activities in Surrounding Soil and Rhizospheric Soil
As shown in Table 2, at the seed filling stage, S-UE in rhizospheric soil of Z106 was significantly lower than that of ZH10; however, S-NR in rhizospheric soil of Z106 was significantly higher than that of ZH10 and Z106G. At the maturing stage, the three treatments were significantly different from each other in terms of S-SC in rhizospheric soil. In the surrounding soil, ZH10, S-UE, and S-NiR were at their highest during the seed filling stage, and S-NR and S-SC were at their highest during the maturing stage. For Z106, S-UE was highest at the flowering stage, while S-NR, S-NiR, and S-SC were highest during the maturing stage. For Z106G, S-UE reached the maximum at the seed filling stage, while S-NR, S-NiR, and S-SC reached their maximum at the maturing stage. In rhizospheric soil, for ZH10, S-UE and S-NiR reached their maximum at the seed filling stage, while S-NR and S-SC reached their maximum at the maturing stage. For Z106, S-UE peaked during the flowering stage, while S-NR, S-NiR, and S-SC peaked during the maturing stage. For Z106G, S-UE and S-NiR peaked during the seed filling stage, while S-NR and S-SC peaked during the maturing stage.

Culturable Nitrogen-Fixing Bacteria and Nodulation Effect in Rhizospheric Soil
In Supplementary Figure S3A, the amount of culturable nitrogen-fixing bacteria in bulk soil of ZH10, Z106, and Z106G at pre-planting, or rhizospheric soil during flowering, seed filling, and maturing stages, did not differ significantly. During the soybean growth stages, the number of culturable nitrogen-fixing bacteria of the three soybean treatments peaked at the seed filling stage, and subsequently declined at the maturing stage.
As shown in Supplementary Figure S3B,C, the number and total weight of nodules per plant showed no significant differences between ZH10, Z106, and Z106G at the flowering, seed filling, and maturing stages. With the abscission and withering of nodules, the number and weight of nodules increased to their peak value during the seed filling stage, and then dropped to their lowest value during the maturation stage.

Alpha Diversity of Rhizosphere Bacterial Communities
For alpha diversity, sobs, Chao1, and ACE reflect community richness, Shannon and Simpson indices reflect community diversity, while coverage reflects community coverage ( Figure 1). The Student's t-test results reveal that, except for Simpson index in the BS samples of ZH10 and Z106G, Z106, and Z106G at the seed filling stage, no significant difference in alpha diversity was observed for different groups of the BS or RS samples between the three treatments at each of the three stages (Supplementary Table S3). Furthermore, there were some significant differences between the BS and RS samples of each treatment at each of the three stages, as well as the three stages for each sampling compartment of each treatment (Supplementary Table S3).

Beta Diversity of Rhizosphere Bacterial Communities
To investigate the similarity or differences in community structure between various samples, the beta diversity analysis was performed using PCoA based on the WUF distance matrix (Figure 2), PCA based on OTU (Supplementary Figure S4), and PLS-DA (Supplementary Figure S5). Bacterial communities in the BS or RS samples of ZH10, Z106, and Z106G were clustered into one group without obvious distinction at the flowering, seed filling, or maturing stage, indicating no marked differences between the three treatments at each of the three stages. The statistical ANOSIM and Adonis based on the Bray-Curtis, WUF, and UUF distance metrics were then performed. The results show that the bacterial communities of ZH10, Z106, and Z106G in the same compartment at a single stage were not significantly different (Supplementary Table S4).

Taxonomic Analysis of Rhizosphere Bacterial Communities
During the taxonomic analysis, more than 290 classified soil bacteria species were discovered. At a taxonomic level, the barplot results show the dominant taxa in each sample as well as the relative abundance of each dominant species in all samples (Supplementary Figure S6). Then, using community heatmap analysis, we examined the relative abundances of major bacterial taxa (Supplementary Figure S7). Afterwards, the relative abundances of some major bacterial taxa were analyzed at the phylum, class, order, and family levels in the different compartments of the three treatments at each of the three stages by using one-way ANOVA (Supplementary Figure S8). The results show that the relative abundance of most taxa in BS or RS samples showed no significant differences between ZH10 and Z106, or Z106 and Z106G at the flowering, seed filling, or maturing stage. Furthermore, the comparison of the relative abundance of main nitrogen-fixing bacterial genera revealed that no genus was substantially different between ZH10 and Z106, or Z106 and Z106G at any of the three stages ( Figure 3). Since the major bacteria shown in the community heatmap at the species level were almost unclassified, the ternary plot analysis was then used to compare the species composition of three groups (Figure 4). The obtained results show that some bacterial taxa such as Ensifer (Sinorhizobium) fredii of the family Rhizobiaceae had a distinct proportional contribution in RS samples between ZH10, Z106, and Z106G. E. fredii was less abundant in ZH10 and Z106G than in Z106 at the seed filling and maturing stage, and its proportional contribution in ZH10 was clearly lower than other two treatments. At the seedling stage, the proportional contribution of Sphingobacterium multivorum, a member of the Sphingobacteriaceae family, was higher in Z106 than in ZH10, but lower in Z106G. In addition, the proportional contribution of unclassified Bradyrhizobium (genus level) at the seed filling stage, and unclassified Bradyrhizobium (genus level), and Micrococcacese (family level) at the maturing stage differed amongst the three treatments. The Circos plot was then produced to see the proportionate distribution of the nitrogen-fixing bacteria E. fredii in all 63 samples (Supplementary Figure S9), and the results are consistent with the ternary plot.

Functional Composition of Microbial Community in Soil Samples
To predict the functional composition of the microbial community in the samples, the software PICRUSt2 was utilized. As shown in Supplementary Figure S10, the functional OTUs related to amino acid transport and metabolism, translation, ribosomal structure and biogenesis, energy production and conversion, and cell wall/membrane/envelope biogenesis were expressed in higher abundances (>5,000,000) in soil samples. We focused on the relative abundance of COG functional classification related to nitrogen fixation (Supplementary Figure S11), which included COG1348, COG1433, COG2710, COG4656, COG5420, COG5456, and COG5554. The results reveal that there were no significant changes between ZH10 and Z106, or between Z106 and Z106G in the same compartment at a single stage for these COGs.

Discussion
The overall framework of the current study was derived from earlier research that investigated the impact of EPSPS/GAT dual transgenic glyphosate-tolerant soybean line with or without glyphosate application on root-associated microbial communities [13]. Unlike the previous study in 2015 [13], we sampled the surrounding soil and rhizospheric soil of a transgenic soybean line and its recipient cultivar during the flowering, seed filling, and maturing stages, and the V3-V4 hypervariable regions of 16S rDNA amplicons were sequenced, which may allow us to observe dominant species in soil samples more accurately. In this study, we looked at the differences between ZH10 and Z106, or between Z106 and Z106G, to assess the potential impact of planting the EPSPS/GAT transgenic soybean line and using a glyphosate-based herbicide on soil ecology and biosafety in modern sustainable agriculture. Furthermore, since the recipient line ZH10 was unable to tolerate glyphosate, we did not treat the recipient soybean with glyphosate.
For the analysis of the basic physicochemical properties of soils and plants, the pH value, carbon content, and nitrogen content in rhizospheric soil showed some significant differences between the three treatments, as well as some key enzymes of the nitrogen cycle and carbon cycle in rhizospheric soil at the three stages. The results reveal that the planting of a dual-gene transgenic soybean line and glyphosate application may have an impact on the rhizospheric soil. However, as shown in Supplementary Figure S3, neither the planting of dual-gene transgenic soybean lines nor the application of glyphosate had an effect on the number of culturable nitrogen-fixing bacteria or the nodulation effect in rhizospheric soil between the three soybean treatments at each stage. This indicates that the above-mentioned considerable variation was insufficient to influence soybean's ability to fix nitrogen. Additionally, the nitrogen content of each treatment's plant did not alter during the crop stages, which is generally variable during soybean growth [36]. This outcome could be due to a variety of external factors or the characteristics of particular soybean lines, which will be investigated further in our future study.
As proved by Yang et al. and Han et al., Bradyrhizobium and Ensifer are the two main genera of soybean microsymbionts in China [37,38]. In our study, the relative abundances of these two genera were not significantly different between the three treatments at each of the three stages, which was consistent with no differences in nodulation effect in the rhizospheric soil in this study. These results indicate that neither planting the EPSPS/GAT transgenic soybean line nor using glyphosate-based herbicide significantly changed the structure of rhizosphere microbial communities, which is in accordance with earlier research [13,39].
According to the ternary plot analysis, the proportional contribution of several soil bacterial taxa was variable. As shown in Figure 4, the proportional contribution of E. fredii was higher in RS samples of Z106 than that of ZH10 and Z106G at the seed filling and maturing stages. E. fredii is reported as a common fast-growing agricultural rhizobia, originally isolated from Chinese soil, with a diverse host range and biological nitrogenfixing capacity [31,40]. In alkaline soils, E. fredii can form nitrogen-fixing root nodules in 79 different legume genera, including soybean [41,42]. At present, the nitrogen fixation capabilities of this rhizobium are gradually being studied and improved, and will bring great benefits to agriculture by reducing the reliance on fertilizer applications [43].
According to the Circos plot results (Supplementary Figure S9), the proportional contribution of E. fredii in the rhizospheric soil increased first and then decreased, peaking at the seed filling stage, which is consistent with the results of culturable nitrogen-fixing bacteria and the nodulation effect in the rhizospheric soil. The proportional contribution of E. fredii in rhizospheric soil in the three soybean treatments at the three growth stages was highest in Z106, followed by Z106G, and was lowest in ZH10. Such results were not observed in our prior study of the 16S V4 region sequencing with three sampling compartments: surrounding soil, rhizospheric soil, and root [13], most likely due to deeper amplicon sequencing and a reduction in sampling compartments. This phenomenon implies that the enrichment of this rhizobia was aided by planting the EPSPS/GAT dual transgenic line and hindered by glyphosate application in rhizospheric soil during the soybean growth stages, and that these changes were regulated by gene expression, root exudation, and the potential biological toxicity of glyphosate herbicide [44,45]. Interestingly, the proportional contribution of E. fredii in Z106G was lower than Z106, but still higher than that of ZH10; the reason for this needs to be investigated further and confirmed by more advanced sequencing technologies, such as full-length 16S rRNA amplicon or metagenome sequencing, which can avoid relic DNA interference [46].
Previous research has revealed that rhizosphere microbes can enrich specific functional genes [47]. We compared the COG function classification and its relative abundance using PICRUSt2. As shown in Supplementary Figure S11, the seven COGs were found to be directly related to nitrogen fixation, including COG1348, which is related to the key enzymatic reactions in nitrogen fixation catalyzed by the nitrogenase complex, which has two components, the iron protein and the molybdenum-iron protein (by similarity), COG1433 related to dinitrogenase iron-molybdenum cofactor biosynthesis protein, COG2710 related to nitrogenase, COG4656 related to the requirement for nitrogen fixation and that may be part of a membrane complex that functions as an intermediate in the electron transport to nitrogenase (by similarity), COG5420 related to nitrogen fixation, COG5456 related to nitrogen fixation protein FixH, and COG5554 related to nitrogen fixation protein. The relative abundances of these COGs showed that either planting the EPSPS/GAT transgenic soybean line or employing glyphosate had no impact on these functional genes.
In conclusion, our results demonstrate that a single planting season of an EPSPS/GAT dual transgenic glyphosate-tolerant soybean line with or without glyphosate application had no significant impact on bacterial communities and the nodulation effect in the rhizosphere. However, the physicochemical properties of rhizospheric soil and the proportional contribution of the nitrogen-fixing bacteria E. fredii responded differently to planting the EPSPS/GAT dual-gene transgenic line and glyphosate application. In general, the current study will aid future researchers to develop field studies to evaluate the biosafety of transgenic soybeans and herbicide application.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.3 390/horticulturae7100374/s1, Figure S1: Rank-abundance curves of all 63 samples, Figure S2: Pan-Core analysis of all samples, Figure S3: A, culturable nitrogen-fixing bacteria among the rhizospheric soil of ZH10, Z106 and Z106G at pre-planting, flowering, seed filling and maturing stages; B & C, number and weight of nodules at flowering, seed filling and maturing stages, Figure S4: PCA based on OTU abundance of bacterial communities, Figure S5: PLS-DA on OTU level at flowering (A), seedfilling (B) and maturing (C) stages, Figure S6: Community bar plot analysis of all samples, Figure S7: Community heatmap at Phylum (A), Class (B), Order (C), Family (D), Genus (E) and Species (F) levels, Figure S8: The relative abundance of the major bacterial taxa in surrounding soil and rhizospheric soil at Phylum (A,B), Class (C,D), Order (E,F), or Family (G,H) levels, Figure S9: Circos plot for Ensifer fredii at Species level, Figure S10: COG functional prediction, Figure S11: Statistics of COG function classification related to nitrogen fixation, Table S1: Names of all 63 soil samples, Table S2: Summary of sequencing data of all samples, Table S3: p-values of six indices for alpha diversity of the different compartment samples between ZH10, Z106, and Z106G according to Student's t-test, Table S4: ANOSIM and Adonis of the different compartment samples between ZH10, Z106, and Z106G based on the Bray-Curtis, Weighted Unifrac and Unweighted Unifrac distance metrics.