Analysis of a native whitefly transcriptome and its sequence divergence with two invasive whitefly species

Illumina sequencing and assembly of Asia II 3 transcriptome

Functional annotation of Asia II 3 transcriptome

For functional annotation, distinct sequences were searched against the non-redundant (nr) NCBI nucleotide database and a total of 16,596 genes returned an above cut-off BLAST result representing about 31.6% of all distinct sequences (Additional file 2). This proportion is similar to the 20% to 40% of annotated sequences from a traditional Sanger sequenced EST library [32]. To determine the possible functions of assembled Asia II 3 genes, Gene Ontology (GO) assignments were used to classify the distinct sequences. Based on sequence homology, 4,819 sequences could be categorized into functional groups under the “Molecular function”, “Biological process” and “Cellular component” divisions (Additional file 3). The functions of genes cover various biological processes and genes participate in “Cellular process” and “Metabolic process” are the most highly represented. Next, we compared the GO classification of the Asia II 3 transcriptome with that of MEAM1 and MED transcriptomes, respectively [29, 30] and found that the distributions of gene functions from these three species are similar (Additional file 3). These results suggest: i) the functions of genes from Asia II 3, MEAM1 and MED are highly conserved; ii) there is no bias in the construction of the libraries from these B. tabaci species.

Analysis of Asia II 3 gene expression

The level of Asia II 3 gene expression was analyzed based on the number of Reads Per Kilobase per Million mapped reads (RPKM); and a list of all the genes and expression levels are shown in Additional file 4. To our knowledge, this is the first global analysis of gene expression level in B. tabaci. Twenty annotated genes with very high expression levels (RPKM> 1,250) were found. Many of the genes with significant expression levels were involved in cell structure (e.g. actin and tubulin), ribosome (e.g. 60S ribosomal protein, 40S ribosomal protein) and energy metabolism (e.g. ATP synthase and glyceraldehyde-3-phosphate dehydrogenase) (Table 2). This finding is not surprising, as these genes are essential for the survival of an organism. Next, we grouped Asia II 3 genes into three categories based on their levels of expression. Roughly, 14% of the genes were highly expressed (RPKM> 50), 25% of the transcripts were in medium level (50>RPKM>20) and 61% of the genes had a RPKM<20 (Figure 1A). Conversely, when calculating the number of reads for each gene, it was found that the small fraction of highly expressed genes constituted approximately 61% of the sequenced reads, whereas only 13% of the reads were from the genes with low expression levels (Figure 1B).

Gene ID	Number of readsa	RPKMb	Swissprot annotation
BT_ZHJ1_ZJU_Unigene4448	13476	3433.5	Actin-5, muscle-specific
BT_ZHJ1_ZJU_Unigene5935	8392	3212.7	Tubulin alpha-1 chain
BT_ZHJ1_ZJU_Unigene44594	6567	2854.2	60S ribosomal protein L18a
BT_ZHJ1_ZJU_Unigene40759	4952	2613.5	ADP,ATP carrier protein 2
BT_ZHJ1_ZJU_Unigene33425	3191	2213.8	Elongation factor 1-alpha 2
BT_ZHJ1_ZJU_Unigene36010	3122	1946.6	Glyceraldehyde-3-phosphate dehydrogenase
BT_ZHJ1_ZJU_Unigene40277	3388	1820.6	ATP synthase subunit beta, mitochondrial
BT_ZHJ1_ZJU_Unigene33344	2299	1595.0	Troponin T, skeletal muscle
BT_ZHJ1_ZJU_Unigene42468	3239	1584.9	Elongation factor 1-alpha
BT_ZHJ1_ZJU_Unigene38168	2689	1552.2	40S ribosomal protein S12
BT_ZHJ1_ZJU_Unigene44051	3336	1498.4	Vitellogenin-A1
BT_ZHJ1_ZJU_Unigene43506	3232	1497.2	40S ribosomal protein S11
BT_ZHJ1_ZJU_Unigene49876	5853	1473.5	60S ribosomal protein L15
BT_ZHJ1_ZJU_Unigene43457	3079	1430.8	40S ribosomal protein S20
BT_ZHJ1_ZJU_Unigene46015	3456	1358.3	60S ribosomal protein L6
BT_ZHJ1_ZJU_Unigene37661	2284	1344.7	Ferritin, middle subunit
BT_ZHJ1_ZJU_Unigene38243	2286	1319.6	Paramyosin, short form
BT_ZHJ1_ZJU_Unigene45492	3074	1258.3	ATP synthase subunit alpha, mitochondrial
BT_ZHJ1_ZJU_Unigene44482	2871	1255.2	Vitellogenin
BT_ZHJ1_ZJU_Unigene34665	1906	1251.8	40S ribosomal protein S2

