SOURCES OF TISSUE-SPECIFIC TRANSCRIPTOME DATA

The raw data of B. dorsalis fat body, testis, and MAG/ED transcriptomes were sequenced and can be downloaded from the Sequence Read Archive of the National Center of Biotechnology Information (NCBI). The accession numbers of the datasets were SRR1026844 (Yang et al. 2014), SRR1032039 (Wei et al. 2015b), and SRR1168415 (Wei et al. 2016). All the samples were collected from the same laboratorial strain under the same condition. The fat body sample was collected from different developmental stages, whereas testis and MAG/EG samples were collected from male adults.

TRANSCRIPTOME RE-ASSEMBLY AND PROTEIN PREDICTION

The raw data of the MAG transcriptome were re-assembled together with the fat body and testis transcriptomes. Overlapping sequences of specific lengths with no base calls (N) were combined in contigs to form longer fragments. De novo re-assembly of these transcriptomes was carried out with the short read assembling program Trinity based on clean reads (Grabherr et al. 2011). The Trinity program consists of 3 independent software modules: Inchworm, Chrysalis, and Butterfly, which were applied sequentially to process large volumes of RNA-seq reads. All of this process was conducted in the same way as in our previous transcriptome sequencing (Wei et al. 2015b).

The resulting sequences were distinct singletons, and then the sequence clustering software Tgicl was used to splice and remove the redundancy to obtain long, non-redundant unigenes ( http://sourceforge.net/projects/tgicl/files/tgicl%20v2.1/). Then, all unigenes were divided into 2 classes, resulting in 1 cluster with the prefix CL followed by cluster ID with similarity more than 70%, and the 2nd cluster with the prefix Unigene followed by cluster ID. In the final step, Blastx alignment was performed between unigenes and proteins in databases such as NCBI nr, Swiss-Prot, and assigned to Kyoto encyclopedia of genes and genomes (KEGG) and clusters of orthologous groups (COG) with an E-value of 10−5, and the best alignment results were used to determine sequence direction of unigenes. The protein annotation was retrieved with the highest sequence similarity to unigenes. If results from different databases were conflicting, a priority order of NCBI nr, Swiss-Prot, KEGG, and COG was followed. If a unigene did not align to any of the above databases, the software of ESTScan was used to determine its sequence direction (Iseli et al. 1999). The direction from 5′ to 3′ ends was provided for unigenes with confirmed directions; for those without direction, sequences from the assembly software were provided.

FUNCTIONAL ANNOTATION OF UNIGENES

Gene annotation was conducted after re-assembly, which provided the necessary information of gene expression and function, including protein functional annotation, COG functional annotation, and gene ontology (GO) functional annotation. Unigene sequences were compared using Blastx to proteins in databases, namely NCBI nr, Swiss-Prot, KEGG, and COG with an E-value <10−5. Using enzyme commission (EC) number terms, biochemical pathway information was collected by downloading relevant maps from the KEGG database, which contains a systematic analysis of inner cell metabolic pathways and functions of individual gene products (Kanehisa & Goto 2000). COG is a database in which orthologous gene products are classified. The functions of unigenes assembled from the 3 tissues were predicted and classified by aligning to COG database. Homology searches were then carried out by query of the NCBI nr database using the Blastx algorithm. After annotation, the Blast2Go program was used to obtain GO annotations (Conesa et al. 2005), and the software WEGO (Ye et al. 2006) was used to perform GO functional classification to understand the distribution of gene functions at the macro level.

DIFFERENTIAL GENE EXPRESSION ANALYSIS

To determine specific expression of genes, expression profiles were calculated using the reads per kilobase mapped (FPKM) method (Mortazavi et al. 2008), which does not include the length of different genes and sequencing level as parameters during calculation of gene expression. The false discovery rate (FDR) control statistical method was used to test multiple hypotheses to correct for P-values (P < 0.01) (Benjamini & Yekutieli 2001). The smaller the FDR, the larger the difference in the expression level between the 2 samples. To identify these differential genes expressed in each tissue in the current analysis, sequences were evaluated with FDR ≤⃒ 0.001. Genes with FPKM < 0.1 were filtered as no expression in tissue. The absolute value of log2Ratio ≥⃒ 1 was used to judge the significance of differences in the initial gene expression analysis. After the analysis, we identified many genes that were only expressed or detectable in MAG/ED. Such genes were thereafter aligned to GO, COG, and KEGG assignments.

VALIDATION OF THE EXPRESSION OF MAG-SPECIFIC GENES

To validate the expression of MAG-specific genes, different tagmata and tissues including head (5 males), thorax (5 males), abdomen (5 males), midgut (20 males), fat body (50 males), Malpighian tubules (50 males), testis (50 males), and MAG/ED (100 males) were cut off or dissected from 3-d-old adult B. dorsalis males. Total RNA was isolated from each sample by using the TRIzol reagent (Invitrogen, Carlsbad, California). After DNA digestion with DNase (Promega, Madison, Wisconsin), total RNA (1 µg for each) was reverse transcribed into 1st-strand cDNA by using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China). The quantitative real-time polymerase chain reaction (qRT-PCR) included 5.0 µL GoTaq qPCR Master Mix (Promega, Madison, Wisconsin), 3.5 µL nuclease-free water, 0.5 µL template cDNA, and 0.5 µL of each primer (10 µM). The relative expression was determined by the StepOne Plus Real-Time PCR System (Life Technologies, Woodlands, Singapore). The thermal program was 1 cycle of 95 °C for 2 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 30 s. At the end of each run, a melting curve was generated to rule out the possibility of primer-dimer formation. The sequences of the specific primer sets were designed using online software Primer3 (version 4.0.0) (Table S1). A fragment of the ribosomal protein subunit 3 from B. dorsalis was also amplified as the inner reference (Wei et al. 2015b). Three biological replicates were performed for each sample. The relative expression levels were calculated using the 2−DDCt method (Livak & Schmittgen 2001). Data were analyzed statistically by 1-way analysis of variance (ANOVA) for expression comparisons in the SPSS 16.0 software (SPSS Inc., Chicago, Illinois), and a value of P < 0.05 was considered to be statistically significant.

