The influence of environmental factors on communities of arbuscular mycorrhizal fungi associated with Chenopodium ambrosioides revealed by MiSeq sequencing investigation

Arbuscular mycorrhizal fungi (AMF) affect multiple ecosystem functions and processes, the assemblages of which vary across ecosystems. However, the influences of environmental factors on AMF communities which may shape these communities are still largely unknown. In this study, AMF communities from roots and rhizosphere soils of Chenopodium ambrosioides in different natural soils were investigated. The root habitat showed significantly smaller numbers of OTUs and lower community richness compared to the rhizosphere soil habitat. Most OTUs in the root habitat were shared by the soil habitat from the same sampling site, indicating that rhizosphere soils represent a pool of AMF species, a fraction of which is recruited by plants. Most of the AMF in root habitats were Glomeraceae, suggesting recruitment preferences of AMF by plants. The relative contributions of environmental factors to explain variations in AMF community composition and phylogenetic structure were assessed. The results revealed soil properties predominantly explained the variation, followed by geographic and climate parameters which explained a small fraction independently, while the host plant showed few explanations. Overall, our results indicated that soil and root habitats as well as soil characters, especially pH, nitrogen and micronutrients (Zn and Cu) affected AMF communities significantly.

Overall sequencing information and taxonomic richness. A total of 1,465,650 paired-end sequences were obtained from the 36 libraries using the AMV4.5NF/AMDGR primer set. These sequences were overlapped to obtain high-quality tag sequences, and the dominant length distribution was 200− 250 bp (> 93%), with an average length of 218 bp. The details of the tag sequences obtained from each of the 36 samples are provided in Supplementary Table S2.
The taxonomic distributions of the obtained sequences in rhizosphere soil and root samples are shown in Supplementary Fig. S1. Fungal sequences other than Glomeromycota were detected. However, the majority of fungal sequences amplified with AMV4.5NF/AMDGR belonged to Glomeromycota (705,323 sequences, corresponding to 60.1% of the total). The second and third largest groups were Basidiomycota (15.2%) and Chytridiomycota (15.1%), respectively. The lowest numbers of Glomeromycota sequences were obtained for root and soil samples of XC (R-XC and S-XC), while the percentage of Glomeromycota sequences obtained from other samples was > 40%. The Glomeromycota sequences were extracted for further analyses. Specifically, pyrosequencing of XC3 and WH3 failed because few Glomeromycota sequences were obtained, and so these two samples were excluded from further analyses. AMF community richness and diversity. All of the rarefaction curves tended to reach saturation ( Supplementary Fig. S2), revealing that the data volume of sequenced reads was sufficient to detect the majority of sequence types. Marked variation in total OTU number among the samples was suggested by the rarefaction curve. The samples from rhizosphere soil habitats had a larger number of OTUs (20.5-55.5) compared to those from root habitats (4.0-19.7) ( Table 2). Consistently, the richness indices, Chao1 and ACE, showed that the samples from root habitats had the lowest number of AMF OTUs, while the samples from rhizosphere soil habitats had the highest number, with significant differences between habitats. However, no significant differences of Shannon diversity were detected between the root and rhizosphere soil habitats.
AMF community composition. Significant variation in AMF community composition at the genus level was detected among both root and soil samples (Fig. 1). Sequences that could be classified were affiliated with nine AMF genera, while those that could not be classified were assigned as others. Glomus and Rhizophagus were the most abundant genera in all root samples, but their relative levels differed. A greater number of genera were present in samples from rhizosphere soil habitats; the most abundant genera were Glomus, Funneliformis, and Claroideoglomus.
A PCoA based on the OTU composition is shown in Fig. 2a. Of the variation in the AMF communities, 26.6% and 13.3% could be explained by the first and second principal components, respectively. The samples from root and soil habitats at the same site clustered together (Fig. 2a). This was confirmed by the hierarchical cluster analysis, which showed clusters of samples from the same site (Fig. 2b).
Indicator species analyses were used to identify specific OTUs associated with the rhizosphere soil or root habitat. Sixteen OTUs were found in rhizosphere soil samples. Among them, 11 OTUs were specific for the rhizosphere soil habitat (Table 3). However, no OTUs associated with the root habitat were detected.  In total, 210 OTUs were found in the AMF communities based on 97% species identity (Fig. 3). The phylogenetic placement of OTUs in different habitats was determined using the phylogenetic tree (Fig. 3). In general, the OTUs in rhizosphere soils were evenly scattered in the phylogenetic tree. However, most of AMF detected in root habitats were Glomeraceae while few AMF belonging to Acaulosporacea, Diversisporaceae or Gigasporaceae were found. Claroideoglomeraceae, Paraglomeraceae were identified in the roots of DN, GP and TL. Although the assemblages of AMF communities were markedly different, most of the OTUs in root samples were shared by the corresponding soil samples from the same site (Fig. 4).
Relationship between AMF community structure and environmental factors. Relationships of AMF richness (indicated by the ACE index) and phylogenetic diversity (revealed by the Faith's index) with plant, soil and climate variables were analyzed (Supplementary Table S3). We found that the soil NH4+ -N level correlated negatively with the AMF richness and phylogenetic diversity (Kendall test, P < 0.05). RDA, the multivariate analysis based on constrained ordination, was applied to analyze the influence of soil properties on the   Table 3. OTUs associated with the rhizosphere soil habitat. A is the probability that the surveyed site belongs to a given environment, given the fact that the species has been found; B is the probability of finding the species in sites belonging to a given environment. Asterisks refer to OTUs specific for rhizosphere soil habitats.
AMF community composition in C. ambrosioides rhizospheres (Fig. 5). The first two RDA components explained 31.7% of the total variation. The results showed that total N (P = 0.007), pH (P = 0.029), Zn (P = 0.033), and Cu (P = 0.047) explained the AMF community in soil habitats most and had a significant correlation between each variable and the ordination scores ( Fig. 5).
To further analyze the relative importance of soil characters, plant parameters, geographic distance, and climate properties in predicting AMF community composition and phylogenetic structure, a stepwise distance-based redundancy analysis (db-RDA) was carried out. In contrast to traditional RDA, an advantage of db-RDA is that it enables the inclusion of environmental factors and the testing of their interaction using non-Euclidean distance matrices. Totally, 71.6% and 60.8% of the variations of the soil AMF community composition and phylogenetic structure were explained by the whole set of the selected variables respectively (Fig. 6). A very large fraction of variations of AMF community composition (62.8%) and phylogenetic structure (41.0%) could be explained by soil properties, followed by geographic and climate variables (31.5% for community composition and 16.5% for phylogenetic structure respectively). The plant variables explained a relative small proportion of the variations in both community composition and phylogenetic structure.
Permutational multivariate analysis of variance (PerMANOVA) was also carried out to examine the relative importance of each single environmental factor to the AMF community composition and phylogenetic structure (Supplementary Table S4). pH showed weak associations with the AMF community composition and phylogenetic structure, while the soil available Cu content revealed significant associations with them. Besides, the PCNM vector (principal components of neighbor matrices, representing the spatial factors) and MAP (mean annual precipitation) showed significant associations to the AMF community composition. For the AMF community phylogenetic structure, soil total Cu, Zn, and Cd contents as well as MAT (mean annual temperature) were significantly correlated.

Discussion
It is commonly accepted that more than 80% of terrestrial plants are colonized by AMF, with which they form associations 29,30 . The extraradical hyphal networks produced by AMF can alter plant physiology, enhance plant nutrient absorption and translocation, and increase plant resistance to HMs [31][32][33] . Root habitats showed a significantly smaller number of AMF OTUs compared to rhizosphere soil habitats, which was confirmed by the lower    Chao1 and ACE richness indices of root samples. Similar results have previously been reported 18,34 . Besides, the PCoA based on the OTU composition showed that the samples from root and soil habitats at the same site clustered together, which was consistent with the hierarchical cluster analysis results, indicating similarities in the AMF community composition among root and soil samples from the same site and variation among samples from different sites. The indicator species analyses revealed 16 OTUs associated with rhizosphere soils while no OTUs with the root habitat were detected. Besides, most OTUs were shared by roots and rhizosphere soils from the same site, while few OTUs were shared among sites. These results indicated that AMF in rhizosphere soils could be considered to represent a pool of species, a fraction of which is recruited by plants 35,36 .
Xu et al. detected a strong influence of the plant community on AMF communities in soil, indicating preferences in plant-AMF associations 37 . These preferences have been observed on both local and global scale systems 36,38,39 , which might result from that host plants exhibited preferential allocation of photosynthates to more beneficial AMF partners 40 . Consistently, Saks et al. showed that root-colonizing AMF represent a phylogenetically clustered subset of AMF available in soil 41 . The conclusion is also proved in this study by the result that most of the AMF in root habitats were Glomeraceae. These results indicated that there might be niche preferences among AMF affecting AMF community composition 34,42 . Different to the study of AMF communities by Xu et al. which covered a wide range of vegetation types 37 , all soil samples were collected from rhizosphere of a single plant with similar size and age to reduce the plant factors affecting the AMF community in this study. Similar sampling strategy was also used for the study of AMF communities in semiarid Mediterranean soils 34 . We did not detect significant influences of the single plant on AMF communities in soils, indicating that the sampling strategy successfully reduce the biotic factors affecting the AMF distribution. Furthermore, by comparing influences of a single plant on AMF communities in roots and rhizosphere soils, we suggest that the effect of AMF communities in soils by plant communities might result from the different composition of plant communities in which each single plant shows the specific recruitment preference of AMF partners.
Glomus and Rhizophagus were the most abundant genera in root samples. Yang et al. reported that almost all of the sequences found in the roots of Elsholtzia splendens were Glomus species 43 . Long et al. revealed that although Ambispora, Kuklospora, and Glomus dominated in the rhizosphere soils of Phytolacca americana, only Glomus was detected in the roots 44 . Our results are in agreement with these previous studies. Glomus has been found to be dominant in various habitats, such as HM-polluted soils 45 , grassland soils 46 , and agricultural soils 47 . The Rhizophagus group is also a general AMF group that has been found in diverse host species and environments [48][49][50] . However, almost all of the fungi detected in R-DN and R-XC (root samples from DN and XC) were Glomus, the relative levels of which were much higher than those in soil habitats. On the contrary, the relative level of Glomus in R-GP (root samples from GP) was markedly lower than that in S-GP (soil samples from GP), while Rhizophagus was significantly more abundant in R-GP than S-GP. The Rhizophagus level was also higher in the root samples from TL and WH compared to the corresponding soil samples, while no obvious recruitment was detected for R-SS. These results indicated that the AMF recruitment preferences by C. ambrosioides differed among sampling sites. Several explanations for the differences in the AMF present in the roots and rhizosphere soils of the same plant have been proposed, including different AMF life strategies, differential sporulation dynamics, and seasonal changes in the AMF community [50][51][52] .
It has been shown that geographical distance influences AMF community at large spatial scales, especially in global-scale studies 37 . However, at local or landscape scales, soil abiotic factors are the key driver in shaping AMF community composition 10,11 . In our study, the geographic and climate parameters could independently explain a small fraction (< 10%) of variances of AMF community composition and phylogenetic structure, while a large fraction of variances (> 32%) could be explained independently by soil variables. These results revealed the complexity of factors regulating AMF communities, and the distribution patterns of AMF communities could not be completely explained by soil heterogeneity.
Despite high levels of Mn in soils of DN, GP and XC, no significant differences in the species diversity and richness of AMF communities were detected between these three samples and others. However, the richness indices of the R-GP and R-XC were markedly lower than those of the other samples, revealing that Mn affects the AMF community to some extent. These results are similar to Wei et al., which reported that root colonization and AMF diversity were negatively correlated with soil Mn concentration 17 . These results suggest also that other environmental factors may affect the AMF community and the effects of Mn on the AMF community may be confounded by those of other factors. Indeed, various factors influence the AMF community 41 . For example, Yang et al. showed that HM contamination is not the only soil parameter that affects the AMF community, and soil pH, Pb, Zn, Cd, and OM levels also have a great influence 18 ; Wei et al. concluded that changes in AMF diversity and colonization are not solely attributed to soil Mn concentration, while soil properties, especially N concentration, were also closely related to it 53 .
In this study, pH was found to be a significant factor influencing AMF communities. Similarly, some studies concluded that pH is the key environmental factor shaping AMF communities 54 . Soil acidity is one of the most important drivers of microbial communities, particularly for AMF 11,55 . Soil pH can directly change the physiological status of indigenous AMF, alter their ecological niches, and it could also indirectly influences the AMF community by regulating soil nutrient bioavailability, and impacting the mobility and sorption of metals. Several studies have reported strong effects of soil Zn and Cu levels on AMF abundance and diversity in soils 18,34,56 . These nutrients play important roles in plant metabolic processes, and their uptake is influenced by AMF 57,58 . It has been shown that total N is related to changes in the composition of soil microbial communities 59

Conclusion
The AMF communities of the roots and rhizosphere soils of C. ambrosioides in six areas were investigated. Both habitats (root or rhizosphere soil) and soil properties are the key environmental factors affecting AMF communities with C. ambrosioides. Total N, pH, Zn, and Cu levels play significant roles in triggering AMF populations. Overall, C. ambrosioides rhizospheric AMF communities and their influencing factors were illustrated, which contributed a better understanding of the AMF community shaping and related ecological factors.

Methods
Root and soil sampling. Samples were collected at six sites in Yunnan, Anhui, and Guangxi Provinces, China. Two sites in Yunnan province were located in Dounan Town (DN, 23°36′ N, 103°41′ E), Yanshan County, and Xincheng Town (XC, 23°38′ N, 102°27′ E), Shiping County. Three sites were located in Anhui Province, including Tongling (TL, 30°49′ N, 117°51′ E), Wuhu (WH, 31°19′ N, 118°26′ E) and Susong (SS, 30°9′ N, 116°10′ E) Counties. Another site was located in Guiping County (GP, 23°31′ N, 110°16′ E), Guanxi province. A single plant, C. ambrosioides, was investigated to reduce the number of biotic factors affecting the AMF community. All samples were collected in August 2015. Three healthy, similarly sized plants were randomly selected and collected at each of the six sites. Root systems were carefully excavated. Using clean tweezers and a brush, the soils bound to the surface of roots were carefully removed. The removed soils were defined as rhizosphere soil samples for DNA extraction. Then the roots were wrapped in tissue paper and stored in sterile Ziploc bags containing silica gel. Topsoil samples (0-15 cm depth) were collected from beneath each plant in at least three directions for chemical property analysis. To remove aboveground plant materials, roots, and stones, the mixed and homogenized soil samples were passed through a 2 mm sieve. After packing in sterile Ziploc bags, these soil samples were transported to the laboratory and stored at 4 °C to analyze physicochemical properties or at − 80 °C for DNA extraction.
Soil physicochemical properties. Soils were air-dried at room temperature, ground into a powder, and passed through a 0.15 mm nylon sieve. The soil pH was determined using a glass electrode with soil suspended in 0.01 M CaCl 2 (soil: solution ratio, 1:5). The nitrate N (NO 3 − -N), ammonia N (NH 4 + -N), and soil organic matter (OM) levels in the soil samples were determined according to our previous study 59 . Available K in soil was extracted with ammonium acetate 64 , and then determined by flame photometry (AAS novAA 400, Analytik Jena AG, Jena, Germany). Using sodium bicarbonate, the available P in soil was extracted and measured using the molybdenum blue method according to Watanabe and Olsen 65 . The total N level in soil was determined using the micro-Kjeldahl method 66 , and measured using a continuous flow analyzer (AA3, Bran + Luebbe, Hamburg, Germany). Soil samples for analyzing the levels of total P, total K, and HM (Mn, Zn, Cd, Cu, and Pb) were dried at 105 °C for 6 h, passed through a 0.15 mm nylon sieve, and then digested in a mixture of HNO 3 and HClO 4 (4:1, v/v) 67  Illumina MiSeq sequencing. To reduce potential early-round PCR errors, three independent PCR products for each sample were combined to construct a PCR amplicon library. The PCR products were subjected to agarose gel electrophoresis and purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer's protocol. Then the amplified DNA was quantified using QuantiFluor ™ -ST (Promega, Madison, WI, USA). The purified amplicons were paired-end (PE) sequenced on an Illumina MiSeq platform (Shanghai Biozeron Biotechnology Co., Ltd, Shanghai, China) according to standard protocols. In total, 36 sequencing libraries were constructed and independently sequenced.
Processing of sequencing data. The raw Illumina sequencing data were quality-filtered using the Trimmomatic software 69 . The reads were truncated at any site that received an average quality score < 20 over a 50 bp sliding window, and the truncated reads shorter than 50 bp were discarded. Then, PE reads were assembled Scientific RepoRts | 7:45134 | DOI: 10.1038/srep45134 according to their overlap sequence with a minimum overlap length of 10 bp, discarding reads that could not be assembled. Sequences that contained more than one ambiguous character or two nucleotide mismatches in the primers, and those with a mismatch ratio within the overlap region of more than 0.2 were removed. The clean sequences were analyzed using the QIIME pipeline 70 . Chimeric sequences were identified and removed using UCHIME 71 . Operational taxonomic unit (OTU) grouping was performed using UPARSE at 97% similarity 72 , and then the representative sequences obtained for each OTU were assigned to taxonomic data using the RDP classifier at a 70% threshold 73 . The rarefaction curves, indices of ACE and Chao1, and Shannon diversity were analyzed using the Mothur software 74 . The Bray-Curtis distances were calculated using the QIIME pipeline 70 , and a principal coordinate analysis (PCoA) was performed using R (http://www.r-project.org/) based on the Bray-Curtis similarities. A Venn diagram based on unique and shared OTUs was produced using R to characterize the differences and similarities among the AMF communities. The indicator species analyses were performed to test whether there were specific OTUs associated with the rhizosphere soil habitat or root habitat, and the indicator value index was used to measure the associations 75 . The analyses were performed using the indcspecies package implemented in R with a permutation test (999 permutations) 76 . To examine the relationships between AMF communities and soil properties, redundancy analysis (RDA) was applied. This analysis was conducted with CANOCO for Windows 77 with 999 permutations of Monte Carlo permutation tests.

Phylogenetic analysis.
A neighbor-joining tree containing the type sequences of all OTUs was constructed using MEGA v6.0 with 1000 replicates 78 . All sequences were aligned, concatenated, and manually adjusted using Geneious Pro v4.8.3 (http://www.geneious.com/). The best-fit model for the datasets was selected using jModel-Test v2 79 .
Statistical analyses. PCNM vectors were calculated using the 'pcnm' function of the 'vegan' package with the R language. Prior to the calculation, GPS coordinates were converted to UTM coordinates in kilometers. PCNM vectors were used as explanatory spatial variables for db-RDA 80 . The weighted Unifrac distance matrix was used as a measure to determine the AMF phylogenetic structure, which was calculated using 'unifracs' function of 'GUniFrac' package with the R language. db-RDA 81 , which is known as constrained analysis of principal coordinates, was used to investigate the variations of AMF communities that were attributable to environmental factors. By a stepwise db-RDA, the contributions of environmental factors to the variation of soil AMF communities were summarized using Bray-Curtis dissimilarity and Unifrac distance matrices separately. The forward-selection for the environmental factors in three groups (soil, geographic and climate, and plant variables) was performed independently with the adjR2thresh stopping criterion 82 , and then the contribution of each of the groups or the combined groups was determined. The amount of variances explained by the individual and combined groups was tested using Monte Carlo permutation tests (999 permutations). The stepwise db-RDA was performed using 'capscale' function of 'vegan' package with the R language. PerMANOVA for the relationship between AMF community composition and each environmental variable was carried out separately based on Bray-Curtis dissimilarity distance and weighted Unifrac distance matrices using 'adonis' function of 'vegan' package with the R language. Tukey's and Kendall test was used for multiple comparisons (P < 0.05), which were performed using SPSS version 19.0 (SPSS Inc., Chicago, IL, USA).