ORIGINAL RESEARCH article

Front. Genet., 29 May 2019

Sec. Livestock Genomics

Volume 10 - 2019 | https://doi.org/10.3389/fgene.2019.00522

Transcriptome Profiling Insights the Feature of Sex Reversal Induced by High Temperature in Tongue Sole Cynoglossus semilaevis

  • 1. Key Laboratory of Marine Genetics and Breeding, Ministry of Education, Ocean University of China, Qingdao, China

  • 2. Laboratory for Marine Fisheries Science and Food Production Processes, Qingdao National Laboratory for Marine Science and Technology, Qingdao, China

Abstract

Sex reversal induced by temperature change is a common feature in fish. Usually, the sex ratio shift occurs when temperature deviates too much from normal during embryogenesis or sex differentiation stages. Despite decades of work, the mechanism of how temperature functions during early development and sex reversal remains mysterious. In this study, we used Chinese tongue sole as a model to identify features from gonad transcriptomic and epigenetic mechanisms involved in temperature induced masculinization. Some of genetic females reversed to pseudomales after high temperature treatment which caused the sex ratio imbalance. RNA-seq data showed that the expression profiles of females and males were significantly different, and set of genes showed sexually dimorphic expression. The general transcriptomic feature of pesudomales was similar with males, but the genes involved in spermatogenesis and energy metabolism were differentially expressed. In gonads, the methylation level of cyp19a1a promoter was higher in females than in males and pseudomales. Furthermore, high-temperature treatment increased the cyp19a1a promoter methylation levels of females. We observed a significant negative correlation between methylation levels and expression of cyp19ala. In vitro study showed that CpG within the cAMP response element (CRE) of the cyp19a1a promoter was hypermethylated, and DNA methylation decreased the basal and forskolin-induced activities of cyp19a1a promoter. These results suggested that epigenetic change, i.e., DNA methylation, which regulate the expression of cyp19a1a might be the mechanism for the temperature induced masculinization in tongue sole. It may be a common mechanism in teleost that can be induced sex reversal by temperature.

Introduction

The types of sex determination are diversified in teleost. Three main types of primary sex determination have been described in gonochoristic species: genotypic sex determination (GSD), temperature-dependent sex determination (TSD) and a combination of both (GSD ++ TSD) (; Yamamoto et al., 2014). The sex of fish had strong uncertainty in the development process. In addition to genetic information, environmental factors could influence the sex determination, such as temperature. Apart from fish, the temperature irreversibly determining gonadal sex has been well established in reptiles and amphibians (; ). Since firstly described in Menidia menidia (), this phenomenon had been widely observed in fish, which showed that sex ratio would become unbalanced if the fish experienced high temperature during thermosensitive period (TSP). The imbalance of sex ratio was caused by sex reversal. Usually, it can be divided into three types: (1) high temperature has positive correlation with the proportion of males. (2) High temperature induces females, and has a negative correlation with males. (3) Both low temperature and high temperature increase the proportion of males (; ; ).

To verify the molecular mechanism of temperature effects (TE), a series of exploration was carried out. Steroid hormone, glucocorticoid, and epigenetic modification have been reported to be related to sex reversal and played critical roles during sex differentiation in TSD (; ; ; Yoshinaga et al., 2010; ; , ; ; ; Zhang et al., 2013). Besides, it was discovered that intron retention of JARID2 and JMJD3 genes in Pogona vitticeps could mediate sex-reversed females (). Androgen-to-estrogen ratio determined whether an undifferentiated gonad differentiated into a testis or ovary in non-mammalian vertebrates (). The regulation of steroid ratio depended on the activity of gonadal aromatase, the product of cyp19a1a, which converts androgens into estrogens irreversibly (). In reptiles, up-regulating or down-regulating cyp19a1a could alter gonad phenotype. The expression level of gonadal cyp19a1a was associated with TSD in Trachemys scripta and Alligator mississippiensis (; ). In teleost, it has been confirmed that high temperature induced masculinization is related to the methylation level of cyp19a1a promoter in Dicentrarchus labrax. Methylation modification in the promoter region could suppress the binding of transcription factors to the corresponding sites (SF-1, FOXL2, and CREB) resulting in the change of expression (; Zhang et al., 2013). Similar conclusions were observed in Oreochromis niloticus and Oncorhynchus mykiss (Valdivia et al., 2014; Wang et al., 2017). Meanwhile, FOXL2 and SOX9, which showed dimorphic DNA methylation patterning were also considered as the candidate genes in A. mississippiensis and Paralichthys olivaceus (; ). Other factors have also been suggested to play a role in GSD + TE, such as heat shock proteins (HSPs), transient receptor potential channels (TRPs), cold inducible RNA binding proteins (CIRBPs), and microRNAs (; ; ; ; ).