ASSEMBLY AND ANALYSIS

Re-assembly of transcriptome data produced 31,782 unigenes with a mean length of 922 bp, and the total length of nucleotides was 29,300,417 bp (Table 1). The results indicated that the re-assembly sequences were longer than those in each transcriptome. A greater number of distinct clusters (7,426) suggested more isoforms were assembled in this study because more raw data were used in the assembly program. Of these unigenes, 3,655 (11.50%) unigenes were longer than 2,000 bp, whereas 1,257 unigenes in fat body, 2,438 unigenes in testis and 1,705 unigenes in MAG tissues were longer than 2,000 bp (Figure S1).

Table 1.

Assembly summary of transcriptomes of fat body, testis, and male accessory glands and ejaculatory duct (MAG/ED) from Bactrocera dorsalis.

ANNOTATION OF PREDICTED PROTEINS

After Blastx against the NCBI nr, nt, Swiss-Prot, KEGG, COG, and GO databases, 21,306 distinct sequences (67.04%) were matched to known genes encoding functional proteins (Table 1). Most of these sequences (20,180) were annotated in the NCBI nr database, whereas 12,456 were annotated in NCBI nt, 15,657 in Swiss-Prot, 13,782 in KEGG, 7,168 in COG, and 13,577 in GO. All the sequences annotated to NABI nr database were thereafter deposited at DDBJ/EMBL/GenBank with the Transcriptome Shotgun Assembly project accession number of GEYS0000000. Taking annotation in the nr database for the following analysis, there were almost 83.42% (16,835 out of 20,180) unigenes that showed strong similarity to 36 Drosophila species with 12.99% of the sequences matching homologous sequences from D. virilis Sturtevant, followed by D. willistoni Sturtevant (12.05%), D. melanogaster Meigen (11.73%), D. mojavensis Patterson (10.86%), D. ananassae Doleschall (8.19%), D. grimshawi Oldenberg (7.82%), and D. pseudoobscura pseudoobscura Frolova (7.22%) (Figure S2A). Only 306 sequences (1.52%) were annotated to the genus Bactrocera. The most important reason was the availability of the genome sequences of many Drosophila species, but the genome sequence of B. dorsalis had not been released when we performing the annotation. For the predicted proteins, over half of the hits had over 60% similarity with known proteins in the nr database (Figure S2B).

FUNCTIONAL ANNOTATION OF UNIGENES

Annotations of the unigenes provided information on their potential functions in each tissue of B. dorsalis. All of the re-assembled sequences were aligned to the COG database for functional prediction and classification. In total, 7,168 genes were assigned to 25 functional categories (Figure S3). Similar with each specific tissue transcriptome, the most common category was “general function prediction only,” which included 2,856 distinct sequences. The top 3 smallest groups were “RNA processing and modification” (1.20%), “extracellular structures” (1.00%), and “nuclear structure” (0.04%). For GO functional analysis, 13,577 sequences were categorized into 59 functional groups of 3 ontologies: biological process, cellular component, and molecular (Figure S4). There were 28 functional groups that contained more than 1,000 unigenes, including 17 biological processes, 8 cellular components, and 3 molecular functions. The most dominant groups in 3 ontologies were “cellular process” (9,880 unigenes) in biological process, “cell part” (7,839 unigenes) in cellular component, and “binding” (7,200 unigenes) in molecular function. KEGG annotation is usually used to categorize gene functions emphasizing biochemical pathways (Kanehisa & Goto 2000). In total, 13,782 sequences were remapped to 256 predicted pathways, of which 46 contained over 200 unigenes. Remarkably, the most abundant pathway was “metabolic pathways” including 1,919 sequences (Table S2).

