Comparative transcriptome analysis of salivary glands of two populations of rice brown planthopper, Nilaparvata lugens, that differ in virulence.

BACKGROUND: The brown planthopper (BPH), Nilaparvata lugens (Stål), a destructive rice pest in Asia, can quickly overcome rice resistance by evolving new virulent populations. Herbivore saliva plays an important role in plant-herbivore interactions, including in plant defense and herbivore virulence. However, thus far little is known about BPH saliva at the molecular level, especially its role in virulence and BPH-rice interaction. METHODOLOGY/PRINCIPAL
FINDINGS: Using cDNA amplification in combination with Illumina short-read sequencing technology, we sequenced the salivary-gland transcriptomes of two BPH populations with different virulence; the populations were derived from rice variety TN1 (TN1 population) and Mudgo (M population). In total, 37,666 and 38,451 unigenes were generated from the salivary glands of these populations, respectively. When combined, a total of 43,312 unigenes were obtained, about 18 times more than the number of expressed sequence tags previously identified from these glands. Gene ontology annotations and KEGG orthology classifications indicated that genes related to metabolism, binding and transport were significantly active in the salivary glands. A total of 352 genes were predicted to encode secretory proteins, and some might play important roles in BPH feeding and BPH-rice interactions. Comparative analysis of the transcriptomes of the two populations revealed that the genes related to 'metabolism,' 'digestion and absorption,' and 'salivary secretion' might be associated with virulence. Moreover, 67 genes encoding putative secreted proteins were differentially expressed between the two populations, suggesting these genes may contribute to the change in virulence.
CONCLUSIONS/SIGNIFICANCE: This study was the first to compare the salivary-gland transcriptomes of two BPH populations having different virulence traits and to find genes that may be related to this difference. Our data provide a rich molecular resource for future functional studies on salivary glands and will be useful for elucidating the molecular mechanisms underlying BPH feeding and virulence differences.

Introduction

Insect herbivore saliva contains digestive enzymes, such as alkaline phosphatase, esterase, amylase and β-glucosidase, as well as other components, such as elicitors that induce plant defense, effectors that inhibit plant defense, and proteins related to pathogen transmission [1]–[3]. Some studies have also found a relationship between saliva components and herbivore virulence [4]. Therefore, herbivore saliva, the first substance to come into chemical contact with the plant, plays important roles in both food ingestion and interactions between herbivores and their host plants [1]–[5]. Characterizing herbivore saliva will provide new insights into plant–herbivore interactions, including induced plant defense and herbivore virulence.

To characterize herbivore saliva, the transcriptome and/or proteome of the salivary glands and/or saliva of several herbivore species – mostly hemipterans such as rice brown planthopper (BPH; Nilaparvata lugens (Stål)); Hemiptera: Delphacidae) [6], [7], pea aphid (Acyrthosiphon pisum; Hemiptera: Aphididae) [8], [9], [10], green peach aphid (Myzus persicae; Hemiptera: Aphididae) [11], [12], whitefly (Bemisia tabaci (Gennadius); Hemiptera: Aleyrodidae) [13], and potato leafhopper (Empoasca fabae (Harris); Hemiptera: Cicadellidae) [14] – have been analyzed. These studies found several hundred proteins in the saliva [4], [6]–[14]. However, whether differences in salivary components exist between herbivore populations with different virulence traits and what functions these components have remain largely unanswered questions.

BPH, one of the most destructive insect pests of the rice plant (Oryza sativa L.) in Asia, causes substantial losses of rice yield every year by sucking phloem sap and transmitting plant viruses, such as the rice ragged stunt virus and the rice grassy stunt virus [15]. The cultivation of resistant rice varieties is an important control measure for the BPH. However, the BPH rapidly overcomes rice resistance by evolving new virulent populations [16]. BPH virulence strains generally correspond to particular resistance genes in rice. For example, rice varieties TN1 (a susceptible variety) and Mudgo (carrying the resistance gene bph1) host virulent BPH strains 1 and 2, respectively [16]. Although some minor differences in morphology and in the composition of bacterial symbionts among virulence genotypes have been reported, the mechanisms underlying changes in BPH virulence are not clear [16]–[19].