^aThe total number of reads mapped to each gene.

^bGene expression level was determined by calculating the number of reads for each gene and then normalizing to RPKM (reads per kilobase per million mapped reads).

Identification and analysis of B. tabaci protein families

To compare the transcriptomes of Asia II 3, MEAM1 and MED, protein families for the three species were built. A total of 10,217 protein families were identified based on sequence similarities (Figure 2). The number of protein families for MED transcriptome was 8,405, which is much higher than that of Asia II 3 (4,938) and MEAM1 (6,174). In this comparison, both diversification and conservation of the B. tabaci transcriptomes were obvious at the protein family level. A total of 3,023 protein families were conserved among the transcriptomes of Asia II 3, MEAM1 and MED which represent the core cellular and physiological proteins (hereafter “Core protein families”) common to the three species. The total number of protein families found in only one species was 3,940 (2,756+582+602), which might be responsible for the differences and unique features of each species (Figure 2). However, as these transcriptomes are incomplete, the number of species-specific protein families is likely to change when more genes are sequenced. Among these specific protein families, the majority was from MED (2,756). This is probably due to the higher sequencing amount of MED (3G) compared to Asia II 3 (1G) and MEAM1 (1G) (Additional file 1) [29, 30]. To reveal the common characteristics of the whitefly transcriptomes, a hypergeometric test was implemented to identify enriched GO terms in the “Core protein families”. A total of 18 protein families were significantly enriched (Table 3). Interestingly, 3 of the 18 enriched GO terms were related to nucleotide binding. This phenomenon is consistent with the previous finding in the analysis of nematode transcriptomes [33] and demonstrates the importance and conservation of nucleotide binding proteins in different species. The other major GO terms enriched in the “Core protein families” were related to amino acid transporter, protein folding, proteolysis and peptidase, suggesting the critical roles of protein transportation and metabolism among the three B. tabaci species (Table 3).

Identification of the orthologous genes between Asia II 3 and MEAM1, Asia II 3 and MED

Possible orthologous genes between the transcriptomes of Asia II 3 and MEAM1were identified according to the diagram in Figure 3 (left). Between Asia II 3 and MEAM1, a total of 20,929 pairs of possible orthologs were identified. Among these sequence pairs, 3,308 pairs of sequences could be unambiguously mapped to the same protein in Swissprot database, suggesting strongly that they are orthologous genes (Figure 3). Based on the coding information of Swissprot hits, the untranslated region (UTR) of each gene pair was predicted (210 pairs of 5'UTR and 337 pairs of 3'UTR) and 2,966 pairs of orthologous coding sequences (CDS) were obtained (Figure 3 and Additional file 5). The features of the 2,966 orthologous genes are listed in Table 4. Similarly, 22,415 pairs of potential orthologs were identified between Asia II 3 and MED (Figure 3, right). After the Swissprot annotation, 102 pairs of 5’UTRs, 102 pairs of 3’ UTRs and 2,529 pairs of CDS were identified between the transcriptomes of Asia II 3 and MED (Table 4). The average divergence for the 2,966 orthologous genes between Asia II 3 and MEAM1 is 1.73%. This difference is comparable to the average divergence between Asia II 3 and MED (1.84%) and much higher than the divergence between MEAM1 and MED (0.83%) (Table 4) [30]. Genome sequences are very useful for constructing phylogenetic trees with high resolution and accuracy. Using all the orthologous genes among the three species, a neighbor-joining tree was reconstructed (Figure 4). The phylogenetic distance between MEAM1 and MED is less than that of the invasive whiteflies to the native Asia II 3. This is not surprising as both MEAM1 and MED originated from the Middle East Asia Minor and Mediterranean Basin regions whereas Asia II 3 is a native species in China [14, 34].

The sequence divergence between Asia II 3 and MEAM1

The sequence divergence between Asia II 3 and MED

For Asia II 3 and MED, the average differences between the 5’UTRs, CDS and the 3’UTRs are 3.57%, 1.84% and 3.13%, respectively (Additional file 6). The sequence divergences are almost the same between Asia II 3 and MEAM1 and nearly twice of the divergence between MEAM1 and MED which is 1.66%, 0.83% and 1.43% (Figure 5) [30]. Likewise, in the 4d and nd regions, the same phenomenon was observed (Figure 5 and Additional file 6). It confirmed the previous analysis that the phylogenetic distance between MEAM1 and MED is shorter than that of the invasive whiteflies to the native Asia II 3 (Figure 4) [2, 16]. These results indicate that despite high-sequence identity in their coding sequences, the Asia II 3, MEAM1 and MED B. tabaci species have diverged substantially between their transcriptomes. The levels of sequence divergence provide further supports to the previous proposition that Asia II 3, MEAM1 and MED are different species [25, 26, 37].

Analysis of synonymous and non-synonymous sites

The nonsynonymous (Ka) and synonymous (Ks) substitution rates have been widely used to measure the intensity of gene evolution. To identify genes undergoing purifying or positive selections, we estimated Ka and Ks rates of orthologous gene pairs. Among the 2,966 pairs of CDS between Asia II 3 and MEAM1, both a Ka and a Ks rate could be calculated for 1,373 orthologs (Additional file 7). The 1,373 sequence pairs had mean values of Ka, Ks and Ka/Ks of 0.0094, 0.0623 and 0.198. These ratios are similar to the average Ka/Ks ratio of MEAM1 & MED (0.225), Ka/Ks ratio of rat & mouse (0.19) and human & chimpanzee (0.22) [30, 38, 39]. The distribution of the Ka/Ks ratio showed that the majority of genes (98.2%, 1348/1373) have Ka/Ks ratios less than 1, indicating the strong purifying selection for these genes. In addition, Fisher's exact test reports that nearly 56% of the genes are statistically significant (P<0.01) (Additional file 7). In this set of data, 25 orthologous gene pairs had a Ka/Ks value > 1 suggesting that strong positive selection acts on these genes. Among the sequences with Ka/Ks values > 1, a number of genes are involved in protein metabolism, such as peptide deformylase, cathepsin, cysteine proteinase and metalloendopeptidase, suggesting this process is critical for the differentiation of the B. tabaci species complex. Between Asia II 3 and MED, a total of 1,221 orthologs had Ka and Ks rates and the mean values of Ka, Ks and Ka/Ks were 0.0092, 0.0729 and 0.201. The list of Asia II 3 and MED genes with Ka/Ks and functional annotation is presented in Additional file 8. The distribution of Ka/Ks between Asia II 3 and MED is similar to that of Asia II 3 and MEAM1, in which 7.9% of the genes had a high Ka/Ks>0.5, 44.8% of the genes were highly conserved (Ka/Ks<0.1) and 25 genes had a Ka/Ks >1. Even though a number of genes under positive selection have been identified in our study, the simple Ka/Ks calculations are quite conservative and may fail to detect positive selection even when it exists [40]. Further studies using more sophisticated site and branch specific models are needed for estimating Ka/Ks.

Sequences with very weak amino-acid similarity

Analysis of SNP

To further understand the mechanism of divergence, we analyzed the potential SNP sites in CDS of the orthologous gene pairs between Asia II 3 and MEAM1. For Asia II 3, a total of 138 SNPs were identified within the 1433.63 kb aligned regions, about 1 every ten thousands bp. The complete list of SNPs with annotation can be found in Additional file 10. Of the 138 SNPs, 96 (69.6%) were synonymous and 42 (30.4%) were non-synonymous. This percentage of SNPs is much lower than those obtained in previous analyses in other insects [41, 42]. The possible reason is that our B. tabaci populations were established from a pair of B. tabaci and the samples were collected within five generations. For MEAM1, a slightly higher number of SNPs (248, about 1.7 every ten thousands bp) were found in the orthologous gene pairs with 196 synonymous and 52 non-synonymous (Additional file 10). Compared with the average divergence between Asia II 3 and MEAM1 in CDS (1.83%), the percentage of SNPs is more than 100 times lower. Next, the potential SNP sites in the CDS of orthologous gene pairs between Asia II 3 and MED were analyzed. Our results showed that a total of 56 SNPs in Asia II 3 and 627 SNPs in MED were identified within the 1,071 kb aligned regions (Additional file 11). The large variation between the numbers of SNPs in Asia II 3 and MED is probably due to the difference in the sequencing among of Asia II 3 (1G) and MED (3G) (Additional file 1). Some of the SNPs in Asia II 3 might have been filtered out because only SNP sites with the minimum read depth of 10 reads were selected. Thses results are consistent with previous findings, in which overall number of SNPs decreases at lower coverage levels [43]. While the SNPs we have identified here are suitable for future research, more rigorous statistical tests are required to confirm the current results as well as to detect specific codons undergoing adaptive changes. In addition, further studies of the SNPs on population samples are warranted as our data were generated from inbred lab colonies.

