Original Research ARTICLE
Transcriptome of Gonads From High Temperature Induced Sex Reversal During Sex Determination and Differentiation in Chinese Tongue Sole, Cynoglossus semilaevis
- 1Key Lab of Sustainable Development of Marine Fisheries, Ministry of Agriculture, Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Qingdao, China
- 2Laboratory for Marine Fisheries Science and Food Production Processes, Pilot National Laboratory for Marine Science and Technology (Qingdao), Qingdao, China
- 3College of Fisheries and Life Science, Shanghai Ocean University, Shanghai, China
- 4Laizhou Mingbo Aquatic Co., Ltd., Laizhou, China
- 5BGI-Shenzhen, Shenzhen, China
- 6Institut de Ciències del Mar (ICM), Consejo Superior de Investigaciones Científicas (CSIC), Barcelona, Spain
The sex of Chinese tongue sole (Cynoglossus semilaevis) is determined by both genetic sex determination (GSD) and environmental sex determination (ESD), making it an ideal model to study the relationship between sex-determination and temperature. In the present study, transcriptomes of undifferentiated gonads from genetic females and males, as well as differentiated gonads from males, females, and pseudomales under high and normal temperature treatments were generated for comparative transcriptomic analysis. A mean of 68.24 M high-quality clean reads was obtained for each library. Differentially expressed genes (DEGs) between different sexes and environmental treatments were identified, revealing that the heat shock protein gene family was involved in the high temperature induced sex reversal. The Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that were enriched in pseudomale and genetic female comparison included neuroactive ligand-receptor interaction, cortisol synthesis and secretion, and steroid hormone biosynthesis. Furthermore, weighted gene co-expression network analyses were conducted on all samples, and two modules were positive correlated with pseudomale under high temperature. An illustrated protein-protein interaction map of the module identified a hub gene, hsc70. These findings provide insights into the genetic network that is involved in sex determination and sexual differentiation, and improve our understanding of genes involved in sex reversal under high temperature.
Sexually dimorphic traits play a dominant role in evolution and behavior of animals, and genes regulating sex-determination and sexual development have long been of central interest in developmental biology. In vertebrates, sex determination mechanisms are broadly divided into two major categories: genetic sex determination (GSD) and environmental sex determination (ESD) (Marshall Graves, 2008). In GSD, the primary sex is determined during fertilization and regulated by heritable genetic components (Budd et al., 2015). In contrast, organisms possess ESD have undifferentiated gonads until they reach a sensitive period perceiving the environmental factors to determine the sex (Shao et al., 2014). Interestingly, some GSD species, mainly fish and reptiles, can change their primary sex without changing the genotype when exposed to environmental factors, a process known as environmental sex reversal (ESR) (Sarre et al., 2011; Shao et al., 2014; Holleley et al., 2015). Temperature is the most broadly studied factor that induces ESR. For example, Nile tilapia (Oreochromis niloticus) exposed to elevated temperature can override genetic sex determination, and the females can be reversed to produce pseudomales (Baroiller et al., 2009). Blue tilapia (Oreochromis aureus), Japanese flounder (Paralichthys olivaceus), European sea bass (Dicentrarchus labrax), turbot (Scophthalmus maximus) and Atlantic silverside (Menidia menidia) also undergo temperature induced ESR (Conover and Fleisher, 1986; Desprez and Mélard, 1998; Yamaguchi et al., 2007; Diaz and Piferrer, 2015; Robledo et al., 2015; Diaz and Piferrer, 2017).
Previous studies have found that epigenetic changes are involved in regulating ESR. In the European sea bass and Nile tilapia, juveniles exposed to high temperatures show increased methylation in the cyp19a1a promoter, which correlated to the masculinization of these individuals (Navarro-Martin et al., 2011; Wang et al., 2017). In red-eared slider turtle (Trachemys scripta elegans) (ESD), high temperature lead KDM6B, the histone H3 lysine 27 demethylase, to promote the transcription of the male sex-determining gene dmrt1. Therefore, a direct genetic link between epigenetic mechanisms and temperature-dependent sex determination has been shown in a turtle species (Ge et al., 2018).
Development of sexual traits is a complex process that involves a network of genes. Genes related to sex-determination (e.g. cyp19a1a and amh) and epigenetic regulators appear to be downstream factors in the ESR related pathways. Therefore, the identification of more upstream genes responding to the change in temperature is critical to further our understanding of ESR. Transcriptomic analyses can serve as a powerful tool to identify potential upstream regulators. The gonad transcriptome of Nile tilapia revealed numerous genes involved in temperature dependent ESR (Sun et al., 2018). In sea bass and pejerrey (Odontesthes bonariensis), transcriptomes of the developing gonads identified pathways involved in sex differentiation affected by temperature (Fernandino et al., 2011; Diaz and Piferrer, 2015). These reports were either in mature fish or did not test for fish that underwent ESR. Thus, genes involved in early development and temperature-induced ESR remains elusive.
Chinese tongue sole (Cynoglossus semilaevis) is an economically important marine flatfish with ZZ/ZW sex chromosome system (GSD) and the genome has been sequenced (Chen et al., 2014). Genotypic females of the Chinese tongue sole can be transformed into pseudomales when exposed to high temperatures during the crucial stage of sex differentiation, and epigenetic regulation is important during this process (Chen et al., 2014; Shao et al., 2014). Therefore, C. semilaevis is an ideal model to illuminate the relationship between sex-determination and temperature during ESR. Previous studies identified genes related to sexual dimorphism through transcriptome analysis of adult C. semilaevis (Guo et al., 2016; Wang et al., 2018; Liu et al., 2019). Gene expression patterns in the gonads of C. semilaevis before and after ESR remain unclear, partly due to the difficulty in distinguishing pseudomales at early stage. Excitedly, recent studies have uncovered sex determining gene dmrt1, and new single nucleotide polymorphism (SNP) loci associated with pseudomales, allowing for the detection of male, female, and pseudomale C. semilaevis (Cui et al., 2017; Jiang and Li, 2017; Cui et al., 2018).
In the present study, transcriptomic analysis was conducted on undifferentiated gonads from genetic females and males, plus differentiated gonads from males, females, and pseudomales under high temperature and normal temperature treatment. A large number of differentially expressed genes (DEGs) between female, male, and pseudomales under high and normal temperature were obtained. DEGs were functionally analyzed via GO term and KEGG pathway enrichment. The heat shock protein (hsp) gene family and several pathways, such as steroid hormone biosynthesis, cortisol synthesis and secretion, and neuroactive ligand-receptor interaction, were differentially expressed. Weighted gene co-expression network analysis (WGCNA) was employed to identify pseudomale-related genetic modules and hsc70 was identified as hub gene regulates the formation of pseudomales. These results provide new insights into the genes involved in sex determination and ESR in C. semilaevis.
Materials and Methods
Fish Culture and High Temperature Induction
C. semilaevis were kept in Laizhou Mingbo Aquatic Co., Ltd., in Yantai, China. Approximately 300 larvae from the same family were fed twice a day with commercial pellets and were reared in filtered seawater at 22°C. For the undifferentiated group (30 dpf_F and 30dpf_M), twenty larvae were sampled at 30-day post fertilization (dpf). The remaining larvae were divided into two equal groups and cultured in two 300 L tanks. The control group was cultured at 22°C (C22_F, C22_M, C22_P), and the high temperature group was cultured at 28°C (C28_F, C28_M, C28_P) until 3-month post fertilization (mpf) when sex differentiation was completed. Gonads were collected from the 30 dpf and 3 mpf C. semilaevis individuals for RNA extraction and immediately stored in liquid nitrogen. The remaining part of each fish was stored in ethanol to identify the genetic sex.
Genetic Sex Identification of C. Semilaevis
Genomic DNA was extracted from small pieces of tail fin using the phenol-chloroform method, and genetic sex was identified by polymerase chain reaction (PCR) against sex-specific simple sequence repeat (SSR) marker with the primers described previously (Liu et al., 2014). Pseudomales at 3 mpf were identified by PCR to detect the two SNP loci associated with sex reversal as previously described (Jiang and Li, 2017; Cui et al., 2018). The relative expression level of dmrt1 was measured to further confirm male, female, and pseudomale individuals by real-time quantitative reverse transcription PCR (qRT-PCR) (Cui et al., 2017).
Total RNA Extraction and RNA Sequencing
Total RNA was isolated from each gonad sample using the Trizol reagent (Invitrogen, USA) according to the manufacturer’s instructions. The qualities of the RNA samples were measured with the Agilent 2100RNA 6000 Nano kit (Agilent, USA). Three biological replicate samples for each group were generated. Thus, 24 samples of 30 dpf_F1-3, 30 dpf_M1-3, C22_F1-3, C22_M1-3, C22_P1-3, C28_F1-3, C28_M1-3, and C28_P1-3 were acquired. RNA sequencing was conducted with a BGISEQ-500 by PTM Biolabs (Hangzhou, China).
Transcriptome Quality Control and Gene Annotation
The sequencing reads were filtered using the SOAPnuke v1.4.0 (Chen et al., 2018) and trimmomatic (Bolger et al., 2014) software using default parameters to remove adapters, low quality reads with more than 5% unknown nucleotides, and reads with more than 20% low-quality bases (Q value <10). The filtered clean data were subsequently mapped to the C. semilaevis genome (NCBI Cse_v1.0) via HISAT2 v2.1.0 (Kim et al., 2015) and aligned with the reference transcript sequence using Bowtie2 v2.2.5 (Langmead and Salzberg, 2012) for removing rRNA sequences and annotating genes.
Analysis of Differentially Expressed Genes
Abundances of transcripts were quantified by RNA-Seq by Expectation Maximization (RSEM) (Li and Dewey, 2011), and gene expression level was further normalized using the Fragments Per Kilobase Million (FPKM) method to eliminate the influence of different gene lengths and amount of sequencing data on the calculation of gene expression. Differential expressed genes (DEGs) were identified using DEGseq and DEseq2 (Wang et al., 2010; Love et al., 2014). Transcripts with a fold change ≥2 and adjusted p (q value) < 0.001 were considered to be significant DEGs. DEGs were then subjected to an enrichment analysis for GO function and KEGG pathways using Phyper in R (Johnson et al, 1992), with p value < 0.05 as threshold. The cluster analysis of DEGs was processed using Pheatmap in R (https://CRAN.R-project.org/package=pheatmap).
Gene Co-Expression Network Analysis
To further understand the relationships between DEGs, a network analysis based on gene-to-gene correlations was performed using the R package WGCNA V1.48 (Langfelder and Horvath, 2008). Expression correlation coefficients of genes were calculated to search a suitable soft threshold to build gene networks using a scale-free topology model. Subsequently, gene modules with similar expression patterns were identified based on the resulting gene cluster dendrogram and using the dynamic tree cut method (GeneFrac Threshold 0.5). GO term and KEGG pathway enrichment analyses of the annotated genes were performed for each module and p value < 0.05 was defined as the threshold for significance. The protein–protein interaction network was predicted using DIAMOND (Buchfink et al., 2015) to map the genes to the STRING database, and score ≥300 was set as the threshold.
Validation of RNA-Seq Data by qRT-PCR
Nine candidate genes related to high temperature induced sex-reversal from the female, pseudomale, and male between the heat versus control temperature were validated with qRT-PCR. Primers (Table S1) were designed based on sequences from the National Center for Biotechnology Information (NCBI) database. β-actin gene was used as the internal control (Wang et al., 2018). One microgram total RNAs for high-throughput transcriptome sequencing was reverse transcribed into cDNA with the PrimeScript™ RT reagent Kit with gDNA Eraser (Takara, Japan). Then, qRT-PCR was performed using QuantiNova™ SYBR Green PCR Kit (Qiagen, Germany) in 20-µl reactions, containing 10 µl 2 × SYBR Green PCR Master Mix, 2 µl QN ROX Reference Dye, 0.7 µM forward primer, 0.7 µM reverse primer, and 1 µl cDNA. The cycling program was carried out at 95°C for 2 min, followed by 40 cycles of 95°C for 5 s and 60°C for 10 s; this was followed by a melting curve analysis in an ABI StepOnePlus Real-Time PCR system (Applied Biosystems, USA). Reactions were performed in triplicate. The relative expression fold changes of nine genes in female, pseudomale, and male gonads under normal and high temperature were analyzed using the 2−ΔΔCt method.
Phenotypic and Genetic Sex Identification of C. Semilaevis
The average weight of larvae sampled at 30 dpf was 0.080 ± 0.004 g, with the average length of 2.513 ± 0.041 cm. At 3 mpf, the average weight of normal temperature group was 1.312 ± 0.080 g, with the average length of 5.689 ± 0.136 cm. The average weight of high temperature group was 1.916 ± 0.127 g, and the average length of 6.911 ± 0.206 cm. The genetic sex of 30 dpf and 3 mpf female and male were distinguished using PCR of SSR marker. The genetic male had Z-specific derived band, and the genetic female had both W-specific and Z-specific derived bands. The genetic females that carried the T allele of Cyn_Z_6676874 and the A allele of Cyn_Z_8564889 may change into pseudomales. In addition, the relative expression of dmrt1 was high in males, low in females, and moderate in pseudomales (Figure S1).
Gonadal Transcriptome of C. Semilaevis
A total of 24 libraries, named C22_F 1-3, C22_M 1-3, C22_P 1-3, C28_F 1-3, C28_M 1-3, C28_P 1-3, 30dpf_F 1-3, 30dpf_M 1-3, were constructed using total RNA from gonads. A mean of 68.24 M filtered clean reads with a Q20 above 96.81% was obtained from each library (Table S2). Of the clean reads, 87.82–90.80% mapped to the C. semilaevis genome (NCBI Cse_v1.0) (Table S3), and a total of 25,022 known genes were identified in the transcriptome based on the C. semilaevis genome. The transcriptome data has been uploaded to NCBI Sequence Read Archive (SRA). The accession number was PRJNA576366.
Analysis of Differentially Expressed Genes
Using an adjusted p value of less than 0.001 and an absolute fold change greater than 2 as the threshold, over 1500 DEGs were obtained for each comparison, and the number of DEGs were highly variable, depending on the sex and environmental temperature treatment (Figure 1).
Figure 1 Number of differentially expressed genes. DEGs identified from each comparison groups. Red and blue indicate up-regulated genes and down-regulated genes, respectively.
GO and KEGG Enrichment Analyses of DEGs between Females and Males
To explore differences in gene expression between genetic females and males before sex differentiation, GO and KEGG enrichment analyses of DEGs between 30dpf_F and 30dpf_M was conducted. In the gonadal primordium region, 13 GO terms were significantly enriched with q < 0.05. The enriched GO terms were mainly related to muscle development, including muscle organ development, striated muscle cell development and differentiation (Figure S2A). Subsequently, KEGG enrichment analysis revealed that 17 pathways were significantly enriched in the gonadal primordium region (q < 0.05), including complement and coagulation cascades, calcium signaling pathway, and apelin signaling pathway (Figure S2B).
Genes that were differentially expressed in phenotypic females and males after sex differentiation were evaluated by identifying DEGs in C22_F vs C22_M. There were 10 GO significantly enriched (q < 0.05). Including DNA integration, microtubule-based process, DNA conformation change, nucleosome assembly, and protein-DNA complex assembly (Figure S2C). KEGG enrichment analysis identified 32 significantly enriched pathways in the differentiated gonads (p < 0.05). Among these, nine pathways were significantly enriched (q < 0.05), including pancreatic secretion, calcium signaling pathway, circadian entrainment, and oocyte meiosis (Figure S2D).
Gene expression changes in individuals exposed to high temperature conditions were identified by comparing C28_F vs C28_M. There were 221 significantly enriched GO terms (p < 0.05), including DNA binding, aminoglycan metabolic process, ion gated channel activity, DNA metabolic process, and temperature homeostasis (Figure S2E and Table S4). In addition, there were 20 enriched KEGG pathways (q < 0.05), including neuroactive ligand-receptor interaction, cell cycle, p53 signaling pathway, and ovarian steroidogenesis (Figure S2F). In addition, circadian entrainment, oocyte meiosis, and calcium signaling pathway were enriched, and these three pathways were also enriched in the in the C22_F vs C22_M comparison.
DEGs that were identified from the 30dpf_F vs 30dpf_M comparison were mainly enriched in genes involved in muscle development, therefore we focused on the remaining shared 943 DEGs that were identified from the C22_F vs C22_M and C28_F vs C28_M comparisons, which included sex-related genes, such as figla, foxl2, amh, and dmrt1. DEGs were further analyzed using GO annotations to identify their potential functions. In biological process, cellular process, single-organism process, and biological regulation were the most abundant GO function items, while binding, catalytic activity, and transporter activity were most abundant in molecular function. In cellular component, GO terms for cell, cell part, and membrane were enriched. In KEGG enrichment analysis, 18 KEGG pathways were significantly enriched (p < 0.05), including cell cycle, oocyte meiosis, hedgehog signaling pathway, calcium signaling pathway, neuroactive ligand-receptor interaction, and pancreatic secretion (Figure 2).
Figure 2 Analysis of DEGs between females and males. (A) Number of DEGs between females and males shown by Venn diagram. (B) GO classifications of 943 DEGs belong to C22_F vs C22_M and C28_F vs C28_M excluding 30dpf_F vs 30dpf_M. (C) Enriched KEGG pathways of the 943 DEGs (p < 0.05).
The GO and KEGG Enrichment Analysis of DEGs Between Fish Exposed to High and Normal Temperatures
To explore the gonadal-specific effects of temperature, DEGs were identified for females, males and pseudomales under high and normal temperature treatments. There were 534 DEGs shared in females, males and pseudomales when exposed to high temperature. DEGs were most enriched in the biological process functional items including cellular process, single-organism process, and metabolic process. Binding, catalytic activity, and transporter activity were most abundant in molecular function. In cellular component, cell, cell part, and membrane, were enriched. In KEGG enrichment analysis, 25 KEGG pathways were significantly enriched (p < 0.05), including nitrogen metabolism, glucagon signaling pathway, phototransduction, adipocytokine signaling pathway, FoxO signaling pathway, and hippo signaling pathway (Figure 3).
Figure 3 Analysis of DEGs in females, males, and pseudomales between high temperature and normal temperature. (A) Number of DEGs between high temperature and normal temperature shown by Venn diagram. (B) GO classifications of 534 DEGs in females, males and pseudomales between high temperature and normal temperature. (C) Enriched KEGG pathways of the 534 DEGs (p < 0.05).
Analysis of DEGs Between Pseudomale Formed Under High Temperature and Normal Temperature
Pseudomales develop from genetic females. To identify the changes in gene expression in pseudomales under high and normal temperature conditions, 7735 DEGs in both C22_P vs 30dpf_F and C28_P vs 30dpf_F comparisons were identified (Figure 4A). These DEGs were further screened by their relative expression pattern in C22_P and C28_P, and then divided into two categories (Figure 4B and Table S5). In the first category, heat shock protein family, including hspa13 (heat shock protein family A (Hsp70) member 13), hspa5 (heat shock protein family A (Hsp70) member 5), hsp90b1 (heat shock protein 90 beta family member 1), hspd1 (heat shock protein family D (Hsp60) member 1), and hspb9 (heat shock protein family B (small) member 9) were up-regulated in C28_P compared to C22_P. The male-related genes dmrt1 and sox6 were also up-regulated in C28_P compared to C22_P. Other genes that were up-regulated in the high temperature treated pseudomales included the vasopressin V2 receptor (LOC103384583), which is involved in signal transduction, the bag3 and BAG family molecular chaperone regulator 4-like (LOC103398065) which were co-chaperones for HSP70 and HSC70 chaperone proteins, and cyp11b-like, which is responsible for glucocorticoid and mineralocorticoid biosynthesis. KEGG enrichment analysis of this group of genes found 24 enriched pathways (q < 0.05). These pathways included steroid hormone biosynthesis, metabolism of xenobiotics by cytochrome P450, and carbon metabolism (Figure 4C). The second category included hsd17b3-like (testosterone 17-beta-dehydrogenase 3-like), fshr (follicle stimulating hormone receptor), aromatase-like, female-related foxl2 and wnt5b, and 11β-hydroxysteroid dehydrogenases hsd11b1-like, which were down-regulated in C28_P compared to C22_P. These DEGs had 30 significantly enriched KEGG pathways (q < 0.05), including neuroactive ligand-receptor interaction, cortisol synthesis and secretion, progesterone-mediated oocyte maturation, and estrogen signaling pathway (Figure 4D and Table S6)
Figure 4 Gene expression patterns in pseudomale formed under high temperature and normal temperature. (A) Number of DEGs shown by Venn diagram. (B) Heatmaps showing the expression of DEGs found both in C22_P vs 30dpf_F and C28_P vs 30dpf_F (red: higher expression, blue: lower expression). (C) Enriched KEGG pathways of DEGs up-regulated in C28_P compared to C22_P (q < 0.05). (D) Enriched KEGG pathways of DEGs down-regulated in C28_P compared to C22_P (q < 0.05), graph showing top 20 significant terms.
Weighted Gene Correlation Network Analysis and Functional Annotation of Pseudomale-Related Modules
To identify the genetic networks involved in the formation of pseudomales in C. semilaevis, a weighted gene coexpression network analysis (WGCNA) was conducted on all samples. Based on the expression values of each gene, we ultimately obtained six modules, and genes within each module had highly interconnected expression patterns (Figure 5A). According to the correlation coefficients between the module and the samples, we found that the blue and yellow modules showed the highest and significant correlation with the pseudomale gonad (Pearson value >0.9). We performed KEGG enrichment analysis on the two modules, and found that the blue module had 140 genes that were significantly enriched in pathways related to signal transduction, including cellular senescence, gap junction, and calcium signaling pathway. This suggests that genes in the blue module are likely involved in sensing external high temperature and activating signals (Figure 5B). In the yellow module, 83 genes were significantly enriched in 19 pathways (p < 0.05). The enriched KEGGs included Wnt signaling pathway, hippo signaling pathway, glycosphingolipid biosynthesis-globo and isoglobo series, adherens junction, and melanogenesis (Figure 5D). The co-expression network of the blue and yellow modules was further explored by protein-protein interaction analysis. The hub gene in the blue module was leucine-rich repeat-containing protein 18 (Figure 5C). The hub genes in the yellow module were heat shock cognate 70 kDa protein-like (XP_008307673.1), elongation factor 1-alpha 1 (XP_020457373.1), and eukaryotic initiation factor 4A-II-like (XP_008332345.1) (Figure 5E).
Figure 5 WGCNA and functional analysis of pseudomale-related modules. (A) The relationship between 6 modules across all samples by weighted gene coexpression network analysis (WGCNA). The color bar indicates correlation value from low (blue) to high (red). (B, D) Enriched KEGG pathways of genes in the blue module (B) and yellow module (D) (q < 0.05). (C, E) The network relationship in the blue module (C) and yellow module (E).
Validation of Transcriptome Sequencing by qRT-PCR
To validate the transcriptome sequencing data, the fold changes in 9 genes including figla, foxl2, wnt5b, hspa13, hsc70-like, hsd11b1-like, hsd17b3-like, aromatase-like, and cyp11b-like were measured using qRT-PCR. Similar up- or down-regulation patterns of these genes were observed in qRT-PCR and RNA-seq results (Figure S3), although few genes showed differentiation between the values. These differentiations may be due to the different calculated methods used between RNA-seq and qPCR. In general, the validation analysis by qRT-PCR indicated the accuracy and reliability of the RNA-seq results.
In the last decade, effort has been made to identify the molecular mechanism that underlies ESR in fish. High throughput sequencing, microarray profiles and genome-wide methylation analysis in zebrafish, European sea bass, and Nile tilapia have identified the role of epigenetic inheritance in ESR (Diaz and Piferrer, 2015; Sun et al., 2016; Ribas et al., 2017). A recent study in T. scripta found that epigenetic factors activate sex-related genes during male sex determination under high temperature (Ge et al., 2018). However, the relationship between high temperature and the epigenetic mechanisms, and the factors that convert external high temperature signals into internal signals to then alter the sex of genetic females to phenotypic males largely remains unclear.
To further identify genes involved in ESR, we exposed C. semilaevis to high temperature conditions during early sex differentiation from 30 dpf to 3 mpf. We focused on identifying gene expression changes in the gonad of males, females and pseudomales at the undifferentiated and post-differentiated stages using transcriptomic analysis. Before sex differentiation, genetic males and females showed differences in genes related to muscle development, suggesting that the sex-related dimorphism in C. semilaevis may be traced back to the pre-sex differentiation stage. After sex differentiation, females had up-regulation of female-specific genes figla and foxl2 and down-regulation of male-specific genes dmrt1 and amh, both in normal and high temperature conditions. DEGs between female and male were enriched in KEGG pathways cell cycle, oocyte meiosis, hedgehog signaling pathway, calcium signaling pathway, neuroactive ligand-receptor interaction, and pancreatic secretion. In the previous study of adult gonads of C. semilaevis, DEGs between sexes were also enriched in KEGG pathways included cell cycle, oocyte meiosis, and calcium signaling, whereas neuroactive ligand-receptor interaction and pancreatic secretion pathways were not enriched (Wang et al., 2018). Thus, these two pathways may be important during early sex differentiation. When comparing gonads of fish exposed to normal or high temperatures, 534 DEGs were found in all sexes of C. semilaevis, suggesting that there is likely a common genetic mechanism to cope with exposure to high temperature. These DEGs were enriched in the KEGG pathways nitrogen metabolism, glucagon signaling pathway, protein digestion and absorption, phototransduction, adipocytokine signaling pathway, and FoxO signaling pathway. Therefore, exposure to high temperature influenced metabolism, signal transduction, cell growth and proliferation. The results were consistent with previous studies in Oncorhynchus mykiss gairdneri and Acanthochromis polyacanthus (Narum and Campbell, 2015; Veilleux et al., 2015).
In order to understand the effect of high temperature on the formation of pseudomales during sex differentiation, we compared the gene expression profiles of 30 dpf genetic females and 3 mpf pseudomales under high or normal temperature treatment. After screening by their relative expression pattern in C22_P and C28_P, the DEGs were divided into two categories: (1) 304 DEGs that were up-regulated in high temperature compared to normal temperature, and (2) 188 DEGs that were down-regulated in high temperature compared to normal temperature. The first category included several members of the heat shock protein family, including hsp70 and hsp90 sub-types. Expressions of HSP70 and HSP90 are known to change in response to changes in water temperature in several species (Palmisano et al., 2000; Kazumi, 2002; Currie SMoyes and Tufts, 2010; Fan et al., 2014). In addition, hsp70 is associated with gonad development in tilapia, where it is highly expressed in the testis during gonad differentiation, along with sox9a, gata4, and lhx6 (Tao et al., 2018). Thus, Hsp70 is a candidate sensor for high temperatures during sex differentiation for future studies. The second gene category contained hsd17b3-like, which was down-regulated under high-temperature treatment. hsd17b3 is involved in androgen metabolism through oxidation-reduction reaction (Carre et al., 1993). Reduced expression of hsd17b3 may facilitate the accumulation of testosterone, thus promoting the formation of pseudomales. These DEGs of the two categories were enriched in several KEGG pathways, including neuroactive ligand-receptor interaction, cortisol synthesis and secretion, and steroid hormone biosynthesis. Within the neuroactive ligand-receptor interaction pathway, the vasopressin receptor gene [vasopressin V2 receptor (LOC103384583)] was significantly up-regulated in the pseudomale under high temperature, while fshr (follicle stimulating hormone receptor) was significantly down-regulated. Vasopressin is produced and released by the pituitary gland, and mediates peripheral and central physiological functions that are important for osmoregulation and reproduction. Receptors of vasopressin are up-regulated in the Thalassoma bifasciatum testis (Liu et al., 2015; Zimmermann-Peruzatto et al., 2015). fshr is expressed in the theca and granulosa cells of the ovary, and sertoli cells in the testis. It belongs to the G-protein-coupled membrane receptors superfamily and regulate gametogenesis and production of sexual steroids in teleost gonads (Shinoda et al., 2010). The neuroactive ligand-receptor interaction pathway was also enriched in Gymnocypris przewalskii brain under high temperature, where neuroactive ligand-receptor genes exhibited sex- and tissue-specific expression patterns (Tian et al., 2019). The expression of these receptors in the gonad is likely induced by neuropeptides through the hypothalamic-pituitary-gonadal (HPG) axis, which then switches on the expression of genes involved in steroid hormone synthesis.
The steroid hormone biosynthesis pathway plays a critical role in sex differentiation across vertebrates (Devlin and Nagahama, 2002). In addition to hsd17b3-like (LOC103398902), aromatase-like (LOC103378914) was also found in the steroid hormone biosynthesis pathway. In teleosts, aromatase activity is correlated with sex steroid levels, converting testosterone into estrogen and thus may play a role in homeostasis of the accumulated androgen.
Cortisol is a well-known hormone that is associated with stress and regulates many physiological processes. Cortisol levels are elevated in fish exposed to high temperature or high-density conditions in tilapia, zebrafish, and brook trout (Salvelinus fontinalis) (Delaney et al., 2005; Chadwick et al., 2015; Nadirah et al., 2017; Valdivieso et al., 2019). In addition, cortisol is related to masculinization in some fish species such as pejerrey, medaka and the Japanese flounder (Hattori et al., 2009; Hayashi et al., 2010; Yamaguchi et al., 2010). In our study, the cortisol synthesis and secretion pathway was enriched in pseudomales under the high-temperature treatment. cyp11b-like was up-regulated and hsd11b1-like was down-regulated in C28_P compared to C22_P. CYP11B can catalyze the conversion of testosterone to 11β-hydroxytestosterone (precursor of 11-ketotestosterone, 11-KT), a major androgen in fish that promotes spermatogenesis (Nadirah et al., 2017). The down-regulation of 11β-hydroxysteroid dehydrogenases hsd11b1-like can indicate a shift from cortisone-cortisol regeneration (Morgan et al., 2014). These expression patterns suggested that cortisol might cross-talk with steroidogenesis to mediate sex reversal under high temperature, which coincide with the study in pejerrey and medaka (Fernandino et al., 2012; Fermandino et al., 2013; Castañeda Cortés et al., 2019).
WGCNA analysis gave rise to two modules, although found in one biological replicate, correlated to pseudomale formation under high temperature. These two modules were enriched in KEGG pathways that were involved in signal transduction and cell developmental processes. In particular, protein-protein network analysis of the yellow module in WGCNA revealed hsc70-like (heat shock cognate 70 kDa protein-like) as the hub gene for pseudomale formation. HSC70 is a regulator of HSF1 (heat-shock factor 1), which plays an essential role in mediating the appropriate cellular response to physiological stresses (Ahn et al., 2005). In humans, hsc70 is also a mediator of the CaM-dependent nuclear import of male-specific SRY (Kaur et al., 2013). A previous study showed that heat shock protein 70 kDa (Hsp70/Hsc70) family proteins interact with steroid hormone receptors through BAG (Bcl-2 associated athanogene), and that Bag1L can enhance the transcriptional activity of the androgen receptor (Knee et al., 2001). In our transcriptome data, bag3 (the Bcl2 associated athanogene 3) and transcript LOC103398065 (BAG family molecular chaperone regulator 4-like) were up-regulated in C28_P compared to C22_P, indicating the linkage between of the heat shock protein family and steroid hormones.
Based on these findings, we propose a mechanism for how the heat signal causes sex reversal. Members of heat shock protein family may detect external high temperature. Then a series of processes including neuroactive ligand-receptor interaction pathway, steroid hormone biosynthesis pathway, and cortisol synthesis and secretion pathway were activated. This leads to changes in the relative amount of sex steroid hormones (androgen and estrogen), promoting the development of pseudomales in C. semilaevis.
The present study provides important insights into ESR induced pseudomale formation in C. semilaevis. Future studies incorporating detailed functional survey of these networks and hub genes will further improve our understanding of sex reversal under external high temperature in fish.
Data Availability Statement
The transcriptome data used in this study has been uploaded to NCBI Sequence Read Archive (SRA) with accession number PRJNA576366.
The collection and handling of C. semilaevis were approved by the Animal Care and Use Committee at the Chinese Academy of Fishery Sciences, and all the experimental procedures were performed in accordance with the guidelines for the Care and Use of Laboratory Animals at the Chinese Academy of Fishery Sciences.
QW analyzed the transcriptome sequencing data and wrote the manuscript. KL and WL reared the fish. BF and ZZ completed fish sampling and SNP locus sequencing. RW finished the qRT-PCR of dmrt1 gene. LT participated in the fish sampling. QL and FP revised the manuscript. CS was in charge of the guidance and revision of this manuscript.
This work was supported by National Nature Science Foundation of China (31722058, 31802275, and 31472269), National Key R&D Program of China (2018YFD0900301), the AoShan Talents Cultivation Program Supported by Qingdao National Laboratory for Marine Science and Technology (2017ASTCP-ES06), the Taishan Scholar Project Fund of Shandong of China to CS, the National Ten-Thousands Talents Special Support Program to CS, and the International Scientific Partnership Program ISPP at King Saud University (No. 0050).
Conflict of Interest
Author WL was employed by company Laizhou Mingbo Aquatic Co., Ltd. (Laizhou, China). The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.01128/full#supplementary-material
Table S1 | qRT-PCR primers used for validating RNA-seq data.
Table S2 | Reads filtered information of each library.
Table S3 | Reads mapping information from each library.
Table S4 | Significantly enriched GO term of DEGs in C28_F vs C28_M.
Table S5 | DEGs between pseudomale formed under high temperature and normal temperature.
Table S6 | Significantly enriched KEGG terms of DEGs down-regulated in C28_P vs C22_P.
Ahn, S. G., Kim, S. A., Yoon, J. H., Vacratsis, P. (2005). Heat-shock cognate 70 is required for the activation of heat-shock factor 1 in mammalian cells. Biochem. J. 392, 145–152. doi: 10.1042/BJ20050412
Baroiller, J. F., D'cotta, H., Bezault, E., Wessels, S., Hoerstgen-Schwark, G. (2009). Tilapia sex determination: Where temperature and genetics meet. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 153, 30–38. doi: 10.1016/j.cbpa.2008.11.018
Carre, J. L., Quemener, E., Amet, Y., Simon, B., Berthou, F., Mangin, P., et al. (1993). Characterization and solubilization of testosterone 17 beta-hydroxysteroid dehydrogenase in human hyperplastic prostate. J. Steroid Biochem. Mol. Biol. 46, 265–267. doi: 10.1016/0960-0760(93)90303-E
Castañeda Cortés, D. C., Arias Padilla, L. F., Langlois, V. S., Somoza, G. M., Fernandino, J. I. (2019). The central nervous system acts as a transducer of stress-induced masculinization through corticotropin-releasing hormone B. Development 146, 172866. doi: 10.1242/dev.172866
Chadwick, J. G., Nislow, K. H., Mccormick, S. D. (2015). Thermal onset of cellular and endocrine stress responses correspond to ecological limits in brook trout, an iconic cold-water fish. Conserv. Physiol. 3, 017. doi: 10.1093/conphys/cov017
Chen, S., Zhang, G., Shao, C., Huang, Q., Liu, G., Zhang, P., et al. (2014). Whole-genome sequence of a flatfish provides insights into ZW sex chromosome evolution and adaptation to a benthic lifestyle. Nat. Genet. 46, 253–260. doi: 10.1038/ng.2890
Chen, Y. X., Chen, Y. S., Shi, C. M., Huang, Z. B., Zhang, Y., Li, S. K., et al. (2018). SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience 7, 120. doi: 10.1093/gigascience/gix120
Cui, Y., Wang, W., Ma, L., Jie, J., Zhang, Y., Wang, H., et al. (2018). New locus reveals the genetic architecture of sex reversal in the Chinese tongue sole (Cynoglossus semilaevis). Heredity (Edinb). 121, 319–326. doi: 10.1038/s41437-018-0126-6
Cui, Z., Liu, Y., Wang, W., Wang, Q., Zhang, N., Lin, F., et al. (2017). Genome editing reveals dmrt1 as an essential male sex-determining gene in Chinese tongue sole (Cynoglossus semilaevis). Sci. Rep. 7, 42213. doi: 10.1038/srep42213
Currie Smoyes, C. D., Tufts, B. L. (2010). The effects of heat shock and acclimation temperature on hsp70 and hsp30 mRNA expression in rainbow trout: in vivo and in vitro comparisons. J. Fish Biol. 56, 398–408. doi: 10.1111/j.1095-8649.2000.tb02114.x
Devlin, R. H., Nagahama, Y. (2002). Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences. Aquaculture 208, 191–364. doi: 10.1016/S0044-8486(02)00057-1
Diaz, N., Piferrer, F. (2015). Lasting effects of early exposure to temperature on the gonadal transcriptome at the time of sex differentiation in the European sea bass, a fish with mixed genetic and environmental sex determination. BMC Genomics 16, 679. doi: 10.1186/s12864-015-1862-0
Diaz, N., Piferrer, F. (2017). Estrogen exposure overrides the masculinizing effect of elevated temperature by a downregulation of the key genes implicated in sexual differentiation in a fish with mixed genetic and environmental sex determination. BMC Genomics 18, 973. doi: 10.1186/s12864-017-4345-7
Fan, Z., You, F., Wang, L., Weng, S., Wu, Z., Hu, J., et al. (2014). Gonadal transcriptome analysis of male and female olive flounder (Paralichthys olivaceus). BioMed. Res. Int. 2014, 291067. doi: 10.1155/2014/291067
Fernandino, J. I., Hattori, R. S., Kishii, A., Strüssmann, C. A., Somoza, G. M. (2012). The cortisol and androgen pathways cross talk in high temperature-induced masculinization: The 11β-hydroxysteroid dehydrogenase as a key enzyme. Endocrinology 153, 6003–6011. doi: 10.1210/en.2012-1517
Fernandino, J. I., Hattori, R. S., Moreno Acosta, O. D., Strüssmann, C. A., Somoza, G. M. (2013). Environmental stress-induced testis differentiation: Androgen as a by-product of cortisol inactivation. Gen. Comp. Endocr. 192, 36–44. doi: 10.1016/j.ygcen.2013.05.024
Fernandino, J. I., Popesku, J. T., Paul-Prasanth, B., Xiong, H., Hattori, R. S., Oura, M., et al. (2011). Analysis of sexually dimorphic expression of genes at early gonadogenesis of pejerrey Odontesthes bonariensis using a heterologous microarray. Sex Dev. 5, 89–101. doi: 10.1159/000324423
Ge, C., Ye, J., Weber, C., Sun, W., Zhang, H., Zhou, Y., et al. (2018). The histone demethylase KDM6B regulates temperature-dependent sex determination in a turtle species. Science 360, 645–648. doi: 10.1126/science.aap8328
Guo, L., Wang, Y., Liang, S., Lin, G., Chen, S., Yang, G. (2016). Tissue-overlapping response of half-smooth tongue sole (Cynoglossus semilaevis) to thermostressing based on transcriptome profiles. Gene 586, 97–104. doi: 10.1016/j.gene.2016.04.020
Hattori, R. S., Fernandino, J. I., Kishii, A., Kimura, H., Kinno, T., Oura, M., et al. (2009). Cortisol-induced masculinization: does thermal stress affect gonadal fate in pejerrey, a teleost fish with temperature-dependent sex determination? PloS One 4, e6548. doi: 10.1371/journal.pone.0006548
Hayashi, Y., Kobira, H., Yamaguchi, T., Shiraishi, E., Yazawa, T., Hirai, T., et al. (2010). High temperature causes masculinization of genetically female medaka by elevation of cortisol. Mol. Reprod. Dev. 77, 679–686. doi: 10.1002/mrd.21203
Holleley, C. E., O'meally, D., Sarre, S. D., Marshall Graves, J. A., Ezaz, T., Matsubara, K., et al. (2015). Sex reversal triggers the rapid transition from genetic to temperature-dependent sex. Nature 523, 79. doi: 10.1038/nature14574
Kaur, G., Lieu, K. G., Jans, D. A. (2013). 70-kDa heat shock cognate protein hsc70 mediates calmodulin-dependent nuclear import of the sex-determining factor SRY. J. Biol. Chem. 288, 4148–4157. doi: 10.1074/jbc.M112.436741
Kazumi, N. (2002). The 70-kDa heat shock protein response in two intertidal sculpins, Oligocottus maculosus and O. snyderi : relationship of hsp70 and thermal tolerance. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 133, 79–94. doi: 10.1016/S1095-6433(02)00115-0
Knee, D. A., Froesch, B. A., Nuber, U., Takayama, S., Reed, J. C. (2001). Structure-function analysis of Bag1 proteins. Effects on androgen receptor transcriptional activity. J. Biol. Chem. 276, 12718–12724. doi: 10.1074/jbc.M010841200
Liu, H., Lamm, M. S., Rutherford, K., Black, M. A., Godwin, J. R., Gemmell, N. J. (2015). Large-scale transcriptome sequencing reveals novel expression patterns for key sex-related genes in a sex-changing fish. Biol. Sex Differ. 6, 26. doi: 10.1186/s13293-015-0044-8
Liu, J. X., Liu, X. B., Jin, C. F., Du, X. X., He, Y., Zhang, Q. Q. (2019). Transcriptome profiling insights the feature of sex reversal induced by high temperature in tongue sole Cynoglossus semilaevis. Front. Genet. 10, 522. doi: 10.3389/fgene.2019.00522
Liu, Y., Chen, S. L., Gao, F. T., Meng, L., Qiao-Mu, H., Song, W. T., et al. (2014). SCAR-transformation of sex-specific SSR marker and its application in half-smooth tongue sole (Cynoglossus semiliaevis). J. Agr. Biotechnol. 65, S22. doi: 10.3969/j.issn.1674-7968.2014.06.015
Morgan, S. A., Mccabe, E. L., Gathercole, L. L., Hassan-Smith, Z. K., Larner, D. P., Bujalska, I. J., et al. (2014). 11β-HSD1 is the major regulator of the tissue-specific effects of circulating glucocorticoid excess. Proc. Natl. Acad. Sci. U.S.A. 111, E2482. doi: 10.1073/pnas.1323681111
Nadirah, M., Ramly, H. R., Abdul Manaf, M. T., Razzak, L. A., Najiah, M. (2017). High water temperature impairs physiological responses in red hybrid tilapia: effects on cortisol and its regulation. AACL Bioflux 10, 1297–1308.
Navarro-Martin, L., Vinas, J., Ribas, L., Diaz, N., Gutierrez, A., Di Croce, L., et al. (2011). DNA methylation of the gonadal aromatase (cyp19a) promoter is involved in temperature-dependent sex ratio shifts in the European sea bass. PloS Genet. 7, e1002447. doi: 10.1371/journal.pgen.1002447
Palmisano, A. N., Winton, J. R., Dickhoff, W. W. (2000). Tissue-specific induction of hsp90 mRNA and plasma cortisol response in Chinook salmon following heat shock, seawater challenge, and handling challenge. Mar. Biotechnol. 2, 329–338. doi: 10.1007/s101260000005
Ribas, L., Liew, W. C., Diaz, N., Sreenivasan, R., Orban, L., Piferrer, F. (2017). Heat-induced masculinization in domesticated zebrafish is family-specific and yields a set of different gonadal transcriptomes. Proc. Natl. Acad. Sci. U.S.A. 114, E941–E950. doi: 10.1073/pnas.1609411114
Robledo, D., Ribas, L., Cal, R., Sanchez, L., Piferrer, F., Martinez, P., et al. (2015). Gene expression analysis at the onset of sex differentiation in turbot (Scophthalmus maximus). BMC Genomics 16, 973. doi: 10.1186/s12864-015-2142-8
Sarre, S. D., Ezaz, T., Georges, A. (2011). Transitions between sex-determining systems in reptiles and amphibians. Annu. Rev. Genomics Hum. Genet. 12, 391–406. doi: 10.1146/annurev-genom-082410-101518
Shinoda, T., Miranda, L. A., Okuma, K., Hattori, R. S., Fernandino, J. I., Yoshizaki, G., et al. (2010). Molecular cloning and expression analysis of Fshr and Lhr in relation to Fshb and Lhb subunits during the period of temperature-dependent sex determination in pejerrey Odontesthes bonariensis. Mol. Reprod. Dev. 77, 521–532. doi: 10.1002/mrd.21179
Sun, L. X., Teng, J., Zhao, Y., Li, N., Wang, H., Ji, X. S. (2018). Gonad transcriptome analysis of high-temperature-treated females and high-temperature-induced sex-reversed neomales in Nile tilapia. Int. J. Mol. Sci. 19, 689. doi: 10.3390/ijms19030689
Sun, L. X., Wang, Y. Y., Zhao, Y., Wang, H., Li, N., Ji, X. S. (2016). Global DNA methylation changes in Nile tilapia gonads during high temperature-induced masculinization. PloS One 11, e0158483. doi: 10.1371/journal.pone.0158483
Tao, W., Chen, J., Tan, D., Yang, J., Sun, L., Wei, J., et al. (2018). Transcriptome display during tilapia sex determination and differentiation as revealed by RNA-Seq analysis. BMC Genomics 19, 363. doi: 10.1186/s12864-018-4756-0
Tian, F., Liu, S., Shi, J., Qi, H., Zhao, K., Xie, B. (2019). Transcriptomic profiling reveals molecular regulation of seasonal reproduction in Tibetan highland fish, Gymnocypris przewalskii. BMC Genomics 20, 2. doi: 10.1186/s12864-018-5358-6
Valdivieso, A., Ribas, L., Piferrer, F. (2019). Ovarian transcriptomic signatures of zebrafish females resistant to different environmental perturbations. J. Exp. Zool. B Mol. Dev. Evol. 332, 55–68. doi: 10.1002/jez.b.22848
Veilleux, H. D., Ryu, T., Donelson, J. M., Van Herwerden, L., Seridi, L., Ghosheh, Y., et al. (2015). Molecular processes of transgenerational acclimation to a warming ocean. Nat. Clim. Change 5, 1074. doi: 10.1038/nclimate2724
Wang, L. K., Feng, Z. X., Wang, X., Wang, X. W., Zhang, X. G. (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. doi: 10.1093/bioinformatics/btp612
Wang, N., Wang, R., Wang, R., Chen, S. (2018). Transcriptomics analysis revealing candidate networks and genes for the body size sexual dimorphism of Chinese tongue sole (Cynoglossus semilaevis). Funct. Integr. Genomics 18, 327–339. doi: 10.1007/s10142-018-0595-y
Wang, Y. Y., Sun, L. X., Zhu, J. J., Zhao, Y., Wang, H., Liu, H. J., et al. (2017). Epigenetic control of cyp19a1a expression is critical for high temperature induced Nile tilapia masculinization. J. Therm. Biol. 69, 76–84. doi: 10.1016/j.jtherbio.2017.06.006
Yamaguchi, T., Yamaguchi, S., Hirai, T., Kitano, T. (2007). Follicle-stimulating hormone signaling and Foxl2 are involved in transcriptional regulation of aromatase gene during gonadal sex differentiation in Japanese flounder, Paralichthys olivaceus. Biochem. Biophys. Res. Commun. 359, 935–940. doi: 10.1016/j.bbrc.2007.05.208
Yamaguchi, T., Yoshinaga, N., Yazawa, T., Gen, K., Kitano, T. (2010). Cortisol is involved in temperature-dependent sex determination in the Japanese flounder. Endocrinology. 151, 3900–3908. doi: 10.1210/en.2010-0228
Keywords: transcriptome, high temperature, sex differentiation, environmental sex reversal, Cynoglossus semilaevis
Citation: Wang Q, Liu K, Feng B, Zhang Z, Wang R, Tang L, Li W, Li Q, Piferrer F and Shao C (2019) Transcriptome of Gonads From High Temperature Induced Sex Reversal During Sex Determination and Differentiation in Chinese Tongue Sole, Cynoglossus semilaevis. Front. Genet. 10:1128. doi: 10.3389/fgene.2019.01128
Received: 30 May 2019; Accepted: 17 October 2019;
Published: 22 November 2019.
Edited by:Peng Xu, Xiamen University, China
Reviewed by:Hongyu Ma, Shantou University, China
Jie Mei, Huazhong Agricultural University, China
Gustavo M. Somoza, CONICET Institute of Biotechnological Research (IIB-INTECH), Argentina
Copyright © 2019 Wang, Liu, Feng, Zhang, Wang, Tang, Li, Li, Piferrer and Shao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Changwei Shao, firstname.lastname@example.org