Wild cicada flowers in China are mainly distributed in the four provinces Anhui, Jiangsu, Sichuan, and Zhejiang. In this study, wild cicada flower and soil were collected from Anji County in Zhejiang Province, China. Our pre-survey revealed that the average density of natural cicada flower in Anji is 0.5 fruiting bodies per square meter, with the highest density being 5 fruiting bodies per square meter. Field investigations revealed that the peaks of wild cicada flowers usually occurred annually from June to August.

Wild cicada flower and soil samples were sampled by the diagonal line five-point method in mid-July 2020. During the sampling, weeds on the ground were removed first, after which the soil was cut into V-shaped pits by using a plastic shovel, and soil slices 5–10 cm deep and 10 cm wide were collected. The soil samples were divided into the cicada flower group (collected from the soil around the fruiting bodies) and the null cicada flower group (collected from a locality where no cicada flower was observed for 10 years). The samples were stored in sterile plastic bags and sent to the laboratory in an ice box. According to the morphological and anatomical characteristics, the wild cicada flowers were divided into two parts: (1) coremium, the aboveground part; and (2) sclerotium, the underground part (a cicada pupa filled with the mycelia of C. cicadae) (Chunyu et al., 2019). The processing and storage method of the collected wild cicada flowers were according to Xia et al. (2015). The fresh samples were stored at −20°C until DNA was extracted. Each soil sample was divided into two parts, one of which was used for DNA extraction and the other for analysis of soil physicochemical properties.

Thirteen physicochemical properties of soil were analyzed, including pH, moisture content (MC), total phosphorus (TP), total nitrogen (TN), available phosphorus (AP), available potassium (AK), organic matter (OM), and magnesium (Mg), aluminum (Al), copper (Cu), iron (Fe), manganese (Mn), and zinc (Zn).

The soil pH (1:2.5, soil:water) was assayed potentiometrically according to HJ 962-2018, China. The MC was determined using a gravimetric method by weighing samples before and after oven-drying at 105°C for 24 h. The determination of TP and AP content was according to LY/T 1232-2015, China. TP content was assessed using the NaOH melting method, which is the Mo–Sb anticolorimetric method; AP was determined using the sodium bicarbonate method. According to LY/T 1228-2015, China, TN was determined using the Kjeldahl method. The determination of AK was according to LY/T 1234-2015, China; AK was determined by atomic adsorption spectrophotometry after the soil was extracted with 1 mol/L of ammonium acetate, pH 7. OM in the soil was assayed using the potassium dichromate titrimetric method in an electric sand bath according to LY/T 1237-1999, China. Mg was measured by atomic absorption spectrophotometer after digested with perchloric acid and hydrofluoric acid according to NY/T 296-1995, China. The determination of Mn, Cu, and Zn was according to HJ 803-2016, China. Mn, Cu, and Zn were detected by an inductively coupled plasma mass spectrometer after the soil was digested with HCl–HNO₃ in a microwave digestion instrument. The determination of Fe and Al was according to LY/T 1253-1999, China. Fe was determined using the phenanthroline spectrophotometric method, and Al was determined using KF replacement–EDTA volumetric method.

Before DNA extraction, each wild cicada flower sample was ground into a powder in liquid nitrogen. Total genome DNA from soil and wild cicada flower sample was extracted using the HiPure Soil DNA Kit B according to the manufacturer’s protocol. The DNA concentration was monitored using a Qubit3.0 Fluorometer. The V3 and V4 hypervariable regions of the prokaryotic 16S rRNA gene were amplified using the forward primer “CCTACGGRRBGCASCAGKVRVGAAT” and reverse primer “GGACTACNVGGGTWTCTAATCC.” The ITS2 region was amplified using the forward primer “GTGAATCA TCGARTC” and reverse primer “TCCTCCGCTTATTGAT.”

