Integrative analysis reveals pathways associated with sex reversal in Cynoglossus semilaevis

Ethics statement

All procedures involving the handling and treatment of fish used in this study were conducted with the approval of the animal care and use committee of Chinse Academy of Fishery Sciences. All animal procedures were carried out according to the guide for the care and use of laboratory animals and the animal welfare in China.

Tissue materials and sex identification

In July 2018, 100 C. semilaevis of 120 days post-hatch (dph) were randomly selected from Tianjin Haisheng Aquaculture Company and their fins and gonads sampled. In our study, experimental fish are taken from the same pond, the temperature of which was 22 °C–23 °C, which will exclude all environmental impacts. Fins were stored in 100% ethanol and the DNA was extracted using the TIANamp marine animal DNA kit following the manufacturer’s protocol (TIANGEN^® Biotech Co., Ltd., Beijing, China). Gonads were frozen in liquid nitrogen then stored at −80 °C and RNA extracted using the TIANamp marine animal RNA kit following the manufacturer’s protocol (TIANGEN^® Biotech Co., Ltd., Beijing, China). The genetic sex was then determined with fin DNA using a previously described sex-specific marker (Liu et al., 2014). Phenotypic sex was determined by sequencing the Cyn_Z_6676874 and Cyn_Z_8564889 loci and verified by tissue sectioning of gonads (Figs. S1 and S2). The details of tissue sectioning are described in our previous study (Jiang & Li, 2017). Finally, gonad tissues from 4 pseudomales and 4 females were selected for transcriptomic and proteomic analysis.

RNA-Seq analysis

