Phenotypic variation and genomic variation in insect virulence traits reveal patterns of intraspecific diversity in a locust‐specific fungal pathogen

Intraspecific pathogen diversity is crucial for understanding the evolution and maintenance of adaptation in host–pathogen interactions. Traits associated with virulence are often a significant source of variation directly impacted by local selection pressures. The specialist fungal entomopathogen, Metarhizium acridum, has been widely implemented as a biological control agent of locust pests in tropical regions of the world. However, few studies have accounted for natural intraspecific phenotypic and genetic variation. Here, we examine the diversity of nine isolates of M. acridum spanning the known geographic distribution, in terms of (1) virulence towards two locust species, (2) growth rates on three diverse nutrient sources, and (3) comparative genomics to uncover genomic variability. Significant variability in patterns of virulence and growth was shown among the isolates, suggesting intraspecific ecological specialization. Different patterns of virulence were shown between the two locust species, indicative of potential host preference. Additionally, a high level of diversity among M. acridum isolates was observed, revealing increased variation in subtilisin‐like proteases from the Pr1 family. These results culminate in the first in‐depth analysis regarding multiple facets of natural variation in M. acridum, offering opportunities to understand critical evolutionary drivers of intraspecific diversity in pathogens.

high-virulence isolates (Gu et al., 2020;Valero-Jiménez et al., 2016;Zhang et al., 2020).Such isolate-specific knowledge on variation in virulence patterns of pathogens is important for understanding the evolutionary ecology of host-pathogen interactions.
The genus Metarhizium Sorokīn (Ascomycota: Sordariomycetes) is composed of globally distributed fungi with important functions from both ecological and economic perspectives.Although species within this group can switch between various fungal lifestyles, growing as saprophytes or root symbionts, they are primarily known for their role in the regulation of insect populations as entomopathogens, across the spectrum from generalist to specialist pathogens (St. Leger & Wang, 2020).The exploitation of this trait in agricultural systems has led to the widespread use of several species as biological control agents, including examples such as spittlebug control on sugarcane crops in Brazil (Li et al., 2010), protection from locust swarms in several African countries (Douthwaite et al., 2001), and genetic modification of strains to target mosquito vectors of malaria (Lovett et al., 2019).In order to improve the effectiveness of Metarhizium in an agricultural context, it is important to consider variable genotypes that generate locally adapted biological and pathogenic diversity.
Similar to other species within this genus, M. acridum has a cosmopolitan distribution, albeit one that is restricted to sub-tropical and tropical regions (Mongkolsamrit et al., 2020), likely due to coevolution with its orthopteran hosts (Brancini et al., 2019).The infection process is comparable to that of other cuticle-penetrating fungal entomopathogens, following the standard six stages of conidial adhesion, germination, appressorium formation, penetration, colonization of haemolymph, and sporulation (Aw & Hue, 2017).However, specialization of M. acridum to hosts within Orthoptera is considered to have led to specific patterns of genetic and pathogenic adaptation.For example, unlike generalist Metarhizium species, which typically kill hosts of many insect orders quickly via toxins, M. acridum causes a systemic infection of host tissues along with increased evasion of the host immune system, eliciting a more gradual killing of their orthopteran host (Kershaw et al., 1999;Li & Xia, 2022).In addition, M. acridum has been shown to be more tolerant to physical stressors than other Metarhizium species, such as to ultraviolet light from solar radiation (Brancini et al., 2019;Fernandes et al., 2010).In terms of genomic differences, M. acridum has more rapidly evolving genes than generalist species, along with a decreased repertoire of G protein-coupled receptor-related proteins and proteolytic enzymes, which is in concordance with a more specific response to environmental triggers and thus a more limited host range (Hu et al., 2014).
Additionally, non-specialist species in the genus Metarhizium are haploid and are composed of a small number of clonal lineages with low rates of recombination (Rehner, 2020).In contrast, M. acridum along with another specialist, M. album, show signs of frequent and recurring sexual reproduction based on fungal-specific genomic signatures.For example, specialists have fewer heterokaryon incompatibility (HI) proteins, which work to inhibit vegetative reproduction cycles by rejecting hyphal fusion between incompatible individuals.Specialists also show evidence for repeat-induced point (RIP) mutations.In fungi, RIP is a genome defence mechanism that helps remove deleterious single nucleotide point mutations, which is only active during meiosis and leaves a distinct base-pair frequency signature in the genome (Selker, 2002).The presence of a RIP-signal in fungal genomes and inferred sexual recombination has been shown to be associated with rapid sequence diversification and specialization (Kemen & Jones, 2012;Ma et al., 2010).Overall, M. acridum has diverged significantly in phenotypic traits and genomic patterns from its generalist counterparts, but the intraspecific diversity of this taxon is still largely unknown, despite its widespread distribution.In consideration of the diversity of its orthopteran hosts, with over 26 000 extant species that have adapted to multiple habitat types, it is conceivable that populations of M. acridum may have adapted to co-occurring host species (Song, 2018).For such host-pathogen coadaptation to occur, we would expect phenotypic differences in virulence and diversity in functional genomic regions associated with virulence.
The focus on these two isolates has had notable benefits for pest management, considering that both are mass produced for largescale locust control (Douthwaite et al., 2001;Peng et al., 2008) and have been whole-genome sequenced (Gao et al., 2011;Nielsen et al., 2021).Within the genus Metarhizium, intraspecific diversity has mostly been analysed with multiple fungal barcode loci and amplified fragment length polymorphism (AFLP) data, which has shown non-clonal diversity within isolates of the broad-host-range species Metarhizium flavoviride W. Gams & Rozsypal found in Denmark (Keyser et al., 2015), as well as phylogenetically distinct clades of both Metarhizium brunneum Petch (Steinwender et al., 2014) and Metarhizium robertsii J.F. Bisch., S.A. Rehner, and Humber (Kepler et al., 2015).Based on a limited number of dominant AFLP markers and single nuclear gene sequences, there is evidence of intraspecific genetic variation and population structure of globally distributed M. acridum isolates (Fernandes et al., 2010).The aim of the current study was to investigate whether phenotypic and genome-wide genotypic variation among globally collected M. acridum isolates is in concordance with a pattern of local adaptation of M. acridum isolates.Specifically, we asked: (1) does phenotypic variation in virulence and in vitro growth rate differ between M. acridum isolates, (2) do M. acridum isolates differ in virulence between orthopteran host species and growth rate on different nutrient sources, and (3) is genetic diversity between isolates structured in genomic regions harbouring genes associated with virulence and/or nutrient acquisition (Figure 1).Using a set of nine M. acridum isolates from four continents representing most of the known worldwide geographic distribution, we evaluated our three questions concerning isolate variation in virulence, growth rate, and genotypes in three steps.First, we measured the virulence of all nine M. acridum isolates against two different species of locusts (Schistocerca gregaria Forsskål, 1775, and Locusta migratoria migratoria Linnaeus, 1758).Second, we measured the growth rate on three different in vitro media composed of different nutrient sources.Third, we used whole-genome sequencing and comparative genomics to identify variable genomic regions among the nine M. acridum isolates.