PCR for bacterial-specific fragments was performed using 3 min of denaturation at 94°C; 24 cycles of 5 s at 94°C, 90 s for annealing at 57°C and 10 s for elongation at 72°C; and a final extension at 72°C for 5 min. For fungal-specific fragments, the PCR procedure included 5 min of denaturation at 94°C; 25 cycles of 30 s at 94°C, 30 s for annealing at 57°C, and 30 s for elongation at 72°C; and a final extension at 72°C for 5 min. The resulting PCR products were evaluated by 1.5% agarose gel electrophoresis. The PCR reaction mixture (25 μl) contained 2.5 μl of TransStart Buffer, 2 μl of dNTPs, 1 μl of each primer, and 20 ng of template DNA.

Next-generation sequencing (NGS) library generation and Illumina MiSeq sequencing were conducted by GENEWIZ, Inc. (Tianjin, China). The concentrations of DNA libraries were validated using a Qubit 3.0 Fluorometer. To normalize the library to 10 nM, DNA libraries were multiplexed and loaded onto an Illumina MiSeq instrument according to the manufacturer’s instructions (Illumina, San Diego, CA, United States). Sequencing was performed in PE250/300 paired-end mod; image analysis and base calling were conducted using the MiSeq Control Software (MCS) embedded in the MiSeq instrument.

The QIIME data analysis package was used for 16S rRNA and ITS rRNA data analysis. The forward and reverse reads were joined by Quantitative Insights Into Microbial Ecology (QIIME) software (Caporaso et al., 2010) and assigned to samples based on barcode and truncated by cutting off the barcode and primer sequence. Quality filtering on joined sequences was performed, and sequences that did not fulfill the following criteria were discarded: sequence length < 200 bp, no ambiguous bases, mean quality score ≥ 20. Then the sequences were compared with the reference database [Ribosomal Database Program (RDP) Gold database] using UCHIME algorithm to detect chimeric sequence, and then the chimeric sequences were removed. The effective sequences were used in the final analysis. Sequences were grouped into operational taxonomic units (OTUs) using the clustering program VSEARCH (1.9.6) against the Silva 132 database and the UNITE ITS database¹ pre-clustered at 97% sequence identity. The RDP classifier uses the Silva (Pruesse et al., 2007) 132 and the UNITE (Koljalg et al., 2005) ITS database, which has taxonomic categories predicted to the species level.

The coremium was placed into a 15-ml centrifuge tube with 6 ml of 0.9% NaCl solution and shaken to liberate the fungal spores. Aliquots comprising 200 μl of serial dilutions of each sample (10^–1, 10^–2, 10^–3, and 10^–4) were spread on plates containing rose bengal agar (RBA), and the dishes were placed in a 26°C incubator for culture. When the colony diameter reached 2–5 mm, a small amount of mycelium was removed and purified 2–4 times on the separation medium. After it was confirmed that there were no other microorganisms, the strains were transferred to potato-dextrose agar (PDA) slants, incubated at 26°C for 7 days, and then stored at 4°C.

The mycelia of the fungal strain growing on PDA were collected, and the fungal genomic DNA was extracted using the E.Z.N.A.^® Fungal DNA Mini Kit. The ITS rRNA gene was amplified by PCR using the universal primers ITS1F (5′-TCCGTAGGTGAACCTGCGG-3′) and ITS4R (5′-TCCTCCGCTTATTGATATGC-3′). The amplified products were separated by 1% agarose gel electrophoresis and sequenced by GENEWIZ, Inc. (Tianjin, China). The sequences were searched against the National Center for Biotechnology Information (NCBI) database using the Basic Local Alignment Search Tool (BLAST), and the selected sequences were aligned using CLUSTALW. Molecular Evolutionary Genetics Analysis (MEGA-X) software was used to construct the phylogenetic tree based on the neighbor-joining algorithm.

