Identification and expression pattern of chemosensory genes in the transcriptome of Propsilocerus akamusi

Insect and tissue collection

Males of P. akamusi midges were collected from Tianjin, China, between October 2017 and November 2017 and then transported live to the lab. All samples were transferred within a few hours of eclosion. The antennae, heads (excluding antennae), and legs were excised from 1500 male specimens under microscope for transcriptome sequencing. The obtained tissues were then immersed in RNAlater (Ambion, AM7020) and stored at −20 °C until processing.

Olfactometer bioassay

The emerging P. akamusi males tend to congregate into cloud-like swarms which usually appear in the early evening. They seek to mate when females are around or just pass through, thus starting the life cycle. However, we still face deficient knowledge of the sex differences in mate preferences and wonder whether males will take the first step to attract females or not.

A mate selection behavioral test was therefore conducted using a horizontal Y-tube olfactometer consisting of a central glass tube with the length of 115 mm as well as two lateral glass arms with the length of 75 mm. The diameter for both the tube and arms was 22 mm. The test was performed in a dark room with a small incandescent lamp working as the only light source over the olfactometer. One lateral arm of the olfactometer was selected as the experimental arm and the other one is the control. Either the female subjects or the male ones were firstly placed on a piece of filter paper (10 ×10 mm) and this paper was then carefully inserted into the midpoint of the experimental arm. Another piece of paper without any subjects was inserted into the equivalent point in the control arm. After the olfactometer was connected to an assembled shunting device, the air sampler was switched on to pump airflow for 30s and then a test individual with the opposite sex was placed at the entrance of the central glass tube. If this test one might prefer to enter into the experimental arm and stayed there for at least 1min, the chosen arm could be considered as the one with proper pheromones for the insects. However, the midge was allowed to acclimate and make a preference within 5 min. Otherwise, we assumed that no response was given by the insect to the ones with the opposite sex and the trial was ended then. Besides, to prevent light from interfering with the taxis of P. akamusi, the experimental and control arms were swapped after completing the behavioral assays of three individuals. Fifty insects were tested for one trial with triplicates.

RNA preparation and cDNA library construction

The antennae, heads, and legs of 100 male insects were subjected to RNA preparation using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. After the treatment of DNase I, the mRNA molecules were harvested with Oligo (dT) from the total RNA. The long material of mRNA was then shortened into something compatible for further sequencing with the supplementation of Thermomixer buffer. With TransScript First-Strand cDNA Synthesis SuperMix (TransGen, Beijing, China), the fragmented mRNA next worked as the template for single-strand cDNA which was subsequently synthesized to double strand cDNA. Purification, end reparation, adaptor-ligation and size selection were carried out for final PCR amplification. During the QC (quality control) steps, Agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System were used for quantification and qualification of the sample library. The library ended up being sequenced using Illumina HiSeq 4000.

De novo transcriptome assembly

After sequencing, a preprocessing of the raw data was performed and consisted of a series of steps, including the trimming of adapter sequences, the removal of reads containing more than 5% ambiguous bases and the elimination of low quality basis. Here the low quality sequences referred to a situation where a single read comprised over 20% of bases with Q score less than 15. The remaining reads generated by these filtering procedures were defined as high quality clean reads and were consequently used for de novo assembly by Trinity v3.0 program (Eisen et al., 1998). This transcript assembly was then clustered into unigenes with the assistance of Tgicl (version: v2.0.6) (De Hoon et al., 2004).

Functional annotation of unigenes as well as prediction of coding sequences and simple sequence repeats

