Meta-Analysis and Evaluation by Insect-Mediated Baiting Reveal Different Patterns of Hypocrealean Entomopathogenic Fungi in the Soils From Two Regions of China

A survey was carried out on forest soils and grassland soils from Hebei and Sichuan provinces using Tenebrio molitor larvae as a bait, and high-throughput DNA sequencing (HTS) of the fungal internal transcribed spacer-2 ribosomal DNA was used to monitor the natural distribution of three leading hypocrealean families of insect fungal pathogens (Clavicipitaceae, Cordycipitaceae, and Ophiocordycipitaceae). The occurrence of insect fungal pathogens in soil samples from 98 different sites was compared. The use of insect bait indicated that entomopathogenic fungi of the genus Metarhizium were predominant, followed by Beauveria and Isaria. Molecular characterization using the Mz_FG543 intergenic region revealed that the Metarhizium species pool was phylogenetically composed of three closely related species as follows; Metarhizium pingshaense (n = 74), Metarhizium robertsii (n = 51), and Metarhizium brunneum (n = 26), as well as one isolate which clustered with Metarhizium flavoviride. Nine potentially new phylogenetic species were delimited within the M. flavoviride species complex by sequencing of the 5′ elongation factor-1 alpha region locus. The Beauveria (n = 64) and Isaria (n = 5) isolates were characterized via sequence analyses of the Bloc region. An intergenic spacer phylogeny of the Beauveria isolate assemblage revealed the phylogenetic species within the Beauveria bassiana clade. Interestingly, the individuals of M. pingshaense (n = 18) and M. brunneum (n = 12) exhibited the presence of both mating types in Sichuan Province. Similarly, for the Beauveria isolates, reproductive mode assays demonstrated that all four B. bassiana subclades possessed bipolar outcrossing mating systems. Of these, 19 isolates contained two mating types, and the rest were fixed for single mating types, revealing opportunities for intra-lineage heterothallic mating. The HTS results showed a significantly higher occurrence of the Clavicipitaceae family and the Metarhizium genus in the soil samples. The Venn diagram showed Metarhizium anisopliae (senso lato), Isaria farinose, and B. bassiana as frequently abundant fungal pathogen operational taxonomic units (core) across sampling sites, while the baiting method showed that the genus of Isaria was isolated locally. The Mantel test verified that community dissimilarity increased significantly with geographical distance, suggesting that geographical coordinates are possible factors that influence the insect fungal pathogen community composition in the studied sites. This study is the first to highlight the usefulness of utilizing soil baiting and deep sequencing to investigate the population dynamics of entomopathogens in soil.

A survey was carried out on forest soils and grassland soils from Hebei and Sichuan provinces using Tenebrio molitor larvae as a bait, and high-throughput DNA sequencing (HTS) of the fungal internal transcribed spacer-2 ribosomal DNA was used to monitor the natural distribution of three leading hypocrealean families of insect fungal pathogens (Clavicipitaceae, Cordycipitaceae, and Ophiocordycipitaceae). The occurrence of insect fungal pathogens in soil samples from 98 different sites was compared. The use of insect bait indicated that entomopathogenic fungi of the genus Metarhizium were predominant, followed by Beauveria and Isaria. Molecular characterization using the Mz_FG543 intergenic region revealed that the Metarhizium species pool was phylogenetically composed of three closely related species as follows; Metarhizium pingshaense (n = 74), Metarhizium robertsii (n = 51), and Metarhizium brunneum (n = 26), as well as one isolate which clustered with Metarhizium flavoviride. Nine potentially new phylogenetic species were delimited within the M. flavoviride species complex by sequencing of the 5 elongation factor-1 alpha region locus. The Beauveria (n = 64) and Isaria (n = 5) isolates were characterized via sequence analyses of the Bloc region. An intergenic spacer phylogeny of the Beauveria isolate assemblage revealed the phylogenetic species within the Beauveria bassiana clade. Interestingly, the individuals of M. pingshaense (n = 18) and M. brunneum (n = 12) exhibited the presence of both mating types in Sichuan Province. Similarly, for the Beauveria isolates, reproductive mode assays demonstrated that all four B. bassiana subclades possessed bipolar outcrossing mating systems. Of these, 19 isolates contained two mating types, and the rest were fixed for single mating types, revealing opportunities for intra-lineage heterothallic mating. The HTS results showed a significantly higher occurrence of the Clavicipitaceae family and the Metarhizium genus in the soil samples. The Venn diagram showed Metarhizium anisopliae (senso lato), Isaria farinose, and B. bassiana as frequently abundant fungal pathogen operational taxonomic units (core) across sampling sites, while the baiting method showed that the genus of Isaria was isolated locally. The Mantel test verified that

