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-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 rRNA-based 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.

Materials and methods

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.

Fig 1

The distributions of field-collected mosquitoes in this study.

A total of 3514 samples of mosquitoes’ whole bodies (convenience samples) were collected in 25 provinces of China, including Aedes albopictus (n = 349), Anopheles sinensis (n = 88), Armigeres subalbatus (n = 502), Culex pipiens (n = 2495) and Cx. tritaeniorhynchus (n = 80). Different colored dots represented different species of mosquitoes, including Ae. albopictus (brown), An. sinensis (blue), Ar. subalbatus (red), Cx. pipiens (green) and Cx. tritaeniorhynchus (purple). The numbers in the region of each province indicated the positive/ total numbers of mosquitoes.

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.

Table 1

Primer pairs used in this study.

Gene	Primer Sequence (5’- 3’)	Annealing Temperature	Amplicon Length (bp)	Reference
COI	GGTCAACAAATCATAAAGATATTGG	51°C	~650	[16]
	TAAACTTCAGGGTGACCAAAAAATCA
wsp	TGGTCCAATAAGTGATGAAGAAAC	53°C	~610	[11]
	AAAAATTAAACGCTACTCCA
16S rRNA	CATACCTATTCGAAGGGATAG	60°C	~438	[19]
	AGCTTCGAGTGAAACCAATTC
gatB	GAKTTAAAYCGYGCAGGBGTT	54°C	~471	[12]
	TGGYAAYTCRGGYAAAGATGA
coxA	TTGGRGCRATYAACTTTATAG	54°C	~487	[12]
	CTAAAGACTTTKACRCCAGT
hcpA	GAAATARCAGTTGCTGCAAA	54°C	~515	[12]
	GAAAGTYRAGCAAGYTCTG
ftsZ	ATYATGGARCATATAAARGATAG	54°C	~524	[12]
	TCRAGYAATGGATTRGATAT
fbpA	GCTGCTCCRCTTGGYWTGAT	59°C	~509	[12]
	CCRCCAGARAAAAYYACTATTC

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 (ddH2O) [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 ddH2O. 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]. 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.

Results

Mosquito species and gender identified

The identification of species and genders were performed on most mosquitoes (77.0%, 2706/3514) employed in this study, revealing the presence of five species of mosquitoes, including Ae. albopictus (n = 341, female = 258, male = 83), An. sinensis (n = 88, female = 62, male = 26), Ar. subalbatus (n = 502, female = 236, male = 266), Cx. pipiens (n = 1752, female = 992, male = 760) and Cx. tritaeniorhynchus (n = 82, female = 72, male = 10) (Table 2) with the accession numbers of OK465203~OK465358. Species identification was not performed on the samples collected in Jiangsu province from September 2013 to March 2014, because the limited number of researchers at that time.

Table 2

The presence and positive rates of Wolbachia in mosquitoes collected from 24 provinces of China.

Mosquito species	Province	Coordinate	Positive of Wolbachia
			Female	Male	Total
Ae. albopictus	Anhui	31.5°N, 117.2°E	100% (2/2)	NAa	100% (2/2)
	Chongqing	29.3°N, 106.3°E	100% (10/10)	NA	100% (10/10)
	Fujian	26.1°N, 119.2°E	100% (1/1)	NA	100% (1/1)
	Guangdong	23.1°N, 113.1°E	80.0% (8/10)	100% (3/3)	84.6% (11/13)
	Guangxi	22.5°N, 108.2°E	100% (17/17)	NA	100% (17/17)
	Hebei	38.0°N, 114.3°E	100% (30/30)	100% (1/1)	100% (31/31)
	Henan	34.5°N, 113.4°E	100% (15/15)	NA	100% (15/15)
	Hubei	30.4°N, 114.2°E	100% (10/10)	100% (3/3)	100% (13/13)
	Jiangsu	32.0°N, 118.5°E	92.6% (25/27)	95.0% (19/20)	93.6% (44/47)
	Jiangxi	28.4°N, 115.6°E	94.7% (18/19)	NA	94.7% (18/19)
	Liaoning	41.5°N, 123.3°E	40.0% (2/5)	100% (3/3)	62.5% (5/8)
	Shaanxi	34.2°N, 108.6°E	100% (3/3)	NA	100% (3/3)
	Shandong	36.4°N, 117.0°E	100% (55/55)	100% (42/42)	100% (97/97)
	Shanxi	37.5°N, 112.3°E	100% (1/1)	NA	100% (1/1)
	Sichuan	30.4°N, 104.0°E	100% (25/25)	NA	100% (25/25)
	Tianjing	39.0°N, 117.1°E	95.0% (19/20)	100% (7/7)	96.3% (26/27)
	Zhejiang	30.2°N, 120.1°E	100% (8/8)	100% (4/4)	100% (12/12)
	Total		96.5% (249/258)	98.8% (82/83)	97.1% (331/341)
An. sinensis	Beijing	39.6°N, 116.2°E	0% (0/2)	NA	0% (0/2)
	Fujian	26.1°N, 119.2°E	0% (0/14)	0% (0/8)	0% (0/22)
	Guizhou	26.4°N, 106.4°E	0% (0/15)	0% (0/10)	0% (0/25)
	Heilongjiang	45.4°N, 126.4°E	0% (0/3)	NA	0% (0/3)
	Jilin	43.5°N, 125.2°E	0% (0/5)	NA	0% (0/5)
	Liaoning	41.5°N, 123.3°E	0% (0/15)	0% (0/6)	0% (0/21)
	Shandong	36.4°N, 117.0°E	0% (0/2)	NA	0% (0/2)
	Shanghai	31.1°N, 121.3°E	0% (0/3)	NA	0% (0/3)
	Yunnan	25.0°N, 102.4°E	0% (0/3)	NA	0% (0/3)
	Zhejiang	30.2°N, 120.1°E	NA	0% (0/2)	0% (0/2)
	Total		0% (0/62)	0% (0/26)	0% (0/88)
Ar. subalbatus	Chongqing	29.4°N, 106.3°E	97.6% (40/41)	79.0% (64/66)	97.2% (104/107)
	Fujian	26.1°N, 119.2°E	100% (77/77)	100% (38/38)	100% (115/115)
	Guangdong	23.1°N, 113.1°E	100% (5/5)	100% (28/28)	100% (33/33)
	Guizhou	26.4°N, 106.4°E	100% (16/16)	100% (9/9)	100% (25/25)
	Hebei	38.0°N, 114.3°E	100% (2/2)	NA	100% (2/2)
	Hubei	30.4°N, 114.2°E	97.1% (33/34)	95.0% (19/20)	96.3% (52/54)
	Hunan	28.1°N, 112.6°E	100% (2/2)	NA	100% (2/2)
	Jiangxi	28.4°N, 115.6°E	NA	100% (1/1)	100% (1/1)
	Liaoning	41.5°N, 123.3°E	63.6% (21/33)	100% (35/35)	82.4% (56/68)
	Shandong	36.4°N, 117.0°E	100% (4/4)	100% (1/1)	100% (5/5)
	Zhejiang	30.2°N, 120.1°E	90.9% (20/22)	97.1% (66/68)	95.6% (86/90)
	Total		93.2% (220/236)	98.1% (261/266)	95.8% (481/502)
