Wolbachia and mosquitoes: Exploring transmission modes and coevolutionary dynamics in Shandong Province, China

Vector-borne diseases leave a large footprint on global health. Notable culprits include West Nile virus (WNV), St. Louis encephalitis virus (SLEV), and Japanese encephalitis virus (JEV), all transmitted by Culex mosquitoes. Chemical insecticides have been widely used to reduce the spread of mosquito-borne diseases. Still, mosquitoes are becoming more and more resistant to most chemical insecticides which cause particular harm to the ecology. Wolbachia belongs to the family Ehrlichiaceae in the order Rickettsiales and is a matrilineally inherited endosymbiont present in 60% of insects in nature. Wolbachia is capable of inducing a wide range of reproductive abnormalities in its hosts, such as cytoplasmic incompatibility, and can alter mosquito resistance to pathogen infection. Wolbachia has been proposed as a biological alternative to chemical vector control, and specific research progress and effectiveness have been achieved. Despite the importance of Wolbachia, this strategy has not been tested in Culex pipiens pallens, the most prevalent mosquito species in Shandong Province, China. Little is known about how the mass release of Wolbachia-infected mosquitoes may impact the genetic structure of Culex pipiens pallens, and how the symbiotic bacterium Wolbachia interacts with mitochondria during host mosquito transmission. Based on the population genetic structure of Culex pipiens pallens in Shandong Province, this study investigated the infection rate and infection type of Wolbachia in Shandong Province and jointly analysed the evolutionary relationship between the host mosquito and the symbiotic bacterium Wolbachia. Our study showed that Wolbachia naturally infected by Culex pipiens pallens in Shandong Province was less homologous to Wolbachia infected by Aedes albopictus released from mosquito factory in Guangzhou. Our results also show that Culex pipiens pallens is undergoing demographic expansion in Shandong Province. The overall Wolbachia infection rate of Culex pipiens pallens was 92.8%, and a total of 15 WSP haplotypes were detected. We found that the genetic diversity of Wolbachia was low in Culex pipiens pallens from Shandong Province, and the mosquitoes were infected only with type B Wolbachia. Visualizing the relationship between Culex pipiens pallens and Wolbachia using a tanglegram revealed patterns of widespread associations. A specific coevolutionary relationship exists between the host mosquito and Wolbachia. Knowledge of this mosquito–Wolbachia relationship will provide essential scientific information required for Wolbachia-based vector control approaches in Shandong Province and will lead to a better understanding of the diversity and evolution of Wolbachia for its utility as a biocontrol agent.