INTRODUCTION
Entomopathogenic fungi (EFs) are a group of fungi that are pathogens of insect species. They include a broad array of diverse fungal species based on morphological characterization, phylogenetic evaluation, and ecological distribution. EFs can be found in five of the eight fungal phyla (Araújo and Hughes, 2016). Ascomycota, with approximately 64,000 known species, forms the largest group in Kingdom Fungi (Kirk et al., 2008). Its species (mostly Hypocreales) infect 13 orders of insects, the most globally diverse and widespread animals (Boomsma et al., 2014;Araújo and Hughes, 2016). The Hypocreales fungi encompass important genera of soilborne fungal pathogen taxa, and the most commonly studied genera of invertebrate pathogens are distributed among three families: Clavicipitaceae (e.g., Metarhizium and Pochonia), Cordycipitaceae (e.g., Cordyceps, Beauveria, and Isaria), and Ophiocordycipitaceae (e.g., Ophiocordyceps and Tolypocladium, formerly Elaphocordyceps) . Communities of EFs exist in many different environments Keyser et al., 2015), and an explicit and comprehensive explanation of their structure is beneficial to address fundamental biological questions of utilizing them as microbial control agents (Behie and Bidochka, 2013). The most widely used approach for the isolation of EFs in soils is baiting with susceptible insects (Zimmermann, 1986). The wax moth larvae, Galleria mellonella (Lepidoptera: Pyralidae) (Vänninen, 1996;Chandler et al., 1997;Bidochka et al., 1998;Klingen et al., 2002;Quesada-Moraga et al., 2007;Sun et al., 2008;Goble et al., 2010;Meyling et al., 2011;Hernández-Domínguez and Guzmán-Franco, 2017;Inglis et al., 2019), and mealworm, Tenebrio molitor (Coleoptera: Tenebrionidae) (Meyling and Eilenberg, 2006;Steinwender et al., 2014;Keyser et al., 2015;Korosi et al., 2019), are the most used bait species for the general sampling of fungal entomopathogens. The outcomes of the above-mentioned surveys of naturally occurring taxa of insect-associated fungi in soil habitats are presented in Table 1. Although most studies are from subtropical or temperate areas (Hajek and Shapiro-Ilan, 2018), baiting with insects has been applied to obtain hypocrealean insect pathogens from soil samples in the tropical regions (Rezende et al., 2015) and the Arctic (Meyling et al., 2012). If part of the pathogen population is temporarily non-infective and/or non-active against the insect host at the time that the baiting occurs, the frequency of pathogens can be underestimated by this method. This issue is raised as one of the apparent biases of this detection technique (Hajek and Shapiro-Ilan, 2018). In contrast, Vega and Kaya (2012) determined this point as the main advantage of baiting on soil samples because the insect baiting method selectively isolates those pathogens which are biologically active. Some fungal species contain a narrow host range. For instance, some Metarhizium species such as Metarhizium album and Metarhizium acridum possess narrow-host-range pathogens of arthropods, and a standard bait method cannot detect them (Wang et al., 2016). This phenomenon might be a reasonable explanation for those studies that used different host insects as bait, resulting in some fungal species that are rarely baited with common host insects (Klingen et al., 2002;Goble et al., 2010). Furthermore, a pathogen should outcompete with a variety of other organisms and/or free-living pathogens for insects. This competition may affect the outcome of the interaction between pathogens and whole insect cadavers produced per sample. Thus, baiting with insects might underestimate the distribution patterns of soilborne insect-associated pathogens in diversity studies (Vega and Kaya, 2012;Hajek and Shapiro-Ilan, 2018). Fungal pathogen distribution patterns can be affected by many environmental factors, including abiotic and biotic factors (Jaronski, 2007(Jaronski, , 2010Hajek and Shapiro-Ilan, 2018), and it is difficult to draw a picture at the global scale or design general predictable patterns for EF distribution in the ecosystems. On the one hand, EF natural occurrence studies rely mostly on using the baiting method to identify insectassociated fungal pathogens on soil samples. On the other hand, various molecular techniques have been incorporated in the ecological study of invertebrate pathogens and the diseases that they cause (Hajek and Shapiro-Ilan, 2018). One of these advanced techniques is high-throughput sequencing (HTS), also named next-generation sequencing. This technique is becoming the go-to avenue for scientists studying invertebrate pathogen ecology (Malacrinò et al., 2017;Hajek and Shapiro-Ilan, 2018), especially for monitoring particular fungal species in ecosystems (Hirsch et al., 2013).
The current research aims to discover the natural occurrence and the diversity of predominant insect-associated fungal pathogens by comparing the insect-susceptible bait method (T. molitor larvae) with HTS, a technique that has surprisingly not been previously examined to study the diversity of hypocrealean EFs broadly. The most abundant genus obtained by the baiting method was Metarhizium, followed by Beauveria. Molecular characterization was conducted to distinguish phylogenetic subgroups within fungal isolates, which were obtained by using the baiting method. HTS was conducted so as to target the fungal entomopathogen communities of the soil to show their taxonomic profile and their relative abundance in comparison to isolates obtained through insect baiting. After that, we turn our attention on screening all obtained Metarhizium and Beauveria using diagnostic mating-type PCR assays.

Ethics Statement
There was not any specific authorization required for the field sampling. The collection areas were not privately owned or protected in any way. Soil sampling was conducted from forest and grassland sites and did not include plant material from imperiled or protected species.

Study Site and Soil Sampling
We collected soil samples at the locations shown on the map (Figure 1), and the soil sampling sites were assigned coordinates via Global Positioning System (GPS). Two main geographical areas were studied (Hebei and Sichuan provinces). We sampled during the peak vegetation growing season (July-August) in 2018 to target maximum microbial activity. Soil collections were conducted from Hebei Province's Taihang Mountains in close proximity to Shijiazhuang City (23 soil samples), from now on referred to as SH (38 • 13 .406 N, 113 • 36 .265 E) and from Saihanba National Forest Park (42 • 14 .084 N, 117 • 08'.124 E) (from now on referred to as FIGURE 1 | Map showing the geographical sites in China (A) at the Saihanba National Forest Park (B), in the Taihang Mountains, and (C) in the area throughout Sichuan to Hubei provinces (D). At each sampled locality, five lines, with a distance of 100 m between each line pair, were selected, and five soil sub-samples, within 50 m 2 with approximately 10 m distance between each sub-sample, were collected (E). A map of China and the sampling sites were generated via Google Earth Pro ver. 7.1.7.2600 (Google Inc.) for Windows. The sampling procedure design and the editing of the photos were performed using Xara Designer Pro X ver. 15.1.0.53605 (PANTONE R LLC 2016, Xara Group Ltd.) for Windows. SNFP). The Taihang Mountains are located on the eastern side of China's secondary step zone of the Bohai Bay basin and draw a physical boundary, and these also represent an important geographic barrier for the Loess Plateau and the north China meadows. The average elevation is 1,000-2,000 m. The highest elevation (2,882 m) is in the northern part of the mountain at Xiaowutai. The annual average temperature is around 10 • C, and the annual rainfall is approximately 600 mm. SNFP expands across 185,000 acres of forest in Chengde. It has an average elevation of 1,400 m and two important specialties of forest and grassland as well as more than 10 plain lakes. The original ecosystem is very well preserved and maintained. It is large and the world's biggest human-made forest park zone; its northern edge borders China's Inner Mongolia autonomous area. The plant canopy is composed of trees, including conifer, Scots pine, birch, larch, and spruce. The sampling at SNFP was conducted at three ecotypes: forest (from now on referred to as F, with 10 samples), grassland (from now on referred to as G, with 10 samples), and the boundaries of forest and grassland (from now on referred to as B, with eight samples). Because of unfavorable weather conditions during sampling, the number of boundary samples was eventually decreased to eight (Figure 1). The second sampling was conducted in the area scattered across a gradient stretching 1,350 km west to east, in southwestern China, situated in both Sichuan and Hubei provinces, a famous biodiversity hotspot known in China and worldwide (Güneralp et al., 2015). In this area, the dominant vegetation canopy is conifer trees such as larch (Larix spp.), fir (Abies spp.), spruce (Picea spp.), and hemlock (Tsuga spp.). The climate is temperate, with an annual precipitation of 800-1,000 mm, with the most rainfall occurring between May and September, and a mean annual temperate of 6-10 • C 1 .
In this region, sampling was performed to provide a comprehensive analysis of insect fungal pathogen diversity at the 47 collection sites along a 1,350-km stretch from Kangding (29 • 49 .516 N, 102 • 54 .276 E, elevation 3,500 m above sea level) to Ping Ke Xian (31 • 06 .962 N, 113 • 15 .18 E, elevation 100 m above sea level). We collected 14 samples from Kangding to Ya'an region (from now on referred to as Y) and seven samples from the Meishan region (from now on referred to as M). Sampling continued from Meishan to Ping Ke Xian (from now on referred to as K, with 26 samples). A sampling at Y and M collection sites was conducted by selecting sites that were 10 km apart from each other. The sampling sites at K location were located 50 km apart from each other (1,350 km). The leading information of the sampling sites is summarized in Supplementary Table S1. For each location, soil sampling was conducted with the identification of five lines, with the distance between each line at 100 m approximately. From each of the five lines, five soil sub-samples at 10 m distance of each site were taken from topsoil to 30-50 cm in the subsoil with a diameter of 30 cm by stepping nearly 50 m from the start point toward the end of the line. The samples were collected using clean shovels and then placing the five soil sub-samples into a plastic bowl to homogenize independently, resulting in one soil sample. As much as possible, the large roots and gravel larger than 0.25 in. were removed. The shovels were cleaned (rinsed with water and wiped dry) for each location. Then, subsamples of the homogenized material were placed in a pre-labeled ziplock plastic bag (14 × 21 cm) separately, and the five soil samples of each site were held in a large, sealable polyethylene bag and placed in a cooler with freezer packs. The schematic soil sampling is presented in Figure 1E. In total, 490 (98 × 5) soil samples were obtained from a diverse array of sites (98 locations) representative of forest and grassland habitats. After transporting the soil samples to the lab, the samples were divided into two parts. Subsequently, 30 g of each soil sub-sample (5 × 30 g) was added to a sterile plastic bag to mix with 150 g soil manually, resulting in one composite sample, and it was immediately stored at −80 • C for the conduct of a DNA extraction procedure from the soil. The second part was air-dried and ground to pass through a 2-mm-diameter sieve within 24 h to employ in the baiting with the insect.