The effect of temperature on the sex differentiation can be profound and far-reaching, and needs comprehensive studies to fully understand the molecular mechanisms. Chinese tongue sole, Cynoglossus semilaevis, is a GSD + TSD sex determination teleost with ZZ/ZW sex chromosomes (Zhou et al., 2005). Female-specific DNA sequences had been identified in C. semilaevis, which could be used for distinguishing genetic female and male (Wang et al., 2009, 2013). Therefore, C. semilaevis is a unique powerful model to explore molecular events associated with GSD + TSD. In previous study, it was reported that epigenetic modification was involved in sex reversal of C. semilaevis by BS-seq and RNA-seq, and transgenerational epigenetic inheritance was observed in offspring generated by sex reversal individuals (; ). We aimed to filter genes related to sex differentiation, explore the relationship of expression level and methylation modification, and analyze whether methylation could regulate the binding of transcription factor. In this study, the genetic female individuals that inversed to phenotypic male individuals are defined as pseudomales. These pseudomales are distinguished from high temperature treatment groups using female-specific markers. RNA-seq was performed on the gonads of females, males, and pseudomales. The whole expression profiles were investigated, and candidate genes involved in sexual gonad development were identified. The methylation patterns of the putative genes were also analyzed. The interaction of upstream regulatory sequence and the corresponding transcription factors was verified by dual-luciferase reporter system. These findings helped us to understand the genetic epigenetic programing driving vertebrate GSD + TE and provide insight for future investigations aimed at clarifying the mechanisms controlling sex differentiation and sex reversal.

Materials and Methods

Fish Rearing and Temperature Treatment

Fish and embryos were collected from Yellow Sea Aquatic Product Co. Ltd., Shandong, China. Embryos were incubated at 20°C, the natural temperature for C. semilaevis spawning, fertilization and hatching. For this study, a batch of embryos collected from three pairs of parents was used. After hatching the fry were reared at ambient temperature (20–22°C). The juveniles at 25 days post fertilization (dpf) with total length (TL) of 13 ± 2 mm were separated into two groups. One group (n = 3000) was reared at ambient temperature throughout the TSP as control group (low temperature group, LT). The other group (n = 3000) was exposed at 28°C during the entire TSP and as the high-temperature group (high temperature group, HT). The temperature was increased to 28°C at a rate of 0.5°C/day, and then maintained for 100 days, until 125 dpf (Figure 1). Then the water was recovered to ambient temperature to follow the natural fluctuations until the end of the study, when the fish were 300 days old. The proportion of phenotypic males and females was counted by gonad biopsy and section confirmation in LT group and HT group, respectively. From these phenotypic males genetic males and pseudomales were identified using female-specific markers (Wang et al., 2009). The survival rate was also calculated for both of the groups.

FIGURE 1

Sample Collection and Gonadal Histology

At 300 dpf, fish were sacrificed and gonadal samples were collected. For each fish, one gonad was processed for histological identification of phenotypic gender and DNA/RNA extraction. Gonads were fixed in 4% PFA in PBS, embedded in paraffin, cut at 7 μm thickness and stained with haematoxylin-eosin. Meanwhile, the other gonad was snap-frozen in liquid nitrogen and stored at -80°C for RNA-seq analysis. Muscle tissues were collected to extract DNA for individual sexing and methylation analysis. The methylation level of muscle was selected as control.

RNA Isolation, cDNA Library Construction, and Illumina Sequencing

Gonads from nine individuals including three biological replicates of females (FO), males (MT), and pseudomales (PMT) were selected for RNA-seq analysis. Total RNA was extracted using Trizol Reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s protocol, treated with RNase-free DNase I (TaKaRa, Dalian, China) to degrade genomic DNA, and then frozen at -80°C. RNAclean Kit was applied to remove proteins. The quality and quantity were evaluated via 1.5% agarose gel electrophoresis and spectrophotometry using NanoPhotometer Pearl (Implen GmbH, Munich, Germany) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States).

The nine RNA-seq libraries were constructed with Illumina TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, United States) in accordance with the manufacturer’s instruction. Then the libraries were subjected to paired-end sequencing of 125 bp on the Illumina HiSeq 2000.

Data Processing and Bioinformatics Analysis

Raw reads generated from the Illumina sequencing platform were cleaned by removing adaptors and low quality sequences using FastQC. The cleaned reads of each sample were mapped to the reference genome () by TopHat with default parameters (). Then the mapping files were analyzed using Cufflinks to assemble the reads into transcripts for each dataset (). Complete transcripts were obtained by merging the assemblies of all datasets using Cuffmerge.

