Diversity and Co-Occurrence Patterns of Soil Bacterial and Fungal Communities of Chinese Cordyceps Habitats at Shergyla Mountain, Tibet: Implications for the Occurrence

Chinese Cordyceps is a well-known medicinal larva-fungus symbiote distributed in the Qinghai-Tibetan Plateau and adjacent areas. Previous studies have involved its artificial cultivation but commercial cultivation is difficult to perform because the crucial factors triggering the occurrence of Chinese Cordyceps are not quite clear. The occurrence of Chinese Cordyceps is greatly affected by the soil environment, including the soil’s physicochemical and microecological properties. In this study, the effects of these soil properties on the occurrence of Chinese Cordyceps were investigated. The results show that the physicochemical properties, including easily oxidizable organic carbon (EOC), soil organic carbon (SOC), humic acid carbon (HAC), humin carbon (HMC), and pH, might be negatively related to the occurrence of Chinese Cordyceps, and soil water content (SWC) might be positively related. Several soil physicochemical parameters (pH, SOC, HMC, HAC, available potassium (APO), available phosphorus (APH), microbial biomass carbon (MBC), and the ratio of NH4+ to NO3− (NH4+/NO3−)) and microbial properties interact and mix together, which might affect the occurrence of Chinese Cordyceps. Soil microbial community structure was also a possible factor, and a low level of bacterial and fungal diversity was suitable for the occurrence of Chinese Cordyceps. The intra-kingdom network revealed that a closer correlation of the bacterial community might help the occurrence of Chinese Cordyceps, while a closer correlation of the fungal community might suppress it. The inter-kingdom network revealed that the occurrence rate of Chinese Cordyceps might be negatively correlated with the stability of the correlation state of the soil habitat. In conclusion, this study shows that soil physicochemical properties and microbial communities could be greatly related with the occurrence of Chinese Cordyceps. In addition, soil physicochemical properties, the level of bacterial and fungal diversity, and correlations of bacterial and fungal communities should be controlled to a certain level to increase the production of Chinese Cordyceps in artificial cultivation.


Introduction
Ophiocordyceps sinensis, also known as Cordyceps sinensis, is a well-known symbiote of fungus (O. sinensis) and larva (Thitarodes, Hepialidae, Lepidoptera) [1]. It is particularly distributed in the Qinghai-Tibetan Plateau and adjacent high-altitude areas [2] and is usually called Chinese Cordyceps for The germination and growth of a fungus can be suppressed by natural soils to a certain extent, and this phenomenon is referred to as soil fungistasis [14]. The intensity of fungistasis is dependent on the soil physical and chemical properties as well as microbial activity [15,16]. Among them, the soil microbial community and activity are influenced by the physicochemical characteristics [17], and in turn, the soil microbiome also plays an important role in biogeochemical processes, such as nitrogen, phosphorus, and other element cycles, and is a key factor for soil health and productivity [18]. Thus, it is speculated that the soil microenvironment is closely related to the survival of O. sinensis and the occurrence of Chinese Cordyceps.
Furthermore, some soil microbes are also influenced by coexisting soil microbes in complicated interaction systems [19,20]. For example, for arbuscular mycorrhizal fungi, fungal spore germination, mycelial growth, and root colonization can be stimulated by some growth factors The germination and growth of a fungus can be suppressed by natural soils to a certain extent, and this phenomenon is referred to as soil fungistasis [14]. The intensity of fungistasis is dependent on the soil physical and chemical properties as well as microbial activity [15,16]. Among them, the soil microbial community and activity are influenced by the physicochemical characteristics [17], and in turn, the soil microbiome also plays an important role in biogeochemical processes, such as nitrogen, phosphorus, and other element cycles, and is a key factor for soil health and productivity [18]. Thus, it is speculated that the soil microenvironment is closely related to the survival of O. sinensis and the occurrence of Chinese Cordyceps.
Furthermore, some soil microbes are also influenced by coexisting soil microbes in complicated interaction systems [19,20]. For example, for arbuscular mycorrhizal fungi, fungal spore germination, mycelial growth, and root colonization can be stimulated by some growth factors produced by the mycorrhiza helper bacteria [21]. Thus, research on the relationships among microbial species could aid in understanding interspecies interactions and promote the understanding of the niche spaces among community members [22,23], and these investigations might help to clarify the mechanism of the occurrence of Chinese Cordyceps. Yang et al. [24] first attempted to investigate the bacterial communities in the habitats of Chinese Cordyceps on the Tibetan Plateau and identified several unique and shared taxa of different soil samples with or without the presence of Chinese Cordyceps. Xia et al. [25] reported fungal diversity of the soil adhering to the surface of the membrane covering Chinese Cordyceps. However, the comprehensive effects considering soil physicochemical properties and bacterial and fungal communities, especially the network interactions among these properties on the occurrence of Chinese Cordyceps, remain insufficient.
In this study, an analysis of the physicochemical properties and the bacterial and fungal communities based on HiSeq sequencing of 16S rRNA and internal transcribed spacer (ITS) genes of the habitat soil of Chinese Cordyceps was performed. The purpose of this study was to investigate and explore the influences of soil physicochemical properties and microbial properties and the correlations among these properties on the occurrence of Chinese Cordyceps.