When sucking sap, BPHs and other phloem feeders secrete two primary kinds of saliva: coagulable and watery saliva. Coagulable saliva forms salivary sheaths, which help to stabilize the insects' stylets and suppress plant defense responses to components of the watery saliva [1], [4]. Watery saliva, which contains a mixture of amino acids, proteins, and digestive enzymes, assists the movement of the stylets inside the salivary sheath, the digestion of plant material, and the suppression of plant defense responses [1]–[4], [9], [12]. Given the important roles of herbivore saliva in plant-herbivore interactions, BPH saliva is likely to play a central role in the interaction between BPHs and rice, including in the evolution of BPH virulence.

Salivary gland morphology and several salivary proteins of BPHs have been reported [20]–[24]. Moreover, Noda et al. (2008) identified 2383 expressed sequence tags (ESTs) in the salivary glands of the BPH by Sanger’s method [6]. However, these data are insufficient to elucidate the nature of BPH saliva. Moreover, whether BPH populations with different virulence have different saliva components remains unclear. To explore these issues, we compared the salivary-gland transcriptomes of two BPH populations with different virulence that were maintained on one of two rice varieties: TN1 (TN1 population) or Mudgo (M population). Using next-generation high-throughput sequencing technology, a total of 37,666 and 38,451 unigenes were identified in the respective populations, and 3,757 unigenes were differentially expressed between the two populations. Moreover, among the differentially expressed genes, 67 unigenes encoded putative secretory proteins. These findings provide an exciting opportunity to elucidate the components of BPH saliva and suggest a possible correlation between BPH saliva and virulence.

Results and Discussion

Illumina Sequencing and Read Assembly

The cDNA libraries of the salivary glands of BPH TN1 and M populations were sequenced using an Illumina platform, resulting in 38,487,746 and 40,350,780 reads, respectively. After cleaning and quality checks, short sequences were assembled into 97,103 and 101,885 contigs, respectively. Using paired end-joining and gap-filling methods, these contigs were further assembled into 73,510 and 76,195 scaffolds, which were then clustered into 37,666 and 38,451 unigenes, respectively (Table 1). Finally, sequence data from the two libraries were combined to obtain 43,312 unigenes with a mean length of 746 bp. The length distributions of the unigenes were similar between the two populations, suggesting there was no bias in the cDNA library construction (Figure S1). However, some unigenes were found in only one population (data not shown). We believe these differences may have been caused by differential long-term ecological adaptations to the different rice hosts. The total number of unigenes obtained was much higher than the number of BPH salivary gland ESTs identified by Noda et al. [6], indicating the depth of coverage possible using HiSeq technology.

Table 1

Summary statistics for salivary gland transcriptomes of two brown planthopper populations.

	TN1 population	M population
Total number of raw reads	38,487,746	40,350,780
Total number of clean reads	35,948,162	37,583,138
Total length of clean reads	3,235,334,580	3,382,482,420
Total number of contigs	97,103	101,885
Mean length of contigs	268	265
GC percentage	45.42%	45.28%
Total number of scaffolds	73,510	76,195
Total unique sequences	37,666	38,451
Mean length of unique sequences	599	609

Annotation of Salivary-gland Transcripts

For functional annotation, the 43,312 unigenes were searched using BLASTx against the non-redundant (nr) NCBI protein database. A total of 19,771 unigenes (45.65% of all distinct sequences) yielded significant BLAST hits with a cut-off E-value of 1.0E–5 (Table S1). The species distributions of the best matches for each sequence are shown in Figure 1. The highest percentages of unique sequences matched genes of the red flour beetle (Tribolium castaneum; 14.34%), A. pisum (11.48%), and a parasitoid wasp (Nasonia vitripennis; 7.90%). A similar result was also reported by Xue et al. [25], Peng et al. [26] and Bao et al. [27]. Further research should explore why the highest percentage of unique sequences matches genes of T. castaneum, a coleoptera beetle, rather than A. pisum, a hemiptera aphid which is much more closely related to BPH.

Figure 1

Species distribution of the BLASTX results of brown planthopper salivary-gland genes.

This figure shows the species distribution of unigene BLASTX results against the NCBI nr protein database with a cutoff E-value ≤1.0E–5 and the proportions of each species, represented by different colors. Species with proportions greater than 1% are shown.

GO and KEGG Orthology Classifications