Soil DNA Extraction and Amplicon-Based Sequencing
Total genomic DNA was extracted in two replicates from 0.25 to 0.5 g of each composite soil sample (98 soil samples) using DNeasy R PowerSoil R Kit (QIAGEN, Hilden, Germany) to obtain genomic DNA free of PCR-inhibitory phenolic compounds following the manufacturer's instructions. The two DNA extractions per soil sample were pooled before further analysis. Genomic DNA purity and integrity were determined by 1% agarose electrophoresis and nanodrop (Nanodrop 2000, Thermo Scientific, Wilmington, NC, United States). DNA was diluted to 1 http://data.cma.cn/data/online/t/1 1 ng/µl using sterile Mili-Q ddH 2 O, based on the concentration. The amplification of the DNA samples was performed using ITS3-2024F (5 -GCATCGATGAAGAACGCAGC-3 ) and ITS4-2409R (5 -TCCTCCGCTTATTGATATGC-3 ) primers, which is target of 350 bp of the fungal rRNA gene internal transcribed spacer-2 (ITS2) region (Orgiazzi et al., 2012). PCR amplifications were conducted in a 30-µl mixture containing 0.5 µl of bovine serum albumin (TaKaRa, Tokyo, Japan), 2.0 µl of forward and reverse primers (10 mM), 5.0 µl of the DNA sample, 7.5 µl of ddH 2 O, and 15.0 µl of Phusion R High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, United States) with HF Buffer (New England Biolabs, Ipswich, MA, United States). PCR amplification was initiated by 98 • C for 10 s, then 30 cycles of denaturation for 10 s at 98 • C, annealing for 30 s at 50 • C, elongation for 30 s at 72 • C, and final extension for 5 min at 72 • C. Subsequently, the PCR products were evaluated by examining 5 µl of the product on 2% agarose gels. The PCR assays were performed in SimpliAmp TM Thermal Cycler (Life Technologies Holdings Pte., Ltd., Marsiling Industrial Estate RD3, Singapore). The samples with an intense band were purified with GeneJET TM Gel Extraction Kit (Thermo Scientific, Waltman, MA, United States). Sequencing libraries were produced using Ion Plus Fragment Library Kit 48 reactions (Thermo Scientific, South San Francisco, CA, United States) according to the manufacturer's instructions, and index codes were added. The library quality was assessed on the Qubit 2.0 Fluorometer (Thermo Scientific, United States) and Agilent Bio-analyzer 2100 system (Agilent, Santa Clara, CA, United States). The clustering of the index-coded samples was performed on a cBot TM 2 Cluster Generation System according to the manufacturer's instructions. Having generated the cluster, the prepared libraries were sequenced on an Ion S5 TM XL SE400/SE600 2 at Novogene Bioinformatics Technology Co., Ltd. (Tianjin, China), and singleend reads were generated. The raw sequences produced in the current study can be accessed through the National Center for Biotechnology (NCBI) 3 under BioProject ID PRJNA551928, accession SUB5873947.

Operational Taxonomic Unit Cluster and Diversity Analysis
To obtain high-quality clean reads and to improve the accuracy of the contig lengths of the resulting assemblies, quality filtering on the raw tags (splicing sequences) was carried out following particular filtering circumstances that were performed using the Cutadapt ver. 1.9.1 (Martin, 2011) with default settings 4 . The comparison between sequence reads with the reference database (UNITE database 5 ) (Abarenkov et al., 2010) was performed to distinguish and to remove chimeric sequences using UCHIME algorithm (Edgar et al., 2011). The UNITE database ver. 6.0 for QIIME TM ver. 1.9.1 software pipeline 6 (Caporaso et al., 2010) was applied as a source file for operational taxonomic unit (OTU) picking and assigning the classification for the three hypocrealean fungal families, and the rest of the OTUs were deleted (Caporaso et al., 2010;Kõljalg et al., 2013). The sequence from each OTU with the highest frequency was then selected as an illustrative sequence of that OTU. The clusters that comprised two sequences, one sequence was randomly picked, and single read clusters were excluded from the data. The clusters were then taken together into a matrix with cluster absolute and relative abundance versus sample identities. Each OTU was taxonomically identified using UPARSE software ver. 7.0.1001 (Edgar, 2013) via a combined comparison with the curated database UNITE (Abarenkov et al., 2010) and with a batch BLAST search against GenBank (NCBI) (Altschul et al., 1990). The ITS2 homology for defining taxa in GenBank was fixed to at least 97-100% and e-values below e −100 for species, 90-97% and e −90 for genus level, and 80-90% and e −80 at the family level, corresponding to at least 90% of the sequence length (Ottosson et al., 2015).
Taxonomic data were assigned to each representative sequence using the QIIME TM software after de-noising the reads using the QIIME TM software (Reeder and Knight, 2010). Multisequence alignments were performed using MUSCLE software ver. 3.8.31 (Edgar, 2004) to examine the differences of the dominant species in different groups and finally to construct the phylogenetic relationship of different OTUs. The frequency of OTUs was normalized using a sequence standard based on the sample with the fewest sequences. Also, we visualized relative abundance and taxonomic hierarchy using Krona radial space filling (Ondov et al., 2011). The order, family, and genus ranks of each sample were selected for representation in the Krona plots. Since EFs were considered in this study, only Hypocreales order was chosen to generate these plots. The less abundant taxa are listed outside the charts along with their relative abundance. Four indices of alpha diversity, including Chao1 (Chao, 1984), Shannon (1948); Simpson (1949), and Faith's phylogenetic diversity (Faith, 1992), were analyzed based on the complexity of species diversity via observed OTUs for a sample using QIIME TM . To evaluate the differences of samples in species complexity, beta diversity on both weighted and unweighted UniFrac (Lozupone et al., 2007) was determined by QIIME TM . The results were visualized in phylogenetic trees and box plots. The abundance of distinguished OTU with high frequency (>0.1%) between sample locations was delimited using t-test and Kruskal-Wallis test (Kruskal and Wallis, 1952). To recognize the taxa with significant differential abundance, a linear discrimination analysis effect size (LEfSe) method was computed, where different sampling locations were assigned as comparison classes. The non-parametric factorial Kruskal-Wallis sum-rank test was applied to analyze all LEfSe-identified features that were significantly different among sampling sites [P < 0.05, linear discriminant analysis (LDA) ≥ 4] (Segata et al., 2011). Furthermore, to evaluate whether the beta diversity of insect fungal pathogen communities was related to the coordinates of the sampling locations, we performed a distancebased redundancy analysis (dbRDA) [Legendre and Anderson (1999) in the package Vegan ver. 2.3.5 (Oksanen et al., 2013) in R ver. 3.2.1 (R Core Team, 2016)]. We used a permutational analysis of variance to test for the significance of terms in the dbRDA. The Mantel test was applied with a significance testing against 999 permutations in R using Vegan package ver. 2.3.5 (Oksanen et al., 2013). Non-metric multidimensional scaling (NMDS) was applied to visualize community differences. Analysis of similarity (ANOSIM), a non-parametric multivariate analysis, was performed to test the significance of the patterns observed by NMDS. The distance matrix for this analysis was calculated using the Bray-Curtis similarity index. ANOSIM was performed using the Vegan package. The effectiveness of sampling was evaluated with a rarefaction analysis of data subsets using the specaccum function (Oksanen et al., 2013) in R (Supplementary Figure S1).