| Fungal isolates
Nine isolates of M. acridum were obtained from the USDA-ARS Collection of Entomopathogenic Fungal Cultures (Ithaca, NY, USA;  S1).Cultures were grown from single spores of isolates revived from a glycerol stock stored at −80°C.All initial cultures were grown on one-quarter strength Sabouraud dextrose agar (SDAY/4: 2.5 g L −1 peptone, 10 g L −1 dextrose, 2.5 g L −1 yeast extract, 20 g L −1 agar) for 14 days in darkness at 23 ± 1°C.Conidia were harvested with a 0.05% Triton X solution to make conidia suspensions that were either used for DNA extraction or virulence assays.The concentration of each suspension was estimated by counting in a haemocytometer (Fuchs-Rosenthal).
Both types of insects were acquired from a commercial supplier (Monis, Viborg, Denmark) as second-instar nymphs and maintained at 25 ± 1°C in 30×30×60 cm glass cages containing 100-300 individuals at a time, with a photoperiod of 14 h light, 10 h dark, while being supplied with fresh lettuce.Juveniles at the third instar were used for experiments, with 10 insects per treatment.In each bioassay, F I G U R E 1 Schematic of analyses for phenotypic and genotypic characterization of nine M. acridum isolates.CNVs, copy number variants; Locust, minimal media with whole crushed locust tissue added; SDAY, Sabouraud dextrose agar with yeast; SNPs, single nucleotide polymorphisms; Trehalose, minimal media with trehalose as carbon source.
3 μL of conidia suspension with a concentration of 1 × 10 6 conidia/ mL was applied beneath the pronotum.Exposed insects were placed individually in 450-mL plastic cups with a ventilated lid and supplied daily with fresh lettuce.All conidial suspensions had a germination rate between 90% and 100% at the start of each bioassay.The exposed locusts were maintained in an incubator at 30 ± 1°C, with a photoperiod of 14 h light, 10 h dark, for 10 days, while being checked daily for mortality.Cadavers were surface sterilized by submerging the insect body in 0.2% NaOH for 10 s, followed by two rinses in sterile water for 10 s each.Cadavers were placed in parafilmwrapped, sterile Petri dishes with a diameter of 5.5 cm, to promote fungal sporulation.The presence or absence of emergent green conidia was recorded for each cadaver.
Locust survival was visualized using the Kaplan-Meier survival analysis (Kaplan & Meier, 1958) as implemented in the survival v3.3-1 (Therneau, 2018) and survminer (Kassambara et al., 2019) packages in R. Statistically significant differences in the mortality of locust species infected with the same isolate were determined with logrank tests.Multivariate modelling of the mortality of fungal isolate treatments was performed using a mixed-effects Cox proportional hazards regression (coxme) model in R, with experimental replicates designated as a random factor (Cox, 1972;Therneau, 2022).In this model, the risk of death is expressed as a hazard ratio, where a value of 1.0 indicates two treatments have the same probability of mortality at any given time.Pairwise comparisons in the coxme model were performed by setting each fungal treatment as the reference and calculating the hazard ratios.Treatments were considered significantly different if the p-value was less than 0.05.

| Media growth rate experiment
The mycelial growth rate, a common proxy for measuring fitness in filamentous fungi (Pringle & Taylor, 2002;Valero-Jiménez et al., 2017), was measured on three types of media representing different carbon sources.The fungi were grown on: (1) Sabouraud Dextrose Agar (SDAY) medium, composed of 16.25 g Sabouraud Dextrose Agar (Sigma-Aldrich), 2.5 g yeast extract, and 11.25 g agar in 1 L of distilled water; (2) trehalose medium, composed of 40 g trehalose, 10 g peptone, 10 g yeast extract, and 15 g agar in 1 L of distilled water; (3) locust minimal medium, composed of 0.4 g KH 2 PO 4 , 1.8 g Na 2 HPO 4 •H 2 O, 0.6 g MgSO 4 •7H 2 O, 1 g KCl, and 0.7 g KNO 3 supplemented with 15 g of agar and 20 g of freeze-dried, whole adult locusts (L.migratoria migratoria) crushed into a fine powder, in 1 L of distilled water.Colony morphologies and growth rates were examined by inoculating 3 μL of a 1 × 10 7 conidia/mL suspension in the centre of a 9-cm-diameter Petri dish containing 20 mL of media.
Plates were maintained at 23 ± 1°C with a photoperiod of 24 h darkness.Plates were measured for colony radius manually at 3, 6, 9, and 12 days.Experiments were performed in triplicate, with three plates of each medium type in each replicate, resulting in a total of nine observations of growth rate for each isolate on every treatment type.
To test for differences in growth rates of the nine isolates, linear mixed-effect models (LMMs) were constructed with the lmer function in the R package lme4 v1.1-29 (Bates et al., 2015).Growth rate values were checked for homoscedasticity and a normal distribution prior to analysis.LMMs for each type of media were constructed with growth rate of each plate as the response variable, the isolate as the explanatory variable, and replicates considered as random effects.Pairwise comparisons of estimated marginal means were calculated with the R package emmeans v1.8-5 (Lenth et. al., 2019), with all post hoc comparisons adjusted using Tukey's HSD.