Gene Ontology (GO) assignments were used to functionally classify the predicted proteins. A total of 18,500 sequences from TN1 population and 17,118 from M population were categorized into 38 functional groups (Figure 2). GO has three ontologies: ‘biological process,’ ‘molecular function,’ and ‘cellular component.’ The results of GO analysis showed a similar distribution of gene functions between the two populations (Figure 2). Among both populations, the three dominant basic ‘biological process’ categories are ‘cellular process,’ ‘metabolic process,’ and ‘biological regulation.’ These results indicated that the cells in the salivary glands were metabolically active, which correlated well with the biological function of the tissue of interest. In ‘molecular function,’ the three most abundant categories are ‘binding,’ ‘catalytic activity,’ and ‘transporter activity.’ This ‘binding’ category includes sequences annotated to be involved in protein, nucleic acid, ion, cofactor and enzyme binding. Saliva secretion is the main function of the salivary glands, and saliva is composed of water, electrolytes, lipids, amino acids and proteins. Not surprisingly, a large number of genes we found are involved in binding and transport. In ‘cell component’, the top three are ‘cell’, ‘organelle’ and ‘macromolecular complex’. This suggests that cells and organelles are the prominent parts of the salivary glands, and macromolecular complex secretion is also critical in the salivary glands.

Figure 2

Gene ontology (GO) classifications of salivary-gland genes of TN1 and Mudgo (M) populations.

The results are summarized in three main categories: biological process, molecular function, and cellular component. The y-axis shows the number of matching unigenes in a category.

To investigate which biological pathways are active in the salivary glands, all of the sequences were assigned to the reference canonical pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG). A similar distribution of biological pathways for TN1 and M populations was also found by KEGG mapping, and a total of 11449 unigenes from TN1 population and 11502 unigenes from M population were mapped separately to 239 and 240 pathways in total (Tables S2 and S3). The salivary glands are specific organs for salivary macromolecule production and, consequently, have a high level of metabolic activity. This characteristic was implied in a previous study by observing the dense cytoplasm and organized whorls of rough endoplasmic reticulum in the salivary-gland cells [28]. Consequently, among these pathways, “metabolism” is important in the salivary glands (1625 unigenes from TN1 population and 1675 unigenes from M population, Figure 3, in pathways associated with human diseases were excluded). Since there are abundant proteins in the saliva of BPHs, as a major original location of these proteins, the salivary glands should be active in protein synthesis and catabolism. Indeed, the transcriptome contained many sequences involved in the “protein processing in endoplasmic reticulum” pathway (254 genes in TN1 population and 257 genes in M population), which is related to the formation and transport of nascent secretory proteins, and the ‘lysosome’ pathway (209 genes in TN1 population and 218 genes in M population), a major pathway related to protein degradation. The GO annotations and KEGG orthology (KO) classifications indicated that the salivary glands might be active in metabolism, binding and transport, and that a high percentage of genes might be involved in secretory protein processing and secretion.

Figure 3

KEGG pathway distributions of unigenes from salivary glands of TN1 and Mudgo (M) populations.

The top fifteen pathways (excluding human disease-related) with the highest percentages of mapped unigenes are shown.

Unigenes Encoding Putative Secreted Proteins

We characterized the salivary secretome of the BPH by combining the transcripts of the salivary glands of the two populations. Secreted proteins are probably delivered into saliva through the eukaryotic endoplasmic reticulum–Golgi pathway [29]. Proteins secreted through this pathway have an N-terminal signal peptide. Therefore, all amino acid sequences were analyzed for the presence of signal peptides and potential cleavage sites. A total of 464 unique genes encoded a secretory signal peptide. Of these, 112 were predicted to have at least one transmembrane domain besides the signal peptide, indicating that their proteins were likely embedded in cell membranes of the salivary glands. Thus, 352 potential secretory proteins were retained (Table S4), which is comparable to the numbers of putative secreted proteins reported in A. pisum (324) and B. tabaci (295) [10], [13]. Interestingly, the possible functions of some putative secreted proteins were closely related to the known roles of insect saliva, such as digestion and suppressing or eliciting plant defenses.