Field Site Description and Sample Collection
Chinese Cordyceps generally occur in the Qinghai-Tibetan Plateau with the altitudes of higher than 3,500 m, and in this study, native habitats of Chinese Cordyceps at Shergyla Mountain, Tibet, were selected as the study region. From 2006 to 2016, a pre-survey was conducted through field investigation on the density of Thitarodes larvae and Chinese Cordyceps, and the occurrence rate of Chinese Cordyceps of each study area was accordingly assessed. Field investigation revealed that the peaks of activity, feeding, growth, development, and population density of Thitarodes larvae usually occurred in June to August every year, especially around mid-July, which is coincidently the time for the eruption of ascospores from mature stroma of Chinese Cordyceps. Furthermore, during these months, the Thitarodes larvae preferred to be situated at the soil depth of 10-20 cm below the ground.
For this study, three sites (each with the area of about 6.6 hectares) with different occurrence rates of Chinese Cordyceps were ultimately selected: these sites were named A, B, and C, and the detailed information is presented in Table 1. Briefly, site A had a high density of Thitarodes larvae and Chinese Cordyceps (occurrence rate of Chinese Cordyceps: 10.0%), site B had a high density of Thitarodes larvae and low density of Chinese Cordyceps (occurrence rate of Chinese Cordyceps: 1.4%), and site C had a high density of Thitarodes larvae and no Chinese Cordyceps (occurrence rate of Chinese Cordyceps: 0). The yearly variations of the occurrence rates remained relatively stable at each site, with relative standard deviations (RSD) at site A and B of less than 10%, and null Chinese Cordyceps was observed at site C for ten years. At each site, five soil samples were sampled by the diagonal line five-point method in mid-July 2016. During the sampling, weeds on the ground were removed first, and then the soil was cut into v-shaped pits with a shovel. The soil slices of 10-20 cm deep and 10 cm wide were specially sampled. The freshly-collected samples were kept at −20 • C in a portable refrigerator and transported to the laboratory. In the laboratory, roots, plant residues, and stones in the soil samples were removed by sieving through a 2 mm mesh. Each soil sample was divided into 2 parts for DNA extraction and analysis of soil physicochemical properties. The soil samples were stored in −80 • C before the analysis. Total DNA was extracted from 0.5 g of soil using a MO BIO PowerSoil ® DNA Isolation Kit (MO BIO Laboratories Inc., Carlsbad, CA, USA) according to the manufacturer's instructions. Pipetting and DNA purification were performed on the Microlab ® STAR line workstation (Hamilton, Bonaduz, Switzerland) and KingFisher Flex purification system (Thermo Fisher Scientific, Vantaa, Finland). The purified DNA was diluted with 100 µL of RNA-free water (TaKaRa, Dalian, China) and stored at −20 • C for further analysis. The DNA concentrations were measured by a NanoDrop ND−3300 spectrophotometer (NanoDrop Technologies, Thermo Scientific, Wilmington, DE, USA).
The V4 region of 16S rDNA gene and the ITS2 region of the fungal ITS gene were used as the bacterial-specific fragment and fungal-specific fragment respectively, and the primer pairs 515F/806R and ITS3/ITS4 were used to amplify bacterial-specific and fungal-specific fragments, respectively. These primer pairs were modified by a 12-bp barcode sequence at the 5 -end to identify each sample. All amplifications were performed in 50 µL reactions containing 0.5 units of Ex Taq DNA polymerase (TaKaRa, Dalian, China), 10 µL 1× Ex Taq loading buffer (TaKaRa, Dalian, China), 8 µL dNTP mix (TaKaRa, Dalian, China), 2 µL of each primer (10 mM), and 10-100 ng template DNA by an ABI GeneAmp ® 9700 PCR System (Applied Biosystems, Waltham, MA, USA). The PCR conditions for bacterial-specific fragments were as follows: 3 min of initial denaturation step at 95 • C, followed by 35 cycles at 94 • C for 30 s, 55 • C for 1 min, and 72 • C for 1 min, and the extension step at 72 • C for 10 min. The PCR conditions for fungal-specific fragments were as follows: 5 min of initial denaturation at 95 • C, followed by 30 cycles at 95 • C for 30 s, 52 • C for 30 s, and 72 • C for 45 s, and the extension step at 72 • C for 10 min. Each sample was amplified 3 times and the products were mixed. After evaluated by 2% agarose gel, the mixed products were purified by an E.Z.N.A. ® Gel Extraction Kit (Omega Bio-tek, Norcross, GA, USA) and quantified with a QuantiFluor TM-ST fluorometer (Promega, Madison, WI, USA). In the end, the PCR products were pooled at equimolar concentrations and sequenced on an Illumina MiSeq PE300 platform at Ozimeks Biotech Co., Ltd. (Shenzhen, Guangdong, China). The obtained sequence data were uploaded to the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) (PRJNA533894).
The sequence data were further analyzed by QIIME (v. 1.9.0; http://qiime.org/), and the UPARSE pipeline [28] was used for taxonomic assignment with similarities higher than 97%. The taxonomic classification was carried out by the SILVA (v. 119; http://www.arb-silva.de) and UNITE (v. 7.0; http://unite.ut.ee/index.php) databases for bacteria and fungi, respectively. The operational taxonomic unit (OTU) and its tags, which were annotated as chloroplasts or mitochondria (16S amplicons) and cannot be annotated to the kingdom level, were removed, then the OTU taxonomy synthesis information table for the final analysis was generated (Supplementary Table S1). All samples were subsequently subsampled based on the minimum soil microbial sequencing depth in the current study to preclude bias from several sequencing depths.

Statistical Analysis
The statistical analysis of soil physicochemical properties was performed with SPSS (v. 21.0; https://www.ibm.com/analytics/spss-statistics-software). QIIME [29] was used to analyze the Shannon, Simpson, and Chao1 bacterial and fungal diversity indices. Heatmaps were analyzed to compare the top 40 classified genera of bacteria and fungi by the gplot package in R (v. 3.5.3; https://www.r-project.org/). Bray-Curtis, weighted and unweighted unifrac beta diversity indexes were calculated by QIIME software, and accordingly, nonmetric multidimensional scaling (NDMS) analysis was performed and visualized by the vegan package of R software to further illustrate the beta-diversity of soil microbial structure. R was also used to calculate the analysis of similarities (ANOSIM), nonparametric multivariate analysis of variance (MANOVA, or Adonis), and multi-response permutation procedure (MRPP) to compare bacterial and fungal community differences among all soil samples with Bray-Curtis distance and 999 permutations in the R vegan package. The linear discriminant analysis (LDA) effect size (LEfSe) algorithm (http://huttenhower.sph.harvard.edu/galaxy/; last access 23 May 2019) [30] was used to identify the taxa that were present in different abundances between the Chinese Cordyceps group (sites A and B) and null Chinese Cordyceps group (site C). LEfSe emphasizes statistical significance and biological relevance. The effect size threshold of LDA score was set to 3.5 for analysis in this study. To examine the effects of soil physicochemical properties on structuring microbial communities, Mantel tests and redundancy analyses (RDA) were performed, and the results were visualized by the vegan package of R software. To demonstrate the relationships among different microbial species at each site, intra-kingdom network analysis was conducted using the 40 most abundant bacterial and fungal genera, and inter-kingdom network analysis was conducted using different microbial families (including bacterial and fungal families) at each site. Highly significant positive (R > 0.80, false discovery rate (FDR) < 0.05) and negative (R < −0.80, FDR < 0.05) Pearson correlations were screened out and co-occurrence patterns were visualized as networks using Cytoscape version 3.6.0 (https://cytoscape.org/) [31]. The size and color of each node represented the number of connections and taxonomy, respectively.

Soil Physicochemical Properties
The soil physicochemical properties for each site are shown in Table 2

Soil Microbial Diversity
Through the raw sequencing reads, there was a total of 602,677 high-quality 16S rDNA sequences and 256,324 high-quality ITS2 sequences. In addition, 10,762 bacterial OTUs and 3361 fungal OTUs were clustered from these high-quality sequences with a 97% identity threshold. Alpha diversity was applied in analyzing complexity of species diversity for a sample through five indices, including two indices to identify community richness: observed species and Chao1, and three indices to identify community diversity: Shannon, Simpson, and Dominance [32]. The alpha diversity indices of bacterial and fungal communities are shown in Table 3. For bacteria, the diversity (represented by Shannon and Simpson indices) of site C was significantly higher than that of site B (p < 0.05), and the bacterial richness (represented by Chao 1) of site C was significantly higher than that of sites A and B. For fungi, the diversity (represented by Shannon and Simpson indices) of sites B and C was significantly higher than that of site A (p < 0.05). Beta diversity analysis was used to evaluate differences of samples in species complexity. In this study, beta-diversity at the OTU level by NMDS [33] analysis is shown in Figure 2. According to Euclidean distance dissimilarity, the bacterial ( Figure 2a) and fungal (Figure 2b) beta-diversities were different for each site and the samples within each site were clearly grouped together. In addition, significant differences (p < 0.05) were found among all soil samples in aspects of bacterial and fungal communities, as shown in Table 4, corresponding to ANOSIM, Adonis, and MRPP analysis.

Fungi
Site Beta diversity analysis was used to evaluate differences of samples in species complexity. In this study, beta-diversity at the OTU level by NMDS [33] analysis is shown in Figure 2. According to Euclidean distance dissimilarity, the bacterial ( Figure 2a) and fungal (Figure 2b) beta-diversities were different for each site and the samples within each site were clearly grouped together. In addition, significant differences (p < 0.05) were found among all soil samples in aspects of bacterial and fungal communities, as shown in Table 4, corresponding to ANOSIM, Adonis, and MRPP analysis.  The influence from soil physicochemical factors on the soil microbial community was examined by Mantel tests (Table S2) and redundancy analysis (RDA), and the significant relationships (p < 0.05) between soil physicochemical factors and microbial communities are shown in Figure 3. In Figure 3, soil physicochemical factors are represented by arrows. The length of the arrow line represents the degree of correlation between the physicochemical factor and microbial communities, and the  The influence from soil physicochemical factors on the soil microbial community was examined by Mantel tests (Table S2) and redundancy analysis (RDA), and the significant relationships (p < 0.05) between soil physicochemical factors and microbial communities are shown in Figure 3. In Figure 3, soil physicochemical factors are represented by arrows. The length of the arrow line represents the degree of correlation between the physicochemical factor and microbial communities, and the projection distance of each soil sample on the arrow line represents the degree of correlation between the physicochemical factor and the sample [34]. Thus, from Figure 3a

Soil Bacterial and Fungal Structure
The relative compositions of soil bacterial and fungal communities at the phylum, class, order, family, and genus level are presented in Supplementary Tables S3 and S4, respectively. The predominant phyla (with relative abundance higher than 0.01%) of bacterial communities were Proteobacteria, Acidobacteria, Verrucomicrobia, Actinobacteria, Planctomycetes, Bacteroidetes, Chloroflexi, Firmicutes, WPS−2, Gemmatimonadetes, WPS−1, Parcubacteria, Armatimonadetes, and Euryarchaeota, occupying more than 90% of the total sequences ( Figure 4a). Five predominant phyla (with relative abundance higher than 0.01%) of fungal communities were identified: Ascomycota, Basidiomycota, Glomeromycota, Zygomycota, and Chytridiomycota, occupying more than 50% of the total sequences ( Figure 4b). Besides these five fungal phyla, the proportion of other fungal phyla was especially enriched in site A.
The 40 most abundant genera of soil bacterial and fungal communities can be seen in Figures 5a and 5b respectively, and the relative abundance of soil microbial community from high to low is represented by red, black, and green. The top genera varied among the sample sites: for bacterial communities, site A had high relative abundance of Gp6, Opitutus, Rhizomicrobium, Ktedonobacter, Aciditerrimonas, Conexibacter, Gaiella, and Flavobacterium, site B had high relative abundance of Blastocatella, Kofleria, Gp4, Gp7, Gp2, Gp1, Candidatus Koribacter, Gp3, Candidatus Solibacter, and

Differential Operational Taxonomic Units (OTUs) Related to the Occurrence of Chinese Cordyceps
In order to discuss the detailed OTUs that might be related to the occurrence of Chinese Cordyceps, the differential OTUs between the Chinese Cordyceps group (sites A and B) and the null Chinese Cordyceps group (site C) were screened out as biomarkers using linear discriminant analysis (LDA) effect size analysis. The LDA score of these biomarkers was illustrated by the histograms (Figure 6a,b), the length of each bar representing the degree of the differences. The taxonomic information of these biomarkers was illustrated by cladograms (Figure 6c,d), with the circles radiating from the center point representing the taxonomic levels from phylum to species. Fourteen bacterial OTUs (mostly belonging to the classes Ktedonobacteria, Solibacteres, Acidobacteria, and Verrucomicrobia) and four fungal OTUs (mostly belonging to Capnodiales, Humicola, and PeltigeraspSW330) presented significantly higher abundance in the soil samples of the Chinese Cordyceps group (sites A and B). Twenty bacterial OTUs (mostly belonging to the classes Methanobacteria, Betaproteobacteria, Chloroflexia, Deltaproteobacteria, Planctomycetia, Pirellula, and Acidobacteria) and 24 fungal OTUs (mostly belonging to the classes Leotiomycetes, Dothideomycetes, Agaricomycetes, Sordariomycetes, Glomeromycetes, and Glomeromycetes) were significantly enriched in the null Chinese Cordyceps group (site C).
In terms of the most concerned Cordyceps-related families [35], Clavicipitaceae, Cordycipitaceae, and Ophiocordycipitaceae were identified but in minor abundance in this study: Cordycipitaceae appeared in all samples, and Ophiocordycipitaceae and Clavicipitaceae presented preference in the Chinese Cordyceps group (sites A and B) (Figure 7a). Among them, Metarhizium, Pochonia, Simplicillium, Elaphocordyceps, Polycephalomyces, Purpureocillium, and Tolypocladium presented preference in the Chinese Cordyceps group (Figure 7b,c). (Figures 6a and 6b), the length of each bar representing the degree of the differences. The taxonomic information of these biomarkers was illustrated by cladograms (Figures 6c and 6d), with the circles radiating from the center point representing the taxonomic levels from phylum to species. Fourteen bacterial OTUs (mostly belonging to the classes Ktedonobacteria, Solibacteres, Acidobacteria, and Verrucomicrobia) and four fungal OTUs (mostly belonging to Capnodiales, Humicola, and PeltigeraspSW330) presented significantly higher abundance in the soil samples of the Chinese Cordyceps group (sites A and B). Twenty bacterial OTUs (mostly belonging to the classes Methanobacteria, Betaproteobacteria, Chloroflexia, Deltaproteobacteria, Planctomycetia, Pirellula, and Acidobacteria) and 24 fungal OTUs (mostly belonging to the classes Leotiomycetes, Dothideomycetes, Agaricomycetes, Sordariomycetes, Glomeromycetes, and Glomeromycetes) were significantly enriched in the null Chinese Cordyceps group (site C).  In terms of the most concerned Cordyceps-related families [35], Clavicipitaceae, Cordycipitaceae, and Ophiocordycipitaceae were identified but in minor abundance in this study: Cordycipitaceae appeared in all samples, and Ophiocordycipitaceae and Clavicipitaceae presented preference in the Chinese Cordyceps group (sites A and B) (Figure 7a). Among them, Metarhizium, Pochonia, Simplicillium, Elaphocordyceps, Polycephalomyces, Purpureocillium, and Tolypocladium presented preference in the Chinese Cordyceps group (Figure 7b,c).

Intra-kingdom Co-Occurrence Analysis
Network analysis is widely performed to explore the interactions of microbial taxon in the complex microbial communities. In this study, network analysis was applied to illustrate the differences of soil bacterial and fungal communities among different sampling sites (Figures 8-11), and the network properties are summarized in Table 5. In Figures 8-11, red and blue lines

Intra-kingdom Co-Occurrence Analysis
Network analysis is widely performed to explore the interactions of microbial taxon in the complex microbial communities. In this study, network analysis was applied to illustrate the differences of soil bacterial and fungal communities among different sampling sites (Figures 8-11), and the network properties are summarized in Table 5. In Figures 8-11, red and blue lines represented that the abundances between the two connected genera (or families) are positively correlated and negatively correlated, respectively.  , , , and represent Ascomycota, Basidiomycota, Glomeromycota, and Zygomycota, respectively.  Based on the top 40 bacterial and fungal genera, the intra-kingdom network analysis of the correlations is shown in Figure 8 and Table 5. For bacterial analysis (Figure 8a,c,e), the co-occurrence analysis showed 16 positive and 4 negative significant correlations, 13 positive and 6 negative significant correlations, and 4 positive and 4 negative significant correlations at sites A, B, and C, respectively. For fungal analysis (Figure 8b,d,f), the co-occurrence analysis showed 6, 7, and 11 significant positive correlations at sites A, B, and C respectively, and only 4 significant negative correlations at site B, with no significant negative correlations observed at sites A and C. The average degree of bacterial communities was 1.143, 1.256, and 0.571, and fungal communities was 0.414, 0.667, and 0.733 respectively, at sites A, B, and C (Table 5). Hence, there were more positive correlations and closer relationships of bacterial taxa at sites A and B, and they were highest for the fungal taxa at site C.

Discussion
Chinese Cordyceps is the outcome of the infection, colonization, and growth of O. sinensis on Thitarodes larvae. Previous studies revealed that the growth of fungi is dependent on physical, chemical, and microbial properties of soil [15,16]. Thus, the soil microenvironment is closely related to the occurrence rate of Chinese Cordyceps. This study revealed that the soil microenvironment is a complex community, and these properties might be correlated or might co-affect the occurrence of Chinese Cordyceps.

Soil Physicochemical Properties
Limited studies have reported that soil physicochemical properties, including available K and pH, could affect the occurrence of Chinese Cordyceps [12,37]. In this study, it was found that there were significant differences in pH, NH4 + -N, NO3 --N, APH, APO, MBC, DOC, HEC, DOC, and HEC between site C and site A (or site B) ( Table 2). These findings reflect that besides K and pH, various soil physicochemical properties could affect the occurrence of Chinese Cordyceps. Nevertheless, the contents of those physicochemical factors did not present progressive increase or decrease of the occurrence rates, i.e., the contents presented A > B > C or A < B < C. For instance: (1) the SWC presented A≈ B > C, (2) the contents of EOC, SOC, HAC and HMC presented A < B ≈ C, (3) pH presented A ≈ B < C, (4) the contents of NH4 + -N and NO3 --N presented C < A < B while (5) the ratio of NH4 + -N/NO3 --N presented B < A ≈ C, and (6) APH, APO, MBC, DOC and HEC presented A < C < B. Therefore, the levels of EOC, SOC, HAC, HMC, and pH might significantly inhibit the occurrence of Chinese Cordyceps, and the levels of SWC might be positively related with the occurrence of Chinese Cordyceps, and the levels of NH4 + -N, NO3 --N, NH4 + -N/NO3 --N, APH, APO, MBC, DOC, and HEC might be related to the occurrence of Chinese Cordyceps, while a detailed relationship could not be clarified based on the current study. In addition, soil is a rather complicated matrix, and all of the physicochemical factors might co-affect Thitarodes larvae, and the occurrence of Chinese Cordyceps

Inter-kingdom Co-Occurrence Analysis and Determination of Hub Taxa
Among all of the bacterial and fungal families detected in the three sites, highly significant positive (R > 0.80, FDR < 0.05) and negative (R < −0.80, FDR < 0.05) Pearson correlations were found among all families detected at the three sites. The corresponding network properties are summarized in Table 5 and occurrence patterns are visualized as a network in Figure 9 (site A), Figure 10  In terms of Cordyceps-related fungal families at site A (Ophiocordycipitaceae, Clavicipitaceae, and Cordycipitaceae), Cordycipitaceae was positively correlated with Parmeliaceae, Venturiaceae, Cladosporiaceae, Geobacteraceae, and Hygrophoraceae, and negatively correlated with Solibacteraceae, at site B, Ophiocordycipitaceae was positively correlated with Minutisphaeraceae and negatively correlated with Syntrophobacteraceae, and Cordycipitaceae was negatively correlated with Acetobacteraceae, Leotiaceae, and Davidiellaceae, at site C, Cordycipitaceae was positively correlated with Vibrisseaceae, Acidothermaceae, and Sphingobacteriaceae, and negatively correlated with Xanthomonadaceae.

Discussion
Chinese Cordyceps is the outcome of the infection, colonization, and growth of O. sinensis on Thitarodes larvae. Previous studies revealed that the growth of fungi is dependent on physical, chemical, and microbial properties of soil [15,16]. Thus, the soil microenvironment is closely related to the occurrence rate of Chinese Cordyceps. This study revealed that the soil microenvironment is a complex community, and these properties might be correlated or might co-affect the occurrence of Chinese Cordyceps.

Soil Physicochemical Properties
Limited studies have reported that soil physicochemical properties, including available K and pH, could affect the occurrence of Chinese Cordyceps [12,37]. In this study, it was found that there were significant differences in pH, NH 4 + -N, NO 3 − -N, APH, APO, MBC, DOC, HEC, DOC, and HEC between site C and site A (or site B) ( Table 2). These findings reflect that besides K and pH, various soil physicochemical properties could affect the occurrence of Chinese Cordyceps. Nevertheless, the contents of those physicochemical factors did not present progressive increase or decrease of the occurrence rates, i.e., the contents presented A > B > C or A < B < C. For instance: (1) the SWC presented A≈ B > C, illustrated that soil physicochemical factors influenced soil bacterial and fungal community structure in different degrees. On the other hand, soil microbial community is also important to the soil ecosystem stability and nutrient transformation, and in turn would affect the soil physicochemical properties [18]. Therefore, soil physicochemical parameters and microbial properties could interact and mix together, which might affect the occurrence of Chinese Cordyceps.

Soil Microbial Structure
Besides the interactions between physicochemical and microbial factors, the interactions between different microbial species might also be related to the occurrence. For instance, soil microbial community structure could interact with fungi and affect the soil fungistasis [14]. On the other hand, some fungi-helper microbiota could produce growth factors to stimulate fungal spore germination, mycelial growth, and host colonization, and reduce environmental stress by alleviating the toxification of antagonistic substances and inhibiting competitors and antagonists [21]. Besides O. sinensis, some other species were proved to be closely associated with Chinese Cordyceps [11] and may be involved in the development of its stromata [40]. Thus, the soil microbial community might be important to the occurrence of Chinese Cordyceps. In the current study, it was found that the bacterial diversity of sites A and B were significantly less than that of site C, and the fungal diversity of site A was significantly less than that of sites B and C ( Table 3). As the occurrence of Chinese Cordyceps was highest in site A, these findings indicate that the decreased microbial diversity was an advantage for the occurrence of Chinese Cordyceps. The increased microbial diversity might result in higher antibiosis and soil fungistasis, which might have negative effects on the development of O. sinensis, and ultimately suppress the occurrence of Chinese Cordyceps.
In terms of the composition of microbial taxa, at the bacterial phylum level, similar patterns were observed among all sampling sites (Figure 4), with Proteobacteria, Acidobacteria, and Verrucomicrobia being the most abundant bacterial phyla, which is in accordance with bacterial research on the habitat soil in Tibet [24]. For fungi, although Ascomycota and Basidiomycota were the most abundant fungal phyla among all sites, a perceived enrichment of other fungal phyla especially existed in site A. Coincidentally, fungal research on the soil adhering to the surface of the membrane covering Chinese Cordyceps also showed a predominance (approximately 90% in relative abundance) of other fungal phyla (unclassified fungi) [25]. Therefore, our study further indicates that some other unclassified fungi might be positively related to the occurrence of Chinese Cordyceps, and their role should be clarified in future research. At the OTU level, NMDS showed that both bacterial and fungal β-diversities were different for each sampling site, and the samples within each site were obviously grouped closely. At a genus level, the top 40 genera of soil bacterial and fungal communities ( Figure 3) differed at different sites. Thus, the soil microbial taxonomic composition was different at both the OTU and genus level among the sites. Beside the top 40 genera above, 20 bacterial and 24 fungal OTUs were significantly enriched in the null Chinese Cordyceps group (site C) based on LDA (with LDA scores higher than 3.5). Among them, most of the bacteria may be responsible for the removal of anthropogenic compounds (Rhodocyclales) [41], autotrophic bacteria (Chloroflexi) [42], or anaerobic bacteria (Methanobacteria) [43]. The presence of those bacteria indicates that the soil habitat at site C might be not friendly to the heterotrophic and anaerobic Chinese Cordyceps. The fungi belonging to Helotiales can decrease during larval development and increase in the survival and fecundity rate of Lepidoptera [44]. Thus, the enriched Helotiales at site C might be a benefit for the Thitarodes insect and indirectly suppress the occurrence of Chinese Cordyceps. Conclusively, the analysis of microbial taxa indicates that the microbial composition in the habitat soil varied among different sampling sites and might be closely related to the occurrence of Chinese Cordyceps.
In the natural environment, O. sinensis first colonizes in alpine plant roots and is transferred to Thitarodes larvae through feeding behavior [45]. The occurrence of Chinese Cordyceps on Thitarodes larvae is dominated by O. sinensis and accompanied by the growth of various microorganisms. Therefore, some scholars have proposed that Chinese Cordyceps should be studied as a unified microbial ecosystem. Various microorganisms have been identified from natural Chinese Cordyceps [11,25], and even sparked controversy regarding the anamorph of Chinese Cordyceps. Zhu et al. found that Paecilomyces hepiali coexisted with Hirsutella sinensis [40] (the most approved anamorph of Chinese Cordyceps [1]) in natural Chinese Cordyceps, and even underwent dynamic changes with the maturation process of Chinese Cordyceps. They proved that a combined infection of Paecilomyces hepiali and Hirsutella sinensis can significantly improve the infection efficiency [46], suggesting that the synergistic infection effect of multiple microbiota may be an important link in the occurrence of Chinese Cordyceps. In this study, Cordyceps-related fungi (Ophiocordycipitaceae, Clavicipitaceae, and Cordycipitaceae) [35] were commonly detected in minor abundance among all detected sample sites (Figure 7). In the soil, most of the entomogenous fungi, including Cordyceps-related fungus and O. sinensis, mostly prefer to colonize in plant roots and derive nutrition from plant sources in the absence of insect hosts [45]. Thus, the detected relative abundance of these three families in this study was extremely low; however, despite that, Ophiocordycipitaceae and Clavicipitaceae showed a preference in sites A and B. These findings indicate that the soil conditions of sites A and B might be better for the survival of Cordyceps-related fungi, and several Cordyceps-related fungal genera, such as Metarhizium, Pochonia, Cordyceps, Engyodontium, Chaunopycnis, Elaphocordyceps, Haptocillium, Polycephalomyces, Purpureocillium, and Tolypocladium, might play a positive role in the occurrence of Chinese Cordyceps. Nevertheless, O. sinensis was absent among all of the soil samples in this study. In order to avoid disturbance of the physicochemical and microbial analysis of the soil, we had removed the roots before analysis, and the O. sinensis might have been removed along with the roots, which would have resulted in its absence. Thus, plants are the essential medium for the occurrence of Chinese Cordyceps, not only for providing food for host Thitarodes larvae, but also for providing the living microhabitat for O. sinensis before encountering the host Thitarodes larvae. These findings indicate that the factors of rhizosphere, such as plant roots and Cordyceps-related fungi, might play a role in the occurrence of Chinese Cordyceps, while their definite role should be proved through an infection trial in future study.

Co-Occurrence Interactions of Soil Microbial Community
Besides the above-mentioned microbial diversity and composition, for each microhabitat, the microbiota establishes a complicated community with interactions among microbial species, including a range of complex positive (commensalism, mutualism) and negative (amensalism, parasitism or predation, and competition) interactions occurring among different microbial species [47]. The network analysis of microbial co-occurrence patterns, especially showing positive or negative correlations, can provide a new perspective to investigate the structure of complex microbial communities, potential microbial interactions, and their ecological roles [48,49]. The occurrence of Chinese Cordyceps is actually the outcome of the infection, colonization, and growth of an entomogenous fungus (O. sinensis) on host insects living in the soil. Thus, the soil ecological status, including the correlations of microbiomes in the habitat [47], would influence the infection process of O. sinensis.
In order to investigate the intra-kingdom correlations among the predominant microbial genera, we carried out an intra-kingdom network analysis based on the positive and negative correlations among the 40 most abundant bacterial and fungal genera of the soil microbial community (Table 5). Co-occurrence patterns show that the network compositions substantially differed among different sites and each network had a distinct set of module hubs and connectors (Figure 8). For the bacterial community, there were more positive correlations and average degrees in the Cordyceps group (sites A and B), while this was reversed in the fungal community, with the most in the null Cordyceps group (site C). The findings indicate that a closer correlation of the bacterial community might help for the colonization of O. sinensis at sites A and B, while antagonism might exist between O. sinensis and the other fungal genera, and these fungal genera were positively correlated and might synergistically suppress the occurrence of Chinese Cordyceps at site C.
Previous studies revealed that dysbiosis in the bacterial kingdom and extensive synergistic networks in the fungal kingdom may enhance colonization by certain fungal families (such as pathogenic fungus) in gut microbiome [50,51]. Thus, the inter-kingdom correlations might also influence the occurrence of Chinese Cordyceps. Inter-kingdom network analysis was performed based on the positive and negative correlations among the different bacterial and fungal families in the soil microbial community. Among all the families, the co-occurrence analysis showed that the number of positive correlations decreased and negative correlations increased for sites A, B, and C, and accordingly, the ratios of positive and negative correlations decreased and were close to 1 at site C (site A: 2.71; site B: 1.54; site C: 1.16), i.e., for the null Cordyceps group, there was little difference between the number of positive and negative correlations. Thus, it can be speculated that the occurrence rate of Chinese Cordyceps might be negatively correlated with the stability of the correlation state (equilibrium between positive and negative correlations) of the soil habitat. Furthermore, the average degree of the Cordyceps group is higher than the null Cordyceps group, indicating that a higher level of inter-kingdom correlations might be helpful for the occurrence of Chinese Cordyceps. In addition, this study found that at site A, the positive correlations with the Cordyceps-related fungal family (Cordycipitaceae) were enriched. This finding indicates that some other microbes might be especially positive for the Cordyceps-related fungal family and might assist in the occurrence of Chinese Cordyceps. Further study is needed to seek out and verify the related microbial species in rhizospheres especially positive correlated with O. sinensis, to ultimately promote the occurrence of Chinese Cordyceps.
Thus, co-occurrence ecological relationships may be important for the occurrence of Chinese Cordyceps, and it is necessary to further prove the speculations proposed in this study and keep the right balance of biological interactions to increase the occurrence rate of Chinese Cordyceps in artificial cultivation.

Conclusions
This study is the first attempt to comprehensively analyze the physicochemical and microbial parameters in the habitats of Chinese Cordyceps with an emphasis on their influence on its occurrence. Soil physicochemical parameters, including EOC, SOC, HAC, HMC, and pH, might be negatively related to the occurrence of Chinese Cordyceps. Several soil physicochemical parameters (pH, NH 4 + -N/NO 3 − -N, SOC, HMC, HAC, APO, APH, and MBC) and microbial properties could interact and mix together, which might affect the occurrence of Chinese Cordyceps. Microbial community analysis revealed that a low level of bacterial and fungal diversity was suitable for the occurrence of Chinese Cordyceps, and soil microbial composition in the habitat soil varied among different sampling sites and might be closely related to the occurrence of Chinese Cordyceps (e.g., some unclassified fungi). Intra-kingdom network revealed that a closer correlation of the bacterial community might help in the occurrence of Chinese Cordyceps, while a closer correlation of the fungal community might suppress it. Inter-kingdom network revealed that the occurrence rate of Chinese Cordyceps might be negatively correlated with the stability of the correlation state (equilibrium between positive and negative correlations) of the soil habitat. Thus, our analysis shows that co-occurring ecological relationships may be important for the occurrence of Chinese Cordyceps, and it is necessary to keep the right balance of biological interactions to increase the occurrence rate in artificial cultivation. Conclusively, both soil physicochemical and microbial parameters in the habitats could be related with the occurrence of Chinese Cordyceps and our study could help in gaining a better understanding of the occurrence of Chinese Cordyceps and provide useful information for conservation and artificial cultivation of this valuable fungus-larva symbiote. However, the structure of the microbial community and the correlations in soil are complicated and vary greatly among different regions and under different stresses [52], thus it is necessary to conduct an extensive and systematic investigation of the habitat microbiome in future studies, and the definite role of the microbial community and the correlations should be confirmed by further experiments. Furthermore, in artificial-cultivated implication, the manipulation of soil microbial communities is difficult to achieve, and the present study also provides a clue that the control of physicochemical factors might aid the realization.
Supplementary Materials: Supplementary materials are available online at http://www.mdpi.com/2076-2607/7/9/ 284/s1. Table S1: OTU for the analysis of fungi and bacteria. Table S2: Result of Mantel test for the test of influence of soil physicochemical properties on soil microbial community. Table S3 and Table S4: Relative compositions of soil bacterial and fungal communities at the phylum, class, order, family, and genus level.