Identification of Differentially Expressed Genes and Functional Enrichment Analysis

All the expressed genes were aligned to databases for homology annotation, including non-redundant protein databases (NR), Swiss-Prot, Gene Ontology (GO), eukaryotic Orthologs Groups (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) by BlastX with e-value of 1e-5 ().

FPKM were used to select the DEGs. The FPKM was calculated by Cuffdiff (Trapnell et al., 2012). To identify the differentially expressed genes (DEGs) among female, male and pseudomale gonads, we set the following standards: genes with an adjusted log2FoldChange ≥ 2 or log2FoldChange ≤-2, and P < 0.01 were considered as DEGs. The DEGs were then enriched by GO terms and KEGG categories using DAVID (). The visualization of global similarities and differences of expression profiles of all individuals was accomplished by principle component analysis (PCA), MA plot and heatmap. These analysis were completed with R package.

qRT-PCR Validation

A total of ten DEGs (Sox9, GATA4, Dmrt1, AMH, HSD11b2, cyp19a1a, esr1, topaz1, GATA6, Sox3) were selected for qRT-PCR validation. Specific primer pairs were designed by IDT. qRT-PCR was performed in a 20 μg solution containing 10 ng of template cDNA and SYBR qPCR SuperMix (Novoprotein, Shanghai, China) by using LightCycler 480 (Roche, Forrentrasse, Switzerland) at 95°C for 5 min pre-incubation, followed by 45 cycles of 95°C for 15 s and 60°C for 45 s. The relative quantities of the target genes expressed as fold variation over GAPDH were calculated using the 2-ΔΔCt comparative Ct method. qRT-PCR data were statistically analyzed using one-way ANOVA followed by LSD test using SPSS 20.0. P < 0.05 indicated statistical significance.

Methylation Levels Measured by Bisulfite-Mediated Genomic Sequencing

Methylation sites were prediction and BSP primers design in promoter were performed by Methprimer. Gonad and muscle tissues of females, males and pesudomales (six individuals each) were used to extract genomic DNA. The DNA samples from the same tissue of the same gender were mixed. The mixed DNA was modified using the EZ DNA Methylation-Gold Kit (ZYMO Research). The primers M-cyp19a1a-Fw1/Rv1 and M-cyp19a1a-Fw2/Rv2 were used for methylation-specific PCR. Eight positive clones were sequenced for each group. Site-specific methylation measurements were analyzed using BiQ-Analyzer.

Cyp19a1a-Luc Reporter Vector Construct and in vitro Methylation

A pGL3-Cyp19a1a-Luc reporter vector was constructed by inserting the cyp19ala promoter fragment into the pGL3-basic vector (Promega, Madison, WI, United States) between SacI and XhoI sites. The promoter was a 1969 bp fragment amplified from genomic DNA with primers pGL3-cyp19a1a-Fw/Rv (Supplementary Table S3). The pGL3-cyp19a1a promoter vector was cytosine-methylated using M. SssI methylase (Thermo Fisher Scientific, MA, United States) (M-cyp19a1a-Luc) according to the manufacturer’s instructions. It could methylate all cytosine residues within the double-stranded dinucleotide recognition sequence. The methylation status of the vector was checked by HhaI, which only digested methylated DNA.

Transfection and Luciferase Reporter Gene Assay

The HEK 293T cell line was used for transfection with unmethylated and methylated plasmids. Before the experiment, a total of 5 × 105 cells were seeded into 24-well plates and cultured for 24 h. Then the plasmids were transfected into HEK 293T cells by LipofectamineTM 3000 Transfection Reagent (Thermo Fisher Scientific, MA, United States) according to the manufacturer’s instructions. At 48 h after transfection, cells were washed with PBS and analyzed for Luc activity using the luciferase assay system (Promega, Madison, WI, United States). Forskolin (5 μM), the activation of cAMP, which binds to CREB site, was added 10 h before the end of cell culture.

Results

Sexual Ratio Changes After High Temperature Treatment

The survival rates and proportion of females and males were counted both in LT group and HT group after treatment. No difference was found between the two groups in survival rate (χ2 = 0.190, P = 0.663). The survival rates were 62.97 and 59.80% in LT group and HT group (Figure 2E). The proportion of females and males as detected by biopsy and gonad tissue sectioning was 56.58 and 43.53% in LT group (Figure 2AC), and 36.51 and 63.49% in HT group (Figure 2D), respectively. The proportion of males was significantly increased more than 20% after treatment with high temperature during TSP (χ2 = 7.624, P = 0.006). These data indicated that masculinization was induced in genetic females following high temperature treatment.

