Integrative taxonomy of Metarhizium anisopliae species complex, based on phylogenomics combined with morphometrics, metabolomics, and virulence data

Metarhizium anisopliae (Clavicipitaceae, Hypocreales) is a globally distributed entomopathogenic fungus, which has been largely studied and used in agriculture for its potent entomopathogenicity. Since its taxonomic establishment as a member of Metarhizium, many closely related taxa have been described with highly similar morphology (cryptic species). A holotype specimen of M. anisopliae is not extant, and the ex-neotype strain (CBS 130.71) does not form a monophyletic clade with other strains, up to now, recognized as M. anisopliae sensu stricto. In this study, we have conducted an integrative taxonomic treatment of M. anisopliae sensu lato by including the ex-neotype strain of M. anisopliae, other unknown strains from our collections identified as M. anisopliae s. lat., as well as other known species that have been previously delimited as closely related but distinct to M. anisopliae. By including whole-genome sequencing, morphometric analysis, LC–MS based metabolomics, and virulence assays, we have demonstrated that M. anisopliae s. str. should also include M. lepidiotae (synonym), and that M. anisopliae s. str. differentiates from the other species of the complex by its metabolome and less severe entomopathogenicity. New taxa, namely M. hybridum, M. neoanisopliae and M. parapingshaense spp. nov., are proposed. The novel taxa proposed here have strong phylogenomics support, corroborated by fine-scale differences in the length/width of conidia/phialides, while the metabolomics and virulence data still largely overlap. We have also demonstrated via population genomics data the existence of local clonal lineages, particularly the one corresponding to the persistence of a biocontrol candidate strain that has been used in the field application for three years. This study showcases the utility of combining various data sources for accurate delimitation of species within an important group of fungal biocontrol agents against pest insects. Supplementary Information The online version contains supplementary material available at 10.1186/s43008-024-00154-9.


