Adult Anopheles mosquitoes were collected during the rainy season in four rural locations (Katana, Rushebeyi, Lwiro and Maziba) in Sud Kivu in the Democratic Republic of the Congo (DRC) in November 2021 and January-February 2024 ( Figure 1). Collections were undertaken using CDC light traps and morphological identification performed at the Centre de Recherche en Sciences Naturelles (CRSN), Lwiro, DRC using identification keys ( Coetzee 2020). Individual adult mosquitoes identified as An. demeilloni were stored in 1.5 mL Eppendorf tubes containing cotton pads and silica gel to prevent microbial contamination and transferred to the University of Warwick (UK) for further analysis. Adult female mosquitoes were subject to dissection to separate abdomens (containing ovaries) from the cephalothorax (containing Salivary glands) to allow a more accurate assessment of Wolbachia and Plasmodium sporozoites prevalence rates. Dissected body parts were individually placed in 96-well extraction plates and homogenised using a Qiagen Tissue Lyser II and Qiagen 5mm stainless beads. DNA was extracted using Qiagen DNeasy Blood and Tissue Extraction kit following manufacturer’s instructions. Extracted DNA was eluted in a volume of 50 μL and concentration quantified using Invitrogen Qubit DNA High Sensitivity Assay kits in combination with an Invitrogen Qubit 4 Fluorometer. Molecular identification of species was performed using a An. demeilloni specific qPCR assay targeting internal transcribed spacer (ITS2) sequence ( Walker et al. 2021). qPCR reactions were run on an Agilent Technologies Strategene Mx3005P in a final reaction volume of 10 μL containing 5 μL of Applied Biosystems™ SYBR™ Select Master Mix (catalogue no. 4472908), 1μL of 10X for each of the forward and reverse primers, 1μL of water and 2 μL of the extracted DNA. Reactions were run at 50°C for 2 minutes, 95°C for 2 minutes, followed by 40 cycles of 95°C for 15 seconds, 59°C for 22 seconds and 72°C for 1 minute. A melting cycle of 95°C for 10 seconds, 65°C for 60 seconds and 97°C for 1 second was run to ensure the correct target was amplified.

The detection and quantification of Wolbachia strain wAnD were undertaken on abdomens and cephalothorax samples using qPCR assays targeting the conserved Wolbachia16S rRNA gene ( Gomes et al. 2017). Each sample and no template controls (NTCs) were run in triplicate in a final reaction volume of 10 uL that consisted of 5 uL of Applied Biosystems™ SYBR™ Select Master Mix (catalogue no. 4472908) with a final concentration of 1 μM of each primer, 1 μL of water and 2 μL of the DNA samples. The qPCR assays were carried out in the Agilent Technologies Strategene Mx3005P under the following conditions: 50°C for 2 minutes, 95°C for 2 minutes, followed by 40 cycles of 95°C for 15 seconds, 57°C for 22 seconds and 72°C for 1 minute. A melting cycle of 95°C for 10 seconds, 65°C for 60 seconds and 97°C for 1 second was run to ensure the correct target was amplified. Wolbachia 16S rRNA gene copies per nanogram of total DNA was calculated from a standard curve of a synthetic oligonucleotide standard (Integrated DNA Technologies) used to calculate 16S rRNA gene copies per ul after ten-fold serial dilutions ( Jeffries et al. 2021). The qPCRs results were analysed using Agilent Technologies software.

