Cordyceps and Beauveria infections drive species-specific microbiome dysbiosis in the mosquito Aedes aegypti

E. Everett1,2, Haley M. Gore1, Swaksha Kallepalli1,2, Kristin R. Duffield1, Lina Flor-Weiler1, John Marino2 and Jose Luis Ramirez1*

1 USDA-ARS, National Center for Agricultural Utilization Research, Crop BioProtection Research Unit, Peoria, IL, United States, 

2 Department of Biology, Bradley University, Peoria, IL, United States

With the rising prevalence of vector-borne diseases and insecticide resistance in mosquitoes, alternative vector control strategies are urgently needed. Fungal entomopathogens offer a promising approach with a decreased likelihood of resistance development in mosquito populations. However, the mechanisms by which each fungus contributes to host mortality remain poorly understood, and the potential role of microbiome disruption as a secondary pathogenic mechanism has received limited attention. We evaluated the impact of four entomopathogenic fungal species (Beauveria bassiana, Cordyceps javanica, C. cateniannulata, and C. amoenerosea) on the microbiome of the yellow fever mosquito (Aedes aegypti) colonized with a defined, field-derived bacterial community. Whole-body bacterial communities were profiled using high throughput 16S rRNA amplicon sequencing, and community structure was assessed through alpha diversity metrics, beta diversity analysis, hierarchical clustering, and linear discriminant analysis effect size (LEfSe). All four fungal species successfully infected the mosquito; however, their effects on the mosquito microbiome were species-specific. Cordyceps javanica and C. cateniannulata reduced community evenness without significantly affecting species richness, a pattern consistent with a dominance-driven dysbiosis rather than broad bacterial loss. Infections by C. amoenerosea significantly increased total bacterial load and drove strong enrichment of the opportunistic genus Pandoraea, suggesting epithelial disruption or immune dysregulation as possible contributing factors. Beta diversity analysis indicated partial community-level restructuring across all fungal infections. B. bassiana showed a distinct genus- level compositional response, with enrichment of core symbiotic taxa and depletion of Chryseobacterium and Kluyvera, which was different from the Enterobacteriaceae-dominated shifts seen across infections with Cordyceps species. LEfSe analysis identified Kluyvera and Burkholderia as the strongest genus-level discriminators of infection state, suggesting potential utility as microbiome-based indicators of successful fungal colonization. These key findings were independently validated using EdgeR and batch-corrected MaAsLin2 analyses, with Pandoraea enrichment under C. amoenerosea and Burkholderia depletion under C. cateniannulata confirmed by both methods. Taken together, these results show that entomopathogenic fungi restructure the Ae. aegypti microbiome in a species-specific manner, inducing community destabilization and opportunistic bacterial enrichment that likely contribute   to the detrimental effects of fungal infection. These results provide a mechanistic insights for the selection and development of fungal biopesticides for mosquito control.

KEYWORDS

entomopathogenic fungi, fungal biopesticides, microbial control, microbiome, mycopesticide

OPEN ACCESS

EDITED BY

Kokouvi Kassegne,Shanghai Jiao Tong University, China

REVIEWED BY

Amr Mohamed, Cairo University, Egypt 

Edward D. Walker, Michigan State University, United States 

Muniaraj Mayilsamy, Vector Control Research Centre (ICMR), India

*CORRESPONDENCE

Jose Luis Ramirez jose.ramirez2@usda.gov

RECEIVED 12 May 2026

REVISED 12 June 2026

ACCEPTED 19 June 2026

PUBLISHED 03 August 2026

CITATION

Everett E, Gore HM, Kallepalli S, Duffield KR, Flor-Weiler L, Marino J and Ramirez JL (2026) Cordyceps and Beauveria infections drive species-specific microbiome dysbiosis in the mosquito Aedes aegypti.

Front. Microbiol. 17:1879658. doi: 10.3389/fmicb.2026.1879658

COPYRIGHT

© 2026 Everett, Gore, Kallepalli, Duffield, Flor-Weiler, Marino and Ramirez. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

1  Introduction