FIGURE 2

RNA Sequencing

The genotypic and phenotypic sex of these individuals were distinguished by molecular marker and tissue section (Figure 2F). A total nine cDNA libraries were sequenced on the Illumina platform, generating 655,677,682 raw reads, encompassing about 30 Gb of sequence. Valid ratio and GC content of each cDNA library were shown in Table 1. Approximately 80.1% of reads exhibited significant hits to the genome. The transcriptome data obtained from the samples has been uploaded to NCBI SRA site, with accession numbers of PRJNA480118 (SAMN09628942, SAMN09628943, SAMN09628989, SAMN09628990, SAMN09628991, SAMN09628992, SAMN09628993, SAMN09628994, and SAMN09628995).

Table 1

SampleRaw readsClean readsErrorQ20Q30GC
Female-1-119757381136152300.0396.8193.4248.92
Female-1-219757381136152300.0395.5691.4248.92
Female-2-119500137132496590.0396.6793.1249.79
Female-2-219500137132496590.0395.8391.8849.79
Female-3-116918995163224680.0396.7091.6148.96
Female-3-216918995163224680.0396.7291.6148.95
Male-1-114492850138532000.0395.5692.0346.69
Male-1-214492850138532000.0395.6192.0046.78
Male-2-115806041121938140.0396.8893.8047.18
Male-2-215806041121938140.0395.7391.9747.14
Male-3-117961443134987870.0396.5493.1146.25
Male-3-217961443134987870.0395.4591.3046.33
Pseudomales-1-119970045130338180.0397.0594.0847.09
Pseudomales-1-219970045130338180.0395.4091.4346.98
Pseudomales-2-121633228127943380.0396.4192.7648.13
Pseudomales-2-221633228127943380.0494.9590.5048.02
Pseudomales-3-118792243125072020.0396.9994.0046.91
Pseudomales-3-218792243125072020.0394.9890.7646.80

Summary statistics of gonad transcriptome sequencing data.

Differential Expression and Functional Enrichment Analysis

Principle component analysis analysis were conducted to detect the global similarities and differences expression profiles among FO, MT, and PMT. It displayed that ovary (FO) replicates clustered closely in a region, and testis (MT and PMT) replicates clustered into another region. The MT and PMT replicates clustered together (Figure 3B). These results demonstrated that the expression patterns of phenotypic females and males was significantly different. However, the expression profiles showed more similarity between males and pseudomales. Although females and pseudomales retained the same genotype, the expression profiles were quite different. Males and pseudomales possessed different sex chromosomes, but the expression patterns were similar.

FIGURE 3

Among the DEGs, 5851 genes were significantly differentially expressed in FO vs. MT. 5611 genes were found differentially expressed in FO vs. PMT. Between MT vs. PMT, only 426 genes were identified as DEGs (Figure 3C and Supplementary Table S1). Regarding the functions of the DEGs, a large number of genes related to gonad development and sex differentiation were identified, which include Dmrt1, Dmrt3, HSD3b1, AMH, HSD3b7, esr1, SOX9, GATA4, GATA6, cyp19a, AMHR2 (Table 1). The heatmap of hierarchical clustering of DEGs was generated to visualize the expression patterns. The profile of phenotypic female was obviously different with all phenotypic male. The expression pattern of pseudomale was prone to that of male (Figure 3A).

After filtration, the DEGs were applied to perform GO analysis and KEGG enrichment. All the DEGs were mapped to GO terms and compared with the background of the whole transcriptome. They were significantly enriched in several GO terms in biological process, cellular component and molecular function (Supplementary Table S2). The results of enrichment were as follow: (1) In DEGs of FO vs. MT, the terms related to sexual differentiation and the regulation of reproduction were enriched, including sperm motility, 3-beta-hydroxy-delta5-steroid dehydrogenase activity and steroid hormone receptor activity. Besides, the terms about immune response were enriched (Figure 4A). (2) In FO vs. PMT, the terms of steroid hormone and helicase activity were detected, such as steroid hormone receptor activity and helicase activity (Figure 4B). They were also involved in reproduction and sexual differentiation and development. (3) In MT vs. PMT, it was found that some terms related to reproduction and the generation and development of testis was detected, comprising of male gamete generation, spermatogenesis, spermatid development, spermatid differentiation, and sterol transport (Figure 4C). Interestingly, the terms about sperm generation and differentiation were detected, including male gamete generation, spermatogenesis and spermatid differentiation and development. Surprisingly, the GO terms about energy metabolism were enriched, including UTP metabolic process, CTP metabolic process, CTP biosynthetic process, GTP metabolic process (Figure 4C and Supplementary Table S2). These terms are involved in meiosis and gamete generation, and may influence sperm activity. Meanwhile, KEGG pathway enrichment analysis was performed. A total of 44 KEGG terms were significantly enriched. The enriched signal pathways were similar in FO vs. MT and FO vs. PMT, including ribosome biogenesis, cell adhesion and metabolism and biosynthesis. Only one signal pathway involved in metabolism was enriched, phosphatidylinositol signaling system (Figure 5).

