Population genomics and evolution of a fungal pathogen after releasing exotic strains to control insect pests for 20 years

Abstract

Entomopathogenic fungi are one of the key regulators of insect populations in nature. Some species such as Beauveria bassiana with a wide host range have been developed as promising alternatives to chemical insecticides for the biocontrol of insect pests. However, the long-term persistence of the released strains, the effect on non-target hosts and local fungal populations remains elusive, but they are considerable concerns with respect to environmental safety. Here we report the temporal features of the Beauveria population genomics and evolution over 20 years after releasing exotic strains to control pine caterpillar pests. We found that the isolates within the biocontrol site were mostly of clonal origins. The released strains could persist in the environment for a long time but with low recovery rates. Similar to the reoccurrence of host jumping by local isolates, the infection of non-target insects by the released strains was evident to endemically occur in association with host seasonality. No obvious dilution effect on local population structure was evident by the releases. However, the population was largely replaced by genetically divergent isolates once per decade but evolved with a pattern of balancing selection and towards expansion through adaptation, non-random outcrossing and isolate migration. This study not only unveils the real-time features of entomopathogenic fungal population genomics and evolution but also provides added values to alleviate the concerns of environmental safety regarding the biocontrol application of mycoinsecticides.

Fig. 1. The sketch map of the isolate collection sites.

The isolates of Beauveria bassiana were collected from the biocontrol site in the Anhui (AH) province and other places across the monsoon (MR) and non-monsoon (NMR) regions in China from 1997 to 2017 (Table S1). A blue arrow points to the infestation of the pine caterpillar and the mycosed cadavers after the application of mycoinsecticide at Magushan Forest Farm, AH province. Pink arrows show the wet/rainy phase of monsoon occurring in the later spring and summer seasons, and purple arrows indicate the dry phase of monsoon occurring in the autumn and winter seasons. The red line separates the NMR and MR regions by the mountain ranges in China. Provinces: XJ Xin-Jiang, TB Tibet, QH Qing-Hai, IM Inner Mongolia, HB He-Bai, NX Ning-Xia, SX Shan-Xi, YN Yuan-Nan, GZ Gui-Zhou, GX Guang-Xi and FJ Fu-Jian. The provinces in the south and north of AH within the monsoon regions are labelled as MR1 and MR2, respectively.

Fig. 2. Relationship among the isolates collected at the biocontrol site over 20 years.

a Phylogenetic analysis. The neighbour-joining tree is generated based on the average pairwise genetic distance of all biallelic SNPs. Isolates were separated into three main clustering groups, i.e., the root group (RG), G1 and G2 groups. The isolates shadowed in the different colours represent the respective collection year. The mating type (MAT) and original insect host of each isolate are indicated in different coloured circles and squares. C1–C24 are the clonal groups that clustered isolates with ANI values > 99.8%. The abbreviations of insect host orders are: COL Coleoptera, LEP Lepidoptera, HEM Hemiptera, HYM Hymenoptera, DIP Diptera and ORT Orthoptera. b Principal component analysis (PCA) of isolates based on all biallelic SNPs. The first two components are presented. c Venn diagram analysis of the clonal and divergent isolates collected in different years.

Fig. 3. Comparative analysis of population structure.

a Ancestry components of two major lineages of AH isolates (as shown in Fig. 2a). b Decay of linkage disequilibrium (LD, expressed in terms of correlation coefficient, r²) as a function of distance averaged over 100 kb for the three populations. Inset shows the declines in LD over the first 100 bp. c Estimation of Tajima’s D of the three populations. Tajima’s D was independently calculated within 5-kb sliding windows across the genome. The median value of each chromosome is shown. d Genome-wide analysis of genetic divergence between populations. The value of F_ST was calculated in 5-kb windows, and the median values shown in parentheses were estimated between two populations. The cut-off lines shown in each panel represent the Wright’s threshold (F_ST = 0.05–0.15) for moderate differentiation. e Derived allele frequency spectra between pairwise populations. The colour of each slot indicates the frequency of each derived allele pair. The scale values along each panel are the exponent of order 2.

Fig. 4. Genetic recombination analysis.

a Variation in recombination rates (RRs) across chromosomes in three populations. b Comparative estimation of RRs along the genome for three populations. The median value of each chromosome was used for plotting each population. c Identification of potential recombinant genomic blocks, which were determined by >1000 consecutive and shared SNPs between the released strains and local opposite mating-type isolates. Each dot represents the number of the conserved SNP blocks present in the MAT1-1 isolates of G1 or G2 lineages shared with the released strain Bb13 or Bb17. The lines in the middle show the median values. d Sexuality induction between the opposite mating-type isolates. Panels 1 and 2 represent the phenotypes of crossings: Bb175 MAT1-1 × MAT1-2, and Bb166 (MAT1-1) × Bb175 (MAT1-2); Panel 3, microscopic observation of the sexual perithecia (PE) and asci (AS) produced on the inducible fruiting bodies. Bar, 50 μm. e PCR verification of the karyotype of the selected isolates. The presence or the absence of the fragment indicates isolates belonging to the MAT1-1 and or MAT1-2 genotype.

Fig. 5. Divergence and gene flow between geographic populations.

a Median genetic differentiation (F_ST) between the pairwise populations. AH, MR1, MR2 and NMR populations are as shown in Fig. 1. b Phylogenetic relationships of all isolates collected from different geographic origins. The Neighbour-joining tree was obtained based on the average pairwise genetic distance of all biallelic SNPs. For simplicity, lineage grouping is determined at a cut-off value of ANI > 99%. The blue branch lines represent the associations with the AH population phylogenies obtained in Fig. 2a. Two released strains Bb13 and Bb17 are highlighted in bold. c Genetic relationships among fungal isolates based on the estimation of individual ancestry. The bar plots in each line show the ancestry proportion over a range of presumed ancestry number, K = 2–6. Abbreviations for geographic origins are as shown in the caption of Fig. 1.