The activated strain was grown on PDA medium at 26°C for 7 days. Spores of C. cicadae were then used to inoculate the fermentation medium and cultured at 26°C and 140 rpm for 3 days. The resulting seed culture was used to inoculate 100 ml of medium (inoculum amount: 5%) and cultured at 26°C and 140 rpm for 5 days. The composition of the seed culture/fermentation medium was extracted by a previous method with some modifications and contained (per liter) the following: 20 g of glucose, 10 g of peptone, 10 g of yeast extract, 2.0 g of KH₂PO₄, 0.5 g of MgSO₄, and 0.5 g of ZnCl₂ (Chunyu et al., 2019).

The mycelia were harvested by centrifugation at 5,900 × g for 30 min and then oven-dried at 60°C. The dried mycelia, coremia, and sclerotia were ground into a powder, and 0.3 g of the powder was ultrasonically extracted with ultrapure water for 90 min. The resulting extract was cleared from debris by centrifugation at 5,900 × g for 15 min. The supernatant was filtered through a 0.22-μm pore-size Minisart^® filter for later use.

Since nucleosides are water-soluble, their presence in the water extracts of wild cicada flowers and artificially cultured C. cicadae mycelia was analyzed by high-performance liquid chromatography (HPLC) on an Agilent ZORBAX Eclipse XDB-C18 column (150 mm × 4.6 mm × 5 μm). The mobile phase consisted of 15% methanol in ultrapure water at a flow rate of 1.0 ml/min. The injection volume was 10 μl, and the detection wavelength was set at 260 nm.

The statistical significance of differences among different groups was determined using the paired-samples t-test and one-way ANOVA for repeated measures in SPSS software version 23.0 (IBM Corp., Armonk, NY, United States). Differences at the level of p < 0.05 were considered statistically significant.

The wild cicada flower mainly grows in bamboo, broad-leaved and evergreen broad-leaved forests, and mixed coniferous and broad-leaved forests. The wild cicada flower we acquired from Anji County in Zhejiang Province was characterized by high quality and yield due to the warm climate, loose soil, and high MC.

In order to disclose the effects of various elements from the soil on the formation of wild cicada flower, the physicochemical properties of soil were determined. As shown in Table 1, compared with the control soil group without cicada flower growth, the cicada flower group exhibited significantly higher levels of MC, TN, and OM, as well as significantly lower levels of pH, TP, AP, AK, Mg, Fe, Mn, Cu, and Zn (p < 0.05).

We performed high-throughput rDNA sequencing of soil and different parts of wild cicada flower, including the aboveground part (coremium) and the underground part (sclerotium), to clarify the external and internal microbial communities of cicada flower and to explore the possible interactions between them. In order to explore the diversities of microorganisms, first, alpha-diversity estimators were used to reveal the richness and diversity of species in the community of the coremia, sclerotia, and soil (Table 2). Then, principal component analysis (PCA) plots demonstrated the marked differences in the composition of the coremia, sclerotia, and soil.

Good’s coverage of samples (>0.99 for all samples) and the rarefaction curves (Supplementary Figure 1) indicated that the majority of microbes were represented in the data and that the sampling depth of sequence dataset in the current study was adequate for our analyses. After 72,875 and 55,554 chimeric sequences were filtered, a total of 746,944 high-quality 16S rDNA sequences and 1,196,194 high-quality ITS rDNA sequences were obtained. With a 97% identity threshold, these high-quality sequences were clustered into 200 bacterial and 346 fungal OTUs. The shared and unique OTUs among the soil, coremia, and sclerotia (Figure 1) were depicted using Venn diagrams drawn by the R software based on OTU. Our results indicated that 100 bacterial and 42 fungal OTUs are shared among the soil, coremia, and sclerotia as shown in Figure 1. Aside from that, 57 bacterial OTUs were only found in the soil, while only 1 and 2 OTUs were exclusive to the coremia and sclerotia, respectively. Among fungi, 196 OTUs were only found in the soil, and 11 were only found in the coremia. By contrast, no OTUs were exclusively found in the sclerotia. The detailed taxonomic information of the OTUs is listed in Supplementary Tables 1–4. In summary, the number of total and unique fungal OTUs in the soil was higher than in the coremia or sclerotia.

