*2.1. 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.

#### *2.2. 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–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).

#### *2.3. 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.

#### *2.4. 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.

#### *2.5. 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).

#### *2.6. 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*, 5*ora*/*V1r* and 57 *olfC*/*V2r* genes.

## *2.7. 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.

#### *2.8. 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.