FIGURE 4

FIGURE 5

Identification of Genes Involved in Sexual Differentiation and Gonad Development

To identify genes involved in reproduction, including gonad development, gametogenesis and steroid biosynthesis in C. semilaevis, three strategies were used. (1) Sex-related genes were retrieved from the enriched GO terms, related to reproduction and steroid. (2) The DEGs were filtered by a set of key words that had been reported in other teleost, including gonad, sex, oocyte, meiosis, steroid, reproduction, and morphogenesis (; ; ). (3) Part genes were chosen from sex-related KEGG pathways. In accordance with these strategies, a set of potential candidate genes were obtained, and qRT-PCR validation were conducted (Table 2). Additionally, the DEGs of MT vs. PMT were analyzed independently, and the genes involved in spermatogenesis, gamete generation and development and energy metabolism were selected (Table 3).

Table 2

GeneFPKMP-valueAnnotation

FemaleMale
FO vs. MT
LOC10338907204.795.00E-05Spermatid perinuclear RNA-binding protein-like, partial
spata2205.945.00E-05Spermatogenesis-associated protein 22
dmrt1062.295.00E-05Doublesex- and mab-3-related transcription factor 1
dmrt3042.995.00E-05Doublesex- and mab-3-related transcription factor 3
HSD3b10.2448.530.00013 beta-hydroxysteroid dehydrogenase
AMH1.20138.895.00E-05Muellerian-inhibiting factor
LOC1033894390.3940.415.00E-05Steroid 17-alpha-hydroxylase
LOC1033841460.2322.015.00E-05EGF-like module-containing mucin-like hormone receptor-like 1
LOC1033815410.3730.155.00E-05Lutropin-choriogonadotropic hormone receptor-like isoform X1
HSD3b70.2518.960.00013 beta-hydroxysteroid dehydrogenase type 7
LOC1033884831.0675.065.00E-05Steroid 21-hydroxylase isoform X1
spef20.062.980.0001Sperm flagellar protein 2
LOC1033811932.3589.685.00E-05Estrogen receptor beta-like isoform X1
esr10.226.485.00E-05Estrogen receptor isoform X2
Sox90.3610.155.00E-05Transcription factor SOX-9
amhr22.1451.655.00E-05Anti-Muellerian hormone type-2 receptor
GATA61.5935.295.00E-05Transcription factor GATA-6
topaz10.7114.885.00E-05Testis- and ovary-specific PAZ domain-containing protein 1
GATA40.5410.335.00E-05Transcription factor GATA-4
srd5a21.6726.285.00E-053-oxo-5-alpha-steroid 4-dehydrogenase 2
ddx173.2226.455.00E-05Probable ATP-dependent RNA helicase DDX17
HSD17b13.4213.935.00E-05Estradiol 17-beta-dehydrogenase 1
LOC1033869022.5510.325.00E-05Oocyte zinc finger protein XlCOF6-like
spata5l145.098.415.00E-05Spermatogenesis-associated protein 5-like protein 1
ebp75.748.355.00E-053-beta-hydroxysteroid-Delta(8), Delta(7)-isomerase
LOC10337789537.954.185.00E-05Oocyte zinc finger protein XlCOF6-like
cyp19a1a5.852.145.00E-05Aromatase-like
FO vs. PMT
LOC10338859905.705.00E-05R-spondin-3-like
LOC10338907204.065.00E-05Spermatid perinuclear RNA-binding protein-like, partial
spata2203.075.00E-05Spermatogenesis-associated protein 22
dmrt1055.045.00E-05Doublesex- and mab-3-related transcription factor 1
dmrt3064.615.00E-05Doublesex- and mab-3-related transcription factor 3
HSD3b10.2451.350.00013 beta-hydroxysteroid dehydrogenase
AHM1.20146.55.00E-05Muellerian-inhibiting factor
LOC1033884831.06121.785.00E-05Steroid 21-hydroxylase isoform X1
HSD3b70.2525.855.00E-053 beta-hydroxysteroid dehydrogenase type 7
LOC1033894390.3936.535.00E-05Steroid 17-alpha-hydroxylase/17,20 lyase-like
LOC1033811932.35110.535.00E-05Estrogen receptor beta-like isoform X1
esr10.229.965.00E-05Estrogen receptor isoform X2
amhr22.1493.255.00E-05Anti-muellerian hormone type-2 receptor
Sox90.368.775.00E-05Transcription factor SOX-9
smox5.3541.325.00E-05Spermine oxidase isoform X2
LOC1033877253.3820.175.00E-05Steroid 17-alpha-hydroxylase/17,20 lyase
meiob123.828.065.00E-05Meiosis-specific with OB domain-containing protein isoform X2

