The developmental transcriptome of the bamboo snout beetle Cyrtotrachelus buqueti and insights into candidate pheromone-binding proteins

Cyrtotrachelus buqueti is an extremely harmful bamboo borer, and the larvae of this pest attack clumping bamboo shoots. Pheromone-binding proteins (PBPs) play an important role in identifying insect sex pheromones, but the C. buqueti genome is not readily available for PBP analysis. Developmental transcriptomes of eggs, larvae from the first instar to the prepupal stage, pupae, and adults (females and males) from emergence to mating were built by RNA sequencing (RNA-Seq) in the present study to establish a sequence background of C. buqueti to help understand PBPs. Approximately 164.8 million clean reads were obtained and annotated into 108,854 transcripts. These were assembled into 24,338, 21,597, 24,798, 21,886, 24,642, and 83,115 unigenes for eggs, larvae, pupae, females, males, and the combined datasets, respectively. Unigenes were annotated against NCBI non-redundant protein sequences, NCBI non-redundant nucleotide sequences, Gene Ontology (GO), Protein family, Clusters of Orthologous Groups of Proteins/ Clusters of Eukaryotic Orthologous Groups (KOG), Swiss-Prot, and KEGG Orthology databases. A total of 17,213 unigenes were annotated into 55 sub-categories belonging to three main GO categories; 10,672 unigenes were classified into 26 functional categories by KOG classification, and 8,063 unigenes were classified into five functional KEGG categories. RSEM software for RNA sequencing showed that 4,816, 3,176, 3,661, 2,898, 4,316, 8,019, 7,273, 5,922, 5,844, and 4,570 genes were differentially expressed between larvae and males, larvae and eggs, larvae and pupae, larvae and females, males and females, males and eggs, males and pupae, females and eggs, females and pupae, and eggs and pupae, respectively. Of these, three were confirmed to be significantly differentially expressed between larvae, females, and males. Furthermore, PBP Cbuq7577_g1 was highly expressed in the antenna of males. A comprehensive sequence resource of a desirable quality was constructed from developmental transcriptomes of C. buqueti eggs, larvae, pupae, and adults. This work enriches the genomic data of C. buqueti, and facilitates our understanding of its metamorphosis, development, and response to environmental change. The identified candidate PBP Cbuq7577_g1 might play a crucial role in identifying sex pheromones, and could be used as a targeted gene to control C. buqueti numbers by disrupting sex pheromone communication.

Introduction system is crucial to sexual communication and reproductive isolation in insects. Finally, differentially expressed candidate PBPs underwent transcriptome data validation.

Insect rearing and collection
Eggs, larvae, and pupae of C. buqueti were collected in July 2015 from in the bamboo planting base of Lushan City, Sichuan Province, China (102˚91 0 N,30˚13 0 E). The field studies did not involve endangered or protected species, and no specific permission was required for the research activity at this location. Pupae were reared in our laboratory at 25˚C ± 1˚C and 70 ± 10% relative humidity, with a 12L: 12D photoperiod. Adults were used in the experiment 3 days after emergence [21]. Female and male adults were placed on ice and were quickly dissected into the head (without antenna), thorax (without thoracic legs), abdomen, antenna, and thoracic legs for qRT-PCR analysis. All samples were immediately frozen in liquid nitrogen and stored at −80˚C until use. Each sample contained eggs, larvae, pupae, males, females and adult tissues from at least five insects, respectively. After mixed sample, three biological replicates were conducted for each treatment.

RNA extraction and sequencing
Eggs, mixed larvae from the first instar to the prepupal stage, pupae, and adults (females and males) from emergence to mating (3 days old) were prepared for RNA extraction. Total RNA was extracted using TRIzol reagent (Qiagen, China Shanghai). The concentration of total RNA was quantified using a spectrophotometer (Implen, Westlake Village, CA), and the RNA integrity was tested using 1.5% agarose gel electrophoresis. After RNA extraction, mRNAs were purified using a Poly A T tract mRNA isolation system and collected using an RNeasy RNA reagent. Mixed mRNAs were fragmented into 300-800 bp pieces using RNA Fragment reagent (Illumina), and the pieces were collected using an RNeasy RNA cleaning kit (Qiagen). Subsequently, RNA fragments were copied into first strand cDNA using MMLV reverse transcriptase (TaKaRa, Dalian Liaoning, Chinese) and random primers. Second strand cDNA synthesis was performed using DNA Polymerase I and RNase H. The Illumina HiSeq2000 system and 100 paired-end reads were used for sequencing. Clean reads were obtained by removing adaptors, low-quality reads, and contaminated reads from raw sequence reads. Statistical analysis of the sequence length was performed to ensure sequence purity.