| DNA extraction and genome sequencing
For DNA extraction, 20 mL of Sabouraud Dextrose Broth (SDB: 2.5 g peptone, 10 g dextrose, 2.5 g yeast extract, 10 g agar in 1 L of distilled water) was inoculated with conidia and incubated on a shaker at 25 ± 1°C at 130 rpm for 2 days.Excess liquid was removed by vacuum filtration, and the resulting fungal mycelia was washed in sterile, de-ionized water, before being flash-frozen in liquid nitrogen.DNA was isolated from all isolates using a CTAB-based extraction method (Möller et al., 1992).Paired-end, 150-bp sequencing of barcoded DNA libraries was performed on the BGI-seq 500 platform (BGI).
Initial quality filtering and adaptor removal was performed by BGI.Paired-end reads from each of the nine M. acridum isolates were subsequently assessed for remaining poor-quality sequence content using FASTQC v0.11.8 (Andrews, 2010).Reads were filtered using fastq-mcf v1.04.636 set with default parameters to trim low-quality terminal ends (≤Q10), discard short sequences (<50 bp), and quality filter for mean Phred scores (Q < 20).Filtered reads were then mapped to haplotype 1 of the M. acridum reference genome for isolate ARSEF 324 (Nielsen et al., 2021), using BWA v0.7.15 (aln algorithm), under default settings.PCR duplicates were subsequently removed from mapped reads using samtools v1.10 (markdup -r; Li et al., 2009), and the average coverage of all reads to the reference assessed using samtools.The coverage of each 1 kb region was determined by merging mapped reads from all nine samples, and regions with excess (>10 000) read depth were removed from variant analysis after visual assessment.Regions with low read depth were retained because subsequent filtering steps use depth of read-mapping information to remove unreliable SNPs in the variant calling pipeline.

| Phylogenetic analysis
In order to confirm the species-level identity of the 9 isolates used in this study, 6 different loci were used to create a multilocus phylogeny (ITS, RPB1, RPB2, SSU, LSU, and tef1).Whole nucleotide sequences of M. anisopliae (CBS 130.71;Mongkolsamrit et al., 2020) loci were used as BLASTn queries against the reference genome used in this study (ARSEF 324).After identifying the location of these loci, the 'intersect' module contained bedtools v2.28.0 was used to extract reads from each isolate that mapped to a particular locus.The 'consensus' module within bcftools v1.9 was subsequently used to create a consensus sequence of each locus for each of the 9 isolates.A maximum-likelihood analysis was performed using RAxML v. 8.2.12 (Stamatakis, 2014), as implemented on the CIPRES portal (Miller et al., 2010), with default settings.The GTR + G model of evolution was used for each partitioned dataset consisting of the 6 distinct loci, with bootstrapping for 1000 pseudoreplicates, to search for the best ML tree.

| Genome-wide diversity in M. acridum
Following read alignment, variants were called using the Genome Analysis Toolkit (GATK v4.2.3.0;DePristo et al., 2011).Within the GATK pipeline, the 'Haplotype Caller' module was employed with the 'ploidy -1' setting to account for haploid genomes.This module was used to call high-quality variants between the reads and the reference assembly.We also called variants from ARSEF 324, an isolate identical to the reference.For only this isolate, we used the setting, 'ploidy -2', as this isolate was known to be a diploid.Genotype likelihoods were computed in gVCF mode at all loci for the 9 samples individually, before joint genotyping across all the samples.Single nucleotide polymorphisms (SNPs) and insertions/deletions (indels) were then filtered separately in GATK using the recommended parameters by GATK developers as follows: Variant Filtration, QD <2.0, FS > 60, MQ < 40, MappingQualityRankSum < −12.5, and ReadPosRankSum < −8.0.To remove SNPs from assembly errors in the reference isolate (ARSEF 324), we removed all regions where ARSEF 324 contains low read depth (<25) or heterozygous SNPs.These likely represent artificial variation resulting from assembly errors in the reference or incorrectly mapped reads from regions only present in haplotype 2 of the reference genome, which was not used for mapping.
An alignment of conserved polymorphisms was obtained by isolating SNPs only present in coding regions of the genome using vcf2fasta.This subset was used to create a NeighborNet phylogenetic network based on uncorrected p-distances as implemented on SplitsTree4 v.4.18.2 (Huson, 1998).A phylogram showing the phylogenetic relationships between the original orthopteran hosts of the isolates that are used in this study is also visualized based on a previous study by Song et al. (2020), for comparison to the phylogenetic network of M. acridum isolates (Figure 2c).For the genomewide clustering analysis, all variants called from the GATK pipeline were used for building a Hamming distance matrix based on pairwise comparisons between the samples.The final distance matrix was visualized in R with the heatmap.2function from the gplots package (Warnes et al., 2005).
In order to further understand intraspecies genetic diversity, we assessed nucleotide diversity (π) and Tajima's D with variants called by GATK using the R package PopGenome v.2.7.5 (Pfeifer et al., 2014).These population genetic statistics were calculated in sliding windows of 1000 bp with no overlap.In quantitative analyses of regions of high nucleotide diversity, a threshold of >3 SD is typically considered sufficient (Baltrušis et al., 2022;Bellis et al., 2016).
However, due to our limited sample size, regions were considered as outliers with high nucleotide diversity if they had at least three consecutive windows with extreme (>3 SD from the mean) values of π.Additionally, regions were conservatively considered to exhibit low-frequency polymorphism if they had Tajima's D values below −2.0.Traditionally, values of Tajima's D below −2.0 for re-sequencing data are significant under the simplest models (Tajima, 1989) and have been shown to be inconsistent with neutrality under multiple demographic models, including exponential expansion, bottlenecks, island structure, and population splitting (Akey et al., 2004;Carlson et al., 2005;Fahrenkrog et al., 2017;Stajich & Hahn, 2005).As a complementary approach, we also assessed the effect of nucleotide variation among the nine isolates by inferring piN/piS ratios for each gene (longest concatenated exons for each predicted gene model) using the Nei-Gojobori algorithm as implemented in SNPGenie (Nelson et al., 2015).In contrast to the traditional dn/ds ratio, piN/piS is used to describe rates of non-synonymous to synonymous mutations in intraspecific populations.We extracted genes that had a piN/piS >1.2, as rates greater than 1 are typically indicative of positive selection.We filtered out genes with piS <0.01, since low values of synonymous mutations can cause unreliable piN/piS estimation.
We also searched for genes that exhibited no non-synonymous mutations.We only retrieved genes that were found in regions with high nucleotide diversity (π; μ >3 SD) or low Tajima's D (<−2.0) as defined above, or with the aforementioned piN/piS values, using bedtools intersect (Quinlan & Hall, 2010).