A homology search of functionally annotated unigenes was performed by BLAST (version: v2.2.23, website: http://blast.ncbi.nlm.nih.gov/Blast.cgi) against NCBI non-redundant protein sequences (NR), non-redundant nucleotide sequences (NT), Clusters of Orthologous Groups of Proteins (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO) and Swiss-Prot databases (Altschul et al., 1990; Conesa et al., 2005). Besides, InterPro analysis was measured by InterProScan5 (version: v5.11–51.0, website: https://code.google.com/p/interproscan/wiki/Introduction) for providing predictive information of gene families and their critical domains (Quevillon et al., 2005). TransDecoder (https://github.com/TransDecoder/TransDecoder/releases) was utilized to recognize candidate coding regions and peptide sequences within the Trinity-created transcript assembly. Unigenes that could not be aligned to any of the databases mentioned above were later predicted by ESTScan (Saldanha, 2004).

Differentially expressed gene analysis

The clean reads were mapped to unigenes using Bowtie2 (Langmead & Salzberg, 2012), and the expression levels of these high-quality unigene sets were computed based on RSEM values (Li & Dewey, 2011). Poisson distribution was applied to accurately describe gene expression variation on those with more than 2-fold changes as well as less than 0.001 false discovery rate [FDR] (Audic & Claverie, 1997). The performance of hierarchical clustering enabled the differentially expressed genes (DEGs) with similar features to be partitioned into distinct clusters and this analysis was done by having the aid of pheatmap function in R. When two or more relatively homogeneous groups were clustered, the intersection and union DEGs between them were performed. The DEGs were classified according to the GO and KEGG annotation results; FDR was calculated for each p-value. In general, the terms with FDR no more than 0.001 were defined as significantly enriched ones.

Identification of chemosensory genes

To identify candidate chemosensory unigenes for P. akamusi, an analysis using the tBLASTn modules was performed with reference to all the publicly available sequences of OBP, OR, CSP, GR, IR, and SNMP from Diptera species. Meanwhile, all the candidates were manually checked by using the BLASTx program.

RT-qPCR analysis

Using previously obtained cDNA from the three body parts of P. akamusi (the antennae, heads excluding antennae, and legs) as templates, OBPs, ORs, CSPs, GRs, IRs, and SNMPs were selected for RT-qPCR analysis. Specific primer pairs were designed by Primer5 according to the transcriptome data (Table S1). The RT-qPCR was conducted on the Roche LightCycler 480 detector (Stratagene, La Jolla, CA, USA) with the following cycling parameters: 94 °C for 30 s, 40 cycles of 94 °C for 5 s, 55 °C for 10 s, and 72 °C for 10 s. The dissociation curves were carried out as a post-qPCR analysis, during which all the components were firstly denatured for 5s at 95 °C, followed by cooling to 60 °C for 1 min. An increase to 95 °C for 30 s then took place, followed by cooling to 50 °C for 30 s. The relative gene expression data were analyzed by normalizing the threshold cycle (Ct) value of each sample to that of endogenous beta-tubulin, which was determined with 2^−ΔΔCt method (Livak & Schmittgen, 2001). The gene amplification from each tissue part of P. akamusi was actually completed in triplicates and the statistical significance was examined with one-way ANOVA test.

Phylogenetic analysis and classification of Chemosensory Genes

A set of chemosensory gene sequences from different Dipteran species were retrieved from the GenBank database, including Drosophila melanogaster, Aedes aegypti, Anopheles gambiae and Culex quinquefasciatus. To visually demonstrate the relationship between candidate chemosensory genes from P. akamusi and other dipteran species, they were first joined with alignments by Muscle 3.8.31 (Edgar, 2004) with default option and then manually refined by BioEdit v7.2.5 (Hall, 1999). The unrooted phylogenetic trees of each family were then generated through running MEGA 7.0 (Kumar, Stecher & Tamura, 2016) and using neighbor joining (NJ) method with Poisson correction, pairwise deletion, and rapid bootstraps (1,000 replicates).

Mate selection behavioral test

To examine whether the male individuals will attempt to mate with female ones, a female insect was placed in the experimental corridor of the Y-tube olfactometer while a male was placed at the open end of the central glass tube. The result showed that the majority of males were relatively active with approximately 58% choosing the female-occupied corridor and 38% choosing the empty side (Fig. 1A). However, this obvious response was not observed when females were given the choice to make a behavioral preference. It turned out that 82% female insects were recorded as non-responders to either corridor and the probability of females selectively orientating to the arm with male subjects was only 8% (Fig. 1B). These results revealed that males tended to seek companions and dominated the mating process probably because of the exposure to airborne pheromones released from the females. We therefore selected male individuals as the main subject for our further investigations and tried to explore how exactly the mutual attraction between the sexes occur during mating.

Transcriptome sequencing and de novo assembly

With the support of high-throughput sequencing platform, a set of raw sequencing reads were constructed from three different tissues of P. akamusi males, including 50.20 Mb from heads (excluding antennae), 51.82 Mb from antennae and 55.06 Mb from legs (Table 1). The trimming and elimination of adaptor-polluted and low-quality basis led to the generation of 4.42Gb, 4.46Gb, and 4.53 Gb clean reads from heads, legs, and antennae, respectively. The data also showed that the assessed intrinsic quality of the clean reads was sufficient for further analysis (Table 2) and the transcript length distribution was identified as well (Figs. S1, S2). In general, 27,609 unigenes were assembled with a total length, mean length, N50, and GC content of 28,312,956, 1,025 bp, 1,832 bp and 38.23%, respectively (Table 3).

Functional annotation and classification of the unigenes

The unigenes were annotated in accordance to seven functional databases and the result of this comparison illustrated that 17,311 unigenes (62.70%) were markedly matched with published proteins in NR database. Similarly, the total of annotated unigenes against known information from NT, Swiss-Prot, COG, KEGG, GO and InterPro was 7,936 (28.74%), 14,614 (52.93%), 7,986 (28.93%), 14,341 (51.94%), 7,501 (27.17%) and 14,752(53.43%), respectively (Table 4).

The homology searching against NR database revealed that the annotated sequences of P. akamusi were partly matched to sequences of Aedes aegypti (2,769 matching hits, 16.00%), followed by Aedes albopictus (1,870 matching hits, 10.80%), Culex quinquefasciatus (1,741 matching hits, 10.06%) and Anopheles gambiae str. (1,311 matching hits, 7.57%) (Fig. 2).

GO analysis, an enrichment tool, rendered 42,572 unigenes to be grouped into three functional categories, i.e., Biological Process (18,955 unigenes), Cellular Component (15,232 unigenes), and Molecular Function (8,385 unigenes). There were 25 subcategories in Biological Process with ‘cellular process’ accounting for 22.2% (4,213 unigenes) followed by ‘metabolic process’ making up 18% (3,411 unigenes). 19 sub-categories were included in the classification of Cellular Component and among which ‘groups of cell’ (3,145 unigenes, 20.6%) and ‘cell parts’ (3,121 unigenes, 20.5%) were the most abundant GO terms. For Molecular Function, sequences were predominately assigned to ‘catalytic activity’ (3,309 unigenes, 39.5%) and ‘binding’ (3, 267 unigenes, 39%) (Fig. 3A).

A more detailed comprehension of gene biochemical function could be further gained via KEGG analysis. The result delineated that a total of 25,377 unigenes were mapped and organized to six functional clusters, including Cellular Process (2,476 unigenes), Environmental Information Processing (2,613 unigenes), Genetic Information Processing (2,531 unigenes), Human Disease (6,293 unigenes), Metabolism (6,405 unigenes), and Organismal System (5,059 unigenes) (Fig. 3B).

Besides, potential functions of these putative unigenes were predicted with the assistance of COG database which was considered as a useful tool for understanding the orthologous relationships of gene products. The retrieved unigene sets were classified under 25 categories, among which the cluster of ‘General function prediction only’ (2,338 unigenes, 19%), ‘Translation, ribosomal structure and biogenesis’ (979 unigenes, 8%) and ‘Posttranslational modification, protein turnover, chaperones’ (941 unigenes, 7.7%) were the most representative three classifications (Fig. 3C). Venn diagram aimed at describing the similarities and differences of unigenes when searched against NR, COG, KEGG, Swissprot and Interpro databases and it turned out that 6,766 unigenes overlapped as the intersection proportion (Fig. 3D).

CDS prediction

A total of 17,522 CDSs were screened out from the 27,609 annotated unigenes. Since ESTScan program was capable of constructing CDSs in the remaining unigenes which were not of particular matches with the aforementioned databases, a total of 2,857 CDSs were predicted from these unannotated unigenes. Altogether, 20,379 CDSs were predicted from unigenes with a total length, a mean length, N50 and GC content of 16,667,634 bp, 817 bp, 1,290 bp and 42%, respectively (Fig. 4, Table 5).

Differentially expressed genes

Differentially expressed genes were examined by comparison among heads, legs, and antennae of male P. akamusi midges (Fig. 5A). The gene expression profiling analysis demonstrated that there were 2,281 up-regulated along with 4,865 down-regulated genes in heads against legs. A total of 3,421 were expressed at a markedly higher level while 5,792 genes were at lower levels in antennae vs. heads. Plus, 2,066 genes were expressed with greater numbers while 7,093 genes with relatively smaller numbers in antennae vs. legs. Shared and exclusively expressed genes among three tissues were shown in Fig. 5B. The number of genes differentially expressed among three tissue types was quite similar to each other with 21,103 in heads, 20,980 in legs, and 21,750 in antennae. Meanwhile, a subset of 5,620 genes that differentially expressed were identified with a tissue-dependent manner, including 2019 specifically expressed in head, 1540 genes expressed in legs and 2071 expressed in antennae.

Candidate odorant binding proteins

Based on functional annotation and tBLASTn results with an E-value of 1E-5 or lower (Table S2), 34 transcripts ranging from 138 to 927 bp were isolated as best candidate OBPs (PakaOBPs) in the P. akamusi transcriptome, and 29 of which contained full-length open reading frames (ORFs). The identified PakaOBP transcripts together with corresponding sequence data of OBPs from A. gambiae, C. quinquefasciatus, D. melanogaster and A. aegypti were used for obtaining phylogenetic inferences, as depicted in Fig. S3.

To determine the transcriptional output of these candidate genes among different tissues, 12 PakaOBPs with RPKM>1.2 were selected for RT-qPCR amplification and the expression profile of each gene differed. Five OBPs (PakaOBP1, 5, 8, 9 and 10) had noticeably higher proportion in antennae than either heads or legs. Meanwhile, seven other OBPs (PakaOBP1, 2, 3, 4, 6, 7, and 12) were expressed with a markedly greater numbers in legs. Of note, we found a unique antennae-specific expression pattern of PakaOBP9 since this particular gene showed almost no evidence of presence in other tissues (Fig. 6).

Candidate odorant receptors

Putative OR genes of P. akamusi (PakaORs) were represented based on their similarities to known insect ORs and tBLASTn results with an E-value of 1E-5 or lower generated a total of 17 PakaORs with length ranging from 312 to 1,671 bp (Table S3). Among them, 16 sequences were available as full-length coding ones. Evolutionary distances were evaluated among ORs from the sampled P. akamusi and four other Dipteran species, as shown in Fig. S4.

RT-qPCR was performed for 6 PakaORs that were relatively abundant in the antennal transcripts with RPKM>1.2. All of these PakaOR genes except for PakaOR4 were highly expressed in antennae compared with heads and legs whereas PakaOR3–6 had a relatively sufficient expression only in heads (Fig. 7). The results of PakaOR expression appeared to be partially consistent with those from the RPKM analysis.

Candidate gustatory receptors

A total of 32 candidate GR genes (PakaGRs) were identified based on the tBLASTn results with an E-value of 1E-5 or lower, with 20 sequences containing a full-length ORF (Table S4). A phylogenetic tree was constructed by using concatenated orthologous sequences derived from GRs in different Dipteran species (Fig. S5).

The relatively abundant transcripts (PakaGR1–10) according to the RPKM analysis were subjected to RT-qPCR detection, aiming to figure out the particular tissues where they preferred to be expressed. It turned out that PakaGR1 and PakaGR2 were mainly expressed in the antennae whereas PakaGR3 had a significant quantity in legs. Meanwhile, a pervasive expression pattern was shown for 10 other PakaGRs, including PakaGR2 and PakaGR3–10 (Fig. 8).

Candidates for other chemosensory genes

Our bioinformatic data facilitated the discovery of 22 transcripts for candidate IR genes (PakaIRs), 6 for candidate CSP genes (PakaCSPs), and 3 for candidate SNMP genes (PakaSNMPs) with the result of tBLASTn being E-value of 1E-5 or lower; of these genes, 19 PakaIRs, 2 PakaCSPs, and 3 PakaSNMPs included full length ORFs (Tables S5, S6). The phylogenetic signals of PakaIRs, PakaCSPs, and PakaSNMPs from P. akamusi with those from other four species were mapped onto the tree, as described in Fig. S8.

RT-qPCR data revealed that 16 PakaIRs were highly expressed in heads and legs (Fig. 9) while CSP1 was distributively expressed at varying levels in all tissues. Additionally, there was an abundance of PakaSNMP1 presented in antennae and legs at transcript level, which was a reliable clue for understanding its role in chemosensory process.