Assembly and functional annotation
Raw sequence data reads in fasta format were firstly processed through in-house Perl scripts [22]. In this step, clean data (clean reads) were obtained by removing reads containing adapter, ploy-N and low-quality reads form raw read data [23,24]. At the same time, Q20, Q30, GCcontent, and sequence duplication level of the clean data were calculated [24]. All downstream analyses were based on clean high-quality data.
A flow chart of transcriptome assembly as described by Grabherret et al. [25] was used in the present analyses. A Perl pipeline as described by Haas et al. [22] was used for analyzing sequence data. As suggested by Haas et al. [22], if multiple sequencing runs are conducted for a single experiment, these reads may be concatenated into two files in the case of paired-end sequencing. The left files (read 1 files) from all samples were pooled into a single large left.fq file, and right files (read 2 files) into a single large right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity (http://trinityrnaseq.github.io) with min_kmer_cov set to two by default and all other parameters set default. The assembled unigenes were annotated by BLASTX and ESTScan against Nr, Nt, Pfam, KOG/COG, Swiss-Prot, KO, and GO databases (E < 10 −5 ), and the best annotations were selected [23,24,26] (S1 Table). Differentially expressed genes were selected by log2 fold change > 1 and q value < 0.005 according the method of DESeq [27]. The nucleotide sequences of each identified PBP gene are listed in S2 Table. Homology analysis A neighbor-joining (NJ) tree was constructed with MEGA version 5.0 [28] and the Jones-Taylor-Thornton model. The olfactory gene sequences of other coleopteran insects were first transcribed into their amino acid sequences using the ORF finder (http://www.ncbi.nlm.nih. gov/gorf/gorf.html). Olfactory genes of other coleopteran species were obtained from the NCBI databases. Bootstrap support values were based on 1000 replicates. All of the candidate olfactory genes were named according to the nomenclature system described previously [29,30].

RT-PCR and qRT-PCR validation
Specific primer pairs were derived from the transcriptome data, and primer pairs for each gene were designed to amplify 100-200 bp products, which were verified by sequencing. A conventional RT-PCR (Bio-Rad S1000, US) analysis was performed for each primer pair using rTaq DNA polymerase (TaKaRa, Dalian, Liaoning, China) before the qRT-PCR (Bio-Rad CFX96, US) analysis to ensure that the correct products were amplified and no primer dimers were present. The qRT-PCR analysis was carried out using an Mx 3000P detection system (Stratagene, La Jolla, CA) as described previously, with thermal cycler parameters of 1 min at 95˚C, then 40 cycles of 30 s at 95˚C, 30 s at 60˚C, and 30 s at 72˚C. GAPDH of C. buqueti (Gen-Bank accession number: SAMN06176790) was used as the housekeeping gene. A standard curve was derived from 10-fold serial dilutions of plasmid containing the target DNA segment to determine the PCR efficiency and for quantifying the amount of target mRNA. All primers tested gave amplification efficiencies of 90%-100%. For each treatment, three biological replicates were conducted. qRT-PCR data were analyzed by the 2 -ΔΔCT method [31]. The primers used in this experiment were designed with Primer premier 5.0 and Oligo 6.0 and are listed in S3 Table. The qRT-PCR data were analyzed and output as PDF files using Graphpad 5.0.

Illumina sequencing and assembly
Clean reads were obtained from raw reads after the removal of those reads with low quality, adapters, duplicated, and ambiguous. This resulted in a total of 31,469,916, 36,773,825, 32,128,345, 33,070,448, and 31,434,121 clean reads in eggs, larvae, pupae, females, and males of C. buqueti, respectively. All clean reads were assembled into transcripts by Trinity software; the longest copy of redundant transcripts was regarded as a unigene [22,24,25]. A total of 108,854 transcripts were obtained and assembled into 83,115 unigenes. Many unigenes were 200-1000 bp in length (Table 1), while approximately 14.7% unigenes exceeded 1000 bp, and 7.2% exceeded 2000 bp (Table 1) Table 2). The number and percentage of unigenes annotated in these databases were counted. The Nr database had the best match against the CP-CE-CL-CF-CM Combined dataset (21,058, 25.33%) (  Table 2) (S1-S15 Texts).

CDS prediction
A total of 21,102 (25.39%) unigenes were predicted via BLASTX with an E-value threshold of 10 −5 in the Nr, and the Swiss-Prot database (Figs 5 and 6). Among these, 14,998 unigenes were in the length of more than 300 bp (Fig 5). Furthermore, 17,703 (21.30%) unigenes were then predicted using ESTScan, which identified 3,415 unigenes of more than 300 bp in length (Fig 6).