Vector-borne diseases represent one of the most significant threats to global public health, with a mortality burden of ∼1.4 million deaths annually (Dorn et al., 2011). Among arthropod vectors, mosquitoes are of particular importance  due  to  their  role in transmitting pathogens that cause diseases such as yellow fever, Zika, West Nile, malaria and dengue (Duguma et al., 2020; Onen et al., 2023). Efforts to control these diseases have been hampered by the limited availability of vaccines or effective treatments (Scolari et al., 2019). Thus, management of vector-borne disease outbreaks still relies heavily on vector control, mainly through the use of synthetic chemical insecticides (Norris, 2004).

However, while synthetic chemical insecticides have historically been the primary tool for mosquito control, their indiscriminate use has led to the rapid development of resistance in all major mosquito vectors (Chaudhry et al., 2019; Meier et al., 2022; Perumal et al., 2024) and, also poses significant environmental risks impacting non-target organisms, including humans and other animals (Kumari et al., 2022; Zhou et al., 2025). These concerns highlight the urgent need for alternative mosquito control strategies. One promising approach is the use of microbial biopesticides, including bacterial and fungal entomopathogens, due to their specificity, reduced environmental impact, and lower potential for resistance development (Thakur et al., 2020; Salimi and Hamedi, 2021).

Currently,  entomopathogenic  fungi  have   been   identified  as a prospective agent  of  mosquito  control,  with  the  potential to kill the  most  important  mosquito  vectors  (Scholte et al., 2007; Mnyone et al., 2009; Perumal et al., 2024). These fungi act by producing conidia that adhere to the  mosquito  cuticle  in  order to penetrate the mosquito exoskeleton and colonize the hemocoel. The mosquito is able to sense the fungal infection at  the cuticle level via the action of pattern recognition receptors (PRRs) binding to pathogen or microbe associated molecular patterns (PAMPs or MAMPs) (Hughes, 2012; Hillyer, 2016; Rodgers et al., 2017; Ramirez et al., 2020) leading to the induction of several immune signaling pathways (Ramirez et al., 2019; Tikhe and Dimopoulos, 2021; Rhodes et al., 2025) and the production  of  antimicrobial  peptides  (Tawidian et al.,  2019; Yan and Hillyer, 2019; Ramirez et al., 2023) that can suppress fungal proliferation. This response leads to a complex reciprocal interaction with direct and indirect effects of the host immune response on the microbiome (Ramirez et al., 2018; Gabrieli et al., 2021; Armitage et al., 2022; Zheng et al., 2023), pH alterations within the cuticle and hemolymph, and oxidative and melanization defenses  that  have  an  effect  on   the   fungal   pathogen   and the host (Santiago-Álvarez et al., 2006; Nakhleh et al., 2017; Smith and Casadevall, 2019). The entomopathogenic fungi  in  turn secretes compounds (enzymes, elicitors, and secondary metabolites) aimed at suppressing the host immune response (Samuels R et al., 2016; Boucias et al., 2018; Qin et al., 2023). These fungal-derived compounds are thought to also be directed at inhibiting the host microbiome, as the fungi compete with resident microbial communities, which in turn can further modulate immune responses and impact host survival (Prompiboon et al., 2008; Boucias et al., 2018). Host death after fungal infection is preceded by the successful suppression of the host microbiota and immune response and is followed by emergence and sporulation on the host cadaver. This emergence is mediated by secondary metabolites produced by the fungi, such as oosporein, which act after host death to reduce competition by the microbiota before sporulation (Nielsen and Hajek, 2006; Fan et al., 2017).

The mosquito microbiome has gained attention for its vital role in host physiology, including effects on larval development, immune function, fecundity, and traits related to pathogen transmission (Hegde et al., 2018; Martinson and Strand, 2021). Mosquito-associated bacteria are influenced by biotic and abiotic factors including the mosquito genome, diet and geographical location  (Muturi et al.,  2018;   Cansado-Utrilla et al.,   2021).  The host microbiota also plays a crucial role in mosquito reproduction  and  development,   influencing   key   processes such as blood meal digestion, reproduction (Coon et al., 2016; Guégan et al., 2018) and successful larval pupation (Wang et al., 2018). During development, certain gut bacteria  residents  produce a  hypoxia  signal  which  promotes  the  transcription  and initiation of processes needed for mosquito growth and development  (Valzania et al.,  2018).  Furthermore,  previous work on axenic (germ-free) mosquitoes has demonstrated significant physiological impairments, including altered gene expression related to nutrient absorption, metabolism, and stress response, alongside developmental failure, thereby highlighting the crucial physiological importance of the mosquito microbiome (Vogel et al., 2017).