DIFFERENTIALLY EXPRESSED GENES (DEGs)

Descriptive and quantitative transcriptome analysis is important to interpret the functional elements of genomes and reveal the molecular constituents of cells and tissues. To identify differentially expressed genes (DEGs) in MAG/ED, fat body, and testis tissues, we compared the transcriptomes from the 3 tissues and identified a number of genes highly presented in each transcriptome. Altogether, 25,054 unigenes were expressed in all the 3 tissues. In this study, transcriptome sequencings were based on a mixture of total RNA from different age/development stages of male flies. Genes with tissue-specific expression may be critical to the function of each tissue. Among all the unigenes, there were 458 unigenes that were expressed specifically in MAG/ED, 385 in fat body, and 2,238 unigenes in testis tissues (Fig. 1). Another 207 unigenes were expressed lowly in the 3 tissues. In testis tissue, mitosis and meiosis take place in the process of spermatogenesis. It is possible that many specific unigenes are involved, explaining the particularly large number of unigenes expressed specifically in testis tissue.

In this study, there were 182 unigenes that were 100-fold higher expressed in MAG/ED than both fat body and testis tissues (Table S3). These unigenes were considered as functional genes highly expressed in MAG/ED. The MAG/ED-specific unigenes (641) were classified into 20 categories (Fig. 2) after COG annotation. Interestingly, the most prominent category was “cell wall/membrane/envelope biogenesis.” Categories of “cell cycle control, cell division, chromosome partitioning” and “transcription” also made up a large proportion. All these MAG/ED-specific highly expressed genes were also classified into 3 main categories and assigned to 45 subcategories by using the GO assignment (Fig. 3). Different from each tissuespecific transcriptome, the predominant groups in biological process were “cellular process,” “single-organism process,” and “multi-organism process.” All the MAG/ED-specific highly expressed genes were thereafter mapped to 6 groups of KEGG pathways including 105 pathways (Fig. 4). Remarkably, “metabolic pathways” was the most important pathway, which contained 32 unigenes. In total, there were 51 unigenes that were involved in immunity, e g., “maturity onset diabetes of the young” involved in endocrine and metabolic diseases, “Vibrio cholerae infection” and “tuberculosis” involved in bacterial infection, “Herpes simplex infection” and “Epstein-Barr virus infection” involved in viral infection.

Fig. 1.

Statistics of sequences expressed specifically in each analyzed tissue of Bactrocera dorsalis.

MINING FOR FUNCTIONAL GENES HIGHLY EXPRESSED IN MAG/ED

To identify MAG/ED-specific genes in this study, we compared the MAG/ED transcriptomes with those of fat body and testis, and identified a number of genes highly expressed in the MAG/ED. Because genes in MAG/ED have a rapid evolution frequency, most of the unigenes were annotated with non-specific functional descriptions; 294 out of 641 unigenes were distinct singleton sequences, and only 53 of them were predicted genes (Table S3). Notable unigenes specific to MAG/ED were circadian clock—controlled protein (Unigene4743), C-type lectin (Unigene12723), vitellogenin-1 (Unigene12065), venom serine protease (Unigene8523), perlucin (Unigene5896), glucose dehydrogenase (Unigene9651), lipase 1 (Unigene5218), and pheromone-binding protein—related protein 5 (Unigene20808). These MAG/ED-specific genes may be involved directly in the male accessory proteins' functions.

The expression profiling of several annotated unigenes (>500 bp) were also validated at transcriptional level by qRT-PCR. The results showed that all of these genes had high expressions in MAG/ED tissue (Fig. 5). Besides, 4 unknown distinct unigenes (>500 bp) with the largest FPKM values were determined in tissues of male adults, and all of the 4 unigenes were highly expressed in MAG/ED tissue (Fig. 6). Other MAG/ED-specific genes sequences of <500 bp or less then 100-fold expressed compared with other tissues were also identified, such as oocyte maturation arresting factor (Unigene21945), antigen 5 precursor (Unigene21744), N-acetylgalactosaminyltransferase (Unigene22691), reverse transcriptase polyprotein (Unigene 21663), and potassium-dependent sodium—calcium exchanger (Unigene21750). These protein-coding genes also had high expression in MAG/ED tissue (Fig. S5). Regulatory mechanism determination of genes expressed differentially or specifically would provide new insights on the regulation of complex processes such as reproduction, development, and evolution.