Fungal Isolation
The insect-susceptible bait method applied in the current study was adjusted from Zimmermann (1986) and Keyser et al. (2015). Enough ddH 2 O was added to the prepared soil samples to moisten them as described before. Approximately 200 g of soil was transferred in a polyethylene container (500 ml), omitting 15 cm of airspace at the top. Ten healthy 5-6 th instar larvae (1-2 cm) of T. molitor L. (Coleoptera: Tenebrionidae) were then transferred to each plastic container. A ventilated lid with oxygenating holes (1-2 mm in diameter) was added on top of each plastic container, and the containers were turned upside down to busy the larvae in a container. Every day the containers were rotated so that the orientation of the cups shifted, driving the T. molitor larvae to move within the soil substrate and enhancing the probability that they would come in touch with insect-associated fungi in the soil. Insect survival was monitored weekly; dead insects were removed, washed with ddH 2 O, and placed in a sterile Petri plate (35 × 10 mm) with a moistened filter paper to determine mycosis. We should mention that, due to the enormous number of soil samples, surface sterilization of the T. molitor larvae did not apply to the fungal isolation procedure (Meyling, 2007). Fungal isolations were made using a sterile toothpick from insect cadavers with mycosis and plated on PDAY media plus 0.6 g streptomycin sulfate (MP Biomedicals, Shanghai, China) and 0.5 g chloramphenicol (MP Biomedicals, Shanghai, China). After 3 weeks of growth, the colonies were evaluated morphologically and diagnosed to the genus level based on the conidial shape, phialides, and hyphae (Humber, 2012). A single-spore isolation method was applied to purify all the obtained fungal isolates using a continuous dilution of the conidia suspension of each isolate (Choi et al., 1999). We considered each infected cadaver as one isolate (Korosi et al., 2019). Finally, a small piece of agar plug of each fungal monosporic isolate was placed into a sterile 2-ml Eppendorf tube containing 20% sterilized glycerol to keep them at −80 • C for further studies. The natural occurrences of the obtained fungal isolates were assessed for homogeneity among all collection sites using Pearson's χ 2 test in GraphPad Prism ver. 8.0.2 (GraphPad Software, Inc.) 7 .
Fungal Genomic DNA Extraction, PCR Amplification, and Sequencing DNA extraction was performed by inoculating the fungal conidia (10 5 conidia/ml) into the wells of a 96-well tissue culture plate (Thermo Fisher Scientific, United States). Each well contained 1 g yeast extract (OXOID, Basingstoke, Hampshire, United Kingdom), 0.5 ml of 1 g peptone (AOBOX, Beijing, China), and 4 g D-glucose anhydrous (Solarbio, Beijing, China) (all weights per liter). The plates were incubated at 27 • C in an incubator for 72 h. The fungal hyphae were collected from the wells by centrifuging 0.5 ml of the re-suspended hyphae in a sterilized micro-centrifuge tube The fungal hyphae were lyophilized using Virtis BenchtopPro Freeze Dryer with Omnitronics TM (Genevac Ltd., Ipswich, United Kingdom). Total genomic DNA for all fungal isolates was extracted by pulverizing lyophilized hyphae with lysis buffer and 0.5 g of 0.5-1-mmdiameter acid-washed, sterilized glass beads (Tiantai Jingong SiLi Glass Beads Co., Ltd., Tiantai, Zhejiang, China) by using Bullet Blender R 24 Gold Tissue Homogenizer (Next Advance, Inc., NY, United States). The rest of the DNA extraction procedure was continued by using the DNeasy R Plant Mini Kit (QIAGEN, Hilden, Germany) following the manufacturer's instructions. Isolated DNA was quantified by using a TU-1950 spectrophotometer (Yima Opto-electrical Technology Com., Ltd., Xi'an, China) and stored at −20 • C until needed. For Metarhizium isolates, PCR amplification and sequencing approaches for the intergenic region, Mz_FG543igs gene, were performed Keyser et al., 2015). The master mix reagents for each PCR reaction consisted of 6.5 µl 2x Easy Taq R PCR SuperMix (TransGen Biotech Co., Ltd., Beijing, China), 2.5 µl of each primer, 11.5 µl Milli-Q H 2 O, and 2 µl DNA sample. The primers used were Mz_FG543igs_1F (5 -ATT CAT TCA GAA CGC CTC CAA-3 ) and MzFG543igs_4R (5 -GGT TGC GAC TCA CAA TCC ATG-3 ). PCR amplification was initiated at 95 • C for 2 min, then 40 cycles of denaturation for 30 s at 95 • C, annealing for 30 s at 62 • C, elongation for 60 s at 72 • C, and a final extension for 15 min at 72 • C. For phylogenetic placement of Beauveria and Isaria isolates, the partial sequence of the nuclear intergenic region Bloc was amplified by PCR using the primer pair B5.1F (5 -CGA CCC GCC AAC TAC TTT GA-3 ) and B3.1R (5 -GTC TTC CAG TAC CAC TAC GCC-3 ) (Rehner et al., 2011). The PCR cycle conditions for Bloc region were as follows: 94 • C for 3 min as initial denaturation, followed by 35 cycles of 94 • C for 30 s, 56 • C for 30 s, and 72 • C for 1 min, and a final extension of 72 • C for 10 min. The PCR assays were performed in a Mastercycler R Pro S (Eppendorf AG, Hamburg, Germany). The PCR products were visualized on 2% agarose gel to ensure strong single bands and purified using Wizard R SV Gel and PCR clean-up kit (Promega Corporation, Madison, WI, United States). All the primers used for amplification were also used for sequencing. Positive amplicons were sequenced in both directions in their entirety to ensure accurate readings at both primers' ends with the diluted PCR primers (10 pmol/µl) using MGISEQ-2000 Genetic Sequencer (BGI-Tech, Beijing, China). The sequences were edited, and the forward and reverse reads were assembled into contigs as implemented in BioEdit Sequence alignment editor ver. 7.2.5 for Windows software (Hall, 1999). All sequences were deposited at NCBI under the accession numbers listed in Supplementary Table S2.

Phylogenetic Analyses
All the fungal sequences were typically examined for initial taxonomic diagnosis through BLAST-based (Altschul et al., 1990) similarity search in the Genbank database. To avoid inappropriate diagnosis, our molecular identification was not completely reliant on the database searches. Therefore, we used phylogenetic analyses to clarify the taxonomic placement of the obtained isolates. Additionally, sequences of known Metarhizium spp. and Beauveria spp. isolates from Rehner et al. (2011) and Kepler et al. (2014) were included in the phylogenetic analyses. In the second step, the fungal sequences were aligned to the reference sequence alignment with MAFFT ver. 7.0 (Katoh and Standley, 2013) with the default settings as implemented on CIPRES Science Gateway ver.3.3 (Miller et al., 2010) 8 . The nucleotide alignments of all datasets were used separately to infer maximum likelihood phylogenies as implemented in RAxML, using the RAxML-HPC2 on XSEDE (Stamatakis, 2014) webserver at CIPRES Portal (Miller et al., 2010). Bootstrap support was determined using rapid bootstrapping with 1,000 replicates. Bootstrap values ≥ 70 were considered as significant (Kepler et al., 2014). The GTRGAMMA model was applied for all datasets. The final dendrograms were depicted using MEGA ver. 7.0.26 for Windows software (Kumar et al., 2016). The final phylogenetic trees were edited by Adobe Illustrator CC, 2019. 9

PCR Diagnosis of Beauveria and Metarhizium Mating Type
All obtained Beauveria and Metarhizium isolates were screened for internal segments of MAT1-1-1 and MAT1-1-2 mating type idiomorphs using PCR. The primer sequences and the PCR cycle conditions for Beauveria and Metarhizium were used according to Meyling et al. (2009) and Kepler et al. (2015), respectively.