Introduction
Metarhizium anisopliae (Clavicipitaceae, Hypocreales) is a globally distributed species of entomopathogenic fungi widely used as a biocontrol agent against insect pests (Lubeck et al. 2008;Cao et al. 2016;Bamisile et al. 2018;Zhao 2023).This species was first discovered by Elias Metchnikoff (Metchnikoff 1879) as a pathogen of Anisoplia austriaca (Coleoptera) larvae in Ukraine (Zimmermann et al. 1995).Before the advent of molecular taxonomy, the nomenclature and classification of M. anisopliae had encountered numerous changes.The species was first described as Entomophthora anisopliae Metchnikoff (Metchnikoff 1879) before being transferred to Metarhizium based on molecular phylogenetics (Bischoff et al. 2009;Kepler et al. 2014).It was at times confused with and placed in Penicillium as P. anisopliae (Vuillemin 1904), and even with Oospora as O. destructor (Delacroix 1893;Aoki 1957;Kawakami 1960).Two major reasons for such confusion were that no original type material of Metchnikoff 's M. anisopliae was available, and molecular techniques were not used until the late twentieth century.Veen (1968) typified M. anisopliae in a French-written thesis, and proposed a neotype culture (CBS 289.67), a proposal that was accepted by Tulloch (1976) in her revision of the genus.ARSEF 7487 was also derived from the original type strain (ex-neotype CBS 289.67), isolated from an Orthopteran insect from Eritrea (part of Ethiopia at that time), and deposited at the culture collection of the United States Department of Agriculture (USDA).This ex-neotype has been used as the reference strain in molecular taxonomic studies of Metarhizium during the twenty-first century (Driver et al. 2000;Bischoff et al. 2009;Kepler et al. 2014).It was thus not isolated from the same insect, nor the same locality as the original specimen of Metchnikoff. Mongkolsamrit et al. (2020) showed that the M. anisopliae strain CBS 130.71, isolated from the root of the cereal crop Avena sativa in Ukraine (the same locality as Metchnikoff 's original specimen) and thus designated as the new ex-neotype culture of M. anisopliae, is phylogenetically distant from ARSEF 7487 and closely related to M. lepidiotae (Bischoff et al. 2009).However, ARSEF 7487 still remains as M. anisopliae, which has left a situation where the ex-neotype culture of M. anisopliae (CBS 130.71) does not form a monophyletic clade with other specimens recognized as the same species, posing a contradiction to the phylogenetic species concept used for delimiting fungal species in the molecular era (Taylor et al. 2000).
Before the neotypification of M. anisopliae by Veen (1968), other Metarhizium species had been discovered, for example M. album (Petch 1931), M. brunneum (Petch 1935), and M. flavoviride (Gams and Rozsypal 1973), while Tulloch (1976) only accepted M. anisopliae and M. flavoviride as valid species.Metarhizium guizhouense and M. pingshaense were subsequently added to the genus in the late twentieth century (Guo et al. 1986).These species were shown to be closely related and many of them were occasionally considered to be only "varieties" within M. anisopliae or M. flavoviride (Driver et al. 2000).With molecular data, M. anisopliae s. lat. is considered to form a monophyletic clade, distinct from M. flavoviride s. lat.(M.flavoviride species complex) (Kepler et al. 2014;Mongkolsamrit et al. 2020).Metarhizium anisopliae s. lat.(M.anisopliae species complex) includes the PARB Clade (M.anisopliae s. str., M. brunneum, M. pingshaense, and M. robertsii), the MGT Clade (M.guizhouense and M. majus), and also M. lepidiotae, M. acridum and M. globosum following Bischoff et al. (2009) and Kepler et al. (2014).Other recent species have been added to this complex mainly based on phylogenetic data, i.e.M. kalasinense (Luangsa-ard et al. 2017), M. gryllidicola and M. phasmatodea (Thanakitpipattana et al. 2020), M. clavatum and M. sulphureum (Mongkolsamrit et al. 2020).These recent studies have demonstrated the cryptic nature of M. anisopliae s. lat.where closely related species have very similar conidial shapes and sizes, whereas the classification of Metarhizium was mainly based on micro-morphological characters of asexual spores in the nineteenth and twentieth century (Sorokin 1883;Tulloch 1976).Despite the discovery of sexual morphs for some species (Sung et al. 2007;Mongkolsamrit et al. 2020), sexual characters (e.g.perithecial size, shape, and embedment; ascospore size and shape) have not been proven to be useful.Therefore, the identification of species within M. anisopliae s. lat.currently relies on molecular phylogenies.
The genus Metarhizium is also well-known for its capacity to infect a broad range of insects (St. Leger and Wang, 2020) and is believed to attack over 200 different host species (Veen 1968).It has been proposed that generalist species (PARB and MGT Clades) evolved from specialists as PARB and MGT Clades form a monophyletic clade branching from the most recent common ancestor shared with M. acridum which is specific to certain locusts and grasshoppers (Hu et al. 2014).The investigation on the variation of virulence and host specificity among intra-specific strains or strains from closely related species has been scarce.Previous studies focused on testing the virulence of a few strains, mainly from the PARB Clade known to infect a broad host range (Shapiro-Ilan 2009;Wang and Feng 2014;Lee 2019).The identification of Metarhizium isolates tested in these studies was not generally done in a taxonomically appropriate manner; either many strains were called "sp." (Accoti et al. 2021) or putative "M.anisopliae" strains were identified as such based only on ITS data with broad morphological recognition (Dong et al. 2016).Therefore, the lack of an evolutionarily and taxonomically welldefined framework hinders a comprehensive understanding of the evolution of virulence as this trait cannot be directly linked to the most accurate and recent molecular phylogenies of this genus.Furthermore, the potential of virulence data to contribute to the taxonomy of Metarhizium has been overlooked.Evaluating virulence of various Metarhizium strains with a clear taxonomic classification will not only improve our comprehension of virulence, but also contribute to an integrative taxonomic elucidation based on virulence data.
It has been demonstrated that the virulence/pathogenicity of Metarhizium species is related to the production of secondary metabolites (Rohlfs and Churchill 2011).Particularly, various destruxins produced by Metarhizium were shown to have insecticidal effects (Pal et al. 2007;Donzelli et al. 2012).Previous chemotaxonomic work has shown more quantitative difference in the production of secondary metabolites between distantly related species than closely related fungal species (Kobmoo et al. 2021;Wongkanoun et al. 2023), in line with the idea of cryptic species being recognizable only by nucleotide data.Recently, Barelli et al. (2022) have shown that the production pattern of destruxins is different between four Metarhizium species, three of which are from M. anisopliae s. lat.(i.e.M. acridum, M. brunneum, and M. robertsii) and one from M. flavoviride.However, there are more species to explore (Mongkolsamrit et al. 2020) and the lack of comprehensive metabolomics data prevents us from testing the hypothesis whether the production pattern of destruxins can be explained by evolutionary history.
Fungal cryptic species are increasingly being discovered thanks to whole-genome sequencing (WGS) data (Passer et al. 2019;Maharachchikumbura et al. 2021;Kobmoo et al. 2021).Combined with other sources of data such as morphological, metabolomics, and ecological data, WGS enables the elucidation of species status among closely related taxa (Seehausen et al. 2014), an approach called "integrative" or "holistic taxonomy" (Grube et al. 2017;Boluda et al. 2019;Stengel 2022).Abundant genomic resources available today for Metarhizium (Pattemore et al. 2014;Prasad et al. 2015;Wang et al. 2016) have been exploited principally for understanding the molecular basis to virulence and pathogenicity (Hu et al. 2014), but less to apprehend natural genetic variations (but see Rizal et al. 2024).Not only has WGS contributed to recent taxonomic discoveries, but this technique can also be used for tracking fungal strains of interest in field applications (Mei et al. 2020;Peng et al. 2021).
In this study, we have used WGS to identify and validate new species status for unknown strains from our old collections along with recent collections from an agricultural area where a candidate strain of Metarhizium (BCC 4849; Wasuwan et al. 2022) was released.
The main objective of this study is to resolve species status within M. anisopliae species complex.We examined (1) genomic diversity, (2) phenotypic variation on conidia and phialides features, (3) variation in secondary metabolite production, and (4) virulence towards an insect species (Spodoptera exigua, Lepidoptera), for species within M. anisopliae s. lat.By including Tulloch's ex-neotype M. anisopliae strain (ARSEF 7487) and the recent ex-neotype culture of M. anisopliae s. str.proposed by Mongkolsamrit et al. (2020) (CBS 130.71), the various sources of data have been combined under an integrative taxonomic framework to propose three new species, namely M. neoanisopliae, M. hybridum, M. parapingshaense, and to synonymize M. lepidiotae with M. anisopliae s. str.Metarhizium neoanisopliae BCC 4849 is shown here to persist in the agricultural area where it was applied to control insect pests.