DEGs associated with sex differentiation and gonad development in FO vs. MT and FO vs. PMT.

Table 3

GeneFPKMP-valueAnnotation

MalePseudomale
Spermatogenesis
LOC1033969970.562.675.00E-05Endothelial lipase-like
LOC1033949591.338.9875.00E-05ATP-binding cassette sub-family G member 4-like
ldlr0.538.2275.00E-05Low-density lipoprotein receptor
spef113.5757.885.00E-05Sperm flagellar protein 1
LOC10338036111.7948.085.00E-05LOW QUALITY PROTEIN: oxysterol-binding protein-related protein 5-like
LOC1033947102.6017.045.00E-05Kelch-like protein 10 isoform X3
tbx32.7713.815.00E-05T-box transcription factor TBX3 isoform X2
LOC1033866553.3919.345.00E-05Kelch-like protein 10
LOC1033806070.3111.755.00E-05MYCBP-associated protein-like
LOC1033901322.169.9655.00E-05Galactosylceramide sulfotransferase-like
LOC1033980254.400.6855.00E-05Spermatid perinuclear RNA-binding protein-like isoform X2
LOC1033911711.7319.685.00E-05Nucleoside diphosphate kinase homolog 5-like
LOC1033970330.411.8955.00E-05Citron Rho-interacting kinase-like, partial
LOC1033953176.1143.275.00E-05DNAJ homolog subfamily B member 13-like isoform X2
Energy metabolism
LOC1033911711.7319.685.00E-05Nucleoside diphosphate kinase homolog 5-like
LOC1033934623.1450.065.00E-05Nucleoside diphosphate kinase A-like
LOC1033929727.1729. 705.00E-05Nucleoside diphosphate kinase, mitochondrial-like isoform X2

DEGs associated with spermatogenesis and energy metabolism in MT vs. PMT.

qRT-PCR Validation

The expression patterns of ten DEGs (Sox9, GATA4, Dmrt1, AMH, HSD11b2, cyp19a1a, esr1, topaz1, GATA6, Sox3) associated with gonad development or steroid biosynthesis were selected for qRT-PCR validation. All the genes displayed consistent expression patterns both in qRT-PCR and RNA-seq (Figure 6). The Pearson correlation coefficient analysis exhibited correlation between qRT-PCR assay and RNA-seq data (R = 0.394, P = 0.031), indicating the accuracy and reliability of RNA-seq.

FIGURE 6

Sex-Specific Methylation Levels of Gonadal cyp19a1a Promoter

The DEGs analysis, qRT-PCR validation as well as the results of previous study (), indicated that cyp19a1a played an essential role in sex differentiation, and sex reversal induced by temperature in C. semilaevis. It was regarded that proper expression of cyp19a1a is essential for maintaining the ratio of androgen and estrogen. The balance might be destroyed by expression changes of cyp19a1a mediated by abnormal environmental temperature. Epigenetic modification is considered as one of the factors that might affect cyp19a1a expression level.

To test our hypothesis, DNA methylation of cyp19a1a promoter in gonad and muscle was examined. The CpG dinucleotides ∼2000 bp upstream of the transcription start site were selected, which had two approximate clusters: 10 CpGs in the distal promoter region (-1857 to -1718, designated as region I) and 6 CpGs in the proximal promoter region (-357 to -220, designated as region II). No difference in the methylation level was detected in the muscle tissue among females, males and pseudomales (Figure 7B). In the gonads, however, significant higher methylation levels were observed in male and pseudomale testis than in the female ovaries (Figure 7A). It was important to notice that high temperature-induced sex reversal from females to pseudomales is accompanied by the significant elevation of methylation level of gonadal cyp19a1a promoter.

FIGURE 7

To investigate if the promoter methylation would regulate the expression of cyp19a1a, qRT-PCR was performed in gonads of three groups. The expression level in females from LT and HT groups was similar, which were significantly higher (P < 0.05) than that in males (both LT and HT groups) as well as pseudomales, No expressional difference was observed between males and pseudomales (P > 0.05) (Figure 7C). Based on the methylation and expression data, we conclude that the expression level of cyp19ala showed highly negative correlation with the promoter methylation levels in gonads.