Wolbachia surface protein ( wsp) gene sequence analysis were carried out using primers wsp81F: 5’-TGGTCCAATAAGTGATGAAGAAAC-3’ and wsp691R: 5’-AAAAATTAAACGCTACTCCA-3’ ( Zhou et al. 1998) in a Bio-Rad T100 Thermal Cycler using previously optimised cycling conditions (Jeffries et al. 2018). The PCR reactions were optimised in a final volume of 20 μL containing 10 μL of 2X Phire Hot Start II PCR Master Mix (Thermo Scientific, USA), 2 μL of 10X for each primer, 4 μL of water and 2 μL of DNA. PCR products were separated and visualised using 2% E-Gel EX agarose gels (Invitrogen) with SYBR safe and an Invitrogen E-Gel iBase Real-Time Transilluminator. PCR products were submitted to Source BioScience (Source BioScience Plc, Nottingham, UK) for PCR reaction clean-up, followed by Sanger sequencing to generate forward reads. Sequencing analysis was carried out in 4peaks ( https://nucleobytes.com/4peaks/) and consisted of sequences being manually checked, edited, and trimmed as required. Sequences were used to perform NCBI BLAST database queries and searches against the Wolbachia PubMLST database ( Jolley et al. 2018) for strain typing. Sequence alignments were constructed in MEGA 12 ( Kumar et al. 2024) ( https://www.megasoftware.net) using the Muscle algorithm to include relevant sequences from the BLAST analysis and Wolbachia MLST database. The phylogenetic tree was generated using the Maximum Likelihood method based on the Tamura-Nei model of nucleotide substitutions. The tree with the highest log likelihood after 1000 bootstrap replicates was shown. The percentage of trees in which the associated taxa clustered together is shown below the branches.

Plasmodium falciparum detection on mosquito samples was performed using a qPCR assay targeting the Pl. falciparum sporozoite cytochrome c oxidase subunit 1 ( Cox1) mitochondrial gene ( Marie et al. 2013). The qPCR reactions were prepared and optimised using 5 mL of Agilent Brilliant III Ultra-Fast SYBR Green Low ROX qPCR Master Mix (catalogue no. 600892) with 1 mL for each primer (5′-TTACATCAGGAATGTTATTGC-3′ and 5′-ATATTGGATCTCCTGCAAAT-3′) with a final concentration of 1mM of each primer, 1mL of water and 2mL of DNA. Each sample and NTCs were run in triplicate in the Agilent Technologies Strategene Mx3005P using the following cycling conditions: 95°C for 3 minutes, followed by 40 cycles of 95°C for 20 seconds, 60°C for 22 seconds and 72°C for 1 minute. A melting cycle of 95°C for 10 seconds, 65°C for 60 seconds and 97°C for 1 second was run to ensure the correct target was amplified.

A sub-sample of DNA from the same individual mosquito for both the abdomen and corresponding cephalothorax underwent Illumina amplicon sequencing of the 16S rRNA gene including variable Wolbachia infection status (based on qPCR). Illumina sequencing targeted the V1-V2 hypervariable regions of the 16S rRNA gene using the universal primers 27F (5’-ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNAGAGTTTGATCMTGGCTCAG-3’) and 338R (5’-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTTGCTGCCTCCCGTAGGAG- 3’). Libraries were prepared and sequenced by the Genomics Facility at the School of Life Sciences, University of Warwick (UK). Amplicon sequencing was achieved using Illumina MiSeq v3 (Illumina, USA) and resulting data received as FASTQ files of the demultiplexed paired end reads. All analysis of resulting reads was carried out in R software version 4.4.0 ( www.R-Project.Org/ ). FASTQ files were then analysed using a DADA2 pipeline ( Callahan et al. 2016) after primer sequences were removed using the Cutadapt package ( Martin 2011). Samples with fewer than 1000 reads per sample were removed from the dataset. Read quality profiles were visualised with a quality score of 20 accepted as a minimum. Reads were then filtered and trimmed based on the observed quality. Paired-end reads were merged and concatenated due to the non-overlapping nature of the V1-V2 hypervariable regions. Taxonomy was assigned using the Silva database (Version 138.1 SSU), and the resulting taxonomic assignments were then visualised and explored using the Phyloseq package ( McMurdie and Holmes 2013).

Comparative statistics were carried out using the JMP 18.0 software (SAS Institute. Inc, North Carolina). General Linear Models were performed assuming a binomial distribution of Wolbachia infection in body parts and across different locations where mosquitoes were sampled. The interaction between both variables was tested but removed from the model if not significant. Wolbachia density across the different body parts was compared using a Mann-Whitney test and a Spearman rank correlation coefficient was performed to compare Wolbachia density between the abdomen and the cephalothorax. Logistic regressions were performed to test the effect Wolbachia infections in the likelihood of Plasmodium infection in mosquitoes. A power analysis was performed using JMP data analysis software ( https://www.jmp.com/en/software/data-analysis-software) to determine the minimum sample size to detect a difference of Plasmodium infection rates between Wolbachia positive and negative samples based on the two-sided test of two independent proportions with a significance level of 0.05 and a power of 0.80. Microbiome diversity was explored by assessing alpha and beta diversity. Alpha diversity was measured using both Simpson and Shannon diversity indices, and beta diversity measured using Bray-Curtis measure of dissimilarity. The difference between the microbiome composition in DNA extracted from abdomen vs cephalothorax was assessed using permutational analysis of variance (PERMANOVA) analysis using the vegan package ( Oksanen 2025). The composition of the microbiome was observed using relative abundance of data agglomerated to the genus level.

A total of 672 mosquitoes morphologically identified as An. demeilloni were analysed from collections in 2021 from four different areas in eastern DRC (Maziba, Katana, Lwiro and Rushebeyi) in addition to a further 96 mosquitoes collected in 2024 from Lwiro. As other Anopheles species are present in DRC including An. funestus which is morphologically similar to An. demeilloni, we first confirmed species using a qPCR assay previously developed to specifically target the ITS2 region of An. demeilloni ( Walker et al. 2021). Using this assay, the accuracy of morphological identification was calculated at 91.7% (704/768). We then determined Wolbachia prevalence rates in An. demeilloni samples using a qPCR assay targeting the conserved Wolbachia 16S rRNA gene. A total of 610 individual An. demeilloni females collected in 2021 were analysed according to body part (cephalothorax and abdomen) to provide a more comprehensive assessment of relative infection rates given mosquito abdomens contain ovaries which are the primary tissues for Wolbachia infection due to maternal transmission. The corresponding cephalothoraxes provided a measure of Wolbachia somatic tissue infection and is often used as a proxy for Salivary glands. Overall, the Wolbachia prevalence rate was significantly higher in abdomens (80%) compared to the cephalothoraxes (17%) (χ ² = 1066.86 Df = 1 P<0.0001) ( Figure 2). There was also a significant difference in Wolbachia prevalence rates across geographical areas, with significantly higher prevalence rates observed in abdomens of individuals collected in Maziba compared to their counterparts in other areas (χ ² = 39.82 Df = 3 P < 0.0001) ( Figure 2). Interestingly, somatic Wolbachia infections in the cephalothoraxes of An. demeilloni were detected in 18% and 2% respectively in Maziba and Lwiro. However, our analysis provided no evidence of Wolbachia infection in the cephalothoraxes of females collected in Katana and Rushebeyi. To provide stronger evidence of somatic infection, we undertook sanger sequencing of the Wolbachia surface protein ( wsp) gene on paired An. demeilloni cephalothoraxes and abdomens shown to be Wolbachia-positive through 16S rRNA qPCR. All our sequences (available at https://osf.io/y6sh8) clustered with existing wAnD wsp sequences, accession numbers in GenBank MW250715.1 and MW250714.1 in BLAST analysis and phylogenetic analysis ( Figure 3A, 3B) indicating no wAnD strain variation was detected.

Quantitative PCR was used targeting a fragment of the conserved Wolbachia 16S rRNA gene with samples run in triplicate alongside standard curves and no template controls (NTCs). Values in the mosaic plot represent the number of individuals in each category based on Wolbachia wAnD strain infection status.

We then assessed Wolbachia densities by compared normalized levels of Wolbachia infection in the different body parts using a synthetic oligonucleotide standard and considering total DNA present in the qPCR reactions. A Mann-Whitney test indicated (as expected) significantly higher Wolbachia densities were present in the abdomens compared to cephalothoraxes ( P < 0.0001) ( Figure 4A). The mean normalised Wolbachia density (copies/ng DNA) in abdomens was 8.23 × 10 ³ (±1.81 ⁴) compared to 3.63 × 10 ² (±6.18 × 10 ²) for cephalothoraxes. This pattern did not significantly vary across years, and the level of infection remained higher in the abdomen 6.79 × 10 ³ (±1.83 ⁴) compared to the cephalothorax 3.25 × 10 ² (±7.42 × 10 ²) ( P < 0.0001) in samples collected in 2024. A Spearman’s rank correlation analysis between Wolbachia densities calculated per nanogram of the total DNA in the cephalothorax and the corresponding abdomens revealed a significant positive correlation ( r _s = 0.2531, P = 0.0046, Figure 4B) indicating a high Wolbachia density in the abdomen could result in high densities in the corresponding cephalothorax. In addition, high density of Wolbachia infection in the abdomen was found in samples with Wolbachia-positives cephalothorax compared to their negative counterparts ( P < 0.0001).

We determined whether malaria parasites reached the infectious sporozoite stage in wild An. demeilloni from our collection locations in the DRC by screening cephalothoraxes (a proxy for Salivary gland stage sporozoite infection). An overall prevalence rate of 1.3% (9/704) was found in all An. demeilloni samples. We then investigated whether the presence of Wolbachia wAnD correlated with the presence of Pl. falciparum in cephalothoraxes using only samples from Maziba and Lwiro that showed positive Wolbachia infections in their cephalothorax. The prevalence rate of Pl. falciparum was not significantly different between Wolbachia-positive individuals compared to Wolbachia-negative individuals ( P = 0.3630). However, interestingly all Pl. falciparum-infected mosquitoes were only found in Wolbachia-negative cephalothoraxes (8/560) ( Figure 5A). A power analysis performed showed a minimum sample size of 1106 mosquitoes would affect the development of Plasmodium sporozoites, from 1.4% in Wolbachia-negative individuals to below the limit of detection in Wolbachia-positive mosquitoes at a level of 0.05 and a power of 0.80. A total of 688 mosquitoes was analysed in our study, which is below the sample size required to detect a significant effect of Wolbachia affecting Pl. falciparum. We then investigated whether the presence of Wolbachia wAnD in the abdomen had any correlation to the presence of Pl. falciparum in the corresponding cephalothoraxes ( Figure 5B). We found significant difference in the prevalence rates of Pl. falciparum between individuals with Wolbachia abdomen positive (4/570) and individuals with Wolbachia abdomen negative (4/118) ( P = 0.0329).

Following Illumina 16S rRNA amplicon sequencing, all samples with fewer than 1000 resulting reads were removed from the dataset, leaving a total of 42 abdomen and 41 cephalothoraxes for further analysis. Wolbachia was identified in abdomens and cephalothoraxes at prevalence rates of 76% and 54% respectively ( Figure 6). As expected, there were no cephalothorax samples positive for Wolbachia for which the corresponding abdomen was not positive. Wolbachia was also found in higher relative abundances in abdomen samples than cephalothorax. For diversity analysis, all Wolbachia-assigned reads within the dataset were removed to assess the diversity within the remaining microbiome community given when present Wolbachia can dominate the microbiome. Overall alpha diversity, assessed using the Shannon’s and Simpson’s diversity indices, showed no significant difference between body parts ( Figure 7A,B). Bray-Curtis dissimilarity analysis was used to assess the beta diversity and there was no distinct grouping suggesting no significant difference in the microbiome composition between mosquito body parts despite a high level of inter-individual differences ( Figure 7C). Following the removal of all Wolbachia-assigned reads within the dataset, PERMANOVA analysis was used to assess the difference of microbiome composition between sample body parts and Wolbachia infection status. There was no significant difference between the composition of the abdomen and cephalothorax microbiome (F(1,82) = 2.136, P = 0.144). There was also no significant difference between Wolbachia-positive and Wolbachia-negative abdomen samples (F(1,40) = 0.138, P = 0.704). However, there was a significant difference in the composition of the microbiome between Wolbachia-positive and Wolbachia-negative cephalothorax samples (F(1,39) = 4.0305, P < 0.05).

Quantitative PCR data, additional microbiome analysis (with Wolbachia reads) and wsp FASTAs from Sanger sequencing are available at https://osf.io/y6sh8

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).