| Gene presence-absence polymorphism calling
To characterize copy number variants among the isolates, we applied CNVpytor v.1.2.1 (Suvakov et al., 2021), a normalized read depth (RD)-based algorithm based on an updated version of CNVnator (Abyzov et al., 2011).BAM files of each isolate mapped to the reference genome (ARSEF 324) were used as input for analysing read depth.Following standard recommendations, we performed the analyses in 500-bp bins.Raw CNV calls were filtered on the following criteria: duplications RD >1.8, deletion RD <0.4,filtered calls with more than half not uniquely mapped reads, and retained calls larger than 1000 bp.We retrieved genes affected by CNVs using bedtools intersect (Quinlan & Hall, 2010) and subsequently performed a gene set enrichment analysis as described below.

| Enrichment analysis of gene ontology (GO) terms
Using the gene annotation for the reference genome ARSEF 324 (Nielsen et al., 2021), genomic regions of interest were used as coordinates in BEDTOOLS v2.30.0 (Quinlan & Hall, 2010).Gene ontology (GO) enrichment analyses were performed using the R package GOstats v.2.62.0 (Falcon & Gentleman, 2007).Original p-values of enriched GO terms were adjusted for multiple testing using the Benjamini-Hochberg (i.e.False Discovery Rate) approach, with the R package stats v4.2.0.Genes with GO terms that exhibited functional enrichment were compared to similar proteins in published Metarhizium genomes using both blastn for further annotation of gene function.Nucleotide alignments with SNPs were translated to protein sequences in Geneious v.2023.0 to determine the number of SNPs, number of amino acid differences, and pairwise identity in proteins.

| Intraspecific genetic diversity
We sequenced 9 isolates of M. acridum with 150-bp paired-end reads using the BGI-seq 500 platform and mapped the reads to haplotype 1 of the diploid reference genome sequence (ARSEF 324), with coverage for each isolate across the entire genome ranging between 41× and 51× (Table 1).Sequence yield and proportions of mapped reads varied between isolates, with the percentage of mapped reads ranging between 92% and 96%.One of the isolates, ARSEF 324, is identical to the diploid reference used for mapping, although most discrepancies represented by SNPs are likely because only one haplotype was used as reference for mapping.All other isolates were presumed to be haploid based on kmer analyses of reads and because the distribution of SNP allele frequencies clustered around 1, whereas the allele frequency of the diploid was clustered around 0.5 (Figure S1).
Variant calling with GATK produced a total of 387 611 SNPs and 51 536 indels among all the isolates following hard filtration of lowquality sites.Of these, 92.8% of SNPs were classified as biallelic, and 15.5% of SNPs were found in coding regions.The vast majority of SNPs (79%) represented transition mutations and were divided similarly between A/G and C/T (mutation types).Within the indels, 28 600 polymorphisms are characterized as deletions, whereas 22 853 are characterized as insertions, relative to the reference genome.
All isolates were confirmed as M. acridum based on phylogenetic analysis with 14 other Metarhizium species (Figure S2).Analyses of the phylogenetic network among isolates were based on the conserved dataset of 60 097 SNPs found in coding regions of the genome (Figure 2d).Genetic distance between the isolates was determined based on 387 611 high-quality SNPs identified with GATK (Figure 3).Both datasets identified two distinct genetics clusters, with isolates from western Africa and Mexico grouping together.Within this group, three of the four western African isolates clustered together with relatively short branches, suggesting that these isolates are more genetically similar to each other than the other individuals within this group.In contrast, isolates from eastern Africa, Madagascar, and Thailand formed a second cluster that was more similar to the reference isolate originating from Australia, with notable genetic divergence of the Tanzanian isolate (Figures 2d and 3).