Sample collection
Unidentified strains of Metarhizium spp.were ordered from the BIOTEC Culture Collection (BCC) (Table 1).These had been collected during numerous excursions around Thailand by the research team at BIOTEC.They were selected for this study based on our in-house taxonomic identification (ITS barcoding) which showed their affiliation to M. anisopliae s. lat.Most were isolated from mycosed insects following Mongkolsamrit et al. (2020); asexual spores (conidia) from the fresh specimens were harvested with a flame-sterilized needle and streaked on a potato dextrose agar (PDA) plate (PDA: freshly diced potato 200 g/L, dextrose 20 g/L, agar 15 g/L), incubated at room temperature for 24 h, then checked for contamination.Uncontaminated germinating spores were transferred to another PDA plate for incubation at 25 °C to get pure cultures before being stored at -80 °C in liquid nitrogen at the BCC.The specimens were dried and deposited at the BIOTEC Bangkok Herbarium (BBH).A few strains were recently isolated by the BIOTEC Biocontrol Technology Research Team from the soil of a fruit orchard (Table 1), where the strain BCC 4849 was used as a biocontrol agent.The soil-borne isolates were obtained using a streak plate method with the following protocol: Soil samples each weighing approx.500 g were collected from a depth of 10-15 cm using a trowel, removing litter.The samples were then placed in plastic bags and transported to the laboratory.Fungal isolation was carried out within two days of collection using a dilution plating technique and cetyltrimethyl ammonium bromide (CTAB) and oatmeal agar (OA) as a basal medium following the procedure described in Abdullah et al. (2015).

DNA extraction and whole-genome sequencing
The cultures were grown on PDA plates.A 1 × 1 cm 2 of agar per strain was cut to be put on another plate for incubation at 25 °C for one to two weeks.DNA extraction was done using the CTAB-based procedure described in Kobmoo et al. (2021).The genomic DNA was checked for purity and quantity with a Nanodrop spectrophotometer (Thermo Fisher) and 1% agar gel electrophoresis.For WGS, approximately 500 ng of genomic DNA was used for library construction following the MGIEasy FS DNA Library Prep Kit (MGI Tech, Shenzhen, China).Paired-end sequencing (150 bp) was performed using the MGISEQ-2000RS Sequencing Flow Cell V3.0 on the DNASEQ-400 according to the manufacturer's protocol.

Population genomics and phylogenomics analyses
MegaBOLT version 1.5.6.11(MGI) was used to extract the sample barcode and to control the quality of the raw read data of each sample.The controlled reads were mapped on the reference genome of M. anisopliae JEF-290 (Lee 2019), using bwa v0.7.17The SNPs were filtered to possess the following features: 1,000 < all-sample depth < 10,000, mapping quality > 40, quality normalized by depth > 20, Fisher strand bias < 10 and strand odd ratio < 3. The individual genotypes with depth < 10 and genotyping quality < 20 were marked as missing data.The filtering was based on the a posteriori distribution of the features.Only the SNPs with no missing data and minor allele frequency (MAF) > 0.05 were retained, resulting in a total of 98,085 SNPs.An initial snapshot of the diversity of Metarhizium spp.included in this study was obtained by an F84 distancebased neighbor-joining tree, inferred directly from the SNPs data.The population structure was inferred with a Bayesian clustering analysis implemented via the software FastStructure (Raj et al. 2014).Genotypes were then classified into multi-locus genotypes (MLG) to represent clonal groups.This classification was done using the R package poppr (Kamvar et al. 2014) following the method of Grünwald and Hoheisel (2006).To construct a phylogenomic tree, we first extracted gene sequences with nucleotides altered by SNPs using FastaAlterna-teReferenceMaker from the GATK toolbox (McKenna et al. 2010).This extraction was based on the gene model associated to the reference genome of M. anisopliae JEF-290 (Lee 2019).Positions with indels were marked as missing data.The identity of the single-copy genes was determined via the web server OrthoVenn2 (Xu et al. 2019).Only single-copy genes with a minimum length of 500 bp and a ratio of SNPs/nucleotide greater than 0.005 were retained.Such criteria ensured finding at least two SNPs for each gene alignment to guarantee a minimum level of phylogenetic signal, resulting in 237 genes.The gene sequences from different samples were aligned using MAFFT (Katoh and Standley 2013).Subsequently, gene-wise alignments were concatenated, resulting in a final alignment of 451,019 bp.The best maximum likelihood (ML)-based phylogenomic tree was constructed using IQ-TREE (Nguyen et al. 2015).Branch supports were determined with 1,000 replicates of ultrafast bootstrap (Hoang et al. 2018).The best sequence evolution model (TVM + I + G4) was inferred using ModelTest-NG (Darriba et al. 2020).The final tree was rooted with the M. flavoviride complex (M.flavoviride and M. frigidum).Additionally, a coalescent-based species tree was also estimated using ASTRAL-III (Zhang et al. 2018).This was done to consider potential conflicts between gene trees, which could introduce a systemic bias in the concatenated gene trees (Liu et al. 2015).

Morphological characterization
The macro-morphological characters and relevant data of the fungi, including the host, were examined under a dissecting microscope (Olympus SZ61).For micro-morphological characterization, phialides and conidia were mounted on a microscope slide with a drop of lactophenol cotton blue solution and measured using a compound microscope (Olympus CX31).The length and the width of conidia and phialides from pure cultures were measured for representative strains from the novel taxa (Table 1).For M. pingshaense and M. parapingshaense, conidia were also collected from dried specimens and measured for length and width.These data were subjected to an analysis of variance (ANOVA) (Kaufmann & Schering 2014) testing for the difference between species with the strains included as a random factor.Thirty phialides and conidia per strain were measured.The fungal cultures were grown on PDA and incubated at 25 °C under a light/dark cycle (L:D = 14:10).The cultures were observed for comparison of essential morphological characters, including the shapes and sizes of conidia, phialides, colony growth, and coloration.
For the descriptions, the color of specimens and cultures incubated on PDA was codified following the Royal Horticultural Society's Colour Chart (Royal 2015).