Among the putative secreted proteins in the BPH, a set of digestive enzymes and hydrolases, including plant cell wall (PCW)-degrading enzymes were found. These putative PCW-degrading enzymes included one β-1,4-endoglucanase (Unigene1860_All), one β-glucosidase (Unigene26172_All), and two β-1,3-glucanases (Unigene10762_All and Unigene23029_All). PCW, a thick, rigid polysaccharide structure comprising an extensive network of celluloses, hemicelluloses and pectins, forms a formidable barrier to herbivores. To overcome this barrier, herbivores have evolved to secrete salivary PCW-degrading enzymes, such as cellulases (consisting of endoglucanases and β-glucosidases) and pectinases [1], [29], [30]. For example, an extract of Homalodisca vitripennis salivary gland contained a variety of enzymes capable of hydrolyzing glycosidic linkages in the polysaccharides of PCWs, including β-1,4-endoglucanase, xyloglucanase and β-1,4-endoxylanase; these enzymes could digest xyloglucans in cell walls to facilitate H. vitripennis stylet penetration and were related to the formation of the salivary sheath [31]. In addition to degrading cellulose cooperating with endoglucanases, salivary β-glucosidases in the B- and D-follicles of BPHs act on phenolic glucosides, such as arbutin and salicin [24], [30]. It has been reported that β-1,3-glucanases secreted by plant pathogens can degrade β-1,3-glucans in callose [32]. Callose has been reported to exist in the PCWs of a wide variety of higher plants and can be induced by either biotic or abiotic stress, including BPH feeding [33], [34]. Callose deposition plays an important role in preventing BPHs from ingesting phloem sap [34]. Thus β-1,3-glucanases in BPH saliva may play a role in countering rice defenses by inhibiting the formation and deposition of callose.

In addition to PCW-degrading enzymes, other digestive enzymes and hydrolases found in putative secreted proteins in BPH salivary glands included 11 trypsin-like proteins or serine proteases (Table S5), 3 glycosyl hydrolases [α-amylase (Unigene9457_All), sucrase (Unigene8043_All) and α-glucosidase (Unigene3598_All)], 2 lipases (Unigene18263_All and Unigene41792_All) and a phospholipase A2 (Unigene23547_All). In the plant phloem sap, there are many proteins (including defense-related proteins), sugars, amino acids, organic acids, ATP and so on [35], [36]. Therefore, these putative secreted proteins, like PCW-degrading enzymes, may function in extra-oral digestion and/or the suppression of defenses. Salivary α-glucosidase in the G-follicle of the principal gland of the BPH, for example, hydrolyzes sucrose and trehalose [24]. α-Amylases aid in pre-digesting starch to maltose by hydrolyzing α-D-(1,4)-glucan bonds; the maltose is then hydrolyzed to glucose by an α-glucosidase [37]–[40]. A recent study showed that a salivary lipase (MdesL1) of Mayetiola destructor was involved in extra-oral digestion and host cell permeability [41]. In another study, phospholipase A2 in the salivary glands of Manduca sexta may have been involved in the hydrolysis of fatty acid moieties from dietary phospholipids and also in killing food-borne bacteria [42]. The functional characterization of these hydrolases and of the digestive enzymes of the BPH may shed light on this species' mode of extra-oral digestion and plant defense suppression.

Of the putative secreted proteins, we also identified 3 odorant-binding proteins (OBPs) and 9 chemosensory proteins (CSPs) (Table S6), and all nine putative CSPs showed high homology to MP10 (E-value ≤1.0E–5), a salivary effector found in aphids [43]. OBPs and CBPs were generally reported to be preferentially expressed in antennae, mouthparts, and other chemosensory organs [44]–[47], and to be involved in perireceptor events of the chemosensory system. Because these proteins transport chemical signals to sensory neuron receptors, they can be said to regulate insect behavior, such as communication among conspecifics, predator detection, and food/host location [45], [48]. Recently, some members of OBPs and CSPs were also found to be broadly expressed in various insect tissues, including salivary glands; in addition, their functions are involved in many areas, such as nymph survival [49], development of embryonic integuments [50], insect leg regeneration [51], and insect-host interactions [43], [52], [53]. CSPs, for example, were the most abundant proteins in Vanessa cardui and V. gonerilla mandibular gland lumens, and they were thought to play an important role in host plant recognition, helping the insect perceive phagostimulants and deterrents, and in the detection of microorganisms on the leaf surface [52]. Some OBPs in mosquito saliva that are predicted to function in olfaction and gustation are also secreted into host cells to manipulate host physiology by scavenging host amines [53]. Of the three identified salivary effectors – C002, Mp42 and Mp10, in aphids – Mp10, homologous to a CSP, induced chlorosis and local cell death in Nicotiana benthamiana and suppressed the oxidative burst induced by flg22 [43]. This suggests that the putative secreted CSPs and OBPs found in BPH salivary glands may also play a role in BPH-rice interactions, an issue that is worthy of future study.