| Virulence bioassays
All nine M. acridum isolates were shown to cause significant mortality to both L. migratoria and S. gregaria in comparison with control treatments.The mortality rates of the different isolates ranged from 70% to 100% in both species of locusts (Figure 4a).Killed insects had a 40%-70% sporulation rate of fungi 3-5 days after mortality, with rates varying across isolates (Table S2).Four of the nine isolates showed significant differences in mortality between the two locust species that were infected with the same fungal isolate.
Hazard ratios were calculated in a pairwise manner using a mixed-effects Cox regression model to account for experimental replicates as a random factor.In L. migratoria, the virulence of the different isolates could be split into three statistically significant groups, with two isolates (ARSEF 3609 and ARSEF 7486) being classified as weakly virulent, whereas the remaining seven isolates were relatively more virulent (Figure 4b; Figure S3, Table S7).In S. gregaria, there were four statistically significant groupings, with three isolates characterized as weakly virulent, whereas the remaining six isolates ranged from medium to high virulence (Figure 4c).The patterns of virulence of the nine isolates between the two locust species are variable, although the two isolates that were weakly virulent against L. migratoria (ARSEF 3609 and ARSEF 7486) were also weakly virulent against S. gregaria.Notably, the reference isolate, ARSEF 324 from Australia, was the single most virulent isolate against S. gregaria, whereas in L. migratoria, it was also highly virulent but grouped together with five other isolates.

| Isolate-specific variation in growth rate
The growth rates of the nine M. acridum isolates were measured on three media that differed in carbon source, as a proxy for fitness differences in varying nutritional environments.In all three media types, a significant correlation was found between isolate and growth rate, indicative of variation in growth rate dependent on the isolate (SDAY: chi-square = 43.9,df = 8, p < 0.0001; trehalose: chi-square = 28.0,df = 8, p < 0.0001; minimal locust: chi-square = 217.0,df = 8, p < 0.0001).In the locust media that was supplemented with crushed L. migratoria, we see that ARSEF 7486 had a statistically significant lower growth rate than the majority of isolates (Figure 5a), and this isolate is also one of the least virulent against L. migratoria.On the trehalose media, we see a more complex pattern of growth rates among isolates (Figure 5b).5c).In a qualitative assessment of growth, visual inspection of the isolates grown on crushed locust media had more uniform phenotypes and consistent sporulation than in SDAY or trehalose media.

| Diversity and neutrality statistics
Pairwise nucleotide diversity (π) and Tajima's D were calculated in sliding windows of 1000 bp across the genotype likelihoods of the nine genetically unique samples of Metarhizium acridum.Using this approach, the mean π across all windows was estimated at 1.49 × 10 −3 (SD = 0.0028) equivalent to one mutation per 671 bases.
Although we did not identify any regions with statistically low nucleotide diversity (3 SD below the mean), 308 windows had extremely high nucleotide diversity (3 SD above the mean).Following the concatenation of regions with at least four consecutive windows, we found a total of 74 genomic regions from 25 different scaffolds exhibiting high nucleotide diversity (Table S3; Figure 6a).For Tajima's D, the mean value across all windows was estimated at 1.49 × 10 −1 (SD = 1.17).We excluded analyses of outlier loci with positive Tajima's D, which is typically indicative of diversifying selection, because these values may be exaggerated by duplicated regions within the genome that can be misinterpreted as polymorphisms.Instead, Tajima's D outlier loci were identified by selecting windows with values less than −2 to find regions inconsistent with neutrality by ways of directional selection (Akey et al., 2004;Fahrenkrog et al., 2017;Stajich, 2004).Here, we find 8 regions that exhibit extremely low values of Tajima's D from 6 different scaffolds (Table S4; Figure 6b).
We further investigated genes that occurred within the regions described above, displaying extremely high values of nucleotide diversity or low values of Tajima's D. We found 56 genes within the 74 genomic regions exhibiting high nucleotide diversity.Enrichment analysis of these gene models revealed 5 molecular functions and 10 biological process GO terms that were significantly enriched (p < 0.01) after adjustment for multiple comparisons (Table 2; Table S8).The GO terms with the lowest p-values in both categories were related to alanyl- Using the full SNP callset, comprising 387 611 bp, we also determined the rate of piN and piS within each gene.We find a total of 72 genes that exhibit no non-synonymous mutations, which resulted in 11 enriched GO terms (Table S6).The GO terms with the lowest p-values in this gene set were largely associated with essential translation and transporter activities.Alternatively, we also retrieved 295 genes with a higher rate of non-synonymous to synonymous mutations (piN/piS >1.2) and found 2 enriched GO terms (Table 2).In this gene set, we find that the enriched GO terms are associated with serine-type hydrolase and peptidase activity (Figure 7).

| Copy number variants
We identified copy number variant (CNV) sites based on read depth following alignment to the reference isolate, ARSEF 324.A total of 640 regions were called, but only regions with an average read depth, across all isolates, less than 0.4 or greater than 2.4 were considered as deletions or duplications respectively.Resultantly, we had 103 regions classified as deletions covering a total of 717 000 bp, and 24 regions classified as duplications covering a total of 307 500 bp.Although the majority of CNVs were intergenic sites, there was a maximum of seven genes occurring within a single duplicated or deleted CNV (duplication, JAGRQK010000011.1: 1 306 501-1 318 500; deletion: JAGRQK010000002.1:269 501-284 500).There was a total of 3 enriched GO terms in deleted regions, with functions primarily related to enzymatic activity, but no significantly enriched terms in CNV duplicated regions (Table 2).

F I G U R E 3 Genetic distance heatmap of all Metarhizium acridum
isolates.Pairwise genetic distance heatmap for all included samples based on the full 387 611 SNP callset, excluding indels.Genomewide estimates of polymorphisms were determined using 1000 bp windows, based on scaffolds from the reference genome, ARSEF 324.

