Prevalence and molecular characterization of Wolbachia in field-collected Aedes albopictus, Anopheles sinensis, Armigeres subalbatus, Culex pipiens and Cx. tritaeniorhynchus in China

Wolbachia are maternally transmitted intracellular bacteria that can naturally and artificially infect arthropods and nematodes. Recently, they were applied to control the spread of mosquito-borne pathogens by causing cytoplasmic incompatibility (CI) between germ cells of females and males. The ability of Wolbachia to induce CI is based on the prevalence and polymorphism of Wolbachia in natural populations of mosquitoes. In this study, we screened the natural infection level and diversity of Wolbachia in field-collected mosquitoes from 25 provinces of China based on partial sequence of Wolbachia surface protein (wsp) gene and multilocus sequence typing (MLST). Among the samples, 2489 mosquitoes were captured from 24 provinces between July and September, 2014 and the remaining 1025 mosquitoes were collected month-by-month in Yangzhou, Jiangsu province between September 2013 and August 2014. Our results showed that the presence of Wolbachia was observed in mosquitoes of Aedes albopictus (97.1%, 331/341), Armigeres subalbatus (95.8%, 481/502), Culex pipiens (87.0%, 1525/1752), Cx. tritaeniorhynchus (17.1%, 14/82), but not Anopheles sinensis (n = 88). Phylogenetic analysis indicated that high polymorphism of wsp and MLST loci was observed in Ae. albopictus mosquitoes, while no or low polymorphisms were in Ar. subalbatus and Cx. pipiens mosquitoes. A total of 12 unique mutations of deduced amino acid were identified in the wsp sequences obtained in this study, including four mutations in Wolbachia supergroup A and eight mutations in supergroup B. This study revealed the prevalence and polymorphism of Wolbachia in mosquitoes in large-scale regions of China and will provide some useful information when performing Wolbachia-based mosquito biocontrol strategies in China.


Introduction
Wolbachia is a genus of Gram-negative intracellular bacteria, belonging to the α-subphylum Proteobacteria, which is the most widely distributed symbiotic bacteria in the world and commonly found in arthropods and some nematodes. It was first reported by Hertig and Wolbach in reproductive tissues of Culex pipiens in 1924 [1]. In mosquitoes and many other arthropods, Wolbachia infection causes sperm-egg incompatibility called cytoplasmic incompatibility (CI) that leads to a reduction in egg-hatch frequencies when a Wolbachia-infected male mates with uninfected females [2]. The consequence of CI provides a reproductive advantage to infected females and contributes to the rapid spread of Wolbachia among the host population [3].
Mosquitoes are the most important vector of numerous arthropod-borne viruses (arboviruses), for example, Zika virus, dengue fever virus, West Nile virus and Japanese encephalitis virus. In order to slow the spread of mosquito-borne diseases, insecticides have been used extensively, but with limited success. At the same time, with the abuse of insecticides, mosquitoes are becoming more and more resistant. In the recent two decades, Wolbachia infection was demonstrated to reduce the potential transmission of mosquito-borne diseases by shortening adult lifespan, affecting mosquito reproduction and interfering with pathogen replication [4][5][6]. In 2013, Bian et al. established a stable infection of Wolbachia strain wAlbB in Anopheles stephensi, which is one of the most important vectors of malaria. The infection of wAlbB conferred resistance in the mosquitoes to the human malaria parasite Plasmodium falciparum [7].
The major lineages of Wolbachia have different host specificity and type of symbiosis, and are currently clustered into 20 supergroups (A-T) according to the order of their description [8,9], while only supergroups A and B have been found in mosquitoes. Initially, 16S rRNAbased PCR and sequencing were used for the detection and molecular characterization of Wolbachia [10], followed by the employment of the more variable Wolbachia surface protein (wsp) gene [11]. Subsequently, a reproducible and portable method, multilocus sequence typing (MLST), was established targeting five alleles (gatB, coxA, hcpA, ftsZ and fbpA) [12]. Although with limitations, it was gradually used for strain differentiation in the community of Wolbachia researchers. Bleidorn et al. recommended that whole-genome typing methods should be applied to the characterization of Wolbachia [8]. Several studies have instigated the prevalence and/or molecular characterization of Wolbachia in certain species of mosquitoes in China, particularly Aedes albopictus [13,14]. However, up to date, more than 45 species of mosquitoes have been identified in 31 provinces, distributed in the east, south, north, central, northeast, southwest and northwest regions of China. Among them, Ae. albopictus, Armigeres subalbatus, Cx. pipiens, Cx. quinquefasciatus and Cx. tritaeniorhynchus were most widely distributed [15]. Our knowledge of the prevalence and diversity of Wolbachia in mosquitoes of different species and geographic regions remains limited, especially in China. Therefore, in this study, we screened the prevalence and diversity of Wolbachia in approximately 3500 field-collected mosquitoes of Ae. albopictus, An. sinensis, Ar. subalbatus, Cx. pipiens and Cx. tritaeniorhynchus from 25 provinces of China with wsp-based real-time PCR and MLST analyses.