Introduction
Culex pipiens pallens, the most abundant Culex mosquito species in northern China, is an essential vector for Bancroftian filariasis, Japanese encephalitis virus (JEV), and potentially West Nile virus(WNV) [1].Culex-mediated disease outbreaks are increasing, posing a significant public health challenge [2].Japanese encephalitis (JE) is a common viral encephalitis in Asia, with an annual incidence of 70,000 cases and 15,000 deaths; China accounts for 50% of the reported JE cases [3].In 2018, the number of WNV infections spiked in 11 European Union/European Economic Area member states, with 1605 human cases, including 166 lethal cases [4].Therefore, understanding the population genetic dynamics and genetic structure of Cx. p. pallens mosquitoes is essential to assessing their ability to transmit diseases and developing practical vector surveillance tools.
Shandong Province is located in eastern China along the lower reaches of the Yellow River, serving as a primary coastal province.The region exhibits a temperate monsoon climate that is often warm and rainy, making it well-suited to mosquito breeding.Historically, mosquitoborne infectious diseases such as malaria, JE, and filariasis have been prevalent in Shandong Province.For example, in the 1960s and 1970s, there were two major malaria outbreaks in Shandong Province, with more than 6 million and 4.6 million malaria cases, respectively [5].In the early days, following the establishment of the People's Republic of China in 1949, more than 5 million people were infected with microfilariae nationwide, and more than 2.5 million symptomatic patients were reported in Shandong Province [6].Climate change may have accelerated the northward expansion of dengue outbreaks in China.Shandong Province is the northernmost focal point of dengue fever cases diagnosed in China [7], underscoring its epidemiological significance.
Mosquito control is essential for preventing vector-borne disease transmission in the human population.Insecticides are the primary weapons for mosquito control, and they have been extensively used to reduce the spread of mosquito-borne diseases, but with limited success.Meanwhile, mosquitoes are becoming increasingly insecticide-resistant due to the negligent and improper use of insecticides [8].In light of these concerns, biological approaches are called upon as alternatives to chemical control [9].For example, Wolbachia can reduce mosquito density and the spread of mosquito-borne diseases.
Wolbachia is a gram-negative symbiotic bacterium which is parasitic in invertebrates and can be transmitted through eggs.It is widely found in arthropods [10].It is estimated that about 65% of insect species and 28% of mosquito species naturally carry Wolbachia [11].The use of symbiotic Wolbachia for insect-borne disease control is primarily based on its ability to induce cytoplasmic incompatibility (CI) and resistance to pathogens.CI is a characteristic of Wolbachia infection that affects insect reproduction; eggs fertilized by a male carrying Wolbachia and laid by a female not carrying Wolbachia, or carrying a different type of Wolbachia, will not hatch [12,13].Wolbachia can control mosquito populations using population suppression and population replacement [14].In population suppression, male mosquitoes carrying Wolbachia are released into areas with uninfected or differently infected mosquitoes.This results in cytoplasmic incompatibility, leading to a decrease in the mosquito population.The targeted mosquito species can be eradicated by introducing many Wolbachia-infected male mosquitoes, leading to long-term infection.Another approach is to introduce female mosquitoes infected with a new Wolbachia strain, which resist the target mosquitoes, causing one-way cytoplasmic incompatibility.By releasing a sufficient number of these female mosquitoes and allowing them to multiply over several generations, Wolbachia can spread to all targeted mosquitoes [15].This would result in replacing the original disease-transmitting mosquito vectors with new disease-resistant mosquito strains, thereby achieving the goal of interrupting disease transmission.
Mosquito vector control technology utilizing the symbiotic bacterium Wolbachia is a promising solution.Successful field trials across China [16], the United States [17], and Australia [18] have demonstrated superior efficiency, environmental protection, and sustainability.As an innovative approach to vector control, Wolbachia has been introduced into Aedes albopictus mosquitoes in Guangzhou, China, through microinjection, significantly suppressing the wild population by over 90%.The application of the insect symbiotic bacterium Wolbachia for mosquito vector control is based on its ability to induce CI.Controlling the Wolbachia type carried by the target population and establishing new Wolbachia mosquito strains are the basis of the application.However, it is unclear whether the mass release of Wolbachia-infected mosquitoes will impact other mosquito species, such as Cx.p. pallens.Understanding the effects of the mass release of Wolbachia-infected mosquitoes on the prevalence of Wolbachia, host population genetic structure, and the dynamics of gene flow patterns is crucial for assessing this population modification strategy's long-term effectiveness and sustainability.Further research and surveillance are necessary to monitor any changes in Wolbachia prevalence and evaluate the impact of this intervention on disease transmission dynamics in Guangzhou and other areas implementing similar strategies.
Preparing incompatible vector insects often requires a lot of laboratory work, time, and technical expertise from the experimenter, who performs tasks such as removing the initially infected Wolbachia from the vector insect population and artificially infecting another type of Wolbachia, repeatedly backcrossing two vector insect populations (e.g., several generations of Wolbachia-infected and uninfected populations), and so on.The application of the resulting incompatible quins depends on the long-term stability of the new insect-Wolbachia symbiosis.Through long-term selection and evolution, natural populations are more stable than artificially created new mosquito populations.Therefore, whether natural populations can be directly used as incompatible insects for insect vector population control is worthy of further study.
Our research aimed to explore the following: i) Investigate Wolbachia's infection rate and distribution pattern in Cx. p. pallens in Shandong Province and construct a phylogenetic tree for typing and tracing; ii) Utilize mitochondrial cytochrome c oxidase subunit I (COI) gene fragments to assess the population genetic structure and spread of Cx. p. pallens and predict disease risk; and iii) Examine the evolutionary relationship between Wolbachia bacteria and the mitochondria of Cx. p. pallens to understand Wolbachia's contribution to shaping genetic diversity in Cx. p. pallens.This research will help to assess the roles of Cx. p. pallens and Wolbachia in disease transmission and develop more effective management strategies for controlling mosquito-borne diseases.