Metabolomics
Strains (Table 1) were grown in 200 mL of yeast malt extract broth (YM broth; malt extract 10 g/L, yeast extract 4 g/L, D-glucose 4 g/L, pH 6.3 before autoclaving) at 23 °C under shaking condition (140 rpm).Each culture was inoculated with 5 pieces of mycelia (7 mm diam) from a fully grown YM agar plate.Due to the growth rate differences among strains, the prolongation of each culture after glucose depletion was set to be half of the time required for that particular strain to reach glucose depletion.This method is justified as exhaustion of the first limiting nutrient influences the growth curve and should be considered for experimental standardization (Vrabl 2019).The glucose content during the culturing was estimated using urine glucose test strips (DIRUI ® , Jilin, China).
The extraction of secondary metabolites was performed according to Phainuphong et al. (2017).The mycelia were separated from the culture broth via vacuum filtration and the mycelial yield was determined.The culture broth was extracted twice with equal amount of ethyl acetate (v/v).The resulting organic phase was filtered through anhydrous sodium sulphate and evaporated on a rotary evaporator to dryness.The mycelia were soaked in acetone and extracted twice under sonication (30 min).The solvent was evaporated to yield an aqueous phase that was extracted twice with ethyl acetate, in a similar manner to the culture broth.The samples (Table 1) were processed in two batches separated by a year due to experimental limitations.
Each sample was analyzed using the instrumental settings and conditions reported previously (Charria-Girón et al. 2023).Raw data were pre-processed with MetaboScape 2022 (Bruker Daltonics, Bremen, Germany) in the retention time range of 0.5-25 min, and the obtained features were dereplicated based on their accurate molecular weight and MS/MS spectra against the compounds reported from Metarhizium spp. in the Natural Product Atlas (NP Atlas) database (Van Santen et al. 2019).For this purpose, MetaboScape performed automatic in silico MS/MS matching based on the InChI-encoded structures using the MetFrag algorithm in the absence of MS/MS reference data (Ruttkies et al. 2016).
The metabolomics analysis was conducted using the data of peak area.As the extracts were prepared in two batches, to remove any systematic bias due to separate batch processes, any feature absent in either batch was discarded and the raw data from the respective batch were normalized within the batch before merging the whole data.A principal component analysis (PCA) was conducted to detect any pattern corresponding to divergence between taxa.The analysis was conducted separately between the data from BE and CE.The data of annotated compounds were extracted to construct a heatmap based on Manhattan distance using the gplots package v.3.0.1 (Warnes et al. 2016) in R (R Core Team 2020).

Virulence assays
The virulence assays were conducted on beet armyworms (Spodoptera exigua, Lepidoptera).To assess interspecific variation in virulence, representative strains from species within PARB and M. anisopliae s. str.(Table 1) were subcultured on PDA at 25 °C for varying periods (1 week-1 month) until sporulation.These strains were selected as they represent the proposed novel species (M.hybridum, M. neoanisopliae, and M. parapingshaense) along with closely related known species (M.anisopliae s. str., M. pingshaense, and M. robertsii).For the assays on beet armyworms, spores from each strain were harvested into 1 ml of sterilized water, counted, and the concentration adjusted to 10 8 spores/ml.Afterwards, 3 μl of the spore suspension was injected into each beet armyworm.Thirty beet armyworms, divided into three replicates of ten individuals, were injected with the spore suspension.The mortality of the insects was monitored daily for one week.Two types of mortality data were considered: (1) unconditional mortality, i.e. any insect found dead was considered as non-survival, and (2) mortality with mycelia, i.e.only insects found dead and covered with fungal mycelia were considered as non-survival.This distinction aimed to disentangle different fungal fitness components as the fungi would be considered as having successfully reproduced only if they managed to develop through the insects' body to produce spores for further generations.The mortality rate was calculated for each replicate as the number of dead insects with or without covering fungal mycelia divided by the total number of insects per replicate (10).
To analyze the mortality rate data, a linear model was fitted using generalized least squares, i.e. gls command from the package nlme in R (R Core Team 2020), accounting for the correlation between the observations within replicates nested within strains.The species and days of observations were considered as interacting fixed effects and were tested using ANOVA.