Dysbiosis of the mosquito microbiota has been linked to altered vector competence, affecting arboviral transmission dynamics (Ramirez et al., 2012; Hegde et al., 2015; Gómez et al., 2022). Additionally, symbiotic bacteria such as Wolbachia can interfere with pathogen replication and influence vectorial competence (Johnson, 2015; Dacey and Chain, 2020). A deeper understanding of how these fungal infections impact the mosquito microbiome is critical for optimizing fungal virulence and enhancing biocontrol strategies. Identifying microbial changes in response to fungal infection can  reveal  key  bacterial  taxa  that  may  play  a  role  in host susceptibility and infection dynamics. Here, we define “dysbiosis” as infection-associated shifts in community evenness and composition relative to the inoculated community. However, mosquito microbiomes are inherently variable, and the concept    of a stable “eubiotic” state in these systems is often contested (Hooks and O’Malley, 2017; Muturi et al., 2017).

While Wei et al. (2017) demonstrated that entomopathogenic fungi interact with  the  gut  microbiota  to  accelerate  mortality  in Anopheles gambiae, that study focused on a single fungal species and the outcome of mortality rather than the specificity    of microbiome restructuring. The present study extends this framework to Aedes aegypti, an important vector with different physiology and vectorial capacity. We also employ four distinct fungal species that for the first time allows direct within-study comparison of how different entomopathogenic fungal species differentially reshape the mosquito microbiome. This comparative analysis is valuable for guiding the selection of fungal biocontrol agents based not only on killing efficiency but on their specific effects on the mosquito microbiome, providing comprehensive insights into diverse host-pathogen interactions. Our results demonstrate that fungal infection significantly impacts microbiome composition and diversity in a strain-dependent manner, offering new insights into host-pathogen-microbiome interactions and guiding the development of future mosquito biocontrol strategies.

2  Materials and methods

2.1  Mosquito rearing

Aedes aegypti mosquitoes (Rockefeller strain) were reared under a 12 h light/dark cycle at 28◦C and 70–80% relative humidity according to standard rearing conditions. Eggs were produced  post blood feeding of adult females using an artificial membrane feeding system and bovine blood (HemoStat Laboratories Inc., Dixon, California, U.S.A.). Larvae were provided a 2:1:1 diet of rabbit food:liver powder:fish flakes and adults were maintained with a 10% sucrose solution in sterile cotton.  Once  emerged, adult mosquitoes were initially fed a 5% sterile sucrose solution supplemented with a defined bacteria set  previously  isolated  from field-caught mosquitoes. This solution thereby mimicked the normal gut microbiota found in natural settings. The bacterial set consisted of seven taxa: Pandoraea sputorum, Strenotrophomonas spp, Leclercia spp, Burkholderia spp, Pseudomonas spp, and two strains of Serratia marcescens. Bacteria were incubated overnight at 27.4◦C, 186 RPM, and concentrations were determined using a hemocytometer. A 5% sterile sucrose solution was used to adjust concentrations to 1 x 108  cells/mL and provided to mosquitoes   ad libitum via sterile cotton. Mosquitoes were maintained on this diet for 5 days until subsequent exposure to fungal pathogens (see below). This experimental design follows the synthetic community (SynCom) approach, where a defined assemblage of field-derived bacteria is used to colonize germ-free or conventionally reared hosts under controlled conditions (Coon et al., 2014; Sommer et al., 2015; Gould et al., 2018). This method offers significant advantages, including enhanced experimental reproducibility and the ability   to attribute observed phenotypes to specific bacterial consortia. Moreover, SynCom approach allows for the incorporation of ecologically relevant taxa directly isolated from the target host population (Valzania et al., 2018).

2.2  Fungal infection assays