However, it was not the same case in the muscle tissue, where the cyp19a1a was only basally expressed. The average methylation levels of cyp19a1a promoter were similar and high in all, regardless of temperature treatment (Figure 7B). Two-way ANOVA analysis showed absolutely no differences among three groups in the cyp19a1a promoter methylation level in terms of temperature treatment (P > 0.05) and sex (P > 0.05). Either, no significant interaction between the two factors was found (P > 0.05).

DNA Methylation Inhibits cAMP-Stimulated cyp19a1a Promoter Activity in vitro

Transcription factor binding sites in cyp19a1a promoter were predicted using MatInspector. Two binding sites for CREB were found in the CpGs in position -1818 and -226, respectively (Figure 8B). In vitro study demonstrated that the methylation could decrease the activity of cyp19a1a promoter. The activity of unmethylated promoter could be significantly induced by forskolin stimulation. In contrast, no significantly change was observed in methylated promoter (Figure 8A).

FIGURE 8

Discussion

Since the initial discovery of vertebrate GSD + TE, the mechanism by which temperature exerts its influence on sex determination has been extensively investigated (; ; ; ; Yatsu et al., 2016). Sex reversal can be induced when the temperature achieved a threshold, and cause sex ratio change. The sex of embryos, larva or juveniles can be reversed completely or partly under a threshold temperature in reptiles and teleost (; ; ; ). The pseudomales have the same chromosome complement with females, but the phenotype is completely different (). Usually, temperature exerts its influence at TSP of embryo, larva or juvenile development, when the individuals remain sexually flexible (; ).

In this study, a teleost, C. semilaevis, sensitive to temperature was used. We demonstrated the significant imbalance of sex ratio and survival rate after high temperature treatment in TSP. The proportion of males was about 20% higher in HT group, indicating masculinization was induced by high temperature. This phenomenon have also been reported in M. menidia, D. labrax, O. niloticus, and O. mykiss (; ; Valdivia et al., 2014; Wang et al., 2017). Interestingly, only part of the females were easy to be induced sex reversal. In a recent study, a SNP (A/T) of FBXL17 had large controlling effect on sex reversal in C. semilaevis, and ZAW genotype would never reverse into phenotypic males, while those with ZTW genotypes would sometimes undergo sex reversal (). Based on these results, we speculated that some mutation might cause the females to be sensitive to temperature, and sex reversal emerged when the temperature exceeded threshold.

To explore the different expression profiling of female, male and pseudomale, RNA-seq was performed. A lot of GO terms involved in reproduction and steroid biosynthesis were identified by DEGs and enrichment analysis. Interestingly, some GO terms related to immune responses were also enriched. Similar results were displayed in Pogona vitticeps, in which the expression levels of prominent immune genes were significantly lower in pseudomales than in females and males. Further, canonical stress-related GO terms were enriched, including defense response, response to biotic stimuli (). It has been known that immune system was intertwined with stress. Meanwhile, evidences showed that stress and sex were connected in vertebrates. In Amphiprion akallopisos and Odontesthes bonariensis, cortisol was considered the regulator of sex change in response to environmental or social stress (; Yoshinaga et al., 2010; , ; ; Todd et al., 2016). In reptiles, POMC and corticosterone-mediated stress was observed in sex-reversed individuals (). In birds and rats, elevated maternal corticosterone and ACTH skewed the sex ratio of offspring (; ). In human, evidence indicated that maternal stress could enhance the circulation of corticosterone and affect neuroendocrine system. These stresses had long-lasting effects on offspring morphology, behavior, physiology, and phenotype, which could cause the imbalance of sex ratio (; ; ; van den Heuvel et al., 2018). According to a series of studies, we speculated that the C. semilaevis larva was stressed by high temperature, and immune response was activated. Then, these responses influenced endocrine system, which caused the up-regulation or down-regulation of cortisol. The biosynthesis and secretion of steroid were interfered, which leaded to sex reversal under the stress of high temperature treatment. Till now, the evidence has not been adequate, so the interaction of stress and endocrine and specific mechanism need further study.