| DISCUSS ION
Within controlled laboratory conditions, we find evidence for intraspecific phenotypic variation in both virulence and metabolic capability among nine isolates of M. acridum.We show that all tested isolates were pathogenic to both L. migratoria and S. gregaria locusts.
However, there is considerable variation between isolates in virulence against the two locust species, with some isolates displaying increased virulence to one species.This result suggests that isolates are better adapted to specific host orthopteran taxa, and we highlight functional regions of the genome with increased diversity that may be correlated with virulence or metabolism.Images depict each isolate grown in the corresponding media for 14 days after inoculation of 3 μL of 1 × 10 7 conidia/mL spore suspensions.Letters indicate different statistically significant growth rates between two isolates (i.e. two isolates with the same letter have similar growth rates, whereas two isolates with different letters have statistically significant growth rates).Isolates in different circles exhibit statistically significant differences.Isolates in overlapping regions are not statistically different from isolates that co-occur in the same circle.
between distantly related isolates is more likely than in generalist counterparts.Nonetheless, we find that total genetic variation does not directly correlate to patterns of virulence or growth rate, which suggests that additional factors influence the structuring of genetic diversity in M. acridum.The countries that these isolates originate from range from tropical rainforest to moist savannah to dry desert biomes.Therefore, besides pathogenicity, isolates may be adapted to climactic conditions, such as humidity, seasonality, and temperature.However, a lack of detailed information regarding collection sites of the sampled isolates impedes the identification of environmental factors that can influence variation.
Although there has been some recent effort regarding the intraspecific ranges of virulence in fungal pathogens, particularly plant pathogens where differences in pathogenicity and disease symptoms have been studied (Elena et al., 2015;Zeng et al., 2018), there are fewer studies focused on entomopathogens.Of the research that has considered intraspecific variation, the focus is largely skewed towards determining whether the cumulative virulence of a fungal species is significantly different in comparison with an alternative species or genus (Liu et al., 2002;Mburu et al., 2009;Paradza et al., 2021).Nevertheless, there is growing evidence that natural variation within entomopathogenic species is evident and important in terms of biological control.isolate capable of inducing nymph mortality (Reingold et al., 2021).
Both of these studies consider the variation within generalist pathogens that are known to infect several orders of insects, where such differences may seem more intuitive as a result of potential adaptation to local insect hosts.In contrast, variability in specialist entomopathogens, such as M. album and M. acridum, has largely been overlooked.
In the present study, the most dramatic variation in virulence between isolates of M. acridum is exhibited in ARSEF 324, which is four times more virulent than both ARSEF 3609 and ARSEF 7486 against S. gregaria (Figure 4c, Hazard Ratio = 0.27).This pattern is, however, host-specific because isolates do not possess the same virulence towards both locust species used in this study, as isolates with high virulence in L. migratoria migratoria were not always highly virulent against S. gregaria, or vice versa.For example, although ARSEF 3341 and 6600 was shown to be highly virulent in L. migratoria migratoria, it was among the lowest in virulence against S. gregaria.None of the two locust species used in virulence assays represent the original hosts of the investigated M. acridum isolates (Figure 2a).Although four isolates from northern Africa have overlapping distributions with S. gregaria, they are not more virulent towards this locust species.Our results do not exclude local host adaptation as a driver of genetic variation.Rather, our results show that virulence in M. acridum is not a fixed trait against all locusts and that some isolates may be better adapted to infect specific orthopteran groups.Even isolates collected from the same host species in Niger (ARSEF 3341,7486) exhibit significantly different virulence in L. migratoria migratoria, suggesting that such differences may even be apparent on relatively small geographic scales.It is important to note that fungal isolates in this study may be adapted to other species among the >26 000 extant species within Orthoptera (Peck et al., 2008;Song, 2018) and not the locust species used for virulence bioassays in this study.This may explain the lack of specificity to S. gregaria by co-occurring isolates.Conversely, previous studies have shown that host resistance may play an important role in pathogenicity, with some species of mosquitos being more susceptible to fungal infection than others (Daoust & Roberts, 1982).Similarly, we find that isolates are generally more virulent against L. migratoria migratoria, which may suggest that S. gregaria has increased defence mechanisms to infection by M. acridum.There is evidence of the efficacy of using indigenous fungal isolates against local lepidopteran pests, such as Tutu absoluta (Tadele & Degaga, 2017) or Helicoverpa armigera (Fite et al., 2020) in Ethiopia.These results have notable ecological implication, providing preliminary evidence for patterns of divergence in virulence between isolates due to local host adaptation and ecological factors, suggesting that natural variation may be used as an exploitable trait when targeting specific locust pests occurring in a specific region.(3) a piN/piS ratios higher than 1.2, indicating that there were higher rates of nonsynonymous to synonymous mutations; and (4) copy number variant regions when averaged across eight isolates, excluding the reference.An average read depth less than 0.4 was indicative of deletions, whereas a read depth greater than 2.4 was indicative of duplications in the genome in comparison with the reference.Gebremariam et al., 2021;Zhang et al., 2011), other studies have failed to find such patterns (Talaei-Hassanloui et al., 2006;Valero-Jiménez et al., 2014).However, our results indicate that testing various nutrient sources can help identify differences in metabolic capability.Further exploration of biologically representative nutrient sources may provide more feasible alternatives for simulating potential virulence as opposed to traditional bioassays.