Four entomopathogenic fungal species were used for infection in this study: Beauveria bassiana (MBC 076), Cordyceps javanica (MBC 439), Cordyceps cateniannulata (MBC 6241), and Cordyceps amoenerosea (MBC 987). Fungi were cultured on Sabouraud dextrose agar and yeast extract (SDAY) medium and incubated at 20–25◦C for 14 days. Spores were harvested the day of infection into 1 mL of soybean oil and vortexed for 10 min before being filtered through cheesecloth. Conidial concentrations were determined using a hemocytometer and adjusted to 1 x 108 conidia/mL for C. javanica, C. cateniannulata, and C. amoenerosea, and 1 x 107 conidia/mL for B. bassiana. The lower concentration used for B. bassiana was based on prior work with these strains and it was necessary to prevent early mortality and ensure sample collection at day 5 post-infection. Female mosquitoes were cold anesthetized 5d post emergence and 50.6 nL of the conidial oil suspension was topically deposited to the ventral surface of the coxal region using a Nanoject III micropipette (Drummond Scientific, Philadelphia, PA, USA). Control mosquitoes received 50.6 nL of soybean oil alone, following the same protocol. Two independent experiments were conducted, with 30 mosquitoes infected per treatment group in each experiment. From each experiment, 10 mosquitoes were individually collected per treatment group,  yielding  a  total  of  20 mosquito samples per treatment for downstream sequencing analysis. This number (n = 10 per experiment, 20 total per treatment) was selected taking into account mosquito mortality and the cost associated with sequencing. Fresh conidial oil suspensions were prepared independently for each experiment. Following infection, mosquitoes were maintained under the standard rearing conditions described above and provided a 10% sucrose solution.

2.3  DNA extractions and 16S rRNA amplicon sequencing

Mosquitoes were CO2 anesthetized 5 days post-infection, collected and stored at −80◦C until DNA extraction. Individual mosquitoes were homogenized in a lysis buffer using a 3.2 mm and 2.4 mm beads in a TissueLyser II (Qiagen, Germany). Whole body mosquito DNA was extracted using the DNeasyQR  Blood and Tissue Kit (QIAGEN, Netherlands) following the manufacturer’s protocol and DNA  concentration was measured using the Qubit    4 Fluorometer (Fisher Scientific, USA). PCR amplification was done on the V3 and V4 regions  of  the  16S  ribosomal  RNA  gene (16S rRNA), yielding an amplicon product of ∼460 bp. Reactions were prepared using 2x KAPA Hifi Hotstart ReadyMix (KAPA Biosystems Willington, MA), 2.5 μL of template DNA, and 0.5 μL of each primer (Supplementary Table S1) and included negative controls. Thermocycling conditions consisted of an initial denaturation step (3 min, 94◦C) followed by 25 cycles consisting of a denaturation step (30s, 95◦C), annealing step (15s, 58◦C), and extension step at (15s, 72◦C), followed by a final extension step (5 min, 72◦C). Library quantification was done using the Qubit     4 Fluorometer (Fisher Scientific, USA) and samples were pooled, normalizing to 4 nM. Libraries were denatured with 0.2M NaOH and then diluted to 8pM. PhiX of the same concentration was added to make up 15% of the final pool and loaded onto the MiSeq V3 reagent cartridge. Libraries were sequenced on Illumina MiSeq with the 600 cycle v3 kit (Illumina) and 300 cycles of paired end reads. Whole-body microbiome profiling is common in mosquito microbiome studies (Hegde et al., 2018; Noskov et al., 2021) and captures bacteria that translocate to the hemocoel during fungal infection, which is relevant to fungal-infections of mosquitoes.

2.4  Quantitative real-time PCR for bacterial and fungal load

To confirm successful colonization and assess its effect on   the overall bacterial community, we quantified fungal load (via 18S rRNA) and total bacterial load (via 16S rRNA) across all treatment groups. Bacterial and fungal loads were determined using quantitative PCR (qPCR) using SYBRQR   Green Reagents (Qiagen) with samples being run on a 384-well block in the QuantStudioTM 6 Flex System. Targets for quantification included 16S rRNA for bacterial load and ITS for fungal load (Supplementary Table S1). qPCR mixtures were prepared, in duplicate, using 5 μL of SYBR Green dye, 0.3 μL of  both  forward  and  reverse  primers,  3.4  μL of deionized water, and 1 μL of the template DNA. Runs included a non-template negative control to test for potential contamination. qPCR was performed using an initial melting cycle (20s, 95◦C), followed by 40 cycles of denaturation (1s, 95◦C) and annealing and extension (20s, 60◦C). Melt curves were generated following the initial PCR with a denaturation step (15s, 95◦C),      a cooling step (1 min, 60◦C), and a final heating at a rate of 0.05 ◦C/s (15s, 95◦C). Expression of these genes was normalized using the ribosomal protein Rsp17 (AAEL004175) as a reference gene (Ramirez et al., 2021), and relative expression was calculated using the ti. ti. CT method (Livak and Schmittgen, 2001). Rps17 was selected as a reference gene because of its demonstrated stability across immune-challenged Ae. aegypti (Ramirez et al., 2019, 2020, 2021). This primer  also  provides  a mosquito-specific signal that normalizes for variation in DNA input across individual samples. This  normalization  yields relative rather than absolute microbial load estimates.

2.5  Sequencing data processing and analysis

Samples for both experiments were sequenced in the same MiSeq run. Raw paired-end reads were demultiplexed based on barcode primer-tags and processed using the DADA2 pipeline, which included quality trimming, filtering, denoising, merging of paired reads, and chimera removal. Taxonomy was assigned against the SILVA reference database (version 138.2) at a 97% taxonomy similarity threshold. ASVs classified as archaea, mitochondria, chloroplasts were removed prior to downstream analysis. Alpha diversity was evaluated using observed ASVs, the Shannon Diversity Index, and the Simpson Diversity Index within the MicrobiomeAnalyst platform.

Beta diversity at the ASV level was analyzed by Principal Coordinate Analysis (PCoA) based on the Bray-Curtis dissimilarity matrix, with group-level differences tested by pairwise PERMANOVA. To assess the contribution of experimental replicate (experiment 1 vs. experiment 2) to whole community composition, we conducted a sequential PERMANOVA (Type-I SS, with 999 permutations) with “experiment” as the first term followed by “treatment” in R using the vegan package (adonis2 function). To visualize compositional differences across treatment groups, Z-score normalization was applied to the 16 most abundant bacterial genera and a heatmap was generated using hierarchical clustering. A five-way Venn diagram was also constructed to assess the number of shared and unique ASVs among control and fungal-infected groups.

Linear discriminant analysis effect size (LEfSe) analysis was performed to identify differentially abundant ASVs and genera among groups. A Kruskal-Wallis alpha value of 0.05  and  an  LDA threshold score of 2.0 were used to select discriminative features. Graphical outputs were used to visualize LDA scores and relative abundances (Supplementary Figures S3, S4). To validate the heatmap and LEfSe results and to identify additional associations, we used two differential abundance analyses for each fungal treatment vs. uninfected control. First, a Multivariable Association with Linear Models (MaAsLin2) (Mallick et al., 2021) was run as the primary  validation,  with  experimental  replicate as a fixed-effect covariate to account for between experiment variation. This analysis was performed at the ASV level using TSS normalization, log transformation,and a linear model (min_prevalence−0.1, BH FDR correction). Any significant ASVs were subsequently mapped to the genus using SILVA 138.2 taxonomy (Supplementary Table S5). This approach was selected as the primary validation because it uses the same TSS normalization as the original heatmap  and  the  relative  abundance  analyses  are conducted at the ASV level (similar to  LEfSe  analysis).  More importantly, it takes into account the experimental batch- effects. Second, EdgeR (Robinson et al., 2010) was run as a complementary analysis using the EdgeR package. ASVs were filtered prior to testing using “filterByExpr” using min.count = 10,  min.total.count  =  15,  min.prop  =  0.7  (Chen et al.,  2016). EdgeR differential abundance analysis for each fungal treatment vs. uninfected controls using the TMM normalization  was applied, dispersions were estimated with “estimateDisp()”, and pairwise comparisons were evaluated using the “exactTest” and    a Benjamini-Hochberg FDR threshold of 0.05 was used. Because the two methods differ in normalization method and the level at which the bacteria are tested, they are not expected to produce identical results.

2.6  Statistical analysis

All statistical  analyses  and  visualizations  were  conducted  in R (version R  4.2.2)  using  packages  phyloseq,  vegan,  ggplot2, pheatmap, and in GraphPad Prism 10. Statistical differences in bacterial and fungal load were evaluated via Kruskal-Wallis test followed by Dunn’s correction for multiple comparisons. Alpha diversity across  groups  were  evaluated using Kruskal-Wallis tests followed by pairwise Wilcoxon rank-sum tests, with false discovery rate (FDR) correction for multiple comparisons.