Differentially expressed genes
A total of 7,273, 3,661, 4,570, and 5,844 genes were differentially expressed between pupae and males, pupae and larvae, pupae and eggs, and pupae and females, respectively, with 1,484 common to pupae, males, larvae, eggs, and females ( Fig 7A). A total of 8,019, 3,176, 4,570, and 5,922 genes were differentially expressed between eggs and males, eggs and larvae, eggs and pupae, and eggs and females, respectively, with 1,110 common to eggs, males, females, pupae, and larvae (Fig 7B). A total of 4,316, 5,922, 5,844, and 2,898 genes were differentially expressed between females and males, females and eggs, females and pupae, and females and larvae, respectively, with 580 common to females, males, eggs, pupae, and larvae ( Fig 7C). A total of 4,316, 8,019, 7,273, and 4,816 genes were differentially expressed between males and females, males and eggs, males and pupae, and males and larvae, respectively, with 2,043 common to males, females, eggs, pupae, and larvae ( Fig 7D). A total of 4,816, 3,176, 3,661, and 2,898 genes were differentially expressed between larvae and males, larvae and eggs, larvae and pupae, and larvae and females, respectively, with 539 common to larvae, males, eggs, pupae, and females ( Fig 7E) (S19-S28 Texts). More genes were shown to be expressed in eggs than in pupae, in females than in pupae, in females than in eggs, in females than in larvae, in larvae than in pupae, in larvae than in eggs, in males than in pupae, in males than in eggs, in males than in females, and in males than in larvae (2, 576, 3,223, 3,080, 1,554, 2,129, 1,737, 4,563, 4,705, 2,706, and 3,266, respectively; Fig 8). Conversely, fewer genes were shown to be expressed in eggs than in pupae, in females than in pupae, in females than in eggs, in females than in larvae, in larvae than in pupae, in larvae than in eggs, in males than in pupae, in males than in eggs, in males than in females, and in males than in larvae (1,994, 2,621, 2,842, 1,344, 1,532, 1,439, 2,710, 3,314, 1,610, and 1,550, respectively; Fig 8).

Phylogenetic analysis of candidate PBP
We constructed a phylogenetic tree comparing Cbuq7577_g1 (Gene name: CbuqPBP1) and the olfactory genes from 28 coleopteran insects (Fig 9) (S29 Text). In this tree, the olfactory genes are well separated with strong bootstrap support. Cryptolaemus montrouzieri CmonOBP2 and Colaphellus bowringi CbowOBP26, along with CbuqPBP1, were located on the same clade.

Expression profiles of pheromone-binding proteins
We identified 19 candidate PBPs through the Nr database (nucleotide sequences are listed in S30 Text). Of these, significant differences in expression profiles were identified in 13 candidate PBPs in male adults and larvae (Table 3), 10 candidate PBPs in female adults and larvae (Table 4), and three candidate PBPs in female and male adults ( Table 5).

Validation of transcriptome data
To validate the transcriptome result, we selected 12 significant differentially expressed genes from Tables 3-5 for quantitative reverse transcriptase-PCR (qRT-PCR) confirmation. Six PBPs transcripts which have demonstrated by RNA-seq to be enriched in larvae were confirmed by qRT-PCR (Fig 10) (S31 Text). Additionally, Cbuq7577_g1 had significantly higher transcriptional level in male than in female and larvae with 5.36 and 85.19 fold exchanges. Moreover, to further explore tissue-and sex-specific expression, we selected Cbuq7577_g1 for Candidate PBP in Cyrtotrachelus buqueti qRT-PCR confirmation. We observed the highest expression of PBP Cbuq7577_g1 in male antennae, followed by male heads, compared with low levels of expression in female antennae and heads (Fig 11).

Overview of transcriptome data
The transcriptome is the complete set of expressed RNA transcripts in one or more cells [32], and its analysis enables the study of gene transcription and the characteristics of transcriptional regulation. High-throughput sequencing technology has been applied to the transcriptome study of many species in class Insecta, such as Phyllotreta striolata [33], D. melanogaster [34], Biston betularia [35], Aedes aegypti [36], Brugia malayi [37], and Bemisia tabaci [38].
In the present study, developmental transcriptomes were established of C. buqueti eggs, mixed-age larvae, pupae, and male and female adults, providing a relatively comprehensive gene pool. The numbers of clean reads from egg, larval, pupal, female and male transcriptomes were 31,469,916, 36,773,825, 32,128,345, 33,070,448, and 31,434,121, respectively. All these clean reads were assembled into 108,854 transcripts by Trinity software. A total of 83,115 Candidate PBP in Cyrtotrachelus buqueti Fig 8. Volcano plot of differentially expressed genes in eggs, larvae, pupae, males, and females. a Volcano plot of differentially expressed genes between CE and CP. b Volcano plot of differentially expressed genes between CF and CP. c Volcano plot of differentially expressed genes between CF and CE. d Volcano plot of differentially expressed genes between CF and CL. e Volcano plot of differentially expressed genes between CL and CP. f Volcano plot of differentially expressed genes between CL and CE. g Volcano plot of differentially expressed genes between CM and CP. h Volcano plot of differentially expressed genes between CM and CE. i Volcano plot of unigenes were annotated by Nr, Nt, KO, Swiss-Prot, PFAM, and KOG/COG. Thousands of differentially expressed genes were identified, which facilitates developmental and evolutionary studies of C. buqueti, and contributes to future work in bamboo snout beetle comparative genomics.
differentially expressed genes between CM and CF. j Volcano plot of differentially expressed genes between CM and CL. Splashes represent different genes. Blue splashes means genes without significant different expression. Red splashes mean significantly upregulated genes. Green splashes mean significantly downregulated genes. CE, CL, CP, CM, and CF represent eggs, larvae, pupae, males, and females of C. buqueti, respectively. https://doi.org/10.1371/journal.pone.0179807.g008    Candidate PBP in Cyrtotrachelus buqueti