Ethics statement
All experimental protocols in this study were reviewed and approved by the Institutional Animal Care and Use Committee of Yangzhou University.

Mosquito collection
Convenience whole-body mosquito samples (n = 3514) were trapped in summer (July and September, 2014) in 24 provinces of China, and in the whole year (September 2013 to August 2014) in Yangzhou, Jiangsu province with hand nets, for an epidemiological survey of Rickettsia (Fig 1) [16,17]. Once captured, the mosquitoes were immediately protected separately in sterile tubes containing 400 μL of DNA/RNA stabilization reagent (Roche, Basel, Switzerland). To prevent cross contamination, disposable gloves were changed before the collection of each mosquito. Then the samples were transported at room temperature to the laboratory for species and gender identification and DNA extraction.

Mosquito species and gender identification
The species and genders of most mosquitoes employed in this study have been identified in our previous study [16]. The remaining samples were identified with standard morphological criteria (https://www.wrbu.si.edu/) with an ordinary microscope (Olympus, Tokyo, Japan). Each mosquito was removed from DNA/RNA stabilization reagent and placed on to a glass slide for morphological identification. In addition, the mosquitoes of each species from each province were subjected to further identification using a genomic approach targeting mitochondrial cytochrome c oxidase subunit I (COI) gene [18] (Table 1). After species and gender identification, the mosquitoes were stored at -80˚C until genomic DNA extraction.

Genomic DNA extraction and molecular detection of Wolbachia
After thawing at room temperature for approximate 15 min, the samples were homogenized individually with Precellys 24 homogenizer (Bertin Technologies, Montigny-le-Bretonneux, France) at 5800 x rpm (two short durations of 15 s with a 15-s break in between) and subjected to DNA extraction with High-Pure PCR Template Preparation Kit (Roche, Basel, Switzerland) following the manufacturer's instructions.
The presence of Wolbachia was screened using SYBR Green real-time PCR targeting a unique and highly conserved fragment of outer surface protein (wsp) gene [11], and the amplicon contained all the four hypervariable regions (HVR1-4) of wsp. In addition, the negative samples for wsp gene were further confirmed with the amplification of 16S rRNA gene [19] ( Table 1). The real-time PCR were performed in a Roche LightCycler 480 II real-time PCR platform (Roche, Basel, Switzerland) with 20-μL volumes composed of 10 μL of 2 x ChamQ Universal SYBR qPCR Master (Vazyme, Nanjing, China), 0.4 μL of forward primer, 0.4 μL of reverse primer, 2 μL of DNA template and 7.2 μL of double-distilled water (ddH 2 O) [20]. Thermal cycling consisted of a 2-min denaturation step at 94˚C followed by 37 cycles of 94˚C for 30 s, 53˚C (wsp) or 60˚C (16S rRNA) for 30 s and 72˚C for 1 min, final elongation at 72˚C for 10 min and final holding at 37˚C. The PCR products of wsp (~610 bp) with different melting temperatures were gel purified using a QIAquick Gel Extraction Kit (Qiagen, Valencia, USA) and sequenced by a commercial sequencing laboratory (Tsingke Biotechnology, Tianjin, China) with sanger dideoxy sequencing. Using the gel-purified PCR products as quantitative standards, the concentrations of DNAs were determined with the Quanti-iT PicoGreen dsDNA Assay Kit (Invitrogen Corporation, Carlsbad, USA). The 10-fold dilutions were performed to give the reaction system containing 100,000, 10,000, 1,000, 100 and 10 copies of the wsp gene of Wolbachia supergroup A or B, respectively. The DNA samples of Rickettsia bellii and R. felis were stored in our laboratory and served as negative controls.

Phylogenetic analysis with wsp and MLST
In order to investigate the molecular characterization of the mosquitoes captured from 25 provinces, one to four positive samples of each mosquito species in each province were selected to Wolbachia polymorphism identification with wsp gene sequencing and MLST targeting five alleles (gatB, coxA, hcpA, ftsZ and fbpA). A total of 109 Wolbachia-positive PCR products amplified from female or male mosquitoes of different species and different geographic regions were gel purified with the QIAquick Gel Extraction Kit (QIAGEN, Dusseldorf, Germany) and sequenced and verified by Sanger dideoxy sequencing (Sangon Biotech, Shanghai, China).
The five MLST loci were amplified with specific primers published in a previous study [12] (Table 1). PCR reaction were performed in a final volume of 25 μL containing 1 μL of DNA, 12.5 μL of 2 x Phanta Max Master Mix (Vazyme, Nanjing, China), 1 μL of forward primer, 1 μL of reverse primer and 9.5 μL of ddH 2 O. Thermal cycling consisted of a 2-min denaturation step at 94˚C, followed by 37 cycles of 94˚C for 30 s, optimal annealing temperature for 45 s (54˚C for gatB, coxA, hcpA, ftsZ and 59˚C for fbpA) and 72˚C for 1.5 min, final elongation at 72˚C for 10 min and final holding at 4˚C. The PCR products (gatB:~471 bp, coxA:~487 bp, hcpA:~515 bp, ftsZ:~524 bp and fbpA:~509 bp) were bi-directional sequenced by a commercial sequencing laboratory (Sangon Biotech, Shanghai, China) with sanger dideoxy method.
The sequences of partial wsp gene and MLST loci of Wolbachia obtained in this study were assembled with DNASTAR Lasergen 15.2 (DNASTAR Inc., Madison, WI) and aligned using CLUSTAL W in MEGA 7.0 (MEGA, Pennsylvania State University, University Park) along with those of Wolbachia strains represented different supergroups obtained from previous studies or GenBank [11,13]. A neighbor-joining phylogenetic tree was constructed using the Tamura-Nei model and the robustness of clusters was assessed by bootstrapping 1000 replicates. Maximum-likelihood phylogenetic analysis was performed to confirm the results [21].

PLOS NEGLECTED TROPICAL DISEASES
The sequences of five MLST loci of each strain were submitted to a public database (PubMLST) to obtain allelic profiles and sequence type (ST).

Statistical analysis
Statistical analyses were performed with the Statistica 7.0 software package (StatSoft Inc., USA). Positive rates of Wolbachia in different sampling time or locations were compared with the Chi-squared or Fisher's exact test. P-value<0.05 indicated a significant difference.

Presence of Wolbachia in mosquitoes from 25 provinces of China
All the mosquitoes captured in 25 provinces were screened for the presence of Wolbachia with wsp-based PCR. The detection limit of this PCR was 10 copies of the wsp gene per reaction for both Wolbachia supergroups A and B, and nonspecific amplifications were not observed (S1 Fig)

Molecular characterization of Wolbachia with wsp and MLST
The sequences of partial Wolbachia wsp gene identified in this study were submitted to Gen-Bank with the accession numbers of MW717996~MW718104. The neighbor-joining and maximum-likelihood phylogenetic trees both showed that 26.6% (29/109) and 73.4% (80/109) of the sequenced samples were clustered in supergroup A and B, respectively (Figs 2 and S2). The  = 4), and two sequences obtained from Ae. albopictus (Fig 2). Submissions to the MLST database of Wolbachia are currently suspended until a new curator can be appointed (https://pubmlst.org/organisms/wolbachia-spp). So for the strains with novel sequence types, we are unavailable to receive the ST codes immediately, and "ST-novel

PLOS NEGLECTED TROPICAL DISEASES
#" was marked in front of these strains instead of the ST codes (Fig 3). There were 12 different sequence types of the Wolbachia strains in this study, including four archived sequence types, ST-2 (n = 1), ST-9 (n = 56), ST-464 (n = 15) and ST-465 (n = 4), and eight novel ones (n = 33). The performance of MLST-based phylogenetic tree showed a high consistency with that constructed with wsp. Based on MLST, the Wolbachia strains identified in this study belonged to supergroups A (n = 30) and B (n = 79). Interestingly, a Wolbachia strain identified in an Ae. albopictus mosquito from Chongqing was clustered into supergroup A and B in the phylogenetic trees of MLST and wsp, respectively (Figs 2 and 3). There were three main clusters in both wsp and MLST-based phylogenetic trees, and it worth noting that most strains in these clusters corresponded to each other in the two phylogenetic trees (Figs 2 and 3).

Discussion
Approximately 3500 species of mosquitoes belonging to 41 genera are widely distributing in the world, of which at least 350 species have been found in China [22][23][24], and Ae. albopictus, An. sinensis, Ar. subalbatus, Cx. pipiens and Cx. tritaeniorhynchus were dominant populations [25,26]. In this study, we investigated the prevalence of Wolbachia among these four genera of mosquitoes captured in 25 provinces of China. Ae. albopictus, the major vector of arboviruses, is the natural host with a prevalence of almost 100% of Wolbachia [13,27,28]. In this study, the  (Fig 1, Tables 2 and 3). A few studies have demonstrated the presence of natural Wolbachia infection in Ar. subalbatus [29][30][31], however, the size and geographical diversity of samples were very limited. In this study, a total of 502 mosquitoes of Ar. subalbatus were collected in ten provinces, locating in eastern (Fujian, Jiangxi, Shandong and Zhejiang provinces), northern (Hebei province), northeastern (Liaoning province), central (Hubei and Hunan provinces), southwestern (Chongqing and Guizhou provinces) and southern (Guangdong province) of China. The overall infected rate of Wolbachia in Ar. subalbatus in this study (95.8%) is similar to the previous epidemiological survey conducted in Sri Lanka

PLOS NEGLECTED TROPICAL DISEASES
(100%) [30], indicating that Ar. subalbatus is the naturally harbor of Wolbachia (Fig 1 and Table 2). Cx. pipiens complex mosquitoes are widely distributed throughout China, including at high altitudes, such as Tibet [16,32,33]. Although considered as an important disease vector, there is little information on the distribution and diversity of Wolbachia in this kind of mosquitoes in China. In this study, samples of Cx. pipiens mosquitoes were identified in 23 out of 25 sampled provinces (except for Guizhou and Hunan provinces). The prevalence of Wolbachia varied greatly among these mosquitoes from different provinces. For example, only one mosquito of Cx. pipiens (1.4%, 1/71) was infected in Guangdong province, while there were 100% positive rates in Chongqing, Henan, Liaoning, Shanxi, Sichuan and Yunnan provinces (Fig 1, Tables 2 and 3). We noticed that Cx. pipiens is not a native mosquito species in Guangdong [15]. Was there a research group that artificially released mosquitoes, resulting in a low rate of Wolbachia infection? Cx. tritaeniorhynchus was only identified in the samples from Guizhou province. Our results showed a prevalence of 17.1% for Wolbachia in natural  Table 2). For decades, Anopheles mosquitoes were considered resistant to Wolbachia. Although have been experimental infected in the laboratory, natural infections of Wolbachia in wild-collected Anopheles mosquitoes have barely been reported [34,35]. Anopheles mosquitoes were always considered to be exempt of Wolbachia until some recent studies showed evidence of natural  [41]. However, in this study, the mosquitoes of An. sinensis sampled from ten provinces in this study were all negative for Wolbachia by the screening of wsp and 16S rRNA genes ( Table 2). The temporal distribution of Wolbachia in arthropods is largely unknown. In this study, mosquitoes were captured at the first day of every month from September 2013 to August PLOS NEGLECTED TROPICAL DISEASES 2014, except for January and February, when the temperature was below zero degrees Celsius in average. Although the information on the species and genders of the captured mosquitoes was only identified and noted from April to August, 2014, it was available for us to find a significant reduction of the infected rates of Wolbachia in Cx. pipiens in July (P<0.01), in both females and males. Previous studies have indicated that the rearing temperature have significant effects on the fertility, frequency and density of mosquitoes and other arthropods infected with Wolbachia, especially wMel and wAlbB [42][43][44][45][46][47]. Further studies are needed to determine the effect on the infection of Wolbachia at more extreme temperature conditions (Table 3).
Although do not perform particularly well in the differentiation of closely related Wolbachia genomes [8], compared with 16S rRNA and cell division protein gene (ftsZ), phylogenetic tree of improved resolution can be constructed with wsp gene, which is evolving at a much faster rate [11]. Molecular phylogeny based on wsp sequences revealed that most Wolbachia infection in Ae. albopictus clustered with wAlbB (supergroup B) (88.9%, 24/27), while 11.1% (3/27) belonged to wAlbA (supergroup A). The genetic variation was observed within both wAlbA and wAlbB strains, including some strains sampled from the same location (Figs 2 and 4). Compared with previous studies, our results suggested that Wolbachia could have invaded and spread throughout populations of these mosquitoes for a long time [48][49][50][51]. The information of the diversity of Wolbachia infection in Ar. subalbatus is quite limited. Nugapola et al. reported the presence of Wolbachia supergroup A in two mosquitoes in Sri Lanka in 2017 [30]. To the best of our knowledge, this was the only published study regarding the diversity of Wolbachia infection in Ar. subalbatus, until we demonstrated the presence of Wolbachia supergroup A in 11 provinces of China. Interestingly, no polymorphism was observed among these strains based on both wsp and MLST loci (Figs 2 and 3). In Cx. pipiens, 98.1% (51/52) of samples were infected with wPip with non-polymorphism based on wsp sequence, except for an individual collected from Heilongjiang province (Fig 2). Further studies are needed to investigate the diversity of wPip with higher-resolution molecular biological assays, for example, whole-genome typing methods, more accurate markers and PCR-restriction fragment length polymorphism (RFLP) [8,52]. Bleidorn and Gerth have provided a characterization of 252 single copy loci for the application of few loci approaches, and some of these loci have the potential to perform better than wsp and MLST. However, some of these loci may not be harbored by mosquitoes or have no polymorphism in mosquitoes. Further studies are still needed to screen accurate markers from these loci. Although Rickettsia were screened in Cx. tritaeniorhynchus with a low prevalence [53], this mosquito specie was not considered to harbor natural Wolbachia infections [54,55]. Interestingly, in this study, Wolbachia-positive samples were detected in Cx. tritaeniorhynchus mosquitoes in Guizhou province and identical with wPip strains (Figs 1-3, S2 and Table 2). Although there were a few outliers, molecular phylogeny based on five MLST loci was generally consistent with that of wsp (Figs 2 and 3).
The present study identified a total of 12 unique mutations of the deduced amino acid (AA) with wsp gene sequencing and aligning, including four mutations in Wolbachia supergroup A and eight mutations in supergroup B, respectively. Two insertions (G 179 and TKDA [187][188][189][190] were identified among the three strains (MW718054, MW718057 and MW718097) of Wolbachia supergroup A. Both the phylogenetic tree and nucleotide/ AA alignment showed that these three strains were mostly closed to the reference sequence (KJ140127) with low polymorphism, although the mosquito vectors were sampled from three far separate provinces (Tianjin, Guangdong and Sichuan provinces) (Figs 1, 2, 4 and 5).