A series of evidences of environmental influences on phenotype plasticity in vertebrate mediated by epigenetic mechanisms, such as DNA methylation and histone deacetylation has been obtained (; ). Epigenetic regulation can inhibit or stimulate gene transcription, which alters gene expression from the same genetic blueprint and thus affects development and differentiation (). In previous studies, whole-genome methylation has been found to be involved in sex-induced by high temperature in C. semilaevis, and methylation modification in sex-reversed males was inherited. Besides, dosage of Z chromosomal region was related to sex reversal in C. semilaevis (; ). However, it was found higher levels of methylation of cyp19a1a and also higher levels in gene expression of cyp19a1a (). In the present study, we found that the methylation level of C. semilaevis cyp19a1a promoter was significantly higher in males than in females. Importantly, the methylation profiles of pseudomales were similar with males, but absolutely different from females, although pseudomales had the same genotype (ZW) with females. Based on the methylation and expression data, we concluded that the expression level of cyp19a1a showed highly negative correlation with the promoter methylation levels in ovaries and testes. In Oryzias latipes, DNA methylation of cyp19a1a promoter was reported to be related to sex differentiation (). The methylation levels were twice in males compared with females in D. labrax gonads (). Besides, the allied discoveries was also observed in O. niloticus and P. olivaceus (; Wang et al., 2017). Cyp19a1a played important roles in sex differentiation by regulating estrogen synthesis. In C. semilaevis, females and pseudomales had the same genetic background (ZW), but different DNA methylation and expression levels of cyp19ala. Epigenetic modification caused by high temperature might transform the topology of DNA and block the binding of transcription factor, which could change the expression of cyp19a1a. In vitro study demonstrated that the methylation of -1818 and -226 sites in cyp19a1a promoter inhibited the binding of transcription factor CREB and suppressed the promoter activity, which could regulate the expression level of cyp19a1a. Thus, our results clearly showed that epigenetic modification, most likely DNA methylation, regulated the expression of gonadal cyp19a1a, which then mediated sex differentiation.

Interestingly, a lot of DEGs between males and pseudomales were enriched to GO terms involved in spermatogenesis, including spermatogenesis, male genitalia development, male gamete generation, spermatid development, and spermatid differentiation. Both males and pseudomales generate sperms, but the process seemed to be significantly different. In males, only Z type sperms were generated, but both Z type and W type sperms were generated theoretically in pseudomales. The generation of different types of sperm might influence spermatogenesis and spermatid differentiation and development. Surprisingly, GO terms related to energy metabolism such as UTP, GTP, and CTP biosynthetic process and metabolic process were also enriched in DEGs between males and pseudomales. Energy metabolism could affect the sperm vitality. The results implied the quality of sperm generated from males and pseudomales might be significantly different. In theory, super female (WW) individuals could be generated by W type sperm fertilized with W type eggs. However, super females were never observed in the larval stage in our lab produced by pseudomales (unpublished data). These lines of evidence suggest that W type sperm generated from pseudomale might have weak vitality. Pseudomales might unable to generate function W type sperms or the WW embryos could not develop normally to larva (Figure 9).

FIGURE 9

Conclusion

In conclusion, we demonstrated that high temperature could induce masculinization in C. semilaevis. The expression patterns of pseudomales was similar to males, but the genes involved in spermatogenesis and energy metabolism were differentially expressed. Besides, high-temperature treatment could change the epigenetic modification of cyp19a1a promoter, leading to DNA methylation level increase in pseudomales, which results in the decrease of cyp19a1a expression. There was a negative correlation between methylation levels and expression of cyp19ala. Thus the epigenetic regulation of cyp19a1a might play an essential role in the sex reversal induced by high temperature in C. semilaevis.

Statements

Ethics statement

This study was carried out in accordance with the recommendations of the Administration of Affairs Concerning Experimental Animals. The protocol was approved by the College of Marine Life, Ocean University of China.

Author contributions

JL, XL, CJ, and XD performed the experiments and analyzed the data. JL and YH prepared the figures and wrote the manuscript. QZ designed the experiments.

Funding

This research was supported by the National Natural Science Foundation of China (31802327), the National Key Research and Development Program of China (2018YFD0900601), and the China Postdoctoral Science Foundation (2017M622282).

Conflict of interest

The 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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00522/full#supplementary-material

TABLE S1

DEGs in FO vs. MT, FO vs. PMT, and MT vs. PMT.

TABLE S2

GO terms in biological process, cellular component, and molecular function levels.

TABLE S3

List of primer sequences used in the study.

References

Summary

Keywords

high-temperature treatment, RNA-seq, cyp19a1a, DNA methylation, Cynoglossus semilaevis

Citation

Liu J, Liu X, Jin C, Du X, He Y and Zhang 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

Received

16 November 2018

Accepted

13 May 2019

Published

29 May 2019

Volume

10 - 2019

Edited by

Paulino Martínez, University of Santiago de Compostela, Spain

Reviewed by

Gustavo M. Somoza, CONICET Institute of Biotechnological Research (IIB-INTECH), Argentina; Laia Ribas, Institute of Marine Sciences (ICM), Spain

Updates

Copyright

*Correspondence: Quanqi Zhang,

This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics