Wolbachia Sequence Typing in Butterflies Using Pyrosequencing

Abstract

Wolbachia is an obligate symbiotic bacteria that is ubiquitous in arthropods, with 25-70% of insect species estimated to be infected. Wolbachia species can interact with their insect hosts in a mutualistic or parasitic manner. Sequence types (ST) of Wolbachia are determined by multilocus sequence typing (MLST) of housekeeping genes. However, there are some limitations to MLST with respect to the generation of clone libraries and the Sanger sequencing method when a host is infected with multiple STs of Wolbachia. To assess the feasibility of massive parallel sequencing, also known as next-generation sequencing, we used pyrosequencing for sequence typing of Wolbachia in butterflies. We collected three species of butterflies (Eurema hecabe, Eurema laeta, and Tongeia fischeri) common to Korea and screened them for Wolbachia STs. We found that T. fischeri was infected with a single ST of Wolbachia, ST41. In contrast, E. hecabe and E. laeta were each infected with two STs of Wolbachia, ST41 and ST40. Our results clearly demonstrate that pyrosequencing-based MLST has a higher sensitivity than cloning and Sanger sequencing methods for the detection of minor alleles. Considering the high prevalence of infection with multiple Wolbachia STs, next-generation sequencing with improved analysis would assist with scaling up approaches to Wolbachia MLST.

Introduction

Wolbachia is an obligate intracellular symbiotic bacterium. It belongs to the order Rickettsiales and is known to infect numerous species of arthropods and nematodes. The bacterium is known to manipulate the reproduction of its insect hosts via cytoplasmic incompatibility, parthenogenesis, male-killing, or feminization [8,18,30,31]. Wolbachia species are classified into 11 (A–K) supergroups [22] based on multilocus sequence typing (MLST) and Wolbachia surface protein (WSP) typing [2]. The standard MLST system for Wolbachia determines a sequence type (ST) based on the combination of alleles for five conserved housekeeping genes (ftsZ, gatB, coxA, hcpA, and fbpA). The WSP typing system determines STs using four hypervariable regions (HVRs) of WSP: HVR1 (aa 52 to 84), HVR2 (85 to 134), HVR3 (135 to 185), and HVR4 (186 to 222) [1].

A host can be infected with a single strain of Wolbachia, or multiple strains, as has been observed in a wide range of insects [13, 20, 21]. Classical MLST can be used for the universal characterization of Wolbachia, but its application is currently limited to hosts infected with a single strain of Wolbachia [2]. To assess the diversity of Wolbachia in insects infected with multiple strains, sequencing of the wsp gene is generally conducted [32]. Amplification of ST-specific regions is also used, but this method has limited application because the design of specific primers is only possible after a strain’s sequence is known. Improved methods for determining STs in insects infected with multiple strains would contribute to understanding Wolbachia diversity.

Butterflies are important plant pollinators; therefore, they are used as model systems in a variety of research fields [4]. Wolbachia has been detected in a wide range of butterfly species [24], with 24 Wolbachia STs recorded in butterflies (http://pubmlst.org/Wolbachia/) [25]. There are 268 known butterfly species in Korea, and five of these are yellow butterfly species (Gonepteryx rhamni, Gonepteryx mahaguru, Colias erate, Eurema hecabe and Eurema laeta) [26]. The yellow butterfly E. hecabe is widely distributed in East Asia and is known to be infected with Wolbachia [8, 15-17]. Japanese populations of E. hecabe are infected with two distinct strains of Wolbachia, wHec1 and wHec2 [17]. Tongeia fischeri is a small butterfly found in Eastern Europe and Northeastern Asia [12];T. fischeri butterflies in west Siberia, Russia are all infected with a single Wolbachia ST300 and personal communication [11]. However, Wolbachia infection in butterflies has not been widely investigated in Korea. Considering the ecological role and abundance of butterflies, the infection status of Wolbachia in butterflies should be investigated. The accurate sequence typing result would provide the genetic framework for tracing the movement of Wolbachia within insect communities, the regional distribution of butterflies, and phenotypic effects on butterflies [2, 7].

Wolbachia strains are impossible to isolate from the host because they are obligate symbionts. Thus, the MLST for Wolbachia has depended on metagenomic cloning and sequencing methods. The public MLST database for Wolbachia is also based on cloning methods with Wolbachia strains found in diverse singly infected host species. Because the metagenome-based sequencing is the standard protocol for Wolbachia MLST, we thought next-generation sequencing (NGS)-based MLST would have merit as a massive screening technique for this particular taxon. Thus, we developed and applied pyrosequencing-based MLST and applied the method to three butterfly species (E. hecabe, E. laeta, and T. fischeri). Our results suggest that pyrosequencing-based MLST could be used for the largescale screening of multiple Wolbachia STs.

 

Materials and Methods

Sample Collection and DNA Extraction

Three butterfly species (E. hecabe, E. laeta, and T. fischeri) were collected from Goheung-gun, Jeonnam, and Han River (Seoul, Korea), between October 2013 and September 2014. One adult butterfly for each species was collected and transported to our laboratory. The body of each butterfly was homogenized, and DNA was extracted from 100 μl of the homogenized sample using a QIAamp DNA Minikit (Qiagen, Germany).

DNA Amplification and Sanger Sequencing

Polymerase chain reaction (PCR) amplification of Wolbachia wsp and conserved housekeeping genes (gatB, coxA, fbpA, ftsZ, and hcpA) was performed using specific primers and protocols described by Baldo et al. [2]. GoTaq Colorless Master Mix (Promega, USA) and PCR Thermal Cycler Dice TP600 (Takara, Japan) were used for the PCR. PCR products were purified with the QIAquick PCR Purification Kit (Qiagen), and were ligated into the pTOP TA V2 vector using the TOPcloner TA core kit (Enzynomics, Korea). Plasmids were transformed into competent Escherichia coli cells; colonies containing recombinant plasmids (10 white colonies) were picked and grown overnight. Recombinant plasmid DNA was extracted from overnight E. coli cultures using a LaboPass Plasmid Miniprep Kit (Cosmo Genetech, Korea). The sizes of inserts were determined by restriction endonuclease digestion of the plasmid DNA using EcoRI (Takara, Japan). Plasmids with inserts of the correct size were subjected to Sanger sequencing by Macrogen Inc. (Korea).

Pyrosequencing

For pyrosequencing, Wolbachia wsp and the five housekeeping genes were amplified using barcoded fusion primers. The forward fusion primer was composed of 454-adaptor, key sequence, linker, and target gene-specific sequences (5’-CCTATCCCCTGTGTGCCTTGGCAGTC-TCAG-AC-target sequence-3’). The reverse fusion primer was composed of 454-adaptor, key sequence, barcode, linker, and target gene-specific sequences (5’-CCATCTCATCCCTGCGTGTCTCCGAC-TCAG-barcode-AC-target sequence-3’), with a unique barcode designed for each subject. The full list of oligonucleotide primer sequences used in the current study is summarized in Table S1. Each PCR was performed in a 50 μl volume and contained 1.25 μl of each primer, 25 μl of Taq DNA polymerase buffer, 21.5 μl of distilled water, and 1 μl of template DNA. Reactions were incubated at 94℃ for 5 min, followed by 37 amplification cycles (94℃ for 30 sec, the optimal annealing temperature for 45 sec, 72℃ for 90 sec), and then a final elongation step at 70℃ for 10 min, with the temperature then held at 4℃. Optimal PCR annealing temperatures were 53℃ for hcpA, 54℃ for gatB and ftsZ, and 55℃ for coxA, fbpA, and wsp. To reduce amplification bias, three independent PCRs were performed and pooled. Amplicons were confirmed by agarose gel electrophoresis and purified using a QIAquick Gel Extraction kit (Qiagen). The DNA concentration was measured using a Picodrop (Bioneer), and equal quantities were mixed to create amplicon pools. Pyrosequencing was performed at ChunLab Inc. (Korea) using a Roche 454 GS Junior platform.

Sequence Analysis and Identification

All processing of pyrosequencing data was performed using Mothur ver. 1.29.2 [27]. The sequencing reads from the different samples were separated by their unique barcodes, and the sequences of the barcode, linker, and PCR primers were trimmed from both ends of the sequencing reads. For gatB, coxA, and fbpA sequences, reads that were <350 bp were filtered out, and for hcpA, ftsZ, and wsp, reads that were <400 bp were filtered out using the screen following command: minlength = 350, maxhomop = 8, maxambig = 0. The pre-filtered reads were aligned using reference sequences obtained from GenBank and the Wolbachia MLST database (http://pubmlst.org/Wolbachia/). During clustering, sequencing errors were ignored by allowing for mismatch errors of up to 2 bp, given the substitution sequence error rate of pyrosequencing (0.5%) [6,10,14], using the pre.cluster command diffs = 2. Chimeras were identified and removed using UCHIME [5] (chimera.uchime). For the five housekeeping genes used in MLST, the top five largest clusters were chosen as representative sequences and subjected to identification. For the wsp gene, reads with an exact match to any HVR (HVR1, 2, 3, or 4) were chosen for identification. The allele numbers of representative sequences were identified according to the Wolbachia MLST database.

Phylogenetic Analysis

To align sequences and construct phylogenetic trees, MEGA5.2 [29] and MrBayes 3.2.3 [9] were used. The wsp gene sequences generated in this study were aligned with homologous sequences deposited in the Wolbachia MLST database. The consensus sequence comprised 495 bp for wsp gene fragments. A maximumlikelihood tree was inferred using the Jones-Yalor-Thorton substitution model and evaluated by 1,000 bootstrap replicates. A Bayesian tree was generated through Markov Chain Monte Carlo methods. The standard deviation was below 0.003 after 5,000,000 generations, and a 10% “burn in” of total samples was conducted.

Sequence Data Availability Statement

Raw data files of pyrosequencing are available in the NCBI Sequence Read Archive (Accession No. SRP 058686). Sequences obtained from Sanger sequencing were deposited in the GenBank database with the accession numbers KP763428–KP763455 and KR006333–KR006338.

 

Results

Sanger Sequencing-Based MLST

Wolbachia infection in the three butterfly species was confirmed by PCR amplification of the wsp gene (514–566 bp). Cloning and Sanger sequencing of the amplicons revealed that T. fischeri was infected with a single Wolbachia ST. In contrast, the yellow butterflies were infected with two Wolbachia STs, with two different wsp sequences seen among 38 clones. The five housekeeping genes were amplified from the same template DNAs, with at least 10 clones for each gene. The resultant gatB, hcpA, and fbpA sequences revealed that both yellow butterfly species were infected with at least two different Wolbachia STs by showing two STs for each gene. However, sequence variability was not consistently observed; for coxA, only one clone was successfully sequenced, thus resulting in a single allele type. For ftsZ, all clone sequences (10 clones for E. hecabe; 11 clones for E. laeta) contained a single allele (ftsZ 36).

Pyrosequencing Reads

The total number of pyrosequencing reads for each sample ranged from 59 to 3,519 (Table 1). After quality filtering, alignment, and the removal of chimeric sequences, the numbers of sequencing reads were reduced by 23%. The pre-filtered reads were grouped into 12–518 clusters, with a cluster cut-off value of a 2 nt difference. The size of the resulting clusters varied from 2–2,063 reads, with the top five largest clusters containing 59% of the pre-filtered reads. In most cases (genes and organisms) analyzed in this study, the increment curve of clusters reached a plateau before the fifth cluster (Fig. 1). In all cases, the allele identification was complete by the top five clusters, and no more allele was found even though the sixth and the rest of the clusters were identified fully. Thus, we considered the other clusters as error sequence-based reads, and the reads of the top five largest clusters were chosen as representative sequences and subjected to further identification.

Table 1.aA cluster cut-off value of a 2 nt difference. bFor gatB, coxA, hcpA, ftsZ, and fbpA, the top five largest clusters were chosen as the representative sequences and subjected to identification. For wsp, reads with an exact match to any hypervariable region (HVR1, 2, 3, or 4) were chosen for identification.

Fig. 1.Plot of the number of clusters as a function of the number of sequencing reads of gatB. The curves grow rapidly at first, but then plateau as only the rarest clusters with a low number of sequencing reads remain.

Wolbachia STs Determined by MLST Allelic Profiles

Allele sequences obtained by Sanger sequencing and pyrosequencing were identified using the Wolbachia MLST database. A summary of the allele IDs for the five housekeeping genes is provided in Table 2. A single ST in T. fischeri was evident from the Sanger sequencing and pyrosequencing results. According to our pyrosequencing results, the top five clusters of the tested genes were identified as an identical allele set and determined as ST41. In contrast, the two yellow butterfly species exhibited two alleles for gatB (39, 38), coxA (14, 38), fbpA (4, 42), ftsZ (36, 35), and hcpA (40, 29) (Table 2). Theoretically, this gene profile could result in 25 combinations, and consequently 32 different STs. However, according to the Wolbachia MLST database, only ST40 and ST41 were reported among the 32 theoretical allele profiles. Eurema hecabe and E. laeta were infected by the same Wolbachia population, and were equally infected with the mixed population of Wolbachia ST41 and ST40.

Table 2.aAlleles detected by pyrosequencing only.

The allelic profiles generated from the pyrosequencing method revealed differences in their composition or proportion when compared with the Sanger sequencing ones. CoxA allele 14 and ftsZ allele 35 were not detected by Sanger sequencing, but were evident through pyrosequencing. The proportion of alleles in a population also differed, depending on the sequencing method employed. As an example, the proportion of gatB in E. hecabe was 53% according to pyrosequencing, but only 20% according to Sanger sequencing.

Wolbachia Phylogeny Based on wsp

Based on our pyrosequencing results, 85% of clusters matched with at least one HVR. The proportion of reads containing HVR1 (0.7%) or HVR2 (8.9%) was lower than that for reads containing HVR3 (59.0%) or HVR4 (73.3%). Since sequencing of the wsp gene was started from HRV4 to HRV1, these differences were likely due to limited read lengths and higher error rates at the 3 end of the sequence. The length of wsp required for haplotype identification is 483 bp; however, the majority of sequences generated by pyrosequencing were shorter than this. For the wsp gene, reads that matched any HVR were chosen for identification. The wsp allele of T. fischeri was identified as wsp-10, which corresponded with the Sanger sequencing results. In the yellow butterfly species, both sequencing methods revealed two wsp alleles, wsp-64 and wsp-10, with unique HVR peptide profiles (10-8-10-8 and 35-35-38-44, respectively) (Table 3). The dominant wsp allele was wsp-10 with complete HVR peptide profiles; wsp-64 contained a unique partial HVR peptide profile.

Table 3.Wsp allelic profiles of Wolbachia strains from yellow butterflies.

A phylogenetic tree based on wsp sequences revealed that the Wolbachia strains detected in this study belonged to Wolbachia supergroup B (Fig. 2). The detected wsp alleles were positioned at different phylogenetic branches within supergroup B, and demonstrated a remarkable genetic distance. The wsp-64 sequences from the two yellow butterflies clustered with those for the orange butterfly (Ariadne merione), and the small brown plant hopper (Laodephax striatellus), and were distinctly separate from wsp-10. The two wsp-10 sequences from yellow butterflies clustered with that for Indian E. hecabe. The wsp-10 sequence of T. fischeri differed slightly from that in the yellow butterflies and belonged to a sister group.

Fig. 2.Bayesian phylogenetic tree based on wsp loci (495 bp). Numbers in parentheses indicate wsp alleles. The numbers above the nodes are posterior probabilities from 5,000,000 generations. The numbers below the nodes are percentages of bootstrap support, from 1,000 resampled datasets of maximum-likelihood and/or neighborjoining tree-inferring methods. The scale bar indicates a 0.02 nt change per position.

 

Discussion

In the current study, we applied NGS technology to conduct MLST and investigate cases of Wolbachia infection involving multiple strains. Pyrosequencing read lengths fell within the size range of genes used in MLST (450–500 bp) [19]; therefore, a large number of trimmed pyrosequencing reads could be successfully assigned to allele number. According to the results of this study, pyrosequencing-based MLST is superior to Sanger-based methods with respect to sensitivity of detection. As an example, minor allele sequences, such as coxA allele 14 and ftsZ allele 35, were only detected by pyrosequencing in this study. However, there were several limitations with the pyrosequencing-based MLST. First, the time and cost effect of NGS-MLST is inferior to Sanger-based MLST. Second, the short sequence length of pyrosequencing reads hindered the complete identification of wsp hypervariable regions. Third, current NGS technologies, such as pyrosequencing, exhibit high error rates (0.49–2.8%) [7] compared with Sanger sequencing (0.0001–1%) [7]. The high error rates could result in over-calculation of STs, especially in cases where multiple strains are present. However, NGS technologies are rapidly evolving, and error rates are likely to decrease while read lengths increase. Improvements in experimental methods and bioinformatics analyses regarding next-generation sequencing-based MLST will allow us to conduct large-scale Wolbachia MLST in the near future.

The coxA of E. hecabe was hard to amplify. We tried various PCR conditions and various cloning vectors, but could get only one Sanger sequence. In line with Sanger sequencing results, only 59 reads of the coxA gene were obtained from pyrosequencing of E. hecabe. The coxA gene of E. hecabe probably had mutation in the primer region and produced a low yield of DNA amplification. This is a usual case that a part of MLST genes is not amplified or sequenced through the traditional Wolbachia MLST scheme [1, 3, 23]. Although the number of pyrosequencing reads were as small as 59 for E. hecabe, it was enough to determine STs in E. hecabe.

Wolbachia infection in E. hecabe and E. laeta was previously reported in India and Japan [24, 28]. Japanese E. hecabe and E. laeta are infected with both ST40 and ST41 [28], the same STs that we identified in both Korean yellow butterflies. In contrast, Indian E. hecabe (ST41) and E. laeta (ST149) are infected with a single Wolbachia ST [24]. This finding implies that Northeastern Asia might have homogeneous yellow butterfly populations. Further investigation of larger geographical regions is required to confirm this. Although two alleles were identified by MLST (ST41 and ST40) and wsp sequence typing (wsp-10 and wsp-64), it is difficult to conclude that only two strains of Wolbachia were present in the yellow butterflies. Certain Wolbachia strains are known to share the same wsp allele despite harboring different STs. For instance, Wolbachia strains sharing wsp-10 but harboring two distinct STs (ST41 and ST157) were reported in four E. hecabe butterflies from India [24]. The two Korean yellow butterflies we investigated were infected with at least two strains of Wolbachia.

We found that the three butterfly species we investigated were infected by a multitude of Wolbachia strains belonging to supergroup B. This is the first report of Wolbachia infection in butterflies from Korea, and would give basic knowledge for researchers on butterflies and insects. However, the number of organisms and species we surveyed were limited. Further investigation of Wolbachia infection in a larger number of butterfly species across many geographical areas will be necessary to obtain an accurate indication of the distribution and diversity of Wolbachia species among Korea insect populations.