Statistic
In addition to phenotypic characterization, we also examined genetic variation of the nine isolates with whole-genome re-sequencing and alignment to the reference genome, ARSEF 324.A limitation of our analysis is because genomic variability is solely based on comparison to a single reference, which itself is an anomaly, since it is the only known diploid of M. acridum.Therefore, assembly errors, duplications, or deletions in the reference may indicate decreased or inflated variation in these genomic regions.To reduce such potential sources of error, we removed sites found to be heterozygous when mapping the re-sequenced reference to itself.After this correction of our SNP data, the results showed an average nucleotide diversity across five species (Kobmoo et al., 2021).In contrast, the highly host-specific Ophiocordyceps sinensis found only in the Tibetan plateau has a value of π ranging from 0.53 × 10 −3 to 1.05 × 10 −3 , which is lower than the diversity of M. acridum (Xia et al., 2017).Although only based on 9 isolates, these results may suggest that nucleotide diversity is correlated with host range or global distribution, although further comparative genomics studies are required to confirm this pattern.
We scanned the genome of M. acridum for regions of exceptional diversity using genome-wide polymorphism data.We focused our analysis on extreme outlier regions displaying high diversity ( ) or gene models with increased rates of non-synonymous mutation (piN/ piS).Enrichment analysis revealed that regions of high nucleotide diversity were characterized by functions related to alanine tRNA ligases.Alignment of the two alanine tRNA ligases (KAG8408506.1,KAG8408507.1) shows up to four protein variations from multiple non-synonymous mutations.Although the molecular effect of such polymorphism is unknown, previous studies have shown modification in tRNA ligases can lead to possible amino acid misincorporation (Magnuson et al., 2016;Mühlhausen et al., 2016).Such plasticity is quite well-tolerated in fungi (Simões et al., 2016;Wang & Pan, 2016), suggesting that translation diversity during stress may provide selective benefits.In addition to this, we also find enrichment of mycotoxin biosynthesis, although the biosynthetic product is unknown entomopathogenic fungi are known to produce a variety of mycotoxins, such as non-ribosomal peptide synthetases (NRPS), which are toxic to insect tissues (Pedrini, 2022).However, in M. acridum, genomic analysis has mostly revealed gene loss events of such enzymes, which is proposed to have played a role in the host specialization of this pathogen (Wang et al., 2012(Wang et al., , 2016)).The two genes that we highlight here show high diversity among our isolates, suggesting that both gene loss and diversification of retained mycotoxins may play a role in the orthopteran host specificity of M. acridum.
We also observed that serine-type endopeptidase activity genes were an enriched category in both regions of high nucleotide diversity and in genes with high rates of non-synonymous mutation.Interestingly, we found two Pr1 subtilisin-like serine peptidases among significantly enriched peptidase genes (KAG8413868.1 and KAG8412168.1).Subtilisin-like serine peptidases are one of the three types of proteases used in the initial cuticle breaching stage of Metarhizium, and a total of 11 isoforms are known from the subtilisin-like serine proteases of the Pr1 family in all species of Metarhizium (Freimoser et al., 2003).These genes are thought to be under strong selective pressure to evolve against rapidly evolving host defence mechanisms involving protein inhibitors (Boucias & Pendland, 1987;St Leger et al., 1996).
The Pr1 family can be split into two classes: class I (bacterial-like), which is only comprised of Pr1C, and class II (proteinase K-like), F I G U R E 7 Overall piN and piS distribution.Filtered set of piN (non-synonymous)/piS (synonymous) analysis and functional annotation of all genes with mutations in the re-sequencing data of 9 M. acridum isolates.The overall piN/piS distribution.The dark green points represent 295 genes with higher rates of nonsynonymous mutations, whereas the red points represent genes with higher rates of synonymous mutations.Grey points indicate genes with roughly equal rates of non-synonymous to synonymous mutations.The labelled, light green points indicate genes with GO terms that were significantly enriched.
which comprises the remaining 10 isoforms.We found high levels of diversity in both Pr1C and Pr1E within isolates of M. acridum, but the 29 SNPs in Pr1C led to 337 amino acid changes, suggestive of frame-shift mutations making the protease non-functional (Table S5).A study regarding Pr1 evolution in insect-pathogenic fungi (Gao et al., 2020) showed that Pr1E is unique to Metarhizium compared to the genus Beauveria and has previously been predicted to have more important roles in the diversification of earlydiverging species in this lineage, such as M. acridum.Although we cannot conclude on the phenotypic output of these subtilisin-like serine proteases, we see that ARSEF 3609 from Thailand had unique SNPs in both subtilisin-like serine proteases and was also one of the least virulent isolates among both locust species used in this experiment.
We also looked for regions that exhibited strong indications of purifying or positive selection based on outlier regions of Tajima's D. Within the few genomic regions with signs of such selection, we only found chitin-binding activity to be significantly enriched.
Chitin is one of the main structural components of the insect cuticle, and a recent study has shown that seven chitin synthase genes were identified in M. acridum, contributing to growth, development, stress tolerance, or virulence (Zhang et al., 2019).Further analysis of the enriched chitinase gene showed that all SNPs were only found within ARSEF 3391, resulting in 36 amino acid changes and a potential frame-shift mutation (Table S5).Read depth information from our CNV analysis showed that this genomic region appears to be duplicated within ARSEF 3391, suggesting that the excess of SNPs that we find are likely due to relaxed selection within a duplicated gene copy of the chitinase.If the duplicated gene would remain functional, duplication events of genomic regions unique to isolates can present an alternative method of diversification that may affect virulence or other phenotypic characteristics.

| CON CLUS ION
The present study reveals notable intraspecific variation within the specialist fungal entomopathogen, Metarhizium acridum, in terms of both phenotypic characteristics and genetic diversity.
The natural variation in virulence indicates that the isolates are more suited to infecting some species of Orthoptera more than others, which is an exploitable trait in terms of targeted management of local locust pests on agricultural crops.However, we do not find correlation between geographic region and isolate virulence in M. acridum, which has also been shown in other studies of entomopathogenic fungi (Daoust & Roberts, 1982;Lappa, 1978).
The results also reveal evident genetic diversity between the isolates suggesting undiscovered population structure among globally distributed isolates.Ultimately, although we do see phenotypic diversity in these fungal pathogens, we cannot determine the evolutionary driver of this variation within the scope of this study, although we highlight genomic regions of interest.