Population structure and clonal diversity
The F84 distance-based neighbor-joining tree constructed from 98,085 SNPs revealed that the unidentified strains clustered with three Metarhizium species, M. anisopliae, M. pingshaense, and M. sulphureum (Fig. 1A).The M. anisopliae clade (blue clade, Fig. 1A) comprises known strains of what had been previously identified as M. anisopliae s. lat.(ARSEF 549) and M. anisopliae s. str.(ARSEF 2080, ARSEF 7487, ARSEF 7450) in the catalog of ARS Collection of Entomopathogenic Fungal Cultures (ARSEF).Surprisingly, the ex-neotype strain of M. anisopliae s. str.(CBS 130.71) did not cluster with this clade but branched as a basal taxon close to M. lepidiotae.Other unknown strains from our collections clustered either with known strains of M. pingshaense (pink clade, Fig. 1A), or M. sulphureum (brown clade, Fig. 1A).Other strains of known species each formed a monophyletic clade, corresponding to their distinctive taxonomic status.
We further conducted Bayesian clustering analyses focusing on the PARB group (M.anisopliae clade, M. pingshaense clade, M. robertsii, and M. brunneum: Fig. 1B, upper panel) and MGT group (M.sulphureum clade, M. guizhouense, M. gryllidicola, M. clavatum, M. phasmatodea, M. kalasinense, and M. majus: Fig. 1B, lower panel).The analyses confirmed the genetic differentiation between different species but revealed also strains/species which might have originated from hybridization; ARSEF 549, clearly showed a genetic mixed ancestry, probably between M. anisopliae s. lat.(= M. neoanisopliae) and M. robertsii or M. brunneum; ARSEF 2107, the type strain of M. brunneum, with the dominant genetic origin from M. robertsii.Metarhizium phasmatodea within the MGT group also exhibited a signal of mixed ancestry, probably between M. sulphureum and M. kalasinense/M.clavatum.There are two distinct genetic groups within the M. pingshaense clade (Fig. 2B, upper panel).
The 70 strains included in this study could be classified into 55 multi-locus genotypes (MLG), i.e. clonal groups (Additional File 1: Table S1).Besides ARSEF 549 with its mixed genetic ancestry, thus constituting on its own a clonal group (MLG 37), the M. anisopliae clade comprised six MLG (23-28).Metarhizium pingshaense clade is composed of eight MLG (29, 30, 33-36, 38, 39).Metarhizium gryllidicola and M. majus were respectively composed of a unique clonal lineage (MLG 47 for M. gryllidicola and 65 for M. majus) while all the remaining species were composed of multiple MLG, each represented by a single strain.Overall, this showed that some species were probably genetically more diverse than others.It is to be noted that the biopesticide strain BCC 4849 is part of MLG 24, represented also by four other strains (SM2387 = BCC96565, SM2402 = BCC96580, SM2398 = BCC96576, SM2400 = BCC96578) isolated from the soil of the same agricultural field.

Phylogenomics species tree
We constructed an ML-based phylogenomic tree from a matrix of 237 genes, partitioned to accommodate the best-fitted model for each gene (Fig. 2).In agreement with the clustering analysis, the ex-neotype culture of M. anisopliae s. str.(CBS 130.71) clustered with the strains known as M. lepidiotae, forming a strongly supported monophyletic clade branching from a deep node at the base of the PARB and MGT clades (Fig. 2), supporting the view that M. lepidiotae should be synonymized with M. anisopliae s. str.Other strains designated as M. anisopliae s. lat.formed a separate strongly supported monophyletic clade (i.e.M. neoanisopliae sp.nov.) while ARSEF 549 (M.hybridum sp.nov.) branched at the base of M. neoanisopliae and was shown above to have a mixed ancestry.These two taxa are thus considered to be different species.Two strongly supported monophyletic clades were observed for the strains considered as M. pingshaense s. lat.(Fig. 3).One clade contained the type strain of M. pingshaense (CBS 257.90) and so is considered M. pingshaense s. str.while the other clade is recognized here as a new species (M.parapingshaense sp.nov.).
The monophyletic clades supporting different species were also recovered in the coalescence-based species tree (Additional File 2: Fig. S1).The inference of the coalescent multi-species tree reconciliates the topologies of the 238 genes into the most probable species tree.This latter method corresponds to the phylogenetic species concepts based on the concordance between different gene trees (Taylor et al. 2000;Maharachchikumbura et al. 2021) and can outperform the concatenated genes-based approach (Liu et al. 2015).Yet, our results using this method confirm the species hypotheses formed with the concatenated-genes tree; all the clades supporting different species were supported with 100% posterior probability in the coalescent multispecies tree (Fig. 2).

Morphological analysis
Based on the measurements from axenic cultures (PDA), we statistically analysed whether clades/strains considered to be distinct but closely related species possessed different morphological features in terms of the width and the length of the phialides and conidia (Fig. 3).We observed a significant difference in the width of conidia between M. neoanisopliae and M. hybridum (F = 31.323,p-value = 0.002) but this was non-significant for the other traits (Additional File 1: Table S2).For the M. pingshaense s. lat.group, there was a significant difference in the length of phialides in culture between M. pingshaense (See figure on previous page.)Fig. 1 Results of population genomics analyses based on 98,085 SNPs.A F84 distance-based neighbor-joining tree showing the identification of unknown strains among three species complexes; colored labels represent strains previously recognized as respective species, M. anisopliae (blue), M. pingshaense (pink) and M. sulphureum (brown).B Bayesian clustering analysis focusing on the PARB group (upper panel), and MGT group (lower panel); MLG Multi-Locus Genotype (clonal group) Fig. 2 A ML-based phylogenomic tree from a concatenated matrix of 237 genes harnessed across the genome.Thick branches are supported by 100% bootstrapped, the red dots represent the nodes with 100% posterior probability following a coalescent-based species tree (Additional File 2: Figure S1) and M. parapingshaense (F = 6.421, p-value = 0.044) but not for the other traits (Additional File 1: Table S2).We also tested whether M. anisopliae s. str.would have different ranges of conidia and phialides width/length to the M. anisopliae s. lat. group (M.hybridum + M. neoanisopliae); surprisingly, none of the traits were found to be significantly different between these two groups (Additional File 1: Table S2).We found that these traits highly overlapped between different Metarhizium groups (Fig. 4).We also obtained measurements from a dry specimen each from M. pingshaense (MY5150) and M. parapingshaense (MY12337) (Additional File 1: Table S2).There was a significant difference in the width (F = 28.583,p-value = 8.676e-07) and the length (F = 96.174,p-value = 2.968e-15) of conidia between these two strains, each representating different species.