As a consequence of co-evolution, BPH saliva may contain some components to counter plant secondary chemicals. In the BPH salivary secretome, we found five putative secreted carboxylesterases (Unigene3085_All, Unigene12962_All, Unigene13873_All, Unigene14160_All, and Unigene35902_All). It has been well documented that carboxylesterases in insect midguts and fat body tissues play important roles in metabolizing xenobiotics (i.e., plant toxic allelochemicals and synthetic chemicals, such as pesticides) [27], [54]–[58]. In plant phloem sap, there are also many defense chemicals [36], [59]. Thus, these putative secreted carboxylesterases in BPH salivary glands might also be related to detoxification.

Interestingly, alkaline phosphatase (ALP) (Unigene1963_All) was identified in the salivary secretome; a ubiquitous enzyme in all organisms, it is characterized by its ability to hydrolyze orthophosphate monoesters at alkaline pH [60]. Evidence suggests that insect ALPs are expressed in various tissues, and are involved in absorption and transport, development, neural and renal function, cuticle sclerotization and xenobiotic tolerance [61]–[64]. ALPs have also been detected in the salivary glands of several insects, and their activity was confirmed in the saliva of wheat aphids (Diuraphis noxia) and B. tabaci Biotype B [65]–[68]. Salivary ALPs have been suggested to be involved in producing the stylet sheath, and thus they may aid in plant feeding [67], [68]. This possible role of salivary ALPs in the BPH will require further investigation.

In addition to aiding in digestion and suppressing plant defense, these putative proteins secreted by salivary glands may also be recognized by plants, thus eliciting plant defense responses. Thus far, some of these enzymes and/or their products, such as lipase and β-glucosidase, have been reported to elicit signaling pathways that activate both direct and indirect plant defense responses [69]–[72]. The lipase activity in Schistocerca gregaria oral secretions, for example, plays a central role in the accumulation of oxylipins, especially 12-oxo-phytodienoic acid [70]. In BPH saliva, β-glucosidase was also found to induce salicylic acid-, H2O2- and ethylene-mediated signaling pathways [71], [72]. It is expected that ongoing efforts will point to more such elicitors in herbivore saliva.

Although many homologs of insect secretory proteins were found in the BPH salivary secretome, some were not identified. For example, the pectinase complex, comprising endopolygalacturonases, rhamnogalacturonan lyases and pectin methylesterases, plays an important role in pectin degradation, which subsequently facilitates the decomposition of (hemi)cellulose and makes the PCWs more susceptible to further breakdown by other enzymes [29], [69], [73]. Pectinases occur mainly in Hemiptera and Coleoptera [69]. However, in the BPH salivary secretome, pectinases were not identified. This suggests that BPH may use unique feeding strategies. Alternatively, these enzymes may have been filtered out during the prediction of secretory proteins because they have incomplete 5′ sequences. In our experiments, for instance, the unigene Unigene1860_All was annotated to encode a cellulase. However, because it lacks the 5′ end, this unigene was filtered out. Therefore, to completely explore the protein composition of BPH saliva, other research approaches, such as the proteome analysis of BPH saliva, should also be used.

Differences in Salivary-gland Gene Expression Profiles between TN1 and M Populations

Gene expression levels can be estimated from Illumina sequencing data based on the number of raw reads [74]. To identify differentially expressed genes in the salivary glands of TN1 and M populations, the number of fragments mapped to each gene was calculated and normalized using fragments per kb per million fragments (FPKM) [74], [75]. We found 3,757 genes with significantly different expression levels (Figure 4, Table S7). Of these, 2,084 were up-regulated in TN1 population relative to M population. The detected fold differences (log2 ratio) ranged from −17 to +16, and 89.88% were up- or down-regulated by from 1.0- to 5.0-fold (Figure 5). To validate these gene expression data, we compared the expression profiles of salivary glands of TN1 and M populations using quantitative real-time polymerase chain reaction (QRT-PCR). Of the 21 randomly selected genes, 20 showed concordant fold differences between the types of analyses (Table S8), indicating that our results were reliable.

Figure 4

Bitmap of differentially expressed genes.

The red and green points represent up- and down-regulated genes in the salivary glands of TN1 population, respectively; the blue points represent genes with no differences in expression based on the criteria of the false discovery rate (FDR) ≤0.001 and an absolute value of the log2 ratio ≥1.

Figure 5

Fold difference distribution of differently expressed genes.

The x-axis shows the fold change (log2 ratio) in gene expression in TN1 population compared to Mudgo (M) population.