Figure
Figure 2a,b; TableS1).Cultures were grown from single spores of
Phylogenomic networks constructed from our high-resolution genomic data revealed two distinct clusters of isolates, with the implication that the five isolates from western Africa and Mexico were more divergent in comparison with the reference originating from Australia.A previous study examining AFLPs in M. acridum identified clustering patterns that reflected geographic origin, and our data support the geographic clustering of isolates into distinct phylogenetic groups(Fernandes et al., 2010).However, the current data suggest that the isolates from eastern Africa are more closely associated with isolates from Thailand and Australia, which was not shown from the AFLP data.When we compare the clustering patterns of the isolates against the phylogenetic relationship of their original orthopteran hosts, we see some evidence for correlation, with the two clusters of M. acridum split between clades of the Acrididae family (Figure2c).Notably, ARSEF 3391 from Tanzania is separate from the rest and originates on a distant host from the Pyrgomorphidae family.These results indicate that further investigation of possible coadaptation between locust hosts and M. acridum populations is needed.The exact cause of genetic similarities within the two clusters of M. acridum in our analysis is unknown, although it is important to note that these results are based on genetic similarity between samples on different continents.Still, the genetic similarity of isolates across continents may be due to either longdistance dispersal, but current knowledge of gene flow at continental distances is understudied in Metarhizium.One study regarding population diversity of Metarhizium generalist isolates showed very little gene flow in populations spanning large distances across the continental western USA, but this may also be due to primarily clonal reproduction in generalist species (Rehner, 2020).Considering the known signatures of sexual recombination in M. acridum, gene flow F I G U R E 5 Growth rates of isolates across different media.Estimated marginal means of the growth rates of all tested isolates of M. acridum in (a) locust minimal, (b) trehalose, or (c) SDAY media.Colours indicate geographic origin, and letters indicate significant groupings in relation to pairwise comparisons of growth rates.

F
Cumulative daily survival of locusts exposed to spores of different Metarhizium acridum isolates.Kaplan-Meier survival curves of (a) Locusta migratoria migratoria (blue) and Schistocerca gregaria (red) individuals inoculated with each of the nine isolates of M. acridum.Experiments were carried out in three replicates of each species, with 10 juvenile locusts on each treatment over 10 days with 5 μL of Triton X (control) or a fungal suspension of 1 × 10 7 conidia/mL.Isolates that exhibited statistically significant differences of mortality between the two locust species are specified with p-values.Isolates are indicated with labels and colours to indicate geographic origin.Pairwise comparisons of all treatments were determined using the coxME model for (b) L. migratoria migratoria and (c) S. gregaria.Values within each box indicate with hazard ratios of the treatment on the y-axis in comparison with the x-axis, with shaded boxes indicating statistical significance (p > 0.05).Significant differences in virulence are summarized in Venn diagrams, with low to high virulence from left to right.
In a study performed by Valero-Jiménez et al. (2014), a significant range in virulence is shown among different isolates of Beauveria bassiana resulting in a difference in effectiveness of up to tenfold against Anopheles coluzzi mosquitos.In another case, isolates of M. brunneum are shown to have altered conidial adhesion dependent on the life stage of aphids, with only one

F
I G U R E 6 Patterns of genome-wide polymorphism in M. acridum.Genome-wide statistics were calculated by (a) average pairwise genome-wide nucleotide diversity, with points highlighted in red indicating windows that exhibit values >3 SD from the mean, and (b) Tajima's D values, with points ± 2 highlighted in red.
Metabolic flexibility is tightly interconnected with the complex infective process of M. acridum.During the cuticle penetration stage, M. acridum employs a wide array of cuticle degrading enzymes, which releases substantial amount of nutrients from the host(Freimoser et al., 2003;He et al., 2012).Multiple metabolic pathways utilizing hydrolysing enzymes and nutrient transportation are activated during this period of rapidly fluctuating nutrient availability.In this regard, mycelial growth rate across different types of media can represent capability for metabolic flexibility and adaptation to specific environments.The three types of media differed in sources of carbon, with growth rate denoting a proxy for fitness across different environmental conditions: (1) SDAY media represented a standard control, as a common media for growing filamentous fungi; (2) trehalose media represented the primary sugar found in insect haemolymph; and (3) minimal media supplemented with crushed L. migratoria migratoria represented an in vitro proxy of the natural host environment.Here, we find some evidence for variation in accordance with virulence, such as in the case of the minimal locust media, where ARSEF 7486 was shown to have a significantly slower growth rate than most isolates, while also being one of the least virulent isolates of L. migratoria migratoria.However, this relationship does not appear consistently, as ARSEF 3609 had a similar growth rate to other isolates, despite TA B L E 2 List of enriched GO terms (p < 0.01) for tests of neutrality and diversity across the genome.Indicated GO terms within molecular function (MF) or biological processes (BP) are found in regions exhibiting: (1) high nucleotide diversity π (μ + 3 SD); (2) low Tajima's D values (<−2.0),indicative of rare polymorphism; (π) of 1.49 × 10 −3 (SD = 0.0028) in M. acridum.Although genomic nucleotide diversity values are unknown in other Metarhizium species, they have been determined for other entomopathogenic fungi in the same order (Hypocreales).The globally distributed, primarily generalist genus Beauveria has a π ranging from 3.0 × 10 −3 to 3.6 × 10 −2

Functional
regions of diversity are associated gene function related to specific subtilisins in the Pr1 family, which could help in understanding the roles of these proteins in terms of localized specialization to specific orthopteran hosts.The evident intraspecific diversity warrants future research on population structures in M. acridum in conjunction with a more in-depth analysis of local orthopteran host ranges.Additionally, we propose that further exploration of intraspecific diversity of subtilisins and duplication events in M. acridum can provide insight into diversification or host specialization in fungal pathogenic systems.This is the first study to consider multiple characteristics of intraspecific variation in M. acridum and shows that further research regarding host-pathogen interactions within Metarhizium can lead to new insights regarding antagonistic reciprocal coevolution.