Transcriptome and expression analysis of Sogatella furcifera (Horváth) (Hemiptera: Delphacidae) in response to cycloxaprid


 Background

The white backed planthopper, Sogatella furcifera (Horváth), causes substantial damage to many crops by direct feeding or transmitting plant viruses, especially southern rice black-streaked dwarf virus (SRBSDV) that has become a great threat to rice production in East Asia. Cycloxaprid, a novel cis-nitromethylene neonicotinoid insecticide, has good industrialization prospects because of its high efficiency against rice planthoppers, including imidacloprid-resistant populations. This chemical would be used extensively in future. However, very little is known about the influence on S. furcifera after cycloxaprid application at the molecular level. We sequenced S. furcifera transcriptome of female adult of S. furcifera using the Illumina sequencing.

Results

By de novo transcriptome assembling and massive parallel pyrosequencing, we constructed two transcriptomes of S. furcifera and profiled the alternation of gene expression in response to cycloxaprid treatment in transcriptional level. Over 157,906,456 nucleotides and 131,601 different unigenes were generated using Illumina technology from both cycloxaprid-treatment and no-treatment S. furcifera. And a total of 38,534 unigenes matched known proteins in at least one database, accounting for 29.28% in total unigenes. The number of Coding DNA Sequence (CDS) were 28,546 and that of the amino acid sequence of coding region were 22,299. A total 15,868 simple sequence repeats (SSRs) were identified. The trinucleotide repeats accounted for 45.1% (7,157) in total SSRs and the (AAG/CTT)n was the most frequent motif. There were 359 differentially expressed genes (DEGs) that might be induced by cycloxaprid. There were 131 genes up-regulated and 228 down-regulated. 22 unigenes may take participate in the resistance to cycloxaprid, such as cytochrome P450, Glutathione-s transferase (GST), Acid phosphatase (ACP), cadherin, etc.

Conclusions

Our study will provide a splendid database for future investigations of the resistance mechanism induced by cycloxaprid and will provide new strategies for pest management and crop protection.


Conclusions
Our study will provide a splendid database for future investigations of the resistance mechanism induced by cycloxaprid and will provide new strategies for pest management and crop protection.

Background
Cycloxaprid, which firstly reported in 2008, patented and developed in China and first named in 2011 [1][2][3]. Its action mode was known as targeting insect nicotinic acetylcholine receptors (nAChRs) and was similar to that of imidacloprid and nitenpyram but cycloxaprid has already been proven more effective than three other neonicotinoids: imidacloprid, thiamethoxam and nitenpyram for the control of piercing and sucking insect pests (aphids, leafhoppers, whiteflies, planthoppers etc), especially Sogatella furcifera in future [4] [5]. As an oxabridged cis-configuration neonicotinoid insecticide, cycloxaprid has been considered a substitute for imidacloprid, especially the management of imidaclopridresistant insects because it has performed well in controlling a broad spectrum of insect pests and has low toxicity for mammals [6] [7].
The ecological , toxicological and physiological perspectives of S. furcifera and other hemipteran insect pests, such as Nilaparvata lugens, Laodelphax striatellus and Bemisia tabaci have been extensively studied [18][19][20][21][22][23][24]. Transcriptomes of Nilaparvata lugens and Laodelphax striatellus were reported using next-generation DNA sequencing techniques [25] [26], and that responsed to SRBSDV [16]and to high and low doses of cycloxaprid [4] had also been studied. Furthermore, the olfactory receptor gene family [27] and the salivary glands [28] of S. furcifera had been sequenced, ect. The possible factors related to cycloxaprid-resistance and related genes by insecticide applications had been conducted. By treated with a low-dose and a high-dose treatment, these findings, based on transcriptome, showed that CYP4, CYP6, and GSTd contributed significantly to insecticide detoxification resistance [4]. However, it is little known to the whole transcripts and unigenes change after treated with cycloxaprid, In this study, hence, we constructed two transcriptomes of S. furcifera and profiled the alternation of gene expression in response to cycloxaprid induction at the transcriptional level. As a whole, 131,601 unigenes have been identified and 359 DEGs have been discovered.

Bioassay
By bioassay, the LC10 and LC50 concentration of cycloxaprid-treatment was established.