The differentially expressed transcripts were assigned to 31 functional groups, such as ‘biological process,’ ‘molecular function,’ and ‘cellular component’; ontology distributions are shown in Figure 6. Among the ‘biological process’ assignments, many genes were assigned to ‘cellular process,’ ‘metabolism process,’ and ‘biological regulation.’ In the ‘molecular function’ category, most genes were related to ‘catalytic activity,’ ‘binding,’ and ‘transporter activity.’ ‘Cell’ formed the main part of the ‘cell component’ category. The differentially expressed genes coding for these physiological processes may affect the saliva composition of these populations; further functional studies must be performed to validate this hypothesis.

Figure 6

Distribution of significantly differentially expressed transcripts in gene ontology (GO) subclasses.

The y-axis shows the number of the matching differentially expressed genes in that category.

To gain insight into the dominant biological pathways of the differentially expressed genes that mapped to KO, a hypergeometric test was performed to explore statistically enriched pathways. In total, 33 enriched pathways (P≤0.05) were identified (Table S9; pathways associated with human disease were excluded). Ten gene sets were correlated to ‘metabolism’ and four to ‘digestion and absorption’ (Table 2). Carbohydrates, especially sucrose, are the dominant chemical component in rice phloem sap, and they are particularly important macronutrients for BPHs. Digestible carbohydrates are used primarily for energy but can also be converted to fat and stored, and their carbon skeletons can contribute to amino acid production. Low levels of carbohydrate metabolism can limit insect growth and performance. Indeed, many gene groups (72 genes) involved in ‘carbohydrate metabolism’ and ‘carbohydrate digestion and absorption’ were enriched, including ‘galactose metabolism’, ‘starch and sucrose metabolism’, and ‘carbohydrate digestion and absorption’. Of these 72 genes, 46 (63.89%) were up-regulated in the salivary glands of M population (Table 2). There are also many differentially expressed genes related to ‘metabolism’ and ‘digestion and absorption’ of other nutrients and substances, such as α-linolenic acid, glycerolipid, amino acids, fat, vitamins, and proteins (Table 2). Interestingly, the ‘salivary secretion’ pathway including 30 genes was also enriched; of these genes, 20 genes (66.67%) were up-regulated in the salivary glands of M population (Table S9). Salivary secretion plays an important role in plant-insect interactions. This analysis indicates that the difference in ‘metabolism’, ‘digestion and absorption’ and ‘salivary secretion’ may contribute to the change in virulence between the two populations. However, more evidence is required to evaluate this hypothesis.

Table 2

Differentially-expressed genes related to ‘metabolism’ and ‘digestion and absorption’ pathways.

KEGG Pathway	P-value	Total1	Up-regulated2	Down-regulated3
etabolism
Alpha-linolenic acid metabolism	6.76E-05	14	8	6
Glycerolipid metabolism	4.89E-04	26	12	14
Amino sugar and nucleotide sugar metabolism	1.98E-03	24	21	3
Starch and sucrose metabolism	3.01E-03	30	17	13
Pyrimidine metabolism	3.15E-03	52	15	37
Retinol metabolism	4.24E-03	15	5	10
Galactose metabolism	4.56E-03	22	14	8
Glycerophospholipid metabolism	5.26E-03	27	18	9
Drug metabolism - cytochrome P450	1.76E-02	12	4	8
Linoleic acid metabolism	4.41E-02	5	4	1
Digestion and absorption
Fat digestion and absorption	8.13E-07	29	13	16
Vitamin digestion and absorption	1.29E-05	25	7	18
Protein digestion and absorption	8.98E-03	39	15	24
Carbohydrate digestion and absorption	1.14E-02	20	15	5

Number of differentially expressed genes in salivary glands belonging to each KEGG pathway.

Number of genes up-regulated in brown planthopper M population relative to TN1 population in each KEGG pathway.

Number of genes down-regulated in brown planthopper M population relative to TN1 population in each KEGG pathway.

Among the genes encoding putative secreted proteins of BPH salivary glands, 67 were found to be differentially expressed between the BPH populations. Of these, 43 genes were up-regulated, while 24 were down-regulated in M population compared to TN1 population (Table S10). The detected fold differences (log2 ratio) ranged from –6.82 to +4.86, and 79.10% of the genes were up- or down-regulated by from 1.0- to 3.0-fold (Table S10). As stated above, many proteins in herbivore saliva play important roles in herbivore feeding and herbivore-plant interactions; thus these differentially expressed putative secreted proteins might mediate BPH virulence.