Virulence
Our statistical analysis revealed a significant difference in both the unconditional mortality (UM) and the mortality with mycelia (MM) between the species within the PARB group (UM: F = 11.140,p-value < 0.0001; MM:  4).Considering the pairwise difference between species, M. anisopliae s. str.were almost always significantly less virulent than all species within the PARB group for both types of mortality (Additional File 1: Table S3).There was no significant difference between neither M. parapingshaense and M. pingshaense (UM: t = -1.588,p-value = 0.567; MM: t = -0.070,p-value = 1.000), nor M. neoanisopliae and M. hybridum (UM: t = 1.148, p-value = 0.755; MM: t = -1.708,p-value = 0.528), respectively.Interestingly, a notable variation of mortality between strains within species was observed (Additional File 3: Fig. S2).

Metabolomics
The PCA on the data of peak area of selected features showed that, for both the broth extract (BE) and cell extract (CE), most of the samples had the distribution of secondary metabolites overlapped without any clear difference between species (Additional File 4: Figure S3).When categorizing different taxa into species complexes based on the phylogenomics results: PARB-M.brunneum, M. hybridum, M. neoanisopliae, M. parapingshaense, and M. pingshaense; Anisopliae Strict (M.anisopliae s. str.);MGT-M.majus, M. guizhouense, M. clavatum, M. gryllidicola, M. kalasinense, M. phasmatodea, and M. sulphureum; Acridum-M.acridum and M. globosum; Flavoviride-M.flavoviride and M. frigidum; no discriminating pattern could be observed between these species complexes (Additional File 5: Figure S4).However, by narrowing the data to only the PARB group, in which we propose three new species and M. anisopliae s. str., we could observe for CE a segregation, to an extent, between M. pingshaense and M. neoanisopliae.However, M. parapingshaense and M. neoanisopliae largely overlapped respectively with M. pingshaense and M. anisopliae s. str.(Fig. 5A).Interestingly, for BE, M. neoanisopliae clearly segregated from M. anisopliae s. str.while M. parapingshaense and M. pingshaense were found to overlap largely between them (Fig. 5B).
Most of the identified known compounds produced by our fungal strains belong to various secondary metabolite families such as cytochalasins, destruxins, metarhizins, subglutinols, and metacytofilins (Fig. 6, Additional File 6).By focusing on the annotated compounds between PARB and M. anisopliae s. str., the Manhattan distancebased heatmap showed that M. anisopliae s. str.tended to closely group, and produce to a lesser extent various destruxins compared to the strains from the PARB group (Fig. 5C).
Notes: Our phylogenomic analyses indicated that the ex-type strain M. lepidiotae ARSEF 7488, and ARSEF 7412, cluster with the ex-neotype M. anisopliae CBS 130.71 from Ukraine.Therefore, this species is synonymized with M. anisopliae s. str.The samples of M. anisopliae s. str.showed greater variability in conidial length and width than M. neoanisopliae sp.nov.and M. hybridum sp.nov., with a significant portion of conidia being longer (> 10 µm) and larger (> 3.5 µm) than M. hybridum and M. neoanisopliae.The virulence of M. (Fig. 8).MycoBank MB 850069.Etymology: The species name is derived from the genetic mixed ancestry inferred from genomic analyses.
Diagnosis: Metarhizium hybridum is highly similar to its sister species M. neoanisopliae sp.nov., but has conidial width statistically smaller than the latter.The conidial width is always less than 3.5 µm, and the conidial length, like its sister species M. neoanisopliae, is always shorter than 8.5 µm (reserved for M. anisopliae).It is also more virulent against Spodoptera exigua than M. anisopliae s. str.Type: Brazil, isolated from unknown source, collection date unknown, D. W. Robert, (BBH50656 -holotype preserved in a metabolically inactive state; ARSEF 549 -culture ex-type).
Diagnosis: Metarhizium neoanisopliae is similar to its sister species M. hybridum sp.nov.but the conidial width is statistically larger than the latter.The conidial width is never below 2 µm, and can be beyond 3.5 µm, but never surpasses 4 µm as in M. anisopliae.The conidial length, as in its sister species M. hybridum, is never longer than 8.5 µm.It is more virulent against Spodoptera exigua than M. anisopliae s. str.
Etymology: The epithet refers to the very close phylogenetic position to Metarhizium pingshaense, a sister species evolving in parallel from a common ancestor.
Diagnosis: Metarhizium parapingshaense is very similar to its sister species, M. pingshaense, but differs in conidial size.The conidia on specimens are statistically larger for M. parapingshaense (mostly > 6 µm in length and frequently > 2.5 µm in width).The phialides produced on PDA are also statistically shorter than in M. pingshaense, with the phialide length never over 12.5 µm (as occurs in the latter).