Ion S5 XL Sequencing and Sequence Analysis
In order to delineate the composition of insect-fungal pathogen communities among different geographical locations, we conducted high-throughput Ion S5 XL paired-end sequencing. The sequencing statistics have been summarized (Supplementary Table S3). A total of 7,790,252 reads were obtained from 98 soil samples; 5,979,723 high-quality reads (49,619 ± 42,532 reads per sample), with an average length of 304 bp, were retained after length filtering as well as removal of chimeric sequences. The validated reads ranged from 73,884 (M) to 80,335 (G). The total count ranged from 2,2430,757 (SH) to 2,436,401 bp (G). The length of the average reads ranged from 302.31 to 305.307 bp. The average GC content percentage ranged from 49.46 to 50.87% (Supplementary Table S3). Finally, the reads were sorted into samples by detecting the barcodes for global fungal databases (UNITE) for the three entomopathogenic hypocrealean family genera, and the rest of the sequences were removed. Consequently, 31,713 sequence reads of insect-fungal pathogens ITS2 (20,427 sequence reads of Clavicipitaceae, 7,885 sequence reads of Cordycipitaceae, and 3,401 sequence reads of Ophiocordycipitaceae) were clustered into 615 OTUs. The K sampling sites were detected with the maximum number of EF OTUs (219) out of 1,681 average fungal OTUs (13.02%), and the lowest number of EF OTUs (14)

Taxonomic Composition and Abundance Analysis
The average relative abundance of EFs from the overall fungal community in terms of the sampling site revealed results in the following order Y > K > SH > M > F > G > B (0.60, 0.53, 0.41, 0.40, 0.087, 0.079, and 0.067%). The summary of the results revealed that insect-fungal pathogen sequences belonged mainly to Clavicipitaceae (46.33%), followed by Cordycipitaceae (37.55%), and Ophiocordycipitaceae (16.12%) (Figure 2A). The Clavicipitaceae family had the highest relative abundance in all the sampling locations in forest ecosystems and in grassland ecosystems except at the SH location (38.30%) (Figure 2B), and Cordycipitaceae (55.70%) prevailed at the SH location (Figure 2A). The most abundant EFs genus was Metarhizium (44.49%), followed by Beauveria (20.61%) and Lecanicillium (11.43%) ( Figure 2C). All sampling locations were dominated by the genus Metarhizium, but its abundance followed the trend of G > B > M > K > F > Y > SH (56.0,50,48.6,45.4,44.0,42.9, and 36.5% of total EF OTUs, respectively) ( Figure 2D). The relative abundance of Metarhizium differed significantly between B and SH (P = 0.0017), between G and SH (P = 0.0038), and between B and Y (P = 0.0259). Metarhizium was the most predominant genus in all the sampling locations except at the SH location, in which Beauveria (42.01%) was found to be prevalent.
The Venn diagram displays the shared OTU counts as well as the unique OTU counts among the collection sites. The K collection site had the highest number of unique OTUs (11), representing 8.46% of the K collection site OTUs. The lowest number of unique OTUs was displayed by the Y collection sites and the M collection sites (one). The SNFP and the SH sampling sites displayed six and four unique OTUs, respectively. Y had the highest number of shared OTUs, with K (six). Overall, only five of the total number of OTUs determined in the different samples were present in all the individual samples, and these may be considered as forming a part of the 'core' soil insect fungal pathogens (Figure 3). The core OTUs were identified as Metarhizium sp., Beauveria sp., Metarhizium anisopliae, Beauveria bassiana, and Isaria farinose.
The frequency of 25 EF species from three fungal families in each sampling location has been shown in the species-level assignment plot (Figure 4). The most abundant EF species was Metarhizium sp. from the G sampling location (39.67%), while Metarhizium flavoviride and Lecanicillium spp. from the K sampling location showed the lowest abundance (0.7%). Several taxa at the species level, including Purpureocillium sp.

Alpha Diversity
Species richness, as indicated by the Chao1 index, differed significantly between the Y and the G sampling locations (P = 0.0491). Maximum chao1 and minimum chao1 indices were presented by the SH (4.717) and the G (2.45) sampling sites, respectively. The Shannon index displayed a significant difference between the G and the M (P = 0.0356) and the K sampling locations (P = 0.  Table S4).

Beta Diversity
In order to determine whether the community compositions were statistically different, we conducted an ANOSIM. The results indicated that the structure of insect-fungal pathogen communities at the K collection site differed from that of the collection sites at F (R = 0.1752; P = 0.013), SH (R = 0.1960; P = 0.001), Y (R = 0.1695; P = 0.003), and G (R = 0.2112; P = 0.004). There were differences between the structure of the fungal pathogen communities at the M collection site and those of the collection site at B (R = 0.3783; P = 0.006), F (R = 0.4846; P = 0.002), G (R = 0.5528; P = 0.003), and Y (R = 0.3139; P = 0.004). Besides the differences mentioned above, marked differences were observed between the G collection site and the Y collection site (R = 0.1701, P = 0.013) and SH (R = 0.1934, P = 0.01) (Supplementary Table S5). In addition, distance boxplots of the weighted-and the unweighted-UniFrac matrices were used to evaluate the diversity between sampling groups. The unweighted UniFrac results indicated that the K sampling site showed the highest variation, while the weighted UniFrac results indicated that the G collection showed the maximum variation (Figures 5A,B). The NMDS plot showed that the sample distribution of OTUs at the various collection sites was relatively concentric (Figure 5C). The weighted UniFrac and the unweighted UniFrac phylogenetic beta diversity yielded substantially different results regarding the structure of insectfungal pathogen diversity. The results of the weighted UniFrac analysis suggested two main clusters, with the SH collection site forming a separate cluster. By contrast, the results of the unweighted UniFrac analysis suggested three main clusters (Figures 5D,E). Grouping the B and the F collection sites together using both the weighted-and the unweighted-phylogenetic beta diversity revealed a relative similarity between the insect-fungal pathogen communities at both collection sites. We compared the collection sites at the family level using t-test and found that the SH collection site differed from the collection sites F (P = 0.040), G (P = 0.006), K (P = 0.019), and M (P < 0.001) for the Cordycipitaceae family. Furthermore, the SH collection site differed from the F collection site (P = 0.026) for the Ophiocordycipitaceae family (Supplementary Figure S3). At the genus level, the relative frequencies of Beauveria, Lecanicillium, Isaria, and Tolypocladium differed at eight, five, two, and one pairwise comparison, respectively. More specifically, the results of the pairwise comparisons showed that the frequency of  Metarhizium marquandii differed at six, five, four, four, and three pairwise comparisons, respectively (Supplementary Figure S5). LDA, coupled with effect size measurement (LEfSe), was used in order to detect differentially abundant EF taxa from different sampling sites that may serve as biomarkers (Figure 6). The results showed that, at the family level, Cordycipitaceae (P < 0.05) was the discriminator in the soil samples of the SH collection site. Moreover, at the genus level and at the species level, Beauveria (P < 0.05) and B. bassiana (P < 0.01) were over-represented at the SH collection site, respectively. When all sites were investigated, Ophiocordyceps and Ophiocordyceps arbuscula of the M collection site emerged as discriminators at the genus level and at the species level, respectively. At the species level, L. psalliotae (P < 0.05) and Tolypocladium cylindrosporum (P < 0.01) were over-represented in the soil samples of the G and the F collection sites, respectively (Figures 6A,B).