RNA isolation
The values of Optical Density (OD260/280) ranges from 1.962 to 2.005, the RNA integrity number (RIN) value of six samples were larger than 7.0, indicating that they were pure and not degraded. So, the Ribonucleic Acid (RNA) quality of all samples meet the requirements of sequencing. (Table 1) Transcriptome Assembly of S. furcifera The S. furcifera transcriptomes were generated from 5 female adults emerged in 48h after treated by cycloxaprid LC10 concentration (HYFA) and 5 no-treatment female adults (CKFA) via the Illumina sequencing. Three replicates were prepared. We then constructed a mass gene database that contains a total of 157,906,456 nucleotides (nt). After eliminating low quality reads from the raw reads, there were 45,934,376 average clean reads remained in cycloxaprid-treatment, which accounted for 96.00% of the raw reads and 42,495,937 average clean reads remained in no-treatment, which accounted for 93.90% of their raw reads. The short reads assembly software Trinity was used to assemble, compare and ligate all clean data so as to get the unigenes from the S. furcifera transcriptome. These reads were assembled into 131,601 unigenes, the mean length is 720 bp. The 67.38% of unigenes ranged from 200-500 bp and 15.28% of unigenes ranged from 500-1000 bp. The 17.33% of unigenes ranged from more than 1000 bp in length. (Figure 1, Table 2 Table 4. A total of 26,204 unigenes returned an above cut-off blast hit to NR database of the National Center for Biotechnology (NCBI), which account to 19.91% of 131,601 unigenes from the S. furcifera transcriptome. The species distribution of the top blastx hits for each unique sequence was shown in Figure 2. The unambiguous assembled sequences revealed that the percentage of matches sequences was 15.3% from zootermopsis nevadensis, followed by Acyrthosiphon pisum (5.2%), Tribolium castaneum (4.7%), and the greatest number of no-matches was 67.8%.
Based on the S. furcifera transcriptome assembly, 27,079 (20.57%) sequences were annotated in the GO database, which were divided into a total of 49 groups in three ontology categories (biological process, cellular component, molecular function). The "biological process" ontology category contains 24 groups. The largest group is "cellular process" with 15,163 unigenes, followed by "metabolic process" with 14,033 unigenes and the smallest group is "rhythmic process" with 30 unigenes. The "cellular component" ontology category contains 17 groups. The largest groups are "cell" and "cell part" with the same numbers, 8,460 unigenes, respectively. And the smallest group is "extracellular matrix component" with only four unigenes. The "molecular function" ontology category contains 8 groups. The largest group is "binding" with 14,136 unigenes, and the smallest group is "transcription factor activity, protein binding" with 314 unigenes (Figure 3).
In order to annotate the detail function of genes, KOG database was used. In total, 15,447 unigenes (11.76%) were annotated and these genes were divided into 26 categories. and a total of 3,104 unigenes (20.06%) were placed into the "General function prediction only" category. Followed by "signal transduction mechanisms" (1,986, 12.83%), "posttranslation modification, protein turnover, chaperones" (1,585, 10.24%), and the smallest category is "Cell motility" with 21 unigenes. There are three unigenes annotated "Unamed protein", accounting for only 0.02% of the functionally annotated unigenes. (Figure 4) To identify the genes that involved in metabolic pathways, a total of 15,058 unigenes (11.44%) were mapped to the KEGG pathway database. These unigenes were divided into 32 pathways. The largest pathway is "signal transduction" with 1,829 unigenes (15.98%).
To validate the transcriptome, we compared the expression profiles of the cycloxapridtreatment and no-treatment using qRT-PCR. We selected 18 unigenes randomly, Primers used in qRT-PCR for validation of DEGs were shown in Table 5. 14 of which demonstrated a concordant direction of change and 4 unigenes were inconsistent between transcriptome and qRT-PCR ( Table 6). The results indicated that our results are reliable.
These results are confirmed with these reports showing that trinucleotide repeats (AAG/CTT)n are the most generous microsatellites in coding ESTs [16][29] [30]. Similarly, Xu(2012) also reported that (AAG)n was the most frequent motif (39.8%) in S. furcifera [16]. However, it had been confirmed by the former work that (AAC)n is the most frequent motif in the small brown planthopper (L. striatellus) [26]. The difference of predicted trinucleotide repeats in the two transcriptomes may give an expression on distinct adaptability of these insects. The numerous potential molecular markers gained from our study would play a particularly important role on genetic mapping, genotyping and parentage analysis of these insect species [16].

Identification of DEGs
To identify DEGs, the numbers of clean reads for each gene were counted, then individual sets of reads were mapped back to the previously assembled transcript and counted as a proxy for gene expression. DESeq [31] was used to identify the differentially expressed transcripts between the two samples (HYFA and CKFA), padj<0.05 (with biological replicates). The numbers of up-regulated genes was 131 and down-regulated was 228, in total 359 DEGs ( Figure 7).
217 DEGs were annotated in the GO database, inluding 65 up-regulated genes, which was approximatelty half of the 131 genes, and 152 down-regulated, which accounted to 66.67% of the 228 genes. These DEGs were assigned into groups according to these genes functions, biological process, cellular component and molecular function, and in the biological process, the largest groups were annotated into the "phosphorylation" and "protein phosphorylation". In the "cellular component" ontology category, the main groups are showed in Figure 8. The largest group of the "molecular function" ontology category is "catalytic activity".
The DEGs were mapped to the KEGG pathway database, which involved in metabolic pathways, These genes were divided into 179 pathways. The largest pathway is "Pathways in cancer" with 59 unigenes. Followed by "Proteoglycans in cancer" with 40 unigenes.
Based on the annotation results of basedata, we identified 3 categories of genes that associated with metabolic resistance, 2 Cytochrome P450s and 2 GSTs and target-site resistance, 1 Cadherin, etc in 359 DEGs. (Table 8)

Conclusion
In this study, we obtained 131,601 unigenes with mean lengths of 720 bp, which are a major genomic resource for clarifying the genes expression induced by cycloxaprid in S. furcifera. A total of 29,326 unigenes matched known proteins in the NCBI database. A total 15,868 simple sequence repeats were identified. 359 differentially expressed genes could be related with the cycloxaprid-induction in S. furcifera. And 22 unigenes may associated with cycloxaprid. such as cytochrome P450, GST, ACP, cadherin, etc.

Discussion
Insect resistance to insecticides generally related to each of three main mechanisms: lower epidermal tissue penetration; metabolic resistance [32], which include the overexpression or detoxification of some critical enzymes, such as cytochrome P450 [4] [33] [34], esterase and GST [35]. And target-site resistance, such as GABA receptor [36], AChE, sodium channel [37] etc.
Cytochrome P450 superfamilies were universal and constituted by many complex members, including 70 families, with 127 subfamilies, of hydrophobic, hemecontaining enzymes. P450s are involved in the biosynthesis of several essential endogenous compounds and the detoxification of many xenobiotics [38] [39]. In insects, P450s are induced by a wide range of inducers [40] and could mediate resistance to all classes of insecticides [41]. And overexpression of P450 genes has been reported to be related with development of insecticide resistance in many insect species, such as Nilaparvata lugens (Stål) to imidaclorprid [33] [42], Laodelphax striatellus (Fallen) to ethiprole [43], Musca domestica to phenobarbital [39], Anopheles gambiae [44] [45] and Meligethes aeneus Fabricius [35] to pyrethroids, etc. More than half of the P450 genes belong to either the CYP4 or CYP6 families [46] [47]. The CYP4 family includes sequences from vertebrates and are involved in hydroxylations of steroids, fatty acids and prostaglandins [48], whereas CYP6 enzymes, insect specific, appear to be involved in the metabolism of xenobiotics and plant products, e.g., pyrethroids and furanocoumarins [49]. In our study, only 2 P450 genes from DEGs were annotated as 1 CYP 4 family and 1 CYP 6 family.
Phosphatases, including acid and alkaline phosphatases(ACP and ALP), are capable of hydrolysis and transphosphorylation. They are important in the metabolism of carbohydrates, phospholipids and nucleotides and have been reported as enzymes significant in resistance to insecticides. ACP belongs to a group of enzymes that hydrolyze phosphomonoesters at acidic pH [63] and occupies a significant position in the detoxification of toxic compounds entering the body, acts as an enzymatic defence against foreign compounds and plays an essential role in maintaining normal physiological functions [64][65][66][67][68]. In the study, only one gene was founded to code acid phosphatases in 359 DEGs.
Cadherins represent a large and complicated family of calcium-dependent, transmembrane glycoproteins and a cytoplasmic tail that binds catenins which support to link the cadherin to the actin cytoskeleton, as well as also function in cellular signaling [69][70][71][72] and are responsible for maintaining the integrity of cell-cell contacts in multicellular organisms. In addition, cadherins include a class of proteins which were known to be related with Bt Cry protein binding and toxicity to Lepidoptera, Diptera and Coleoptera insects [73] [74], and have been identified in other invertebrates as crucial target receptors for Cry toxins [75][76][77]. Plutella xylostella was the first insect species which reported to evolve resistance to Bt in field populations [78] and its high resistance levels to Bt have been confirmed with the loss of binding affinity of Cry toxins to their protein receptors, including cadherin, on the midgut brush border membrane [79] [80]. Vadlamudi et al reported that a cadherin localized in the midgut epithelium of Manduca sexta serves as the receptor for Cry1A toxins of Bacillus thuringiensis [81] [82]. Gahan et al.(2001) showed that the resistance to Cry1Ac toxin in Heliothis virescens is associated to retrotransposon-mediated disruption of a specific cadherin gene. Moreover, confirmed that midgut epithelial cadherins participated directly in the entomopathogenicity of B. thuringiensis [76]. However, there is little evidence that cadherin was related with the resistance to other insecticides beside Cry toxins. In our study, only one gene from DEGs could code cadherin.

Insects and Insecticides
S. furcifera were collected from Huaxi district of Guiyang in Guizhou Province, China in July 2014 and bred down 20 generations under condition as hereinafter. The insects were reared on 10 day old rice seedlings cultured in plastic boxes (34×23.5×20cm) at 25±1°C temperature, 50% to 60% relative humidity, and L16:D8 h photoperiod in constanttemperature incubator [83].

Bioassay
The rice-stem dipping method [22][84][85] was used for the dose-responses of S. furcifera to cycloxaprid. Cycloxaprid was firstly dissolved into mother liquid by acetonitrile plus 10% Triton-100 (m/V; Solarbio Science & Technology Co., Ltd. Beijing, China) as an emulsifer. Then the mother liquid were serially diluted into 5-9 different concentrations with distilled water, respectively. The mortality was recorded 96h later for cycloxaprid.
The nymphs were considered dead if they failed to move when gently prodded with a fine bristle. The mortality data were corrected by the control mortality using Abbott's formula.

Sample collection and preparation
The mother liquid of cycloxaprid was diluted into 1×LC10 concentration and distilled water was used as the control. The rice seedlings were dipped into these solution for 30s, respectively. More than one hundred of the fifth instar nymphs were introduced on these treated rice seedlings when air-dried at room temperature. After 48h treated, emergency adults would be distinguished between male and female and collected in EP tubes, respectively. Then, 48h, 5 female adults (HYFA) were randomly sellected for subsequent RNA isolation, transcriptome sequencing. Three replicates were prepared (plus the control (no-treatment) (CKFA) in total 30 nymphs). All collected bodies were immediately frozen by liquid nitrogen and stored at -80°C until use [83].

RNA isolation, Library construction
Qubit® RNA Assay Kit (Life Technologies, CA, USA) was used to measure RNA concentration. 1% agarose gels was used to monitor RNA degradation and contamination.
The NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) was used to check the RNA purity. The RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) was used to assess RNA integrity. preferentially. Subsequently polymerase chain reaction (PCR) was performed with T100 thermal cycler PCR Instrument (Bio-rad, USA) and then PCR products were purified with AMPure XP system. Finally, library quality (insert size) was assessed on the Agilent 2100 system and effective concentration of the library was analysized with Q-PCR (>2nM).
Illumina Hiseq platform was used to sequence the library and generate the paired-end clean reads. TruSeq PE Cluster Kit v3-cBot-HS (Illumia) was used to perform the clustering of the index-coded samples on a cBot Cluster Generation System. Two transcriptomes of S. furcifera were constructed treatment by cycloxaprid (HYFA) and no-treatment (CKFA) [83].

