Transcriptome Analysis of Zebrafish Olfactory Epithelium Reveal Sexual Differences in Odorant Detection

Animals have evolved a large number of olfactory receptor genes in their genome to detect numerous odorants in their surrounding environments. However, we still know little about whether males and females possess the same abilities to sense odorants, especially in fish. In this study, we used deep RNA sequencing to examine the difference of transcriptome between male and female zebrafish olfactory epithelia. We found that the olfactory transcriptomes between males and females are highly similar. We also found evidence of some genes showing differential expression or alternative splicing, which may be associated with odorant-sensing between sexes. Most chemosensory receptor genes showed evidence of expression in the zebrafish olfactory epithelium, with a higher expression level in males than in females. Taken together, our results provide a comprehensive catalog of the genes mediating olfactory perception and pheromone-evoked behavior in fishes.


Introduction
Olfaction is a sense for detecting environmental odorants that plays essential roles in many aspects of animal activities, such as foraging, migration, prey avoidance and mating [1]. In most terrestrial vertebrates, there are two distinct chemosensory organs: the olfactory epithelium (OE) and the vomeronasal organ (VNO) [2]. In general, the OE is assumed to detect environmental odorants, while the VNO senses pheromones, although some exceptions were also reported in recent studies [3,4]. Different set of chemosensory receptors is thought to be expressed in each organ. For the OE, odorant receptors (ORs) and trace amine-associated receptors (TAARs) are known to be expressed, while for the VNO, vomeronasal receptors type 1 (V1Rs) and type 2 (V2Rs) are considered to be expressed [5]. Therefore, in terrestrial mammals, ORs and TAARs are suggested to detect "ordinary" odorant molecules, while V1Rs and V2Rs are employed to detect pheromones [6][7][8].
However, unlike the terrestrial counterparts, fish have only the OE to sense their olfactory environment, due to an absence of the VNO [9]. Without directly functional studies of the chemosensory receptors in fish, it is unclear whether the different types of chemosensory receptors in fish respond to different classes of odorants-and which family of the chemosensory receptors is used for recognizing odorants or pheromones. Interestingly, there are also some studies using electrophysiological experiments that provide indirect evidence of putative ligands for chemosensory receptors [10][11][12][13][14][15][16].
For example, members of ORs from goldfish can recognize F-prostaglandins, while both ORs and V2Rs can perceive amino acids [17,18].
Several studies have found that sex-specific behaviors can be due to sex differences in responses to external stimuli, including courtship songs, colors and chemosensory cues [19]. Indeed, chemosensory cues in mice, such as pheromone molecules, have been shown to be detected specifically by their vomeronasal receptors [20]. Therefore, sex-specific behaviors maybe initiated by sex-specific pheromones [21][22][23] or by sex-differential expression of chemosensory receptors in response to the same pheromones [24]-or even by sex-specific alternative splicing [25]. For instance, Darcin (the MUP20 encoded by Mup20)-a male-specific pheromone in mouse urine-only attracts females to affect their memory [24]. ESP1 (exocrine-secreted peptide 1) is detected by both male and female mice, but only stimulated female-specific mating behaviors, suggesting the existence of sex-specific neuronal circuits [23,26]. Male fruit flies require the protein encoded by the fruitless (fru) gene to complete courtship, which is produced in different sex-specific isoforms via alternative splicing [25]. However, to date, few studies have reported whether sexual differences in odorant-sensing exists between males and females, especially in fish.
To determine the full repertoire of chemosensory receptors expressed in the zebrafish olfactory epithelium and evaluate the extent of sexual differentiation in odorant-sensing between sexes, we used RNAseq to profile their transcriptomes in male and female zebrafish. We found that a very high percentage of chemosensory receptors are indeed expressed in the zebrafish olfactory epithelium, with a higher expression level in males than females. However, the olfactory transcriptomes between males and females are highly similar, with limited genes showing differentially expressed or alternatively splicing, which may be associated with odorant-sensing between sexes. Collectively, our results provide a comprehensive catalog of the genes mediating olfactory perception and pheromone-evoked behavior in fish.