Sample collection and preparation
Culex p. pallens were collected between July and September 2022 from 11 locations in Shandong Province (Fig 1 and Table 1).Captured mosquitoes were immediately preserved in separate sterile tubes containing 200 μl of RNAprotect Tissue Reagent.The samples were transported to the laboratory at room temperature for species identification and DNA extraction.

Mosquito species identification and DNA extraction
Mosquito species were identified microscopically based on their morphological characteristics (https://www.wrbu.si.edu/) (Olympus, SZX7, Japan) and further confirmed using molecular marker analysis.Molecular identification was achieved with COI and internal transcribed spacer 2 (ITS2) barcoding, using BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to determine sequence percentage identity.DNA was extracted from individual mosquitoes using the cador Pathogen 96 QIAcube HT Kit following the manufacturer's directions, and stored at -80˚C.

PCR identification of Wolbachia infections
Detection of the Wolbachia endosymbiont in mosquitoes was performed using PCR-based molecular approaches using the most common Wolbachia-specific DNA marker (the WSP gene), with the forward primer 81F (5 0 -TGGTCCAATAAGTGATGAAGAAAC-3 0 ) and reverse primer 691R (5 0 -AAAAATTAAACGCTACTCCA-3 0 ) [19].PCR reaction mixtures comprised 25 μl of 2× Phanta Max Master Mix, 1 μl each of 10-μM forward and reverse primers, 2 μl of template DNA, and nuclease-free water to a final volume of 50 μl.PCR conditions were as follows: 95˚C for 2 min; 35 cycles of 95˚C for 30 s, 54˚C for 45 s, and 72˚C for 1 min; and 72˚C for 5 min.PCR products were separated electrophoretically and sequenced at Sangon Biotech (Shanghai, China).

Amplification and sequencing of the COI gene
The COI gene of the mitochondrial genome was amplified using primers LCO 1490 (5 0 -GGTCAACAAATCATAAAGATATTGG-3 0 ) and HCO 2198 (5 0 -TAAACTTCAGGGTGAC-CAAAAAATCA-3 0 ) to determine population genetic structure [20].PCR reaction mixtures comprised 25 μl of 2× Phanta Max Master Mix, 1 μl each of 10-μM forward and reverse primers, 2 μl of template DNA, and nuclease-free water to a final volume of 50 μl.PCR conditions were as follows: 94˚C for 1 min; 5 cycles of 94˚C for 40 s, 45˚C for 40 s, and 72˚C for 1 min; 30 cycles of 94˚C for 40 s, 53˚C for 40 s, and 72˚C for 1 min; and 72˚C for 5 min.The PCR products were sequenced by Sangon Biotech (Shanghai, China).

Data analysis
All sequences were manually aligned, checked and edited using BioEdit version 7.0 and compared with other sequences available in the GenBank database to determine the percentage identity using BLAST.The most similar sequences were downloaded for phylogenetic analysis.Specimens showing more than 99% nucleotide sequence identity with the available species sequences in the database were considered Cx.p. pallens.Based on the mitochondrial COI gene sequences, we used Arlequin v3.5 [21] and DnaSP v6 [22] to calculate the number of segregating sites (s), haplotypes (h), haplotype diversity (Hd), and nucleotide diversity (pi), and neutrality tests [23], namely Tajima's D test and Fu's Fs test, to investigate the genetic diversity of Cx. p. pallens [24].To determine the genetic structure of the mosquito populations, we used analysis of molecular variance (AMOVA) to partition the genetic variation among groups (Fct), populations within groups (Fsc), and populations among groups (Fst) [25].Pairwise Fst and gene flow (Nm) values were calculated for all the populations.In addition, haplotype networks of Cx. p. pallens were created using the haplotypes network (TCS) method in PopART 1.7 to visualise the relationships among populations.Geographical population structure was evaluated using the Bayesian clustering method in STRUCTURE v.2.3, which identifies the most probable number of genetic clusters (K) and assigns individuals to these clusters.Different runs were conducted using different datasets for further clustering.Subsequently, COI and WSP sequences were aligned using ClustalX, and a neighbour-joining phylogenetic tree model based on genetic distance values was created using MEGA X software [26].

Polymorphism of the mitochondrial gene COI sequence
A total of 1782 COI sequences were generated from the 11 populations.The COI sequences were aligned, yielding a total length of 603 bp and 27 variable sites.The overall haplotype diversity (Hd) and nucleotide diversity (Pi) were 0.352 and 0.98 × 10 ˗2 , respectively (Table 2).Among the 11 populations, Cx. p. pallens from Dezhou exhibited the highest haplotype diversity, and that from Yantai exhibited the highest nucleotide diversity (Table 2).In general, Cx. p. pallens from coastal cities exhibited higher COI haplotype and nucleotide diversity than those from other cities.Based on differences in the nucleotide composition of the COI gene, we identified 26 mitochondrial haplotypes in the 11 studied populations, which were denoted as H01-H26 (Fig 3).A haplotype network graph was constructed using the TCS method.The size of each circle is proportional to its corresponding frequency.Different colours indicate different populations.The number of small diagonal lines refer to variable asynchronous numbers.

Haplotype network reconstruction of Cx. p. pallens
The statistical parsimony network was generated with TCS using COI sequences of Cx. p. pallens, identifying 26 haplotypes across the sampled sites (Fig 3).H01 and H06 were the most frequently observed haplotypes.H01 was the central haplotype that was highly connected to the haplotype lines and was the only haplotype found in all the localities in Shandong Province.Nearly all other haplotypes originated from H01 through one or more mutations.The haplotype distribution showed that H01, H05, H06, H08, H11, H13 and H15 were shared among multiple populations, whereas the remaining haplotypes were only observed in a single population.

Population genetic structure
A total of 891 individuals from 11 populations were analysed for population genetic diversity.Pairwise comparisons indicated that Fst values were significantly different from zero, with the lowest value (0.00109) found between Binzhou and Linyi and the highest value (0.54103) between Qingdao and Heze (Table 3).Generally, highly significant pairwise population differentiation was observed between Qingdao and most other populations (Fst > 0.25).However, frequent gene flow occurred in other populations (Nm > 1).The AMOVA showed that the majority of genetic variance occurred within populations (87.91%) (Table 4).The total Fst was 0.13135 (P > 0.05) and Nm was 1.65, reflecting moderate population differentiation.The Fu's Fs and Tajima's D values for the Cx.p. pallens populations in Shandong Province were almost entirely negative, except in Qingdao (Table 2).Tajima's D value (˗1.97) and Fu's Fs value (˗27.58) for the overall population were significant, suggesting a high number of lowfrequency mutations and a demographic expansion of Cx. p. pallens in Shandong Province.Among specific populations, Jining and Binzhou had significantly negative D values, whereas Dongying and Zibo exhibited significant negative Fs values.Tajima's D and Fu's Fs tests revealed that Jining, Binzhou, Dongying, and Zibo were significantly negative, suggesting recent population expansion or selection.Multilocus cluster Bayesian analysis of all 11 population samples indicated the genetic structure among Cx.p. pallens populations (K = 2-11) (S2A Fig) [27].The results show that when K = 4, a larger Delta K value is obtained, and when K = 7, the Delta K value is the largest.Therefore, the optimal K value should be 4-7 (S2B Fig) .Based on the geographical distribution and genetic differentiation within Shandong Province, we divided the 11 populations into five groups: the Northwest Shandong Plain (Dezhou), the Jiaolai Hills (Yantai, Qingdao), the Southwest Shandong Plain (Liaocheng, Heze), the alluvial plain of the Yellow River (Dongying, Binzhou), and the south-central Shandong Plain (Zibo, Linyi, Rizhao, Jining).

Migration and gene flow patterns
LAMARC analysis revealed that historical gene flow rates ranged from 0.88 to 96.02.The highest migration rates were detected among neighbouring populations within each of the five locality groups (Fig 4).Moderate migration levels were observed in the five locality groups.Based on the coalescent analysis, migration was asymmetrical (Fig 4).

Evolutionary relationships between mosquitoes and Wolbachia
We recorded mosquito-Wolbachia associations in wild-caught mosquitoes from Shandong Province.Visualisation of these associations using a tanglegram revealed patterns of broad associations (S3 Fig) .The distance-based quantitative test showed no significant consistency between the mosquito and Wolbachia phylogenies at the global level (ParaFit Global test: Par-aFit Global < 0.001, P = 1).Among the host-endosymbiont links, the associations between H07 and Wol 08, H07 and Wol 12, H07 and Wol 13, and H15 and Wol07 were statistically significant (Table 5).

Detection of Wolbachia infection and its distribution in wild mosquitoes
In this study, we assessed the prevalence of Wolbachia in Cx. p. pallens collected from Shandong Province.The overall infection rate of all tested mosquitoes was 92.8%, indicating that Wolbachia was widespread in Cx. p. pallens.
The WSP molecular marker enabled successful detection of Wolbachia infection across numerous taxa, strain genotyping, and evolutionary comparison between the detected Wolbachia strains.In this study, a low genetic diversity of Wolbachia strains was found in Cx. p. pallens, and natural populations were infected only with type B Wolbachia, not with type A. Studies have shown that animal habitats such as wetlands and forests tend to have relatively high infection rates of mosquito-borne viruses and bacteria [28].However, Wolbachia infection rates of Cx. p. pallens in the Yellow River Delta wetland of Dongying, an important habitat for migratory birds [29], were lower than those in other places.This presents a new opportunity for the application of Wolbachia to mediate the transmission of mosquito-borne diseases by Cx. p. pallens.

The impact of the Guangzhou"mosquito factory" on mosquitoes in other areas
The centrepiece of the "mosquito factory" is a colony of Wolbachia mosquitoes, called the brood stock, from which all future populations of Wolbachia mosquito offspring are bred.
Owing to the role of cytoplasmic incompatibility, the release of Wolbachia-infected male mosquitoes into the wild has become a promising strategy for suppressing wild mosquito populations [30].However, more attention should be paid to the effects of Wolbachia on hosts'  behavioural patterns to evaluate the environmental impact of releasing these Wolbachiainfected insects into the field [31].Notably, a large number of Aedes albopictus infected with Wolbachia were released in Guangzhou.So, will the increase or decrease in the abundance of the target species lead to the disruption of the original ecological balance.It is generally agreed that the success of Wolbachia is due to its ability to switch from one host species to another [32].This raises further questions about its effects on other species.The population of Cx. p. pallens is expanding in Shandong Province.We believe that apart from assessing the geographic regions where climate change, the taxa of hosts and transmission types of pathogens, future research should also consider the potential effects of Wolbachia more carefully.Finally, we sought to inform future research avenues, policy, and practices via the trends and impacts identified herein.
In the present study, we found that the Wolbachia infection rate of Cx. p. pallens in Shandong Province is high, and phylogenetic tree analysis shows that it has low homology with Wolbachia infection of Aedes albopictus in Guangzhou, but high homology with Wolbachia infection of Culex mosquitoes in Guangzhou.Questions about the influence of Wolbachia infection on host behaviour under field conditions, and consequently on the ecosystem, remain to be addressed.The molecular mechanisms by which Wolbachia affect host behaviour also need to be elucidated in more detail.The effects of Wolbachia on host fitness traits may be multidimensional, and a number of host genes, miRNAs and proteins may be modified by Wolbachia infection.

Geographical isolation of Cx. p. pallens populations from Qingdao
In the present study, we analysed the population genetics of Cx. p. pallens collected from different regions of Shandong Province based on the COI gene.Most (87.91%) of the genetic variation occurred within individuals, whereas only 12.09% of the variation was detected among populations.Our data showed significant genetic differentiation and limited gene flow between the Qingdao population and populations from most other cities.In contrast to other cities within Shandong Province, Qingdao is primarily situated on hilly terrain that slopes from east (Laoshan Mountain Group) to west (Jiaozhou Bay), with rolling hills in the north and the Fushan Mountain Range to the south along the Pacific coast.The expansive Pacific Ocean thus forms a natural geographical barrier between Qingdao and other cities within the province.In addition to Qingdao, other cities also have varying geographical environments.Therefore, precise countermeasures should be adopted for mosquito vector prevention and control.

Demographic expansion of Cx. p. pallens in Shandong Province
The demographic expansion of Cx. p. pallens closely corresponded to the COI haplotype network.The haplotype profiles were star-shaped, reflecting their recent appearance and rapid population growth.The neutrality test results were significantly negative, which also supports this phenomenon.A strong link exists between demographic expansion and climatic variability.Meteorological factors, including temperature, humidity, and precipitation, have a significant impact on the number, density, and distribution of disease vectors, as well as the spatial and temporal dynamics, epidemic frequency, and intensity of vector-borne diseases.In addition, cluster analysis revealed that Cx. p. pallens community structure existed spatial and temporal heterogeneity to some extent.Mosquito development and survival and viral replication depend on environmental conditions, particularly climatic conditions.A study has shown that an increase in temperature was observed from 1970 to 2021 in most places in China, and annual change rates varied substantially among different sites, from ˗0.22˚C/year to 0.58˚C/ year [7].Notably, Shandong Province reported an autochthonous dengue case for the first time on 16 August 2017.Ninety-five cases were subsequently reported across the entire province, and 79 indigenous cases occurred in Jining, the northernmost region where local dengue fever cases have been detected [33].Therefore, upgrading surveillance systems for vectors and vector-borne diseases in the context of climate change will strengthen research on risk assessments, predictions, early warnings, control strategies, and intervention measures to effectively cope with the new challenges of climate-sensitive vector-borne diseases.

Evolutionary relationships between mosquitoes and Wolbachia
Symbiosis is a major force of evolutionary change that influences virtually all aspects of biology, from population ecology and evolution to genomics and molecular/biochemical mechanisms of development and reproduction.A broad association pattern between mosquitoes and Wolbachia strains was observed based on the tanglegram (S3 Fig) .A previous study reported that Aedes mosquitoes tended to be associated with Wolbachia supergroup A, whereas other mosquito species, particularly those of the genus Culex, were largely associated with Wolbachia supergroup B. This indicated that closely related Wolbachia strains are likely to establish themselves in related hosts.In the present study, we found that all Wolbachia infections in Cx. p. pallens were associated with Wolbachia supergroup B, which is consistent with the above results.Significant coevolution was detected between Wol08, Wol12, and Wol13 and haplotype H07, and between Wol07 and haplotype H15.These regional variations in mosquito-Wolbachia interactions may represent an ongoing evolutionary process, or infections may occur by chance or be associated with local environments.Therefore, further investigation is warranted.Understanding Wolbachia host specificity has significant implications, particularly for optimising Wolbachia biocontrol strategies.In addition to selecting strains that can effectively limit pathogen replication, strains should also be selected based on their host specificity.
This study has potential limitations, including a) differences in sampling location and environment, which may lead to potential sample size bias and not fully represent the actual genetic diversity of the population; and b) the existence of a specific coevolutionary relationship between the haplotype of some host mosquitoes and the symbiotic Wolbachia.The significance of this relationship for precise control should be studied in the next step.