In order to reveal the richness and diversity of species in the communities among the soil, coremia, and sclerotia, alpha diversity was analyzed using four indices, i.e., ACE, Chao1, Shannon, and Simpson, and was calculated using QIIME. For bacteria, the richness (represented by Chao1 and ACE) and diversity (represented by Shannon and Simpson indices) of communities in the soil was significantly higher than those of the communities in the coremia and sclerotia (p < 0.05). However, there was no significant difference between the coremia and sclerotia (p > 0.05). For fungi, the richness and diversity of communities in the soil were also significantly higher than those in the coremia and sclerotia (p < 0.05), but the fungal richness of the coremia was significantly higher than that of the sclerotia (p < 0.05), and there was no significant difference in fungal diversity between the coremia and sclerotia (p > 0.05). In conclusion, the richness and diversity of bacterial communities in the coremia and sclerotia were very similar, but the richness of fungal communities in the coremia was significantly higher than in the sclerotia. Notably, C. cicadae was dominant in the sclerotia, which may explain their low fungal diversity and richness.

To infer differences of microbial composition between samples, Euclidean distances were calculated for all samples and analyzed in R using the “vegan” package. The results were visualized using PCA. Our results showed that the bacterial (Figure 2A) and fungal communities (Figure 2B) in the soil samples were different from those in the coremia and sclerotia of cicada flower. However, the microbial community in the coremia samples appeared more similar to that of sclerotia samples.

In conclusion, the external microbial environment of cicada flower was richer and more diverse than the internal microbial environment, possibly following the reason that a portion of the microbes was unable to colonize in the internal environment of cicada flowers, which might be inhibited or killed by the fungus during its invasion into cicada nymphs and larvae. The similar microbial compositions of the coremia and sclerotia indicated that some bacteria from the sclerotia might be brought into the coremia during the growth of mycelia within C. cicadae, which was consistent with earlier findings that fungal mycelia can promote the short- or long-distance migration of bacteria in the soil (Yang and van Elsas, 2018).

The microbial composition of different samples was analyzed at the phylum and genus levels. Figure 3 intuitively illustrates the top 30 phyla and top 30 genera of the bacterial and fungal communities. The relative compositions of bacterial and fungal communities at the class, order, and family levels are presented in Supplementary Figures 2–7, respectively. The dominant phyla or genera were defined as the phyla or genera whose relative abundance accounted for more than 1% in the soil, coremia, or sclerotia. At the phylum level, the structure of the bacterial community displayed the same patterns between the coremia and sclerotia (Figure 3A). Proteobacteria, Bacteroidetes, and Firmicutes were the dominant phyla in both the coremia and sclerotia. The relative abundance of Proteobacteria in the coremia and sclerotia was 62.89 and 87.08%, respectively, followed by Bacteroidetes with a relative abundance of 34.23% in the coremia and 9.71% in the sclerotia. Firmicutes accounted for 2.28% in the coremia and 1.41% in the sclerotia. By contrast, the bacterial structure of the soil exhibited significantly different patterns. The dominant phyla in the soil samples were Acidobacteria (50.48%), Proteobacteria (29.01%), Firmicutes (10.20%), Actinobacteria (5.28%), and WPS-2 (2.80%). The fungal community of the coremia, sclerotia, and soil (Figure 3B) was dominated by Ascomycota (sclerotia, 99.97%; coremia, 97.62%, and soil, 44.84%). In addition to Ascomycota, the dominant phyla were mainly Basidiomycota (2.11%) in the coremia, and Basidiomycota (41.77%), Mucoromycota (7.80%), and Mortierellomycota (3.92%) in the soil.