Gonad samples of 4 pseudomales and 4 females were subjected to RNA-Seq using an Illumina HiSeq Xten at the Novogene Biotechnology Co. Ltd (Beijing, China). For RNA-Seq, 3 µg total RNA from each sample was used to enrich mRNA for cDNA library construction. The cDNA library was generated using the NEB Next^® UltraTM RNA Prep Kit for Illumina^® (NEB, USA) following the manufacturer’s protocol. Raw reads were filtered by removing low-quality reads with ambiguous nucleotides and adapter sequences of “AGATCGGAAGAGC” using FastQC software. Clean data were then aligned using Hisat2 (https://ccb.jhu.edu/software/hisat2/index.shtml) and mapped to C. semilaevis reference genome (NCBI: GCF_000523025.1_Cse_v1.0_genomic). Sequences were assembled based on alignment results and used for further analysis. Gene expression levels were further normalized using the fragments per kilobase of transcript per million mapped reads (FPKM) method to eliminate the influence of differences in gene lengths and the amount of sequencing data on calculations of gene expression (Wang et al., 2018). Differentially expressed genes (DEGs) between females and pseudomales were identified with an absolute log2FC ≥ 2 and a false discovery rate-adjusted P (q value) <0.05 using the ballgown package (https://www.rdocumentation.org/packages/ballgown). The raw data of transcriptomes and proteomics obtained in this bioproject (PRJNA586726), have been deposited in Genbank-SRA under the accession numbers (SRR10492824 –SRR10492831).

qRT-PCR verification

To validate the results of high-throughput transcriptome sequencing, 16 DEGs were selected for real-time quantitative reverse transcription PCR (qRT-PCR). All genes used for quantification were selected based on factors such as PCR efficiency, gene expression level, and expression difference level. C. semilaevis β-actin gene was used as an internal reference (Li et al., 2010). Primers for the 16 DEGs were designed using Primer-blast based on the gene sequences from the RNA-Seq datasets (Table 1). Total RNA (1 µg) was reverse transcribed using the TOYOBO Reverse Transcription Kit then qRT-PCR conducted using SYBR^® Green Reagents with 15 µL reactions, each containing 9 µL SYBR^® Green, 0.3 µL forward primer, 0.3 µL reverse primer, 0.6 µL cDNA, and 4.8 µL ddH₂O. The following PCR amplification procedure was used: 95 °C for 2 min, followed by 40 cycles of 95 °C for 15 s, 60 °C for 30 s, and 72 °C for 30 s. All reactions were repeated 3 times (Livak & Schmittgen, 2001) and disassociation curve analysis performed using an ABI 7500 (Applied Biosystems, Foster City, CA, USA) fast real-time PCR system. Relative expression levels for genes of interest in females and pseudomales were then calculated with the 2^−△△Ct (Livak & Schmittgen, 2001) method and the log2 fold change obtained from qRT-PCR and RNA-Seq was compared.

Protein extraction and labeling

Frozen samples of pseudomale and female gonad tissues were ground in liquid nitrogen then the 100 mg of powdered sample transferred into a 1.5 mL tube containing 600 µL phenol extract and protease inhibitor. Samples were then ultrasonicated on ice then mixed with a phenol-tris-HCL saturated solution of equal volume for 30 min at 4 °C. After centrifugation at 15,000 × g for 15 min at 4 °C, the upper phenol layer was transferred to a new tube and proteins precipitated using 0.1 M ammonium acetate-methanol at −20 ° C. Proteins were quantified using the BCA protein concentration assay (Smith et al., 1985). Next, 100 µg of protein was reduced with 120 µL of reducing agent buffer (10 mM DTT, 8 M urea, 100 mM TEAB, pH 8.0) at 60 °C for 1 h. Reduced proteins were then alkylated for 40 min at room temperature in the dark by the addition of IAA to a final concentration of 50 mM. Proteins were then digested with 2 µL of 1 µg/µL trypsin per sample for 12 h at 37 °C (Wisniewski et al., 2009). Digested peptides were collected, mixed with 50 µL of 200 mM TEAB buffer, then centrifuged at 12,000 rpm for 20 min and the resulting lower layer collected and lyophilized. Samples were then resuspended in 200 mM triethylammonium bicarbonate (TEAB) and the peptides labeled using the TMT labeling reagent for 1.5 h at room temperature (Xiao et al., 2016; Guerreiro et al., 2014).

Proteomic analysis

After labeling, samples were mixed in a tube and dried with a vacuum pump. Pooled mixtures were fractionated using an Agilent 1100 HPLC with an ultraviolet (UV) absorption detector (210 nm and 280 nm) and a flow rate of 300 µL/min. Eluted peptides were grouped into 10 final fractions, desalted on an Agilent Zorbax Extend C18 column (2. 1 ×150 mm), and vacuum dried. LC–MS/MS analysis of each fraction was then performed on an EASY-nLC 1200 liquid chromatograph and peptides eluted using a gradient of 2–80% acetonitrile with 0.1% formic acid for 120 min on a Q-Exactive mass spectrometer (Thermo, USA). The MS scanning range was 300–1,600 m/z and the top 10 most intense signals were selected for further MS/MS analysis. The resolutions of the MS and the MS/MS scans were 70,000 and 17,500, respectively, and the dynamic exclusion set to 15 s. Experimental data was analyzed by Proteome Discoverer ™ 2.2 (Thermo, USA) software. The protein database used was Acanthomorphata uniprot and the false discovery rate of peptide identification was below 1%. After mass spectrometry analysis, results were retrieved from the database and relevant proteins separated. Differentially expressed proteins (DEPs) were identified based on fold changes (FC) calculated by comparing pseudomale and female groups and p-values were obtained by t-tests. Proteins that were significantly differentially expressed were defined by FC >1.5 or FC <2/3 and a p-value <0.05.

GO and KEGG analysis

Gene Ontology (GO) functional enrichment analysis was used to determine biological functions significantly associated with DEGs or DEPs, including relevant molecular functions, cell components, and biological processes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was used to aid biological interpretation of DEGs and DEPs. GO functional enrichment analysis for DEGs and DEPs was conducted using the online tool DAVID (https://david.ncifcrf.gov/) and KEGG pathway enrichment analysis performed with the KOBAS 3.0 online tool (http://kobas.cbi.pku.edu.cn/). Since C. semilaevis is a non-model species, genes were mapped to the zebrafish function database for analysis.

Correlations between the transcriptome and proteome

The Spearman correlation coefficient between mRNA and protein expression was calculated for genes differentially expressed in both transcriptome and proteome analysis to avoid masking effects of random variation. DEG sequences were extracted from RNA-Seq data and used as a nucleic acid database for alignment. Amino acid sequences of peptides from DEPs were aligned as queries to the DEGs database using the tblastn (https://blast.ncbi.nlm.nih.gov/Blast.cgi) command, then the genes differentially expressed at both the transcriptome and proteome levels extracted. An alignment result of identity ≥ 80%, e-value ≤ 1, and gap = 0 indicated a perfect matched between mRNA and protein sequences. Spearman correlation coefficients were calculated based on mean fold changes for mRNA and protein expression.

Characteristics of high-throughput sequencing data

A total of 194 Gb of sequence was generated using eight C. semilaevis gonadal tissue libraries. After quality control and adapter removal, 150.7 Gb of clean reads were obtained, 85.14%–89.43% of which were mapped to C. semilaevis genome. The assembled transcriptome consists of 46,317 transcripts belonging to 25,761 genes. A total of 2,172 (FDR < 1%) reliable proteins were identified by combining high-throughput TMT technology with high-resolution LC-MS/MS technology (Table 2).

Differentially expressed genes between females and pseudomales

We identified a total of 1,893 differentially expressed genes between females and pseudomales (Table S1). Of all the differentially expressed genes (DEGs), 785 genes were more highly expressed in females while 1,108 genes were more highly expressed in pseudomales (Fig. 1). These differentially expressed genes included several sex-related genes, such as dmrt1, dmrt3 and foxl2. The expression levels of dmrt1 and dmrt3 in pseudomales were 5.09 FC and 16.82 FC higher than in females, respectively. The expression level of foxl2 in females was 44.21 FC higher than in pseudomales.

GO enrichment analysis indicated significant enrichment of 49 GO terms (p < 0.05). Among them, 32 GO terms were identified from 785 genes highly expressed in females (Fig. 2) and 17 GO terms were identified from 1,108 genes highly expressed in pseudomales (Fig. 3). Of these 32 GO terms highly expressed in females, DNA replication initiation is most significant. In addition, biological processes including wnt signaling pathways, cell development, and oocyte maturation are enriched. Cellular component category extracellular space is enriched with the largest number of differential genes. For molecular functions classification, transforming growth factor beta receptor binding was most significantly enriched. Moreover, the top three categories were growth factor activity, methyltransferase activity, and cytokine activity. Among the 17 GO terms highly expressed in pseudomales, autophagy enrichment was the most significant and cytoplasm enriched with the largest number of differential genes.

KEGG analysis assigned the 1893 DEGs to 30 pathways, in which 19 pathways of highly expressed in females were enriched (Fig. 4), and 11 pathways that are highly expressed in pseudomales are enriched (Fig. 5). Of the 19 pathways, cell cycle was the most significantly enriched. In addition, progesterone-mediated oocyte maturation and oocyte meiosis closely related to female germ cells were enriched. Of the 11 pathways, adipocytokine signaling pathway is the most significant enriched.

The validation of RNA-Seq

Transcriptional regulation revealed by RNA-Seq data was confirmed in a biologically independent experiment using qRT-PCR. A total of 16 genes, including general transcription factor IIIA (gtf3a), somatomedin B and thrombospondin type 1 domain containing (sbspon), cysteine rich transmembrane module containing 1 (cystm1), POC1 centriolar protein A (poc1a), doublesex and mab-3 related transcription factor 1 (dmrt1), doublesex and mab-3 related transcription factor 3 (dmrt3), wingless-type MMTV integration site family, member 9B (wnt9b), forkhead box L2 (foxl2), zona pellucida sperm-binding protein 4-like (LOC103393374, ZP4), zona pellucida sperm-binding protein 4-like (LOC103396071, ZP4), charged multivesicular body protein 1b-like (LOC103391323), LOC103392682, cilia- and flagella-associated protein 157-like (LOC103394689), tubulin alpha-1C chain-like (LOC103391404), NAD(P)H dehydrogenase [quinone] 1-like (LOC103392072), and LOC103379620 were selected for qRT-PCR analysis. 16 DEGs including nine up-regulated genes were highly expressed in pseudomale (dmrt1, dmrt3, sbspon, cystm1, poc1a, LOC103391323, LOC103392682, LOC103394689, LOC103391404) and others were down-regulated in pseudomale. Except for the gene of dmrt3, other genes were more quantitative than RNA-Seq, but similar up-regulation or down-regulation patterns of these genes were observed in qRT-PCR and RNA-Seq results (Fig. 6). Further Pearson correlation analysis of qRT-PCR and RNA-Seq showed that their correlation was ρ = 0.98. This indicates that qRT-PCR can effectively validate the results of RNA-Seq.

Proteomic analysis

A comparative proteome survey was performed of females and pseudomales using the TMT technique to complement the transcriptome experiments. Out of 2172 reliable proteins, 184 up-regulated proteins and 140 down-regulated proteins were identified in pseudomales compared to females (Fig. 7 and Table S2). Enrichment analysis of DEP’s biological process showed that cytoplasmic mRNA processing body assembly, RNA processing and protein polymerization are the top three females with high expression of DEP enrichment (Fig. 8), while regulation of DNA-templated transcription, elongation is the high expression of DEPs in pseudomales enrich the most significant term (Fig. 9). In addition, U6 snRNP and nucleus are two terms in which female and pseudo-male highly expressed DEPs are enriched in cell component categories, respectively. RNA binding and nucleotide binding are the two main categories of high expression DEPs in females and pseudomales under molecular function, respectively.

KEGG pathway analysis of the female highly expressed proteins revealed that RNA degradation was the most significantly pathway, and the following is spliceosome and metabolic pathways (Fig. 10). However, KEGG pathway analysis of the pseudomale highly expressed proteins founded that protein processing in endoplasmic reticulum, PPAR signaling pathway and fatty acid biosynthesis were the top three pathways (Fig. 11).

Integrative analysis of transcriptomic and proteomic results

A total of 52 genes were differentially expressed between females and pseudomales at both the transcriptomic and proteomic levels. Further correlation analysis for these 52 genes revealed that ρ = 0.59 (p < 0.05) for the mRNA–protein correlations. Among these 52 related genes, 46 DEGs (88%) were well matched in their levels of protein abundance change. Additionally, compared to females, 33 genes were up-regulated in pseudomales, based on both transcriptomic and proteomic analysis, including YOD1 deubiquitinase (yod1), ubiquitin recognition factor in ER associated degradation 1 (ufd1), MINDY lysine 48 deubiquitinase 3 (mindy3), ubiquitin like domain containing CTD phosphatase 1 (ublcp1) and ADP ribosylation factor like GTPase 1 (arl1), importin 4 (ipo4), microtubule associated protein 1 light chain 3 alpha (map1lc3a). Transcriptomic and proteomic analysis both indicated the down-regulation of 13 genes, including Rho GDP dissociation inhibitor alpha (arhgdia), neural cell adhesion molecule 1-like (LOC103395384), eukaryotic translation initiation factor 4E nuclear import factor 1 (eif4enif1), aldehyde dehydrogenase 6 family member A1 (aldh6a1), and zygote arrest 1 (zar1). Six genes displayed decoupling patterns between the transcriptome and the proteome. Both LOC103376887 and heat shock 60 protein 1 (hspd1) were up-regulated in proteome analysis while the corresponding DEGs were down-regulated at transcript level. Creatine kinase M-type (LOC103377896) and sarcoplasmic/endoplasmic reticulum calcium ATPase 2-like (LOC103383309) were both down-regulated at the translation level but up-regulated at the transcriptomic level (Table 3). GO function analysis found that arl1, map1lc3a, and yod1 were enriched in cytosolic component while arl1, ipo4, hspd1, arhgdia, and zar1 were enriched in the cytoplasm component.