Sample Collection
All experiments in this research were approved by the Institutional Animal Care and Use Committee of Institute of Hydrobiology, Chinese Academy of Sciences (Approval ID: Y21304501). Adult wild-type zebrafish from AB background were maintained in the zebrafish facilities in the China Zebrafish Resource Center (CZRC)for one week to familiarize with the laboratory environment. Zebrafish were raised together in a mixed sex population. Olfactory epithelium from each individual was dissected out and frozen in liquid nitrogen quickly. Due to their small size, each of the samples was a pool of mRNA from three individuals of the same sex. Three independent biologic replicates for both male and female samples were prepared.

Library Construction and High-Throughput Sequencing
Total RNA from each of the six samples was extracted using the SV Total RNA Isolation System (Promega). We assessed the RNA quality using agarose gel electrophoresis and measured RNA integrity using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Sequencing libraries preparation and high throughput sequencing were generated by Novogene (Beijing, China) following our previous study [27][28][29]. Briefly, mRNA was purified from total RNA with poly-T oligo-attached magnetic beads and fragmented into short pieces. Then, the first-strand cDNA was synthesized using random hexamer primer and second-strand cDNA was then generated. Finally, the paired-end cDNA library was prepared according to the Illumina's protocols and sequenced on Illumina HiSeq 4000 platform (150 bp paired-end) (Illumina, San Diego, CA, USA). The RNA-seq reads were deposited into the National Center for Biotechnology Information (NCBI) Sequence Read Archive database (Accession No. SRP154651, Supplementary Materials Table S1).