Conclusions
Based on the results of population genetic structure and Wolbachia infection, this study has provided valuable insights for controlling Cx. p. pallens in Shandong Province.From this analysis, it is evident that the population of Cx. p. pallens is expanding, increasing the risk of human exposure to vector mosquitoes and infection with arboviruses, and thus the need to strengthen preventive and early warning measures for related diseases.Additionally, the geographical isolation between the Qingdao population and the majority of populations in Shandong Province suggests that more precise control measures should be implemented based on local conditions.Phylogenetic analysis showed low homology between Wolbachia strains infecting Aedes albopictus released from mosquito factories in Guangzhou and Wolbachia strains naturally infecting Cx. p. pallens.This finding underscores the importance of understanding Wolbachia's non-changing host characteristics for future mosquito control strategies.The investigation based on the co-evolutionary relationship between Wolbachia and Cx.p. pallens mitochondrial markers suggests that the host specificity of Wolbachia is crucial in optimizing its biocontrol strategy.Therefore, in addition to considering strains that can effectively limit pathogen replication, the selection of strains should be based on host specificity to implement biological control.This study provides a valuable reference for the scientific and accurate control of mosquito-borne diseases in the future.

Fig 4 .
Fig 4. Bayesian estimates of historical asymmetrical migration between populations of Culex pipiens pallens.The five locality groups are indicated by dotted circles, and arrows indicate the direction of migration.The base layer of the map is from the Resource and Environment Science and Data Center (https://www.resdc.cn/).https://doi.org/10.1371/journal.pntd.0011944.g004