It is worth noting that due to the vast number of undiscovered microorganisms and the limited resolution of the rRNA gene and the ITS region sequence, “Unclassified” made up a sizeable proportion of the OTUs at the genus level. In the bacterial community (Figure 3C), the most abundant genera detected in the coremia were Pedobacter (26.35%), f__Enterobacteriaceae_Unclassified (12.81%), Pandoraea (8.29%), Allorhizobium–Neorhizobium–Pararhizobium–Rhizobium (7.49%), Phyllobacterium (7.17%), Chitinophaga (7.12%), and Stenotrophomonas (3.63%). The most abundant genera detected in the sclerotia were Achromobacter (14.92%), f__Enterobacteriaceae_Unclassified (12.85%), Stenotrophomonas (12.51%), Burkholderia–Caballeronia–Paraburkholderia (11.39%), Allorhizobium–Neorhizobium–Pararhizobium–Rhizobium (9.19%), Serratia (6.11%), Pedobacter (5.43%), Comamonas (5.09%), Pandoraea (3.78%), f__Mitochondria_Unclassified (3.57%), and Chitinophaga (3.42%). In the soil samples, Burkholderia–Caballeronia–Paraburkholderia (5.21%), Bacillus (5.10%), Acidibacter (4.11%), f__Xanthobacteraceae_Unclassified (4.05%), and Candidatus_Solibacter (3.83%) were the most abundant genera. In the fungal community (Figure 3D), Isaria and C. cicadae were unsurprisingly dominant, accounting for 83.48, 53.57, and 16.50% in the sclerotia, coremia, and soil, respectively. Other major genera in soil included f__Clavariaceae_Unclassified (13.96%), Umbelopsis (6.18%), f__Chaetomiaceae_Unclassified (4.03%), Mortierella (3.80%), f__Sordariaceae_Unclassified (3.58%), and Arcopilus (3.42%). Other major genera in the coremia included Apiospora (5.55%) and Cladosporium (3.99%). Fusarium was also highly abundant in the sclerotia, reaching 15.82%.

In order to elucidate the interactions between environmental factors and the soil microbial community, we assessed Pearson’s correlation between soil physicochemical parameters and all phyla, the top 10 abundant genera, all OTUs, and the alpha-diversity indices (Figure 4).

For bacteria (Figure 4A), Chao1, Shannon, and Simpson were significantly negatively correlated with Mg, while ACE was significantly negative correlated with Cu (p < 0.01). All bacterial OTUs were significantly negatively correlated with Mg (p < 0.05). The phyla Acidobacteria, Chloroflexi, and Patescibacteria were significantly negatively correlated with Cu (p < 0.05), while Verrucomicrobia and Patescibacteria were significantly negatively correlated with Mg (p < 0.05). Conversely, Firmicutes and Actinobacteria were significantly positively correlated with Mg (all p < 0.05 at least), and Actinobacteria were significantly positively correlated with Cu (p < 0.05). The dominant genera Acidibacter, Candidatus Solibacter, Occallatibacter, Bryobacter, and Roseiarcus were negatively correlated with most environmental factors. Among them, Acidibacter, Candidatus Solibacter, and Bryobacter were significantly negatively correlated with Cu (all p < 0.05 at least), while Candidatus Solibacter and Occallatibacter were significantly negatively correlated with Zn (all p < 0.05 at least). Additionally, Burkholderia–Caballeronia–Paraburkholderia, Bacillus, Tumebacillus, Candidatus Koribacter, and Streptacidiphilus were negatively correlated with most environmental factors. Among them, Burkholderia–Caballeronia–Paraburkholderia, Bacillus, and Tumebacillus were significantly positively correlated with Mg (all p < 0.05 at least), Candidatus Koribacter was significantly positively correlated with Fe and MC (p < 0.05), and Streptacidiphilus was significantly positively correlated with Cu (p < 0.05).