Analysis of RNA-Seq Data
Raw RNA-seq reads were first filtered to delete primer dimers and low-quality bases (Phred quality score lower 20) using the Trim Galore! program (version 0.3.7) (https://www.bioinformatics. babraham.ac.uk/projects/trim_galore/). We only retained paired-end reads from which either end was longer than 50 bp after trimming for subsequent analyses. High quality paired-end reads from each sample were aligned to the transcripts from zebrafish genome annotated by Ensembl (release 97) [30] using Bowtie2 [31];abundances of transcripts(FPKM, Fragments Per Kilobase Million)were estimated using RSEM program (v 1.3.1) [32]. We filtered the genes with FPKM > 1 in at least half of the six samples as transcriptionally active genes for subsequent analyses. The raw read counts for each transcript estimated by RSEM were extracted and then normalized using TMM method to control for differences in sequencing depth among samples, and the differentially expressed transcripts were identified using the edgeR package [33] using a minimal fold change of 2 and an adjusted p value cutoff of 0.05. Full lists of differential expression genes can be found in Supplementary Materials Table S2.

Gene Ontology Analysis
Overrepresentation of the gene ontology (GO) terms for upregulated genes between male and female zebrafish olfactory epithelium were identified using Gorilla (http://cbl-gorilla.cs.technion.ac. il/) [34], which allows detection functional overrepresentation in a candidate data set against a list of background genes. We set the false discovery rate (FDR) of 0.001 as our cutoff value and conducted separately for each sex.

Characterization of Alternative Splicing Events (ASEs)
ASEs are divided into five broad categories including skipped exon (SE), alternative 5 splice site (A5SS), alternative 3 splice site (A3SS), mutually exclusive exons (MXE) and retained intron (RI) [35]. We used rMATS [36] to detect and count reads that correspond to each of the five types of ASEs. rMATS can identify these ASEs events from a GTF file of annotated transcripts and count the number of reads that correspond to each of the five events described. To identify differential alternative splicing (DAS) events between male and female zebrafish olfactory epithelium, an FDR-adjusted p value less than 0.05 was used as threshold for DAS events (Supplementary Materials Table S6).

Data Mining for the Chemosensory Receptor Repertoire
We extracted all the sequences from annotated and automatically predicted paralogs for or, taar, ora/V1r and olfC/V2r genes from the Ensembl zebrafish genome (GRCz11, release 97). We only considered a gene as a putative chemosensory receptor gene for a given family by checking the candidates position within each chemosensory receptor family clade in a phylogenetic analysis. By using this method, we obtained a total of 170 or, 126 taar, 5ora/V1r and 57 olfC/V2r genes.

Phylogenetic Analysis
All of the coding sequences for the chemosensory receptor genes were obtained from Ensembl. The coding sequences for each of the 4 chemosensory receptor gene families were translated into protein sequences, aligned with the program MUSCLE [37] and then back-reversed to their coding sequences alignment. The ML trees were reconstructed by RAxML (version 8.1.17) [38] under the GTRGAMMMAI substitution model with bootstrap support values determined using 1000 replicates.

Quantitative Real-Time PCR (qRT-PCR)
In order to confirm the differentially expressed genes detected by RNA-seq, we further employed quantitative real-time PCR (qRT-PCR) on a subset of genes that among the significantly DEGs between the male and female zebrafish olfactory epithelium. Primers of these genes were designed using the NCBI primer designing tool (http://www.ncbi.nlm.nih.gov/tools/primer). We synthesized the first strand cDNA from 500 ng of total RNA samples using M-MLV Reverse Transcriptase (Promega, Madison, WI, USA) and diluted 1:10 as amplification template. qRT-PCR was performed in a 10-µL volume using the LightCycler ® 480 SYBR Green I Master on a LightCycler ® 480 II Instrument (Basel, Roche, Switzerland). Thermocycling conditions were 95 • C for 5 min, followed by 45 cycles of 95 • C for 20 s and 58 • C for 25 s, and a melting curve analysis was performed to confirm the primer specificity after amplification. The relative gene expression between male and female zebrafish olfactory epithelium was determined using the comparative CT method [39] and the fold change values were the mean of six biologic replicates from each group.

Transcriptome Data
A total of six cDNA libraries were constructed in this study. In total, we obtained 23.75, 24.44, 25.65 million reads for male olfactory epithelium and 22.96, 21.70, 22.28 million reads for female olfactory epithelium, respectively. From these, 75.46-77.75% of the reads can be mapped to the annotated regions in zebrafish genome, indicating the good quality of reads (Supplementary Materials Table S1). Overall, we detected a total of 18,658 genes with FPKM > 1 in at least half of the six samples, which were defined as robustly expressed genes. The expression estimates for all annotated transcripts in each replicate were provided in Supplementary Materials Table S2.
We first evaluated the variation in gene expression among the three biologic replicate samples for each sex. We found that the correlation values were highly significant between them all ( Figure 1, Pearson correlation coefficients of at least 0.95, p-value < 2.2 × 10 −16 ). Only very small sets of genes are unusually variable among replicates ( Figure 1). In general, the OE transcriptomes from the male and the female zebrafish were highly correlated (Figure 2A, Pearson correlation coefficient, r = 0.89, p < 2.2 × 10 −16 ). We therefore averaged the FPKM values for each gene across each sex. In the olfactory epithelium from both male and female, a few gene are extremely highly expressed. For example, in male 312 most abundant genes account for almost 50% of the fragments and in female 333 most abundant genes account for almost 50% of the fragments obtained from the tissue (Supplementary  Materials Table S3). Moreover, a total of 295 genes were found to be shared between male and female most abundant genes. These results were generally consistent with the patterns found in mouse olfactory transcriptomes [40].

Sex Dimorphism in Zebrafish Olfactory Epithelium Expression Profiles
To assess whether transcriptional differences in the zebrafish olfactory epithelium can account for sex-specific responses to olfactory cues, we examined their sexually dimorphic gene expression profiles. We found that the overall transcriptomes are highly similar between male and female zebrafish olfactory epithelia (Figure 2A,B). Briefly, we detected a total of 713 transcripts showing higher expression level in male olfactory epithelium and 605 transcripts showing higher expression level in female olfactory epithelium, respectively ( Figure 2C). However, only 68 and 53 transcripts remained significantly differential expression between male and female olfactory epithelium after applying a false discovery rate (FDR) threshold of 0.05 ( Figure 2D and Supplementary Materials Table S4). Among these, some genes identified to be differentially expressed between sexes are expected to be involved in response to olfactory cues (Table 1), such as gene OR132-2 (odorant receptor, family H, subfamily 132, member 2), OTX1 (orthodenticle homeobox 1) and OR115-10 (odorant receptor, family F, subfamily 115, member 10). However, whether these differentially expressed olfactory receptor genes can detect sex relative pheromones needs to be verified in further functional experiments.

Sex Dimorphism in Zebrafish Olfactory Epithelium Expression Profiles
To assess whether transcriptional differences in the zebrafish olfactory epithelium can account for sex-specific responses to olfactory cues, we examined their sexually dimorphic gene expression profiles. We found that the overall transcriptomes are highly similar between male and female zebrafish olfactory epithelia (Figure 2A,B). Briefly, we detected a total of 713 transcripts showing higher expression level in male olfactory epithelium and 605 transcripts showing higher expression level in female olfactory epithelium, respectively ( Figure 2C). However, only 68 and 53 transcripts remained significantly differential expression between male and female olfactory epithelium after applying a false discovery rate (FDR) threshold of 0.05 (Figure 2Dand Supplementary Materials Table S4). Among these, some genes identified to be differentially expressed between sexes are expected to be involved in response to olfactory cues (Table 1), such as gene OR132-2 (odorant receptor, family H, subfamily 132, member 2), OTX1 (orthodenticle homeobox 1) and OR115-10 (odorant receptor, family F, subfamily 115, member 10). However, whether these differentially  To gain further insights into the biologic processes that differed between male and female zebrafish olfactory epithelium, we performed a GO enrichment analysis with the upregulated genes detected in male and female zebrafish olfactory epithelium (Supplementary Materials Table S5). We found that gene ontology biologic process terms associated with response to external stimulus were highly overrepresented among the differentially expressed genes (DEGs) with male-biased expression (Figure 3 and Supplementary Materials Table S5), including response to other organism (GO:0051707), chemotaxis (GO:0006935) and taxis (GO:0042330). However, a significantly different pattern was found in the enriched GO terms for female olfactory epithelium upregulated genes. For example, female olfactory epithelium upregulated genes were enriched in ion import, such as potassium ion import (GO:0010107), inward rectifier potassium channel activity (GO:0005242) and ligand-gated channel activity (GO:0022834). Taken together, these results indicated that there are significant differences in the upregulated genes between male and female zebrafish olfactory epithelium.
Genes 2020, 11, x FOR PEER REVIEW 8 of 16 Figure 3. Gene ontology (GO) terms associated with differentially expressed genes between male and female zebrafish olfactory epithelium.

Alternative Splicing between the Two Sexes of Zebrafish Olfactory Epithelium
To examine the influence of sex on alternative splicing regulation in zebrafish olfactory epithelium, alternative splicing events (ASEs) including skipped exon (SE), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exons (MXE) and retained intron (RI) ( Figure 4A) between male and female zebrafish olfactory epithelium were identified using rMATS [36].
Roughly 46% of genes have 2 detectable isoforms in zebrafish genome ( Figure 4B). Among them, we detected a total of 15,193 ASEs, which were distributed in 7585 genes. The most abundant ASEs were skipped exon, accounting for 76.4% of all ASEs, followed by retained intron (7.0%), alternative 3′ splice site (6.5%), mutually exclusive exons (6.3%) and alternative 5′ splice site (3.9%) ( Figure 4C). Using p value cutoff of 0.05, differential ASEs were identified between the two sexes of zebrafish (Supplementary Materials Table S6). A total of 86 significantly differential ASEs were identified between the male and female zebrafish olfactory epithelium, which included 50 SE, 4 A5SS, 17 A3SS, 6 MXE, and, 9 RI (Table 2). For example, maptb, which is expressed in the developing central nervous system [41], was found to generate different isoforms by sex. stxbp5l, which is involved in syntaxin binding, also showed differential isoforms between sexes, which was identified as sexually dimorphic gene between cattle and rat species [42].

Alternative Splicing between the Two Sexes of Zebrafish Olfactory Epithelium
To examine the influence of sex on alternative splicing regulation in zebrafish olfactory epithelium, alternative splicing events (ASEs) including skipped exon (SE), alternative 5 splice site (A5SS), alternative 3 splice site (A3SS), mutually exclusive exons (MXE) and retained intron (RI) ( Figure 4A) between male and female zebrafish olfactory epithelium were identified using rMATS [36].   Roughly 46% of genes have 2 detectable isoforms in zebrafish genome ( Figure 4B). Among them, we detected a total of 15,193 ASEs, which were distributed in 7585 genes. The most abundant ASEs were skipped exon, accounting for 76.4% of all ASEs, followed by retained intron (7.0%), alternative 3 splice site (6.5%), mutually exclusive exons (6.3%) and alternative 5 splice site (3.9%) ( Figure 4C). Using p value cutoff of 0.05, differential ASEs were identified between the two sexes of zebrafish (Supplementary Materials Table S6). A total of 86 significantly differential ASEs were identified between the male and female zebrafish olfactory epithelium, which included 50 SE, 4 A5SS, 17 A3SS, 6 MXE, and, 9 RI (Table 2). For example, maptb, which is expressed in the developing central nervous system [41], was found to generate different isoforms by sex. stxbp5l, which is involved in syntaxin binding, also showed differential isoforms between sexes, which was identified as sexually dimorphic gene between cattle and rat species [42].

The Chemosensory Receptor Repertoires
As in mouse olfactory sensory neurons (OSNs), each of the OSNs in mature zebrafish randomly expressed only one olfactory receptor [43][44][45]. Thus, the expression level of any given olfactory receptor within the olfactory epithelium will be low. Consistent with this prediction, among the 170 or, 126 taar, 5ora/V1r and 57 olfC/V2r genes in the zebrafish genome, we detected 147 out of 170 or, 104 out of 126 taar, 4 out of 5ora/V1r and 57 out of 57 olfC/V2r genes expressed in either male or female olfactory epithelium with FPKM > 1, confirming that the chemosensory receptors are mainly expressed in the olfactory epithelium (Supplementary Materials Table S7). However, almost all (90%) of these chemosensory receptor genes were found to be very lowly expressed in the zebrafish olfactory epithelium transcriptomes (FPKM < 15), which is consistent with the results from mice [46].
To further assess whether there are differences in chemosensory receptor gene expression levels between the two sexes in zebrafish, we compared the expression levels of chemosensory receptor genes between male and female zebrafish olfactory epithelium transcriptome. From this analysis, we made the following interesting observations. First, the relative receptor abundance level vary greatly between each member, suggesting that each member of the chemosensory receptor genes was expressed asymmetrically. Some members have high expression levels, whereas other members expressed with very low levels ( Figure 5). Second, almost all of the chemosensory receptor genes expressed with no differences between male and female zebrafish olfactory epithelium, with only six or genes showing differentially expressed between the two sexes (Supplementary Materials Table S7). Third, all of the 6 or genes ( Table 1) that showed sexual preference were expressed at a higher level in the male olfactory epithelium than in females, which is generally consistent with the results from mice [40,46]. In order to verify the expression results from RNA-seq, we further employed qRT-PCRto measure relative mRNA levels for a total of 12 candidate genes. Our results demonstrated that the expression patterns for these genes were highly consistent between RNA-seq and qRT-PCR ( Figure 5E), suggesting the reliability of our RNA-seq datasets.

Discussion
During the last few decades, chemosensory receptor gene families have received extensive attention across various organisms, including humans [45,[47][48][49], mice [40,46,50], fish [28,[51][52][53] and insects [21,[54][55][56][57]. However, it remains still unclear whether male and female individual possesses the same ability for odorant-sensing. In the present study, we reported the transcriptional profiles of the olfactory epithelium from the male and female zebrafish obtained by RNA-seq. By comparing the expression levels of genes from three different biologic replicate samples for each sex, we obtained highly correlated samples, which suggested that the subsequent differential expression analyses were reliable.
Sexual dimorphisms in behaviors have been widely observed in the olfactory systems in both mice [50,58] and flies [59]. Although previous studies have suggested that sexually dimorphic behaviors between sexes could be influenced by sensory input such as olfactory cues [60], it is unclear whether differences in gene expression in olfactory epithelia can underlie sexually dimorphic behaviors in fish [40]. Our results showed that the overall transcriptional profiles of zebrafish olfactory epithelium between sexes are highly similar. However, we indeed detected a few genes differentially expressed between male and female zebrafish, whose roles may be associated with odorant-sensing. Moreover, the functional enrichment analyses of the differentially expressed genes were also involved in response to external biotic stimulus. Therefore, all these results from our study suggested that differences in gene expression in the olfactory epithelium between sexes in zebrafish may play a role in odorant sensing to some extent.
However, when we focused on the chemosensory receptor genes that are directly binding to the odorants in surrounding environment, our results suggest that the olfactory receptor genes may not all be involved in behavioral differences between males and females. Although almost all of the chemosensory receptor genes can be detected to be expressed in the olfactory epithelium of zebrafish, most them displayed similar expression levels between sexes, with only six or genes showing differentially expressed between the two sexes. Considering the existence of nearly 358 chemosensory receptor genes in zebrafish genome, the six differentially expressed or genes should play minor roles in the dimorphic behaviors between male and female zebrafish [61,62]. Therefore, the sexually dimorphic behavioral responses to odorants in zebrafish are unlikely to be solely accounted for by transcriptional differences at the level of detection [40].
Interestingly, several functional studies have identified several olfactory receptor ligands in fish [10][11][12][13][14][15][16]. For example, Yabukiet al. reported that or114-1 and or114-2 in the group β or genes are the key olfactory receptors for the sex pheromone prostaglandin F 2α to mediate male courtship behavior in zebrafish [16]. However, none of these olfactory receptor genes were found to be significantly differentially expressed in the olfactory epithelium between male and female in our RNA-seq results. It should be noted that differential expression of peripheral odorant receptors is not the only way to generate sex-specific behavior, as higher neuronal circuits could be different between males and females. Therefore, whether the six differentially expressed or genes can play crucial roles in sexually dimorphic behaviors is still uncertain and need to be confirmed by further functional studies in future.

Conclusions
In this study, we performed a transcriptomic analysis on the olfactory epithelia of both male and female zebrafish using high-throughput RNA sequencing. We found the olfactory transcriptomes between males and females are highly similar, with only few genes displaying differentially expressed. Most of chemosensory receptor genes showed evidence of expression in the zebrafish olfactory epithelia, with a higher level of expression in males than in females. Collectively, these results provide a comprehensive catalog of the genes mediating olfactory perception and pheromone-evoked behavior in fish.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/11/6/592/s1, Table S1: Summary of the RNA-seq datasets in olfactory mucosae from male and female zebrafish, Table S2: Expression estimates in the zebrafish olfactory epithelium. A dataset containing the raw and normalized expression values for all the genes in the zebrafish olfactory epithelium, Table S3: 312 and 333 most abundant genes account  for 50% of the fragments obtained from the tissue for male and female, Table S4: Differential expression analysis between male and female samples. Genes with an FDR < 0.05 were considered significant, Table S5: Functional enrichment analysis for differentially expressed genes, Table S6: Differential alternative splicing events between male and female zebrafish olfactory epithelium, Table S7 Expression level of the chemosensory receptor genes.