Pheromone-binding proteins
Previously reported physiological functions of PBP include: binding specificity, transporting specific pheromone molecules, and filtering odorant molecules entering the sensor chamber [39]; acting as a carrier to transport pheromones through the hemolymph to the receptor [40]; forming the PBP-pheromone complex for receptor recognition [41], cascade initiation, and deactivation to restore receptor sensitivity [42]; and protecting pheromones from enzymatic degradation [41]. Identifying the developmental transcriptome of C. buqueti provides an opportunity to understand the physiological function of PBPs. A total of 19 candidate PBPs were identified in the present study. The alignment of CbuqPBP1 and 27 known coleopteran insect olfactory gene sequences, and a phylogenetic tree indicated that CbuqPBP1, CmonOBP2 and CbowOBP26 are on the same clade. Therefore, it was speculated that such genes may have the same ancestral gene, but were differentiated by adaptation to different types of environmental chemical factors during evolution, and perform the same or similar functions among different species. Candidate PBP in Cyrtotrachelus buqueti Thirteen candidate PBPs in male adults and larvae (Table 3), 10 candidate PBPs in female adults and larvae (Table 4), and three candidate PBPs in female and male adults (Table 5) showed significant differences in expression, with Cbuq7577_g1 demonstrating significant differences in expression among larvae, male adults, and female adults. Cbuq7577_g1 showed 49% identity with the OBP of Dendroctonus ponderosae (AGI05175), and 47% identity with the OBP of Lissorhoptrus oryzophilus (AHE13794). AGI05175 and AHE13794 were previously functionally annotated as insect pheromone/OBP domains. Cbuq7577_g1 in C. buqueti showed very low similarity to genes in the NCBI database, which likely reflects the limited research that has been carried out into Curculionidae and Lepidoptera PBPs. To research gene function, a PBP gene of C. buqueti (CbuqPBP1) was cloned in this study for prokaryotic expression. Using N-phenyl-1-naphthylamine as the fluorescent probe in a competitive binding assay, the ability of CbuqPBP1 to bind 12 sex pheromone analogs and three volatiles of Neosinocalamus affinis shoots was examined. These results will be published later.
qRT-PCR results of the present study showed that candidate PBP Cbuq7577_g1 in C. buqueti is expressed in male and female adult antennae, which is consistent with the expression pattern of PBPs in most insects, such as Heliothis virescens [43], Manduca sexta [44], Spodoptera exigua [45], and B. mori [46]. PBP expression in adult females may enable the identification of hydrophobic pheromones in the male or the monitoring of pheromones released by the female, as well as transporting pheromones and general odorant molecules. Although it is generally believed that insect PBPs are only expressed in the antennae, researchers have also documented their expression in the head, feet, wings, and other body parts [46,47]. Zhang et al. [48]found HarmPBP2 of Helicoverpa armigera was expressed in female's maxillary palp, but the highest expression in the antennae. The candidate PBP Cbuq7577_g1 was mainly expressed in antennae (97.07%). Its expression level in male antennae was 14.43 times that in female antennae. Cbuq7577_g1 may play an important role in the identification of odorant molecules, specifically those involved in identifying females from the external environment through C. buqueti antennae.

Conclusions
We constructed a comprehensive, good-quality sequence resource from the developmental transcriptomes of C. buqueti eggs, larvae, pupae, and female and male adults. This resource enriches what is known about C. buqueti genomics, thus facilitating our understanding of metamorphosis, development, and fitness to environmental change. The identified candidate PBP Cbuq7577_g1 might play a crucial role in identifying sex pheromones, and could be used as a target to control C. buqueti as a pest by disrupting sex pheromone communication.