Conclusions

This study provides the first comprehensive evaluation of the salivary gland transcriptome of BPH using high-throughput sequencing and comparative expression analysis between two BPH populations, TN1 and M, with different virulence. In total, we obtained 43,312 unique sequences and 3,757 differentially expressed genes. GO and KO analysis indicate that genes related to metabolism, binding and transport were significantly active in the salivary glands. A total of 352 genes were predicted to encode secretory proteins, and some might play important roles in BPH feeding and BPH-rice interactions. Comparative analysis of the transcriptomes of the two populations revealed that the genes related to ‘metabolism,’ ‘digestion and absorption,’ and ‘salivary secretion’ might be associated with virulence traits. Furthermore, some putative secretory proteins that were differentially expressed in the two populations were identified. These findings provide a valuable resource for future investigations of the composition of BPH saliva and the role of saliva in changes in BPH virulence and interactions between BPH and rice.

Materials and Methods

BPH Cultures, Salivary-gland Collection, and RNA Isolation

Two different virulent populations of BPH, TN1 and M, were maintained on rice varieties TN1 and Mudgo, respectively, for more than 170 generations at 27±1°C and 70±10% relative humidity under a 14/10 h light/dark photoperiod. The original insects were provided by the Chinese National Rice Research Institute (Hangzhou, China). BPH adult females, within 2 days after eclosion, were placed in a Petri dish on ice, and their salivary glands were dissected in phosphate buffer saline (pH 7.2) using forceps. For RNA extraction, the salivary glands of about 150 females were pooled for each population and directly dipped in RNA Lysis Buffer (Promega, Fitchburg, WI, USA). Total RNA was isolated using the SV Total RNA Isolation System (Promega) according to the manufacturer's protocol. The concentration and quality of total RNA were determined by a NanoDrop spectrophotometer (Thermo Fisher, Waltham, MA, USA).

Library Construction of the Salivary Glands and Illumina Sequencing