Natural Distribution and Phylogenetic Diversity Using the Baiting Method
Fungal isolation based on standard soil baiting with the T. molitor larvae resulted in the isolation of 161 Metarhizium, 64 Beauveria, and five Isaria isolates from the soil samples collected at the 98 forest sites and grassland sites ( Table 2). The overall distribution of the Metarhizium and Beauveria isolates by area was found to be significantly heterogeneous (χ 2 = 39.217; df = 4; P < 0.0001) and homogenous (χ 2 = 11.313; df = 4; P = 0.0233), respectively. All fungal isolates were sub-cultured and phylogenetically assigned to the species level using Mz_FG543 and Bloc loci for Metarhizium and Beauveria, respectively (Figures 8, 9). The Bloc region for the five Isaria isolates was amplified successfully and aligned with the Beauveria reference sequences for phylogenetic assignment (Figure 9). A total of four species lineages (Metarhizium pingshaense, Metarhizium robertsii, Metarhizium brunneum, and M. flavoviride) were identified in Metarhizium, and their phylogenetic placement has been shown (Figure 8). Furthermore, nine Metarhizium isolates were obtained from the SH location (SH3 sampling site) ( Table 2). The initial morphological characterizations revealed the identity of these nine isolates as belonging to the M. flavoviride species complex, where Mz_FG543 primers failed to amplify the intergenic spacer region of those isolates. Amplification and sequencing of the 5' elongation factor-1 alpha region (5 -TEF) were accomplished for these isolates according to the procedure of Kepler et al. (2014). The result suggested a potential new sister cluster with the Metarhizium pemphigi clade within the M. flavoviride species complex (Supplementary Figure S7). The sequence types of the 5'-TEF gene fragment were submitted to GenBank (Supplementary Table S2). The maximum likelihood of the phylogeny of Mz_FG543 was topologically congruent and consistent with that of prior analyses (Bischoff et al., 2009;Kepler and Rehner, 2013;Kepler et al., 2014;Rehner and Kepler, 2017). The backbone and the monophyly of the PARBH species complex (M. pingshaense, M. anisopliae, M. robertsii, M. brunneum, and M. humberi) were strongly supported by the phylogenetic assessment of the Mz_FG543 locus with M. robertsii, the sister to a clade containing M. anisopliae, M. pingshaense, and M. brunneum, creating the basal-most clade (Rehner and Kepler, 2017;Luz et al., 2019). M. pingshaense was the most common species, with 74 isolates. It was divided into four major subclades and subsequently referred to as the M. pingshaense subclades 1 (n = 17), 2 (n = 6), 3 (n = 13), and 4 (n = 38) (Figure 8). Sublades 2, 3, and 4 were not accompanied by any M. pingshaense reference sequences. The second most abundant species was found to be M. robertsii, with 51 isolates, and it was subgrouped into three major subclades, subsequently referred to as M. robertsii subclades 1 (n = 1), 2 (n = 14), and 3 (n = 36) (Figure 8).  Subclade 2 was not clustered with the M. robertsii reference sequences. A majority of the terminal and internal branches of both the M. robertsii and the M. pingshaense isolates were unsupported and formed shallow grades. The two remaining species, M. brunneum and M. flavoviride, were composed of 26 isolates and one isolate, respectively. The Metarhizium brunneum clade was divided into three intraspecific phylogenetic subclades, subsequently referred to as M. brunneum subclade 1 (n = 2), 2 (n = 5), and 3 (n = 19). Subclades 1 and 3 were not grouped with any M. brunneum reference sequences. The phylogenetic diversity and the relationships among Isaria and Beauveria isolates were evaluated via the molecular phylogeny of the Bloc locus. All Beauveria isolates were placed within the B. bassiana clade (Rehner et al., 2011) and further subdivided into four terminal subclades with a strong bootstrap support (Figure 9), whereas the majority of the terminal branches were unsupported. All four subclades were grouped with the B. bassiana reference sequences.