Cx. pipiens	Anhui	31.5°N, 117.2°E	98.1% (53/54)	98.4% (62/63)	98.3% (115/117)
	Beijing	39.6°N, 116.2°E	95.2% (60/63)	94.4% (51/54)	94.9% (111/117)
	Chongqing	29.4°N, 106.3°E	NA	100% (3/3)	100% (3/3)
	Fujian	26.1°N, 119.2°E	0% (0/1)	100% (1/1)	50.0% (1/2)
	Gansu	36.0°N, 103.5°E	59.3% (35/59)	100% (14/14)	67.1% (49/73)
	Guangdong	23.1°N, 113.1°E	0% (0/68)	33.3% (1/3)	1.4% (1/71)
	Guangxi	22.5°N, 108.2°E	70.0% (7/10)	100% (2/2)	75.0% (9/12)
	Hebei	38.0°N, 114.3°E	16.7% (1/6)	50.0% (1/2)	25.0% (2/8)
	Heilongjiang	45.4°N, 126.4°E	42.2% (19/45)	98.4% (61/62)	74.8% (80/107)
	Henan	34.5°N, 113.4°E	100% (35/35)	100% (50/50)	100% (85/85)
	Hubei	30.4°N, 114.2°E	81.0% (17/21)	100% (21/21)	90.5% (38/42)
	Jiangsu	32.0°N, 118.5°E	92.6% (113/122)	82.4 (89/107)	88.2% (202/229)
	Jiangxi	28.4°N, 115.6°E	50.0% (10/20)	86.2% (25/29)	71.4% (35/49)
	Jilin	43.5°N, 125.2°E	72.2% (83/115)	100% (1/1)	72.4% (84/116)
	Liaoning	41.5°N, 123.3°E	100% (7/7)	100% (31/31)	100% (38/38)
	Shaanxi	34.2°N, 108.6°E	98.0% (97/99)	88.9% (8/9)	97.2% (105/108)
	Shandong	36.4°N, 117.0°E	98.2% (54/55)	100% (65/65)	99.2% (119/120)
	Shanghai	31.1°N, 121.3°E	96.3% (77/80)	100% (37/37)	97.4% (114/117)
	Shanxi	37.5°N, 112.3°E	100% (17/17)	100% (97/97)	100% (114/114)
	Sichuan	30.4°N, 104.0°E	100% (4/4)	NA	100% (4/4)
	Tianjing	39.0°N, 117.1°E	98.7% (76/77)	83.3% (5/6)	97.6% (81/83)
	Yunnan	25.0°N, 102.4°E	100% (27/27)	100% (94/94)	100% (121/121)
	Zhejiang	30.2°N, 120.1°E	85.7% (6/7)	88.9% (8/9)	87.5% (14/16)
	Total		80.4% (798/992)	95.7% (727/760)	87.0% (1525/1752)
Cx. tritaeniorhynchus	Guizhou	26.4°N, 106.4°E	5.6% (4/72)	100% (10/10)	17.1% (14/82)
Unidentified	Jiangsu	32.0°N, 118.5°E			78.2% (586/749)
Total					83.6% (2937/3514)

a NA: not available

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). Overall, 83.6% of the mosquitoes (2937/3514) were positive for Wolbachia, including 97.1% (331/341) of Ae. albopictus, 95.8% (481/502) of Ar. subalbatus, 87.0% (1525/1752) of Cx. pipiens, 17.1% (14/82) of Cx. tritaeniorhynchus and 78.2% (586/749) of unidentified mosquitoes, while both wsp and 16S rRNA genes of Wolbachia was unavailable to be detected in An. sinensis (n = 88). In detail, the positive rates of Wolbachia were 62.5~100% in Ae. albopictus, 82.4~100% in Ar. subalbatus, 1.4~100% in Cx. pipiens and 17.1% in Cx. tritaeniorhynchus sampled in different provinces of China (Table 2). The positive rate of Wolbachia in Cx. pipiens was significantly lower than that in Ae. albopictus and Ar. subalbatus (P<0.01), and significantly higher than that in Cx. tritaeniorhynchus (P<0.01).

Presence of Wolbachia in mosquitoes in Yangzhou, Jiangsu province

From September 2013 to August 2014, a total of 1025 mosquito samples were collected in a campus in Yangzhou, Jiangsu province. The monthly positive rates of Wolbachia ranged from 50.8% (66/130) to 98.4% (182/185). Compared with March 2014, the mosquito samples collected from September to December (2013), April (2014) and July (2014) harbored significant lower positive rates (P<0.01) of Wolbachia. Among the samples collected from April to August (2014), of which the species and gender information is available, significant higher infection rates were found with female samples of Cx. pipiens in April (P<0.05) (Table 3).

Table 3

The presence and positive rates of Wolbachia in mosquitoes collected in Yangzhou, Jiangsu province between September 2013 and August 2014.

Sampling time	Daily mean Temperature (°C)	Species	Positive rate
			Female	Male	Total
Sep 2013	20~27	NA	NA	NA	50.8% (66/130)
Oct 2013	13~23	NA	NA	NA	76.0% (225/296)
Nov 2013	7~16	NA	NA	NA	82.1% (78/95)
Dec 2013	0~9	NA	NA	NA	81.4% (35/43)
Mar 2014	6~15	NA	NA	NA	98.4% (182/185)
Apr 2014	11~20	Cx. pipiens	97.9% (47/48)	79.5% (35/44)	89.1% (82/92)
May 2014	17~27	Cx. pipiens	95.7% (22/23)	79.5% (23/23)	97.8% (45/46)
Jun 2014	20~28	Cx. pipiens	96.2% (25/26)	95.0% (19/20)	95.7% (44/46)
Jul 2014	23~30	Ae. albopictus	0% (0/1)	NA	0% (0/1)
		Cx. pipiens	76.0% (19/25)	60.0% (12/20)	68.9% (31/45)
		Total	73.1% (19/26)	60.0% (12/20)	67.4% (31/46)
Aug 2014	22~28	Ae. albopictus	96.2% (25/26)	95.0% (19/20)	95.7% (44/46)
Total	0~30				81.2% (832/1025)

a NA: not available

Molecular characterization of Wolbachia with wsp and MLST

The sequences of partial Wolbachia wsp gene identified in this study were submitted to GenBank 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 Wolbachia was available to be detected in four species of mosquitoes, including Ae. albopictus, Ar. subalbatus, Cx pipiens and Cx. tritaeniorhynchus. Only the positive samples of Ae. albopictus shared both supergroups, while the Wolbachia in Ar. subalbatus (supergroup A), Cx. pipiens (supergroup B) and Cx. tritaeniorhynchus (supergroup B) belonged to a single supergroup, respectively. In supergroup A, the partial sequences of wsp gene amplified from 26 Ar. subalbatus (14 females and 12 males) were identical to KJ140131 identified in Ae. albopictus in China. The remaining three strains (MW718054, MW718057 and MW718097) in supergroup A were closest to KJ140127 (Ae. albopictus, China). In supergroup B, the Wolbachia wsp sequences identified in this study were clustered into two sub-clades. The first sub-clade consisted the 81.5% (22/27) of Wolbachia sequences amplified from Ae. albopictus in this study. The second sub-clade consisted of all the Wolbachia wsp sequences obtained from Cx. pipiens (n = 52) and Cx. tritaeniorhynchus (n = 4), and two sequences obtained from Ae. albopictus (Fig 2).

Fig 2

Neighbor-joining phylogenetic tree based on Wolbachia wsp gene partial sequences (~610 bp) obtained in this study and GenBank database.

Strains identified in this study are identified with open circles (○) for supergroup A (in red) and filled circles (●) for supergroup B (in blue). The numbers at the branches show bootstrap support (1000 replicates). The bar at the bottom of the figure denotes distance.

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 #” 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).

Fig 3

Neighbor-joining phylogenetic tree based on five MLST loci (gatB, coxA, hcpA, ftsZ and fbpA) of the strains obtained in this study and GenBank database.

Strains identified in this study are identified with open circles (○) for supergroup A (in red) and filled circles (●) for supergroup B (in blue). The GenBank accession number of wsp sequences of each strain identified in this study was marked on the right side of ST code. The numbers at the branches show bootstrap support (1000 replicates). The bar at the bottom of the figure denotes distance.

Alignment of nucleotide and deduced amino acid sequences of wsp gene

When compared with the closest reference sequences (KJ140127, KJ140131, KJ140130 and KT964225) deposited in GenBank, we found that the sequences obtained in this study had a total of 15 unique mutations in partial wsp gene. Among them, three mutations (C39G, T42G and C/G394T) were identified from supergroup A and the remaining 12 mutations (T272C, T295A, C308T, G376T, A450G, G451T, G493C, T494C, T519C, C590A, C592A and T595C) were included in supergroup B. Alignment of deduced amino acid sequences demonstrated that there were four (R13G, L14V, A/S131V and K175N/S) and eight (V98D, S125I, S150V, R164P, F173L, S196R, T197D/A and V198A) mutations in supergroups A and B, respectively (Figs 4 and 5).

Fig 4

Alignment of partial-length of Wolbachia wsp gene nucleotide (~610 bp) sequences identified in this study (MW717996, MW718054, MW718057, MW718097, MW718000, MW718001, MW718008, MW718027, MW718034, MW718069, MW718076, MW718098, MW718103) together with 4 mostly identical reference sequences (KJ140127, KJ140131, KJ140130 and KT964225) obtained in GenBank database.

Strains belonging to supergroup A and B are shown in red and blue, respectively, while the reference sequences are indicated with underlines. Numbers, read from top to bottom, above the sequences are nucleotide positions on the wsp gene of MW717996. Dots indicate nucleotides identical to the reference sequences. The triangles above nucleotide position numbers indicate sites of mutations in our isolates (supergroup A: ▽ in red, and supergroup B: ▼ in blue) relative to the one or both of the reference sequences.

Fig 5

Alignment of partial-length of Wolbachia wsp gene deduced amino acid (~200 AA) sequences identified in this study (MW717996, MW718054, MW718057, MW718097, MW718000, MW718001, MW718008, MW718027, MW718034, MW718069, MW718076, MW718098, MW718103) together with 4 mostly identical reference sequences (KJ140127, KJ140131, KJ140130 and KT964225) obtained in GenBank database.

Strains belonging to supergroup A and B are shown in red and blue, respectively, while the reference sequences are indicated with underlines. Numbers, read from top to bottom, above the sequences are deduced amino acid positions on the wsp gene of MW717996. Dots indicate nucleotides identical to the reference sequences. The triangles above nucleotide position numbers indicate sites of mutations in our isolates (supergroup A: ▽ in red, and supergroup B: ▼ in blue) relative to the one or both of the reference sequences.

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-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 total infection rate of Wolbachia was 97.1% (331/341) among Ae. albopictus collected from 17 provinces. The prevalence rates were above 94.7% in most of the sampled provinces, except Guangdong and Liaoning provinces, where the infection rate was significantly lower in the females of Ae. albopictus in Liaoning province (P<0.05). To the best of our knowledge, this study revealed the high prevalence of natural Wolbachia infection in Ar. subalbatus (95.8%, 481/502), Cx. pipiens (87.0%, 1525/1752) and Cx. tritaeniorhynchus (17.1%, 14/82) in vast areas of China for the first time (Fig 1, Tables 2 and 3). A few studies have demonstrated the presence of natural Wolbachia infection in Ar. subalbatus [29-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 (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 populations of Cx. tritaeniorhynchus in Guizhou province, China, considerably lower than the positive rates of Ae. albopictus, Ar. subalbatus and Cx. pipiens (Fig 1 and 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 Wolbachia infections in field populations of several Anopheles mosquitoes (An. arabiensis, An. coluzzii, An. demeilloni, An. funestus, An. gambiae and An. moucheti) targeting 16S rRNA [36-40], but not wsp gene. In addition, Walker et al. determined a heavy infection of Wolbachia in the ovaries of An. moucheti that could induce maternal transmission [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 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-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-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 (G179 and TKDA187-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).

The sensitivity and specify of wsp-based PCR in this study.

The 10-fold dilutions were performed to give solutions containing 100,000, 10,000, 1,000, 100 and 10 gene copies per PCR reaction system. The DNA samples of Rickettsia bellii and R. felis and double-distilled water were served as negative controls.

(TIF)

Click here for additional data file.

Maximum-likelihood phylogenetic tree based on Wolbachia wsp gene partial sequences (~610 bp) obtained in this study and GenBank database.

Strains identified in this study are identified with open circles (○) for supergroup A and filled circles (●) for supergroup B. The numbers at the branches show bootstrap support (1000 replicates). The bar at the bottom of the figure denotes distance.

(TIF)

Click here for additional data file.