Insect rearing

Stock culture of the Asia II 3 (mitochondrial cytochrome oxidase I gene GenBank accession no: AJ867556) was maintained on cotton, Gossypium hirsutum (Malvaceae) cv. Zhe-Mian 1793 in a climate chamber (see conditions below). The purity of the culture was monitored using the random amplified polymorphic DNA-PCR technique with the primer H16 (5’-TCTCAGCTGG-3’) [44]. For sample preparation, a pair of virgin adults of B. tabaci Asia II 3 were released onto a cotton plant to oviposit and develop for five generations at 27 ± 1°C, a photoperiod of 14 h light:10 h darkness and 70 ± 10% relative humidity [45]. The same protocols were used to raise the MEAM1 and MED whiteflies for sample collection and subsequent transcriptome data generation.

Sample preparation and RNA isolation

In order to get an overall picture of the Asia II 3 whitefly transcriptome, we collected the samples from four different developing stages: 1) egg & nymph (the eggs are extremely small, therefore a mixture of eggs and first to third instar nymphs were collected as one sample); 2) pupa; 3) female adult and 4) male adult. To ensure that the whitefly adults are in the same developmental stage, only newly emerged adults were collected. Previously, samples from MEAM1 and MED have been collected using the same strategy [29, 30]. Total RNA was isolated from the four samples using SV total RNA isolation system (Promega) according to the manufacturer’s protocol [46]. RNA integrity was confirmed using the 2100 Bioanalyzer (Agilent Technologies) with a minimum RNA integrated number value of 8. Then, equal amount of RNA from egg & nymph, pupa, female adult and male adult were mixed, and mRNA was purified from the mixture using oligo (dT) magnetic beads.

Library preparation and Illumina sequencing

For transcriptome sequencing, a 200 bp cDNA library was prepared using Illumina’s kit as previously described [29]. The library was not normalized because we intend to use the resulting sequence reads to analyze the level of gene expression. The cDNA library was sequenced for both ends on the GAII Illumina sequencing platform (a single lane) at The Beijing Genome Institute (Shenzhen, China). The total sequencing amount was 1G. The raw reads were filtered by removing adaptor sequences, empty reads and low quality sequences (reads with unknown sequences 'N') [30]. Next, the reads were randomly clipped into 21 bp K-mers for assembly using SOAPdenovo software because the 21-mer provided the best result for transcriptome assembly. Small K-mers made the graph very complex; while large K-mers have poor overlap in low sequencing depth regions [31]. For alleles, the nucleotides with the highest frequency were selected. The resultant contigs were joined into scaffolds based on the mate pairs information and the scaffolds were clustered using TGI Clustering tools [47]. Assembled genes were used for subsequent analyses and are referred to as “distinct sequences”.

Data deposit

The data sets of Illumina sequencing are available at the NCBI Short Read Archive (SRA) database with the accession number: SRR062575. The assembled sequences were deposited in the NCBI Transcriptome Shotgun Assembly (TSA) database under the accession number of HP777244 to HP823074 and can be searched using the GeneID listed in Additional file 2, Supporting information.

Functional annotation and gene expression analysis

Distinct sequences were annotated by Blast search against the NCBI nr database with a cut-off E-value of 10^-5. GO annotation was analyzed by Blast2GO software [48]. The GO terms were retrieved from Blast hits with an E-value threshold of 10^-5. Comparisons of the distribution of GO terms among the Asia II, MEAM1 and MED transcriptomes were done using the Web Gene Ontology Annotation Plot (WEGO) [49]. Pathway annotation was performed using Blastall software against the KEGG database. Based on the number of reads for a gene, gene expression levels can be estimated from Illumina sequencing with great accuracy [50, 51]. To estimate the level of gene expression, the number of reads mapped to each distinct sequence was extracted. Since read mapping is sensitive to the size of the target reference sequence and sequencing amount, we adjusted the raw read count by the total number of reads mapped and the length of the gene by calculating Reads Per Kilobase per Million mapped reads (RPKM) [50].