DISCUSSION
To the best of our knowledge, the current study utilized both insect baiting and HTS techniques to provide the most comprehensive analysis of the natural distribution of insectassociated fungi in forest habitats and grassland habitats to  date. Our results demonstrated that Metarhizium was the most dominant genus, occurring most frequently at all collection sites except at the SH location. Beauveria was the most predominant genus at the SH collection site. Species abundance distribution (SAD) is one of the most basic models used in studies associated with population ecology. SAD describes communities that contain a few highly numerous species and many uncommon species. It is considered as a general pattern in ecology (Shade et al., 2018), and insect-associated fungi are not excluded from this comprehensive regulation. Similarly, our HTS results indicated that the community structure of EFs was composed of some highly abundant genera (Metarhizium, Beauveria, and Isaria) and several rare taxa such as Lecanicillium, Tolypocladium, Purpureocillium, Hirsutella, Haptocillium, and Ophiocordyceps (Figures 1C,D). The three leading hypocrealean families are composed of genera and species which regularly colonize multiple hosts, varying from plants, nematodes, and fungi to a broad range of insects (Kepler et al., 2012(Kepler et al., , 2014Araújo and Hughes, 2016). In addition to the commonly observed species M. anisopliae s.l. and B. bassiana, this work detected other invertebrate pathogens within these three families via deep sequencing of the soil samples. For instance, the Pochonia genus comprised of species that not only show pathogenicity to nematodes but also infect mollusk eggs and rotifers (Zare and Gams, 2001) and demonstrate a capability to colonize the soil rhizosphere of plants (Maciá-Vicente et al., 2009;Moonjely and Bidochka, 2019). In the current study, P. chlamydosporia was detected with a frequency of 1.7 and 3% at the SH and the K sampling collection sites, respectively (Figure 4). Other invertebrate pathogens, such as M. marquandii (Paecilomyces marquandii), which was mostly predominant (6.1%) at the K collection site, were also detected (Figure 4). The nematicidal activity of M. marquandii has been reported in a previous study (Genier et al., 2016). However, this taxon is not a member of the core Metarhizium clade generally considered for biocontrol programs (Kepler et al., 2014). The pathogens of nematodes also exist in other taxa, which are classified outside of Clavicipitaceae, such as some species of Hirsutella and Haptocillium (Kepler et al., 2014;Simmons et al., 2015). Hirsutella (asexual morphs) species are classified as a host-specific entomopathogenic soil fungus. These species are physiologically and morphologically distinct from Metarhizium and Beauveria, which have over 700 hosts (Qu et al., 2018). Hirsutella thompsonii is a widely studied species from this genus, which is potentially pathogenic toward mites (Tigano et al., 2006). The determination of OTUs at the species level detected two species of this genus: H. thompsonii and Hirsutella subulata. The former was detected at the SH (1.7%) and at the M (8.5%) collection sites, whereas the latter was detected at the SH (0.8%), the F (2%), and the G (2%) collection sites (Figure 4). The genus Lecanicillium contains both entomogenous and fungicolous species (Berendsen et al., 2010) and is commonly found in the tropical and the temperate regions (Pearce et al., 2009). The Lecanicillium fungicola and L. psalliotae species from this genus were detected with HTS technique and were present in three sampling sites and five sampling sites, respectively, while also differing in prevalence at their respective sites (L. fungicola 45.99% vs. L. psalliotae 51.92%) (Figure 4). It has been suggested that L. fungicola shows potential for pathogenicity against insects (Berendsen et al., 2010). This may be due to the ability of this fungus to infect mushrooms and the close relatedness of this species to a variety of insect pathogens (Amey et al., 2007;Berendsen et al., 2010). Interestingly, L. psalliotae, which was initially recognized as a pathogen of mushrooms (Treschow, 1941), is presently most often discussed in relation to its nematocidal activity (Pirali-Kheirabadi et al., 2007). We should mention that the species-level assignment (Figure 4) must be interpreted with some caution because using ITS2 for soil-fungal entomopathogen composition structures does not provide informative molecular delimitation below the genus level accurately, especially for the closely related species in some genera such as Metarhizium and Beauveria, and that is a straightforward consequence of the dependence of culture-independent metabarcoding approach on databases. For this reason, proper identification should be integrated with additional standard techniques (Tremblay et al., 2019). In this study, Metarhizium and Beuaveria were relatively abundant across all sampling sites. Prior studies that have laid emphasis on the ability of these generalist fungi to interact with plants, suggesting that the type of vegetation canopy may be associated with their widespread distribution in ecosystems (Wyrebek et al., 2011;Behie et al., 2012;Behie and Bidochka, 2013;Kepler et al., 2015;Keyser et al., 2015;McKinnon et al., 2017McKinnon et al., , 2018Tall and Meyling, 2018). Thus, based on the unique relationship of Metarhizium and Beauveria lineages with the plant rhizospheres, the high frequencies with which these are found in the forest habitats and the grassland habitats can be conceivably illustrated by adding an analysis of fungal-plant associations to the study. Bidochka et al. (2005) hypothesized that the global distribution and the abundance models of Metarhizium lineages may be based on two central principles. In temperate regions, fungal genotypes may be linked with habitats characterized by specific abiotic environmental determinants that select for the genotypes that may be present (habitat selection), while in the tropical and in the subtropical regions, fungal genotypes may be associated with particular hosts (host selection). Bidochka et al. (2002) expanded the habitat selection hypothesis for entomopathogenic soil fungi by including the association between B. bassiana distribution and habitat availability. They contended that the congruence between B. bassiana lineages and habitat nature was rationally indistinguishable from that of the M. anisopliae lineages (Bidochka et al., 2001), in which two distinct genetic groups (M. robertsii and M. brunneum) were found to be associated with different localities. Such circumstantial evidence presented by these studies represent a significant paradigm shift, wherein it is habitat selection, not insect host selection, that determines the community composition of insect-pathogenic fungi (Bidochka et al., 2001(Bidochka et al., , 2002. The four Metarhizium species obtained using the baiting method in this study were not consistently distributed across the investigated collection sites. M. pingshaense was the most common species obtained from the soil collection sites (n = 74), followed by M. robertsii (n = 51), M. brunneum (n = 26), and M. flavoviride (n = 1). M. pingshaense was obtained throughout the investigated sites except at the SNFP collection site and was isolated from 45 samples of the total 490 soil samples. By contrast, M. brunneum, which was recovered from 15 soil samples and since it was not isolated from the G, the B, and the SH collection sites, was more narrowly distributed. These results are in line with several investigations from agricultural soils that have confirmed the locality-specific abundance and distribution of the Metarhizium lineages, at least in northern Europe (i.e., Denmark) (Meyling et al., 2011;Steinwender et al., 2014Steinwender et al., , 2015Keyser et al., 2015). The M. robertsii isolates were obtained from all seven collection sites except the boundary between the grassland and the forest regions (B). Surprisingly, our baiting results, based on the infection of mealworm larvae, highlighted that the M. pingshaense lineages had a higher abundance and frequency compared to those of the M. robertsii lineages, although the M. robertsii lineages had a better distribution and a wider occurrence across the sampling sites ( Table 2). The HTS results verified the highest frequency of the Metarhizium lineages from the G (grassland) collection site (56.0%) (Figure 1A), and the intergenic marker phylogeny indicated that all the Metarhzium isolates obtained from this site belonged to the M. robertsii clade (Figure 8). These results might lead us to conclude that the distribution of the M. pingshaense lineages was restricted to forest regions with a higher average annual temperature (Supplementary Table S1), while the M. robertsii lineages can stand the relatively extreme environmental temperatures and may possibly possess better ecological fitness in different environments. In addition to the abiotic theory, M. robertsii could be more closely associated with grassland habitats (Wyrebek et al., 2011). The average annual temperature at SNFP is approximately −1.3 • C, with a prolonged cold winter and a brief growing season (May to September) (Yang et al., 2014). The capability of Metarhizium for cold tolerance was assessed, and Fang and St. Leger (2010) reported that the Metarhizium lineages acquired the cold shock response genes (CRP1 and CRP2) evolutionarily via horizontal gene transfer from soil-dwelling bacteria within the same ecological niche. The low level of the Metarhizium isolation at the SNFP sampling site (Table 2) encouraged us to perform a selective agar medium to evaluate the Metarhizium frequency according to the procedure of Keyser et al. (2015). The results (not shown here) showed that the number of colony-forming units per plate (CFU/plate) was higher at the F sampling sites, whereas at the B and the G sampling sites the number of CFU was lower. These results were not consistent with the HTS and the baiting results, suggesting that different isolation methods resulted in various fungal frequencies. As highlighted by other studies, using a single detection method raises many doubts regarding the abundance of fungal pathogens in the environment (Keller et al., 2003;Goble et al., 2010;Medo and Cagáň, 2011). During the course of the study, we considered these procedural deficiencies, and the baiting approach may not be genuinely quantitative (Jaronski, 2007;Scheepmaker and Butt, 2010). Due to this reason, we performed HTS to provide an additional dimension that facilitated a deeper understanding of the insect-associated fungal distribution. In contrast to the predominance of the M. anisopliae species complex observed in the sampling sites, we reported a potentially new sister clade to the M. pemphigi cluster composed of nine fungal isolates within the M. flavoviride species complex using the sequence of the 5 -TEF gene fragment (Supplementary Figure S7). These isolates were only obtained from the SH collection site, indicating locality-specific abundance and distribution of these isolates, and further studies are needed to understand the potential interactions of this clade with hosts (plant roots and insects) and its above-and below-ground ecology. Our study indicated that both M. brunneum and B. bassiana were likely to be obtained from the site with a relatively high elevation (2,447 m) and from regions with elevations below (300-100 m). Our major finding for M. brunneum in terms of elevation is substantiated by those of Masoudi et al. (2018) (Vänninen, 1996) and 64 • N in Iceland (Oddsdóttir et al., 2010). Our HTS data indicated a negative correlation between M. anisopliae s.l. lineages and geographical coordinates (42 • 32 .347 N, 29 • 55 .457 E) (Figure 7). A few studies have reported that latitude was correlated with B. bassiana distribution, while longitude was correlated with M. anisopliae natural distribution (Bidochka et al., 1998;Quesada-Moraga et al., 2007;Garrido-Jurado et al., 2011). Our HTS results and, the Mantel tests verified that community dissimilarities increased significantly with geographical distance (the Mantel correlation = 0.228, P = 0.001). Thus, we inferred that geographical coordinates may be the factors that possibly influence insect-fungal pathogen community composition in the investigated sites. Data gathered by natural occurrence studies using bait insect species possibly underestimate not only the diversity and frequency of commonly occurring EF species but also that of rare taxa, such as the hypocrealean insect pathogens (Vänninen, 1996;Chandler et al., 1997;Bidochka et al., 1998Bidochka et al., , 2001Klingen et al., 2002;Ali-Shtayeh et al., 2003;Keller et al., 2003;Meyling and Eilenberg, 2006;Quesada-Moraga et al., 2007;Inglis et al., 2008Inglis et al., , 2019Sun et al., 2008;Goble et al., 2010;Medo and Cagáň, 2011;Meyling et al., 2011;Rocha et al., 2013;Wakil et al., 2013;Steinwender et al., 2014;Kepler et al., 2015;Keyser et al., 2015;Rezende et al., 2015;Hernández-Domínguez and Guzmán-Franco, 2017). Some of these studies have reported lineages of Metarhizium, Beauveria, and Isaria as the common pathogenic microorganisms found in soil environments ( Table 1). The results of our HTS proved that M. anisopliae s.l., B. bassiana, and I. farinosa were the core insect-associated fungal genera in the assessed habitats (Figures 3, 4). The results of HTS help us to distinguish the difference between 'locally abundant' OTUs and the 'core' (frequently abundant) diversity of fungal pathogens in the soil. Moreover, the results of the baiting method showed that the Isaria isolates were only obtained from the SH collection sites (SH4, SH8, and SH23) (Figure 9). Furthermore, Krona analysis indicated that the frequencies of Isaria at SH4 (two isolates), SH8 (two isolates), and SH23 (one isolate) sampling sites were 88, 65, and 73%, respectively (Supplementary Figures S10-S12). Also, B. bassiana (four isolates) was obtained from the SH23 sampling site with a frequency of 17% (Supplementary Figure S11). These results indicate that the detection of pathogenic microorganisms using the baiting method may be strongly influenced by two factors, which are fungal frequency and the potential entomopathogenicity of the pathogens. The co-occurrence of B. bassiana and the Isaria sp. was detected at the SH23 sampling site, where the relative frequency of B. bassiana was less than that of the Isaria sp., but the number of infected cadavers with Beauveria was four times higher than that infected with Isaria. This result showed that, compared to the Isaria lineages, the B. bassiana isolates demonstrated superior virulence even when present at a lower frequency (17%) and exhibited this virulence by competing for hosts and by infection of more insects. These two factors may confer a better competitive capability on the pathogen for colonizing insect larvae, which enables it to exclude and/or affect the presence of other pathogens. Based on these data, it may be inferred that the negative interaction between pathogens caused by the cooccurrence of other fungal pathogens in soils may lead to direct competition in which the subdominant taxa may be excluded from the outcome of baiting (competitive exclusion) (Gold et al., 2009). Additionally, it has been suggested that each bait insect species method may highlight the discriminating isolation of specific fungal species/taxa (Klingen et al., 2002;Goble et al., 2010) under a particular environmental condition (Mietkiewski and Tkaczuk, 1998). Klingen et al. (2002) highlighted that the fungal species obtained varied significantly depending on the insect species used as bait. For example, their results showed that T. cylindrosporum and B. bassiana were obtained only using Delia floralis (Diptera: Anthomyiidae) and G. mellonella, respectively (Klingen et al., 2002) (Table 1). Our HTS data detected that T. cylindrosporum followed the same frequency trend of the SNFP sampling location, F > B > G (14 > 2.5 > 2% for all the three families of the EF OTUs, respectively) and detected T. pustulatum only at the K (2.3%) and at the Y (2.8%) collection sites. These results indicate that some Tolypocladium species may possibly have a limited geographical distribution, although Tolypocladium sp. was detected at all locations except at the F sampling site (Figure 4). This may be a conceivable reason for the failure of our insect bait method, which was similar to the methods used by Steinwender et al. (2014); Keyser et al. (2015), and Korosi et al. (2019), to detect other hypocrealean fungi such as Lecanicillium, Hirsutella, and Tolypocladium. Environmental factors can possibly underlay EFs habitat preference. Jaronski (2007) provided a meticulous overview of the soil physicochemical ecology in relation to hypocrealean insect-associated fungi and concluded that the complex interactions of fungal pathogens with the abiotic parameters of soil may exert both negative as well as positive effects on the community structure of soil-inhabiting entomopathogens. In order to develop a comprehensive picture of insect-associated fungal diversity in soils, additional studies that investigate the influence of edaphic parameters (physicochemical properties) on the abundance and the occurrence of EFs in soil environments may be needed.
The current study explored and compared the hypocrealean EF isolates (Metarhizium and Beauveria) that were obtained using mating-type PCR screening that addressed their potential for sexual reproduction and verified their distribution on the phylogenetic tree (Figures 8, 9). Mating-type PCR assays revealed that some isolates of M. pingshaense, M. brunneum, and B. bassiana contained both mating-type genes (homothallic). In a cutting edge global-scale survey, Rehner and Kepler (2017) reported the presence of a MAT1-1-1 or a MAT1-2-1 mating-type idiomorph in the M. anisopliae PARBH species complex, where they speculated that all PARBH species were fundamentally outcrossing (Rehner and Kepler, 2017). Results from the mid-Atlantic region showed that individual M. robertsii and M. brunneum strains possessed a single mating type, either MAT1-1-1 or MAT1-2-1, indicating that these species were heterothallic . However, Pattemore et al. (2014) highlighted the existence of both MAT1-1-1 and MAT1-2-1 idiomorphs in the genome sequence of an Australian isolate of M. anisopliae (BRIP 53293, APNB00000000) and determined that the isolate was homothallic and thus potentially self-fertile. Rehner and Kepler (2017) reported that all the Australian M. anisopliae isolates in their study were of a single mating type, either MAT1-1-1 or MAT1-2-1. Matingtype characterization of our B. bassiana isolates showed that the 19 isolates possessed both idiomorphs. Meyling et al. (2009) reported that most B. bassiana clades within an agroecosystem in Denmark were fixed at a specific mating type, indicating that they retained the potential at least for heterothallic sexual activity, and that one population (Eu-1) contained two mating types, suggesting a potential for sexual reproduction at that study site for this population. Ramos et al. (2017) reported a skewed distribution of the mating-type system of MAT1/MAT2 (146/4) of the B. bassiana populations from Cuba, indicating a limited potential for the recombination of complementary mating specificity (or mating types). Our results of matingtype characterization of the obtained EF isolates indicate a major development, showing the possible sexual reproduction processes within entomopathogen species of the M. anisopliae PARBH complex and the B. bassiana populations inhabiting Chinese soils. The inclusion of the mating-type assignment PCR assays as a part of natural distribution surveys of hypocrealean entomopathogenic fungi may provide useful information related to the spatial distribution of mating-compatible individuals in the environment, thereby enabling the identification of populations in which the EF lineages experience occasional sexual recombination events with profound consequences for their population dynamics and evolutionary trajectories.

CONCLUSION
Despite the abundance of information regarding the natural distribution of widespread EFs, such as Metarhizium, Beauveria, and Isaria lineages in different habitats, obtained via baiting methodology under laboratory conditions, there is a lack of understanding with respect to the frequency profiles of other EFs in soils, irrespective of certain circumstances, such as incubation temperature or the type of insect used as a bait, that may affect the isolation of the entomopathogens. This current study showed that the forests and the grasslands of two geographical regions of China harbored a relatively high frequency of Metarhizium and Beauveria and highlighted the abundance of those fungal pathogen taxa whose frequency or presence has not been frequently reported by other studies using HTS method for soil samples as well as the discovery of a potentially new cluster within the M. flavoviride species complex. The study indicates that the prevalence of insect fungal pathogens in certain habitats should be a greater focus of future efforts to address the factors influencing the locality-specific abundance and distribution of such fungi. This may lead to a better understanding of the dynamics and the distribution of EFs in order to recommend insect-associated fungi as a biocontrol agent.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the raw sequences produced in the current study can be reached through the National Center for Biotechnology (NCBI) under the BioProject ID PRJNA551928, accession SUB5873947.

AUTHOR CONTRIBUTIONS
AM and HW originated the research. AM, MW, XZ, CW, HW, ZQ, and WW obtained the soil samples. AM performed the data analysis, the laboratory experiments, and wrote the manuscript. JL offered assistance in designing the research and improving the quality of the manuscript. All the authors corrected the first draft of the manuscript.

FUNDING
This study was made possible with the assistance of the National Natural Science Foundation of China (31672365) and the National Key Research and Development Project of China (2018YFC0506900). This research was also partially supported by a postdoctoral research grant from Hebei Normal University to AM.