For fungi (Figure 4B), Chao1 and ACE were significantly positively correlated with TP (p < 0.01), while Shannon and Simpson were significantly positively correlated with Zn (p < 0.01) and Fe (p < 0.05). All fungal OTUs were significantly positively correlated with TP (p < 0.05). Ascomycota, the most abundant phylum, was positively correlated with TN, OM, Zn, Fe, MC, Mg, Cu, pH, and TP. By contrast, most of the other phyla presented negative correlations with most environmental factors. Among them, Basidiomycota, Chytridiomycota, and Kickxellomycotina were significantly negatively correlated with Mg, TN, and TP, respectively (all p < 0.05 at least). Mucoromycota and Rozellomycota were significantly positively correlated with Mg (p < 0.001). At the genus level, the dominant genus Isaria was significantly negatively correlated with Fe, MC, and Zn (all p < 0.05 at least). By contrast, there was a significant positive correlation between most other dominant genera and Mg (all p < 0.01 at least).

According to the rDNA sequencing results, Isaria was the most dominant fungal genus in all samples, including the coremia, sclerotia, and soil; and this fungus also produces the medicinally valuable fruiting bodies known in China as cicada flowers. Based on this, isolating a high yield and high activity strain of C. cicadae is a critical step for the large-scale artificial cultivation of this medicinal ingredient. Therefore, we isolated a fungal strain via a culture-dependent method (Figure 5). In order to identify this fungal strain, we amplified the ITS gene sequences by PCR and sequenced the obtained amplicons, and the results were aligned with sequence data in the NCBI database using the BLAST algorithm.

The phylogenetic tree was constructed using Cordyceps-related sequences in GenBank. The genetic distance between the isolated strain AH10-4 and C. cicadae strain CCF (GenBank accession: KP771879) was the closest, with the rDNA sequence identity reaching 100%. After the sequence data were combined with morphological characteristics used to determine the taxonomic status of wild strain, AH10-4 (GenBank accession: MZ676080) was identified as C. cicadae (Figure 6).

Nucleosides are the major bioactive constituents of Cordyceps and as such play important roles in the regulation and modulation of various physiological processes in the body (Meng et al., 2019). In order to test the fermentation ability of the isolated C. cicadae strain AH10-4, the content of nucleosides in the fermented mycelia was quantified and compared with that of the coremia and sclerotia of wild cicada flower using HPLC.

The established analytical method was applied to simultaneously determine the content of adenosine, HEA, and cordycepin (Figure 7A), the representative nucleosides of cicada flower. For wild cicada flower, the adenosine content was 0.224 ± 0.005–0.518 ± 0.022 mg⋅g^–1 in the coremia and 0.317 ± 0.008–0.927 ± 0.037 mg⋅g^–1 in the sclerotia. The content of HEA was 0.268 ± 0.004–1.034 ± 0.024 mg⋅g^–1 in the coremia and 0.190 ± 0.003–0.691 ± 0.015 mg⋅g^–1 in the sclerotia (Figures 7B,C). The content of adenosine and HEA in the coremia and sclerotia of wild cicada flowers are shown in Table 3. Compared with that of the wild cicada flower, the average content of adenosine and HEA in the artificially cultivated mycelia increased greatly and, respectively, reached 3.042 ± 0.235 and 2.357 ± 0.207 mg⋅g^–1 in AH10-4 mycelia (Figure 7D). However, cordycepin was not detected in any of the samples studied here.

Adenosine and HEA, the most important active substances of wild cicada flower, were also found in cultivated mycelia of C. cicadae, and their content was significantly higher than in wild fruiting bodies (p < 0.05). Notably, the C. cicadae strain AH10-4 isolated in this study exhibited higher adenosine and HEA productivity as compared with relative researches (Chunyu et al., 2019; Ke and Lee, 2019; Qu et al., 2019), indicating that it may be a viable substitute for the rare and expensive fruiting bodies found in nature.