Analysis of protein families

To reveal the functional differences among Asia II 3, MEAM1 and MED transcriptomes, we analyzed their protein families. Using Blastx (E-value <10^-5), the translated region of each gene was identified by aligning the sequence to the Swissprot database. The longest translated protein sequence for each gene was then extracted and sequences less than 200 bp were removed. To identify protein families among the three transcriptomes, an all-against-all Blastp was performed for all the translated genes from the three transcriptomes. Blastp results were further analyzed by a Markov Cluster Algorithm (MCL) with an inflation factor of 1.6. The protein families belonging to all the three transcriptomes were referred to as “Core protein families”. Based on the GO annotation of Asia II 3, MEAM1 and MED transcriptomes, we calculated the total number of genes under each GO term in the “Core protein families” and the three transcriptomes. For each GO term, its enrichment in the “Core protein families” was measured using the hypergeometric test with an cut-off p value of 10^-5[33].

Identification of orthologous genes and prediction of coding and untranslated regions

The orthologous genes between Asia II 3 and MEAM1 and those between Asia II 3 and MED were identified respectively according to the previous description using MegaBLAST [30]. Briefly, pairs of sequences that were reciprocally a best hit and with a minimum length of 200 bp were retained as putative orthologs. To remove potential paralogs, only pairs of sequences unambiguously mapped to the same protein in Swissprot database with an E-value <1×10^-5 were selected. The CDS of the orthologous genes were determined by BLASTx the Swissprot database with an E value<1×10^-5. The start codon was determined by examination of the in-frame ATG codon of the aligned reference protein. The stop codon position was determined by examination of in-frame TAA, TAG and TGA motifs present within 30 bp of the stop codon of the reference protein. The 5’UTR and 3’UTR regions were defined based on the position of start codon, stop codon and predicted CDS. To prevent false positive results, only UTR pairs with an E-value<1×10^-30 were selected for further analyses. CDS containing unexpected stop codon(s) and shorter than 150 bp were removed.

Sequence divergence analyses and estimation of substitution rates

The 5’UTR, CDS and 3’UTR regions were separately extracted from each pair of orthologs. The CDS and UTR regions were aligned separately to each other with a MegaBlast algorithm and checked manually for errors. Only the homologous regions of each gene pair were extracted for sequence comparison. Sequence divergence between the homologous regions of each gene pair was calculated by dividing the number of substitutions by the number of base pairs compared. The average divergence between transcriptomes was determined by dividing the total number of substitutions by the total number of base pairs compared. The sequence divergence at nondegenerate (nd), fourfold degenerate (4d), CpG and non-CpG regions was determined respectively according to the previous descriptions [39]. The ratio of transition over transversion (ts/tv) was determined for the 5’UTR, CDS and 3’UTR as well. Using the KaKs Calculator, we also estimated the substitution rates for non-synonymous sites (Ka) and synonymous (Ks) sites (YN method) [52, 53].

Phylogeny of Asia II 3, MEAM1 and MED

The orthologous genes among Asia II 3, MEAM1 and MED were selected for sequence alignment using MUSCLE [54] and the aligned fragments were extracted for phylogenetic analysis using MEGA 5 [55]. The evolutionary history was inferred using the neighbor-joining method [56]. All positions containing gaps and missing data were eliminated and the within population polymorphisms were not included for divergence estimation. The analysis involved a total of 686,101 positions in the final dataset. The evolutionary distances were computed using the Maximum Composite Likelihood method [57] and are in the units of the number of base substitutions per site.

SNP analysis

To reveal the mechanism of divergence between Asia II 3 and MEAM1, we analyzed the potential SNP sites in CDS of the orthologous gene pairs between the invasive and indigenous whitefly species. The orthologous gene pairs of Asia II 3 & MEAM1 and Asia II 3 & MED were subjected to SNP analysis according to the previous description with slight modifications [58]. In short, the Illumina sequencing reads were mapped to the orthologous CDS regions of each gene using TopHat (V1.2.0) with the following parameters: -g1 -r 200 --mate-std-dev 20 -I 10000 [59]. All possible SNP sites with the minimum read depth of 10 reads were then identified by SAMTools (V0.1.13) based on aligned outcomes [60]. The analyses of amino acid mutation and functional annotation were performed by a custom-written algorithm.