Discussion
Hidden Diversity within the Metarhizium anisopliae species complex Our population genomics-based phylogenetic results agree with previous taxonomic treatments of the genus Metarhizium (Bischoff et al. 2009;Kepler et al. 2014;Mongkolsamrit et al. 2020) in the sense that the M. flavoviride complex is a sister clade to the M. anisopliae species complex and that, within the M. anisopliae complex, M. acridum and M. globosum form a sister clade to the rest of the complex.We recognized three new species within the M. anisopliae complex, namely M. neoanisopliae, M. hybridum, and M. parapingshaense.Metarhizium neoanisopliae and M. hybridum are distinguished from M. anisopliae s. str.based on multiple lines of evidence.First, the genomics data segregated M. neoanisopliae and M. hybridum from M. anisopliae s. str.Second, the virulence data showed M. anisopliae s. str. to be less virulent than M. neoanisopliae and M. hybridum in laboratory conditions.Finally, the metabolomics analysis also revealed the differences in secondary metabolite production between the two former taxa and the latter.These support the distinct species status of M. neoanisopliae and M. hybridum from M. anisopliae s. str.
The original description with illustration of M. anisopliae by Metchnikoff (1879) is in Russian, and currently not accessible.The oldest account of the morphology of this species, that we could find, was that of Delacroix (1893), who stated that he had examined specimens of Metchnikoff and described the length of M. anisopliae's conidia as 7-15 µm.Veen (1968) and Tulloch (1976) made a reference to Metchnikoff 's description of M. anisopliae as having conidia of 4.8 µm long and 1.6 µm wide.The examination of putative M. anisopliae strains by Veen (1968) and Tulloch (1976) resulted in the conidia length being respectively at 4.6-11.5 µm and 3.5-9 µm.These measurements are smaller than those of Delacroix ( 1893) and of M. anisopliae s. str.as interpreted in our study, but more within the range of M. neoanisopliae.Overall, the conidial dimension of M. anisopliae s. str.with CBS 170.71 as the neotype (Mongkolsamrit et al. 2020) better fits the account of Delacroix (1893).Taken with the fact that CBS 170.71 came from the same original locality of Metchnikoff 's species, it is justifiable to accept this isolate as the neotype of M. anisopliae s. str.(Mongkolsamrit et al. 2020).This reclassification means that ARSEF 7450 and many strains previously identified as M. anisopliae should be reidentified as M. neoanisopliae.
With CBS 170.71 accepted as the ex-neotype culture of M. anisopliae s. str., M. lepidiotae should be considered an objective synonym of M. anisopliae s. str.This proposition is supported by the clustering of the ex-type strain of M. lepidiotae (ARSEF 7488) with CBS 170.71.Initially described as a variety of M. anisopliae by Driver et al. (2000), M. lepidiotae was later elevated to species rank by Bischoff et al. (2009), primarily based on molecular phylogeny.Notably, the spore dimensions of M. lepidiotae (conidia 7.3-10.6× 3-4.1 µm; Driver et al. 2000) falls within the range of M. anisopliae s. str.
The difference between M. neoanisopliae and M. hybridum principally relies on genomic data.As M. hybridum clustered next to M. neoanisopliae and was represented by a single strain, one might argue that M. hybridum should be considered as M. neoanisopliae.However, M. hybridum clearly demonstrated a genomic signature of mixed ancestry, which is not the case for M. neoanisopliae.Furthermore, the conidial width of these two species is significantly different, those of M. hybridum being slightly narrower than those of M. neoanisopliae.Although this difference is not easily perceptible by eye, it is statistically significant.
The distinction between M. parapingshaense and M. pingshaense is supported by the genomic data.This showed them to form closely related but distinct clades.In addition, morphological data show that the former has statistically larger on-specimen conidia than the latter, although some individual spores have an overlapping size.The virulence between these two species is at the same level and, like M. neoanisopliae and M. hybridum, higher than that of M. anisopliae.Overall, it is well supported that M. anisopliae sensu stricto is different to the novel species proposed in this study.

Evolution of virulence and host specificity of the Metarhizium anisopliae complex
The phylogeny of the M. anisopliae complex shows that M. acridum and M. globosum formed a sister clade to the other species.Metarhizium acridum is well documented as a pathogen specific to Orthoptera (Wang and Leger 2005;Hu et al. 2014) choff et al. 2009;Hu et al. 2014).The close position of M. acridum to the most recent common ancestor (MRCA) of the M. anisopliae complex have led some authors to formulate the generalization that it "evolved from specialists" via "transitional species with intermediate host ranges" (Hu et al. 2014).Without going into the controversy of whether such a statement is accurate, we can accept that the M. anisopliae complex comprises two distinct groups, one formed by M. acridum and M. globosum which are specialists, and the other consisting of ones with intermediate to very broad host ranges.Therefore, there are distinct evolutionary trajectories into different levels of host specificity.This situation is in contrast with the genus Beauveria, in which most of the species were described as broad generalists; many strains, including the types, were found associated with Coleoptera, the most species-rich order of insects, with some species appearing to become more specialized on other arthropod groups such as Orthoptera (B.gryllidicola, B. namnaoensis), Lepidoptera (B. pseudobassiana, B. thailandica) or Acari (B.varroae) (Kobmoo et al. 2022).The specialization to other insects besides Coleoptera appeared independently multiple times in that genus.
As our study only focused on the virulence of Metarhizium spp.towards an insect (Spodoptera exigua), this does not have the power to refute any hypothesis regarding host specificity, but it did show that the variation of virulence could be partly explained by evolutionary history.Metarhizium anisopliae s. str.(syn.M. lepidiotae), which phylogenetically branched directly from the MRCA of the M. anisopliae complex, is clearly less virulent than other species considered to be M. anisopliae s. lat.(PARB group; i.e. M. hybridum, M. neoanisopliae, M. parapingshaense, M. pingshaense, and M. robertsii).It is tempting to hypothesize that the M. anisopliae complex evolved from a moderately entomopathogenic common ancestor to become highly entomopathogenic.The difference of virulence between M. anisopliae s. str.and other species of the complex has an implication in the development of biocontrol strategies as many of the Metarhizium strains used in biocontrol have been putatively identified as M. anisopliae based solely on ITS (Dong et al. 2016;Alam 2019;Ahmed et al. 2020;Qubbaj and Samara 2022) while they might actually represent different species with different virulence patterns.The difference in virulence is probably explained by the PARB group species which produce more specific destruxins and demethyldestruxins as shown by our metabolomic analyses.Destruxins are well known for their insecticidal activity (Pedras et al. 2002;Wang et al. 2004;Pal et al. 2007).Our data warrant further research on specific cyclic peptides produced by highly virulent strains.
Metarhizium is known to be ecologically versatile (St. Leger and Wang 2020), capable of being soil-borne, associated to plant rhizospheres (Vega et al. 2009), and even occasionally pathogenic to vertebrates (Horgan et al. 2022).Entomopathogenicity can be viewed as an ecologically specialized function within the genus.A lesser entomopathogenic MRCA of the M. anisopliae complex would support this hypothesis.However, it is necessary to gain more virulence data across a larger panel of strains from the M. anisopliae complex (e.g.MGT group) in order to accurately construct an actual ancestral reconstruction.