The salivary-gland cDNA library was prepared using a SMARTer™ PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA) and an Advantage 2 PCR Kit (Clontech) following the kit’s instructions. After end repair and adaptor ligation, the products were amplified by PCR and purified using QIAquick PCR Purification Kit (Qiagen, Hilden, Germany). This cDNA library was then sequenced on a Hiseq2000 Illumina sequencing platform in BGI-Shenzhen (Shenzhen, China). Raw reads were generated using a Solexa GA Pipeline 1.6 (Illumina). After the removal of empty and low-quality reads with undetermined nucleotides (‘N’), processed reads were assembled using SOAPdenovo software and clustered with TGI Clustering tools [76], [77]. All raw transcriptome data was deposited in the NCBI short read archive (SRA) with accession number SRX276865 (TN1 population) and SRX276866 (M population). The assembled sequences whose orientations were predicted via blast or ESTScan (recorded in 5′ to 3′) have been deposited in the NCBI transcriptome shotgun assembly (TSA) database with accession number 214044 TSA. The generated unigenes were screened with BLAST (http://blast.ncbi.nlm.nih.gov) and annotated against the nr and SwissProt databases with an E-value cut-off of 1.0E–5. GO and KO annotations of the unigenes were determined using BLAST2GO (http://www.blast2go.org/) and Inter-ProScan software [13], [25].

Secretory Protein Prediction

To translate the unigene nucleotide sequences into amino acid sequences, we first searched all the combined transcriptome unigenes against protein databases using BLASTx (E-value ≤1.0E–5, switch -F F) in the following order – nr, Swiss-Prot, KEGG, and clusters of orthologous groups (COG) – until the sequences had significant hits. The BLAST results were used to extract coding sequences (CDSs). The CDSs of unigenes with no significant BLAST hits were predicted by ESTScan [13]. We then used the SignalP 4.0 Server (http://www.cbs.dtu.dk/services/SignalP/) to predict the presence of signal peptides and cleavage sites in the amino acid sequences. To predict transmembrane domains, we submitted each amino acid sequence with a signal peptide to the TMHMM Server v. 2.0 (http://www.cbs.dtu.dk/services/TMHMM/) [13], [43]. Putative proteins with a signal peptide and 0–1 transmembrane domain (the signal peptide can be a transmembrane domain) were considered to be potential secreted proteins [13].

Analysis of Differential Gene Expression

Salivary gland genes that were differentially expressed between TN1 and M populations were identified using a table of counts constructed with FPKM values, which adjusted the number of fragments by the total number of fragments mapped and the length of the gene [74], [75]. The false discovery rate (FDR) was used to determine threshold P-values in multiple test and analysis. An FDR<0.001 and an absolute value of the log2 ratio >1 provided significance thresholds for gene expression differences.

QRT-PCR Analysis

To confirm the results of FPKM analysis, the expression levels of 21 randomly selected salivary-gland genes were measured in TN1 and M populations by QRT-PCR. Total RNA from each sample (salivary glands of about 120 females per sample) was extracted using the SV Total RNA Isolation System. For each total RNA sample, 0.75 µg of RNA was reverse-transcribed using the PrimeScript RT–PCR Kit (TaKaRa, Otsu, Japan). QRT-PCR was carried out in a CFX96TM Real-Time System (Bio-Rad, Hercules, CA, USA) using iQ™ SYBR Green (Bio-Rad) under the following conditions: 95°C for 3 min, then 40 cycles of 95°C for 5 s, 60°C for 15 s, and melting curve generation from 65 to 95°C. Primers used for QRT–PCR are listed in Table S11. Each gene was analyzed in three biological replications, after which the average threshold cycle (Ct) per sample was calculated. The endogenous actin gene was used for normalization. A no-template control sample (using nuclease-free water) was run to detect contamination and to determine the degree of dimer formation (data not shown). Finally, relative expression levels were calculated as 2–ΔΔCt.

Identification of Statistically Enriched KEGG Pathways

The differentially expressed genes were used for KO enrichment analysis using the hypergeometric test to measure significantly enriched terms. The formula was:

In this equation, N indicates the number of genes with KO annotations and n the number of differentially expressed genes in N. The variables M and m represent the numbers of genes and of differentially expressed genes, respectively, in each KO term. The threshold to determine significant enrichment of the gene sets was corrected to P-value ≤0.05.

Length distribution of unigenes in salivary-gland transcriptomes of brown planthopper populations. The x-axis shows the calculated lengths of unigenes in the salivary-gland library and the y-axis shows the number of unigenes. (A) TN1 population. (B) M population.

(TIF)

Click here for additional data file.

Top BLASTX hits from NCBI non-redundant (nr) database. All unigenes were screened using BLASTX against the NCBI nr database with a cut-off E-value ≤1.0E–5.

(XLSX)

Click here for additional data file.

KEGG orthology (KO) annotation of unigenes from the salivary-gland transcriptome of TN1 population of the brown planthopper (BPH).

(XLSX)

Click here for additional data file.

KEGG orthology (KO) annotation of unigenes from the salivary-gland transcriptome of Mudgo (M) population of the brown planthopper (BPH).

(XLSX)

Click here for additional data file.

Genes encoding putative secreted proteins. Based on the coding sequences of translated unigenes, the presence of a signal peptide was predicted using SignalP v. 4.0, and cleavage position was calculated, resulting in 464 proteins with putative secretion signals. Transmembrane domains in these proteins were predicted with the TMHMM Server v. 2.0. After removing sequences for proteins that were likely embedded in cell membranes, 352 predicted secretory proteins were retained. They were annotated by searching against the nr, Swiss-Prot, and COG databases.

(XLSX)

Click here for additional data file.

Genes encoding putative secreted trypsin-like proteins and serine proteases.

(XLSX)

Click here for additional data file.

Genes encoding putative secreted odorant-binding and chemosensory proteins.

(XLSX)

Click here for additional data file.

Differentially expressed genes indicated by the FPKM statistical analysis. FPKM, FDR-value, fold change (log2 ratio) of gene expression, and best hits against nr (E-value ≤1.0E–5) for all gene pairs are listed in this table.

(XLSX)

Click here for additional data file.

Gene expression data verified by quantitative real-time polymerase chain reaction (QRT-PCR). Twenty-one unigenes were selected for expression-level validation using QRT-PCR.

(XLSX)

Click here for additional data file.

Statistically significantly enriched KEGG pathways. Human disease-associated pathways were excluded from this analysis.

(XLSX)

Click here for additional data file.

Differentially expressed genes encoding putative secreted proteins indicated by FPKM values.

(XLSX)

Click here for additional data file.

Primers used in quantitative real-time polymerase chain reaction (QRT-PCR) for validation of differentially expressed genes.

(XLSX)

Click here for additional data file.