Quality control
Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts.
In this step, excluding reads containing adapter, ploy-N and low quality reads from raw data, clean reads with high quality were obtained so as to subsequent analysis, Q20, Q30, GC-content and sequence duplication level of the clean reads calculated. All the downstream analyses were based on clean data with high quality [83]. Hochberg's approach was used to modify the resulting P values for controlling the false discovery rate. Genes with an adjusted P-value <0.05 was set as the threshold for obviously differential expression.
To verify the accuracy of the transcriptome, the expression levels of 18 genes were conducted by qRT-PCR. Total RNAs from each sample were extracted using Triozol RNA (Life Technology) and treated with DNase I (Invitrogen) according to the manufacturer's protocol. The concentration of each RNA sample was adjusted to 1 mg/mL with nucleasefree water and total RNA was reverse transcribed in a 20ul reaction system using the  Availability of data and materials Data Availability Statement: All data generated or analysed during this study are included in this published article.
The data that support the findings of this study are available from NCBI SRA database(Accession number: SRR4294200 and SRR4294203), but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of NCBI SRA database.

Competing Interests
This manuscript and its authors are not involved in any potential conflicts of interest, including financial interests and relationships and affiliations.        Table 7 The statistic information for special unigenes associated with insecticide resistance    The detail function of genes in KOG database. The x axis shows the groups name in KOG and the percentage of unigenes annotated in a group account to the total unigenes annotated in KOG database.