Persistence of a biocontrol candidate strain in field application
Some strains of M. neoanisopliae were isolated from the soil of a fruit orchard where BCC 4849 had been applied as biocontrol agent against insect pests two years earlier.Some of these strains were inferred to be of the same clonal lineage as BCC 4849.This shows that this candidate strain can persist in the field for several years.Similarly, a few studies have showcased the persistence of biocontrol strains of Beauveria (Mei et al. 2020) and Metarhizium (Peng et al. 2021) in the field many years after their initial applications.These findings, with ours, are encouraging as they show that biocontrol strains could be sustainably used in long-term efforts to reduce the utilization of chemical pesticides and ensure better food security.

Conclusion
We combined multiple lines of evidence in resolving the taxonomy of a species complex.Fungal taxonomy in the twenty-first century has increasingly relied mostly on molecular data and the phylogenetic species concept (Taylor et al. 2000).Particularly for Ascomycota, molecular phylogenetics provides the basis to delineate evolutionary lineages that frequently do not match other characteristics (Dugan and Everhart 2018;Balasundaram et al. 2015;Leavitt et al. 2016;Mongkolsamrit et al. 2018;Kobmoo et al. 2021;Steenwyk 2023).We have demonstrated here that evolutionary lineages based on wholegenome sequencing can be used to seek for differences in other aspects of fungal biology.The metabolomics and virulence data confirmed the difference between M. anisopliae s. str.and the rest of the species complex while the morphometrics analysis revealed fine-scale difference between closely related species.This work has also taken advantage of isolates stored in various institutional collections, allowing the clarification of the phylogenetic positions of ex-type cultures, and demonstrated the importance of integrating a vast panel of specimens across different institutions.Institutional collections are rich sources of undiscovered or misidentified species that can only be uncovered using multiple lines of evidence (molecular, genomics, metabolomics, ecology).Such integrative approach based on inter-disciplinary and inter-institutional collaborations can greatly benefit future taxonomic work of fungi.
polymorphisms data under the VCF format is available at Dryad data repository (10.5061/dryad.q83bk3jr8).The fungal type materials and type strains of the new species proposed in this study were deposited at BIOTEC Bangkok Herbarium (BBH) and BIOTEC Culture Collections (BCC) with the accession numbers provided in Table 1.

Fig. 3
Fig. 3 Notched boxplots representing the distribution of the width and the length (μm) of phialides and conidia between the novel taxa and closely related species from (A) cultures, and (B) stromatal samples.Statistical significance (one-factor ANOVA): * = p < 0.05, *** = p < 0.001

Fig. 4
Fig. 4 Virulence of Metarhizium species within the PARB group (M.anisopliae sensu stricto, M. hybridum, M. neoanisopliae, M. parapingshaense, M. pingshaense, M. robertsii) along 7 days of observation.A The mortality with mycelia; only dead insects covered by fungal mycelia and spores were counted.B All mortality; all dead insects including those manifesting no fungal material on the external surface were counted.The error bars represent the standard errors of respective data points

Fig. 5
Fig. 5 Metabolomic analyses based on peak area data obtained from liquid chromatography-mass spectrometry (LC-MS) for Metarhizium species of the PARB group.A Principal component analysis (PCA) of data from cell extracts.B PCA of data from broth extracts.The ellipsoids delimit the perimeter of the sample distributions.C A heatmap based on Manhattan distance clustering for annotated compounds.

Fig. 7
Fig. 7 Metarhizium anisopliae sensu stricto: CBS 130.71, A colony obverse on PDA after 2 wk.B colony reverse on PDA after 2 wk.C Phialides and conidia on PDA.(D) conidia on PDA.; ARSEF 7412, E colony obverse on PDA after 2 wk.F colony reverse on PDA after 2 wk.G Phialides and conidia on PDA.H conidia on PDA.; ARSEF 7488 (M.lepidiotae culture ex-type), I colony obverse on PDA after 2 wk.J colony reverse on PDA after 2 wk.K Phialides and conidia on PDA.L conidia on PDA.-Scale bars: b − c = 10 mm; d − f = 10 µm

Table 1 .
List of Metarhizium spp.included in this study ( T = ex-type strain, M = included in the morphometric analysis, C = included in the metabolomics analysis, V = included in the virulence assays).Culture collection and fungarium codes: BBH BIOTEC Bangkok Herbarium, BCC BIOTEC Culture Collection, CBS Fungal and yeast collection of WI-KNAW (CBS-KNAW) Culture Collection, ARSEF U.S. Department of Agriculture (USDA) Culture Collection

Table 1 .
(continued) Metarhizium parapingshaense is separated from its sister species M. pingshaense based on genetic segregation, following the phylogenetic species concept, from genome-wide polymorphisms.The conidia on specimens are also statistically larger and longer than the latter.