Comprehensive Transcriptome Analysis of Hypothalamus Reveals Genes Associated With Disorders of Sex Development in Pigs

doi:10.21203/rs.3.rs-45383/v1

Download PDF

Research

Comprehensive Transcriptome Analysis of Hypothalamus Reveals Genes Associated With Disorders of Sex Development in Pigs

https://doi.org/10.21203/rs.3.rs-45383/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 01 Jun, 2021

Read the published version in The Journal of Steroid Biochemistry and Molecular Biology →

Version 1

posted

You are reading this latest preprint version

Background

Disorders of sex development (DSD) is a chronic autoimmune disease characterized by systemic inflammation, pathological osteogenesis, and endocrine abnormality. However, its genetic etiology remains mostly unknown. In addition, little research focuses on the regulation mechanism from the view of transcriptomics in the hypothalamic-pituitary-gonadal axis (HPGA). The hypothalamus is the integrated center of the HPGA mediating neural, hormonal, and environmental stimulus to sex development.

Methods

Three XX-DSD (SRY-negative) pig (DSD) and three NF pigs (five months old, 40 kg ± 5 kg) were selected by external genitalia observation and sex determining region Y gene (SRY) detection. The hypothalamus were sampled for RNA isolation, and the mRNA, lncRNA and miRNA expression profiles were analyzed by sequencing.

Results

A total of 1,258 lncRNAs, 1,086 mRNAs, and 61 microRNAs were found to differentially express in XX-DSD pigs compared with normal female pigs. Many genes in hormone biosynthesis and secretion pathway are significantly up-regulated, and the up-regulation of GNRH1, KISS1 and AVP may be the candidate genes leading the abnormal secretion of GnRH. Next, we predicted the lncRNA-miRNA-mRNA co-expression triplets and constructed three competing endogenous RNA (ceRNA) potentially associated with DSD. Functional enrichment suggested TCONS_00340886, TCONS_00000204 and miR-181a were related to GnRH secretion.

Conclusions

Our research revealed the first transcriptomic profile in the hypothalamus of XX-DSD pigs and provided new insight in coding and non-coding RNAs that may be associated with DSD in pigs.

Animal Science

Pig

Hypothalamus

mRNAs

non-codingRNA

XX-DSD

DSD are rare conditions characterized by inconsistencies in chromosomal, gonadal, and genital sex development [1]. Humans with DSD suffer from heavy economic and mental pressure. Recent studies have shown that it will increase the risk of gonadal germ cell carcinoma [2]. In domestic animals, DSD decline the product quality, leading to serious economic losses to farmers. Pig intersex is part of animal DSD, called 38, XX disorder of sex development (38, XX-DSD). Typical features of this condition, including 38, XX karyotype, the intersex phenotype of genitalia, and the undeveloped testes or ovotestis. However, the pathogenesis and mechanism of 38, XX-DSD in pigs is unclear. XX-DSD usually accompanies by hormone secretion disorder, and hormone level is one of the criteria for DSD diagnosis. Androstenone produced in testes or ovotestes of XX-DSD pigs can affect meat quality [3]. Progesterone in XX-DSD pigs is 5–90 times higher than that of normal sows [4]. Studying the molecular mechanism of DSD hormone regulation is beneficial for understanding the pathogenesis DSD pig and for improving pork quality.

HPGA is involved in the process of mammalian reproduction from fetal development to puberty to sexual maturity [5]. The differentiation and maturation of reproductive organs are mediated by various reproductive hormones or HPGA [6; 7]. Epigenetic regulation in the HPGA has also gained attention as it is related to DSD, and non-coding RNAs (ncRNAs) are the most studied regulators of gene expression [8]. NcRNAs play important roles in promoting mammalian sexual phenotype [9]. Long non-coding RNAs (lncRNAs) may affect the expression and function of genes regulating sex determination and gonad development, for example, forehead box L2 (FOXL2) [10], Dmrt1-related gene (DMR) [11] and PISRT (polled intersex syndrome regulated transcript) [12]. The role of miRNAs in modulating HPGA is also well documented. For example, miR-124 suppresses Sox9 gene in XX gonad during ovary differentiation [13]. The expression of lin28/let-7 is associated with puberty in the hypothalamus [9], and miR-361-3p is involved in regulating FSH (follicle-stimulating hormone) secretion in a pig hypothalamus cell [14] Moreover, interactions between different protein-coding RNAs (mRNAs) and ncRNAs also play critical roles in various biological processes, such as cancer progression [15], ovarian function [16], and HPAG function [17], via acting as competing endogenous RNAs (ceRNAs) [18; 19]

The hypothalamus is the integration center of the HPGA mediating neural, hormonal, and environmental stimulus to sex development. Gonadotropin-releasing hormone (GnRH), an essential regulatory for HPGA, regulates the gonad development and sex hormone secretion [20]. In our previous study, we also found that the number of organelles in GnRH neurons of the hypothalamus increased and fused in the XX-DSD pigs [21]. However, the molecular mechanism was not clear. Here, whole-transcriptome RNA sequencing (RNAseq) of the hypothalamus in the XX-DSD pigs were performed to uncover the gene regulation mechanism at both the mRNAs, lncRNAs and miRNAs levels. Moreover, ceRNA networks were constructed to identify critical genes associated with DSD.

Experimental Design And Sample Collection

In this study, we selected the hypothalamus of DSD and normal female (NF) pigs raised at the Zhongshan pig farm (Guangdong, China). After external genitalia observation and SRY detection[21], three XX-DSD (SRY-negative) pig (DSD) and three NF pigs (five months old, 40 kg ± 5 kg) were selected to slaughter for sample collection. All samples were frozen in liquid nitrogen immediately and stored in a – 80 °C freezer until RNA isolation.

RNA isolation

According to the manufacturer’s recommended protocol, total RNA were extracted from six frozen hypothalamus using TRIzol reagent (Invitrogen, CA, USA). RNA concentration was measured using the Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). The quality of extracted RNA was determined using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library construction and sequencing

A total of 3 µg of RNA per sample was used as the material for RNAs library preparation. Ribosomal RNA was removed from each sample using the epicenter Ribo-zero TM rRNA Removal Kit (Epicenter, USA). According to the manufacturer’s recommendations, cDNA libraries were prepared using the NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina® (NEB, USA). After the purity and quality check, the libraries were sequenced on an Illumina HiSeq 2500 platform and generated 150 bp paired-end reads. Clean data of RNA-seq were obtained by removing reads containing adapter, reads containing poly-N, and low quality reads from raw data. Then, STAR (v2.5.1b) was used to map the cleaned reads to the reference genome (Sscrofa 11.1).

Small non-coding RNA (sncRNA) libraries were generated using the NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) following the manufacturer’s recommendations. The libraries were purified to recover fractions of 140 to 160 nt (the length of small non-coding RNA plus the 3' and 5' adaptors). The quality was measured with the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. Then the libraries were sequenced on an Illumina HiSeq 2500 platform and 50 bp single-end reads were generated. Several criteria were implemented to generate clean miRNA reads, including removing reads containing poly-N, 5’ adapter contaminants, ploy A or T or G or C. Low quality reads from raw data and reads without a 3' adapter or the insert tag were also removed. Furthermore, the cleaned data of miRNA were matched to the reference sequence (Sscrofa 11.1) by Bowtie [22].

LncRNA and miRNA Identification

StringTie 1.3.2d [23] was used to the transcriptome based on the clean reads mapped to the reference genome. Following transcriptome assembly, the assembled transcripts were annotated using the gffcompare program. The following four steps were used to identify the novel lncRNAs based on their characteristics: (1) Transcripts only one exon and shorter than 200 nt, were first discarded. (2) Transcripts that overlapped with the annotation exon in the database were discarded. (3) Novel transcripts without coding potential were selected by CPC2 [24], Pfam [25] and CNCI [26]. (4) The novel lncRNAs were named based on HGNC (The HUGO Gene Nomenclature Committee) standards and the location of its coding gene.

MiRBase20.0 (http://www.mirbase.org/) was used to identify known miRNA, and Mirdeep2 [27] was used to obtain the potential miRNA and draw the secondary structures. Moreover, the miREvo [28] and miRDeep2 [27] were integrated to predict novel miRNA by exploring the secondary structure, dicer cleavage site, and minimum free energy of the small RNA tags not annotated in the previous steps.

Identification of significantly differentially expressed mRNA, lncRNAs and miRNAs

HTSeq v0.6.0 [29] was used to count the reads numbers mapped to each gene. And then, the fragments per kilobase million (FPKM) of each gene was calculated based on the length of the gene and reads count mapped to this gene. Differential expression analysis between DSD and NF group was performed using the DESeq2 R package (1.10.1) [30]. The p.adj < 0.05 & |log2 (foldchange)| ≥ 1.5 were set as the threshold for significant differential expression.

The miRDeep2 quantifier module was used to quantify expression and retrieve counts for the known and novel miRNAs. The miRNA expression levels were estimated by the transcript per million (TPM) [31]. Differential expression analyses across the two groups were also run with DESeq2.

Functional analysis of DEMs

Functional analysis was performed by gene ontology (GO), gene-set enrichment analysis (GSEA) and Kyoto encyclopedia of genes and genomes (KEGG) to determine the biological significance of DEMs by the clusterProfiler R package. The GO term or KEGG pathway with a p < 0.05 is considered significantly enriched.

Protein-protein interaction (PPI) network construction

The PPI network analysis was performed based on the STRING database (https://string-db.org/). The network was visualized by Cytoscape 3.7.1, and the topological features of the network, including degree and betweenness centrality (BC) were analyzed. Hub genes were obtained by screening the degree and BC of each node in the network. The degree of a node is the number of edges connecting to other nodes. The BC is an indicator that measures the influence of a node spreading information in the network. The high BC represents the critical role of a node in communication and information diffusion [32]. To identify the significant modules in the network, Cytoscape plug-in MCODE (Molecular Complex Detection) was conducted with a score > 2.

CeRNA network construction

The CeRNA network was built by the significantly differentially expressed miRNAs (DEMis), lncRNAs (DELns) and mRNAs (DEMs). First, the targets lncRNAs and mRNA of DEMis were predicted by mi-Randa[33]. Then the Pearson correlation coefficient (PCC) was used to filter the paired the miRNA-lncRNA and miRNA-mRNA (PCC <-0.5, p < 0.05). Co-expressed lncRNA–mRNA pairs were identified by the PCC > 0.99 (p < 0.05). According to the selected standard mentioned above, both of lncRNA and mRNA in a co-expressed pair must be the target gene of a common miRNA. The qualified miRNA–mRNA–lncRNA was identified as one co-expression competing triplet. Finally, a ceRNA network was drawn by Cytoscape 3.7.1 with all the predicted competing triplets. The degree of a node was analyzed, which is the number of edges connecting to other nodes.

Confirmation and quantification of miRNAs, lncRNAs and mRNAs by qRT-PCR

To validate the sequencing datasets, we performed the Quantitative real-time PCR for DEMis, DELns and DEMs. U6 was selected as the endogenous reference gene of miRNA, and glyceraldehyde-3-Phosphate dehydrogenase (GAPDH) was selected as the endogenous reference gene of mRNA and lncRNAs. Total RNA for qRT-PCR was collected as described above. According to the manufacturer’s instructions, cDNA from lncRNA or mRNA was synthesized using the M-MLV RTase cDNA Synthesis kit (Takara, Dalian, China).cDNAs of miRNAs were reverse transcribed using One Step PrimeScript miRNA cDNA Synthesis Kit (Takara, Dalian, China). qRTPCR was carried out on the QuantStudio 3 (Applied Biosystems) with SYBR Premix Ex TaqTM II (Takara, Dalian, China) according to the manufacturer’s protocol. The primer sequences are listed in Supplementary Table 1. All reactions were tested in triplicates. Differences in threshold cycles between target genes and housekeeping genes were calculated using the 2^−ΔΔCT method.

Identification of mRNA, lncRNA and miRNA in the DSD hypothalamus

The external genitalia of XX-DSD pig showed atypical female sexual characteristics. The vulva was located below the anus, and the clitoris swelling protruded outside the vulva. After external genitalia observation and molecular identification (Fig. 1), the hypothalamus from DSD and NF groups were sampled for RNA-seq and small RNA-seq analysis. Sequencing of six transcriptome libraries generated a total of 651,894,346 raw paired-end reads, and a total of 630,405,324 clean reads were obtained (Table 1). After mapping the clean reads to the reference genome (Susrofa 11.1), we detected 22,409 mRNAs, 474 annotated lncRNAs, and 13,248 novel lncRNAs. Comparison of the genomic characterizations of the lncRNAs with mRNAs showed that lncRNAs had shorter transcripts length and open reading frame (ORF) lengths than mRNAs. A higher percentage (86.73%) of lncRNAs had 2 to 4 exons (Fig. 2A-C).

Table 1

Summary of sequencing data
	Sample
Groups	Sample	Raw_reads	Clean_reads	Mapping rate (%)	Q30 (%)	Raw_reads	Clean_reads	Mapping rate (%)	Q30 (%)
	H1	110,678,932	107,929,138	93.3	91.56	17,640,048	17,037,451	97.67	96.58
DSD	H2	109,396,514	107,296,926	87.97	91.5	22,671,555	22,196,313	97.05	96.87
	H3	95,137,076	910,320,02	91.06	91.01	19,599,180	19,375,616	97.00	97.01
	HC1	95,595,932	921,263,80	91.02	93.07	18,892,330	18,507,737	96.23	97.73
NF	HC2	124,027,710	118,310,522	85.06	92.36	18,285,778	17,932,141	97.56	97.70
	HC3	117,058,182	113,710,356	84.37	91.77	21,421,570	21,121,242	98.27	97.91

Meanwhile, we obtained a total of 118,510,461 clean reads from 116,170,500 raw reads generated in all the miRNA libraries (Table 1). The lengths of most clean reads were 20–24 nt (Fig. 2D). After mapping and annotation, about 70% of clean reads were classified as miRNAs, including 363 known miRNAs and 245 novel miRNAs (Fig. 2E).

Expression profiling of mRNAs, lncRNAs and miRNAs

Expression levels of mRNA and lncRNA were quantified using FPKM, while TPM was used to determine the expression levels of miRNA. A total of 1,086 DEMs were identified between the DSD and NF pigs (Fig. 3A). Most genes from the hormone biosynthesis and secretion pathway, such as oxytocin/neurophysin I prepropeptide (OXT), GNRH1, agouti related neuropeptide (AGRP), and follicle stimulating hormone subunit beta (FSHβ) were significantly up-regulated in DSD pigs. For ncRNAs, 61 miRNAs (37 up-regulated and 24 down-regulated; Fig. 3B) and 1,258 lncRNAs (755 up-regulated and 503 down-regulated; Fig. 3C) were differentially expressed in DSD pigs, respectively.

Functional analysis was performed to determine the biological significance of the DEMs using GO and KEGG method. Our results showed that 595 GO functions (p < 0.05) (Supplementary Table 2) were significantly annotated. Interestingly, t DEMs enriched in sex development and hormone signaling were almost up-regulated, such as doublesex and mab-3 related transcription factor 1 (DMRT1), estrogen receptor 1 (ESR1), FSHβ, Wnt family member 5A (WNT5A), steroidogenic factor 1 (NR5A1; SF-1) and GNRH1. The expression profiles of DEMs enriched in at least sex development and hormone signaling pathways were listed in Fig. 4A. The top 20 significant pathways were showed in the (Fig. 4B). According to the pathological study before, we focus on the GnRH secretion pathway (P = 0.0004). Further, KEGG-GSEA analysis showed that genes involved in the GnRH secretion pathway were up-regulated (Fig. 4C). Kiss1 metastasis-suppressor (KISS1) gene, the critical factor of mediating the negative feedback of sex steroids upon GnRH release [34], was up-regulated (Fig. 4D).

PPT network construction

The PPI network with 62 DEMs was established based on the candidate GO and KEGG analysis to investigate the relationship of DEMs and their potential role in XX-DSD (Supplementary Fig. 1). The genes with the highest BC ≥ 0.01 [35] and top five percent of highest degree [36; 37] is considered as the hub components. Four hub genes were identified and listed in Table 2, including pro-neuropeptide Y (NPY), GNRH1, AVP and glycoprotein hormones α polypeptide (CGA). We also investigated the modularity feature of the PPI network. Two functional modules, comprising 24 genes, were identified using the MCODE method (Supplementary Table 3). NPY, GNRH1, AVP, and CGA were observed to occur in the same cluster, suggesting these genes had similar functions in the DSD hypothalamus.

Table 2

Hub genes in the PPI network
Symbol	Betweenness centrality	Degree	LogFC
NPY	0.174	24	4.888
GNRH1	0.156	24	12.904
AVP	0.083	22	11.098
CGA	0.041	17	9.643

CeRNA network construction and Function analysis

A total of 117 DEMs and 101 DELns were targeted by one or multiple DEMis. We also found 23,300 DELns-DEMs pairs by the PCC method. Three ceRNA networks were constructed containing a total of 625 co-expression competing triplets. One of the ceRNA networks (CeNET1) was composed of 12 miRNA nodes, 70 mRNA nodes and 44lncRNA nodes (Supplementary Fig. 2A). CeNET2 concluded 16 miRNA, 45 mRNA and 55 lncRNA were included (Supplementary Fig. 2B). In the CeNET3, there are only one miRNA, two mRNAs and two lncRNAs (Supplementary Fig. 2B).

To further validate the potential functional implication of CeNET in the DSD hypothalamus, we performed GO and KEGG analysis of DEMs in the CeNET. The top three significant GO classifications were serine family amino acid metabolic process, odontogenesis and the ensheathment of neurons (Supplementary Table 4). The pathway annotations of DEMs involved in the cAMP signaling pathway, GnRH secretion and ErB signaling pathway (Fig. 5A) were selected for the subnetwork construction (Fig. 5B). In the subnetwork, TCONS_00198346 and TCONS_00157173 were the hub lncRNA, while ssc-miR-181a was the hub miRNA (Supplementary Table 6). Moreover, TCONS_00198346 and TCONS_00157173 shared ssc-miR-181a with calcium/calmodulin dependent protein kinase IV (CAMK4), AP-1 transcription factor subunit (FOS), glutamate ionotropic receptor AMPA type subunit 2 (GRIA2) and P21 activated kinase 6 (PAK6) gene, which enriched in the cAMP signaling pathway and ErbB signaling pathway.

Validation of differentially expressed miRNAs, mRNAs and lncRNAs

To validate the RNA-seq results, the qRT-PCR and the stem-loop qRT-PCR were performed to measure the expression levels of seven DEMs, four DELns and two DEMis. The expression profiles of these RNAs detected by qRT-PCR were consistent with those obtained by sequencing (R² = 0.867), which confirmed the reliability of our sequencing results (Fig. 6A). Moreover, four co-expression competing triplets were selected from the sub-CeNET for qRT-PCR to confirm the expression correlation of miRNA and their targets. The expression profiles of co-expression competing triplets also related (Fig. 6B).

DSD is a congenital genetic disease that is not fully understood. DSD occurs in approximately 0.18% of the human population [38] and 0.75-20% of the porcine population [39]. Previous studies have focused on the genetic mechanism of DSD. Candidate genes were identified, including WNT4 [40], R-spondin1 (RSPO1), and Wilms tumor-associated gene 1 (WT1) [41; 42]. Recently, many researchers have shown that epigenetic participates in various biological processes by regulating gene expression. In DSD patients, higher methylation levels in the CpG box located 3 kbp upstream from the SRY [43]. However, the role of ncRNAs during sex development is unclear. Moreover, hormone disorder is one of the main clinical symptoms of DSD. Nevertheless, the regulation mechanism of HPAG, which plays a vital role in the endocrine system, is not well characterized. Therefore, it is critical to study the function of ncRNAs in DSD and investigate their potential implications for DSD diagnosis.

In our study, we used RNA- and miRNA -sequencing to demonstrate the transcriptional regulation of DSD in the hypothalamus. The hypothalamus contains sexually dimorphic structures responsible for driving sex differences in behavior and physiology, controlling many facets and phases of reproduction [44; 45]. We found that ssc-miR-9 and ssc-miR-9-1 were the highest expressed miRNAs in the hypothalamus form both DSD and NF groups, which is in accordance with studies in the human [46] and goat [17]. A total of 1,086 mRNAs, 1,258 lncRNAs and 61 miRNA differentially expressed in the hypothalamus between DSD and NF pigs. Function enrichment analysis revealed that hormone receptor genes, such as ESR1, prostaglandin F receptor (PTGFR), thyroid stimulating hormone receptor (TSHR) and progesterone receptor (PGR) were significantly up-regulated, suggesting the endocrine system in the XX-DSD pig was disrupted. Moreover, two male-predominant genes, DMRT1 and SF-1, involved in the development and maintenance of the male sex features were also up-regulated [47]. SF-1 is essential for ventromedial hypothalamus (VMH) development, energy homeostasis [48], female biased brain circuitry and behavior [49]. In SF-1 KO mice, the number of immunopositive PGR cells was reduced in the VMH [49]. We thus hypothesized that the high expression of hormone receptors genes is partially related to the high expression of SF-1.

GnRH is the main factor that integrates the central and peripheral regulatory signals to dictate the gonadotropin release during sex development [50]. In the XX-DSD pigs, the number of organelles in the GnRH neurons cell and nerve fibers around cells was increased, leading to the abnormal increment of GnRH transportation [21]. Similarly, the up-regulation of the GnRH secretion pathway and GNRH1 gene was also noticed in the present study, which may stimulate CGA and FSHβ expression [51]. Generally, the secretion of GnRH is regulated by the sex steroid feedback [52]. Sex hormones can inhibit the release of GnRH byγ-aminobutyric acid (GABA) or β-endorphin, but promote the secretion by kisspeptin in the rostroperiventricular continuum of the third ventricle (RP3V) [53; 54]. The AVP in suprachiasmatic nucleus (SCN) neurons can innervate the RP3V kisspeptin neurons, and then orchestrate downstream GnRH secretion [55]. We speculated that the up-regulated of AVP/KISS1/GNRH1 signaling transduction in the XX-DSD pigs may affect the up-regulated of GnRH secretion. Furthermore, many studies also claimed that there is coordination between energy states and reproductive functions [56–58]. It has been proposed that leptin induces GnRH releasing by increases the synthesis and secretion of kisspeptin [59], and play a key roles in the hypothalamus to control the reproductive function and energy balance [60]. In our study, the up-regulated expressions of leptin receptor (LEPR), KISS1 and GNRH1 mRNA may also relate to positive energy balance caused by the abnormal elevation of testosterone disorder in DSD pig [21; 61]. In a positive energy balance, leptin has negative effects on the expression of AgRP and NPY [62], which can repress GnRH secretion in a direct or indirect fashion [57; 63; 64]. Controversially, the expressions of AGRP and NPY also increased in our study. This suggested that the GnRH secretion up-regulation in the XX-DSD was mainly due to the disruption of the endocrine regulation system, and furthering decelerated the HPAG disorder. Combing the results of PPI network, we consider that AVP, KISS1, GNRH1 are the marker genes in the XX-DSD pig hypothalamus.

ceRNA networks were constructed to identify candidate coding and non-coding RNAs associated with pig XX-DSD. By integrating the RNA-seq and miRNA-seq data, three ceRNA networks were established. DEMs in the ceRNA networks enriched in the cAMP signaling pathway, GnRH secretion and ErB signaling pathway were selected to construct a sub-CeNET. The GnRH signaling pathway continues to signal through the cAMP signaling pathway. The ErB signaling pathway [65] and the cAMP signaling pathway [66] positively regulates GnRH secretion. Ssc-miR-181a is the hub miRNA in the sub-CeNET, and it is reported that miR-181a was participated the ovarian follicle development [67]. Here ssc-miR-181a was predicted to target PAK6. PKA6 interacts with multiple binding partners including sex-steroid receptors [68], and co-expressed and inhibit AR signaling [69; 70]. Decreased expression of ssc-miR-181a may account for the higher PKA6 protein content in XX-DSD pig, leading disruption of the sex-steroid feedback of GnRH secretion. Ssc-miR-181a also targeted CAMK4, FOS and GRIA2, suggesting that ssc-miR-181a may participate in the cAMP signaling pathway.

A growing number of studies have indicated that lncRNAs regulating gene expression to function in a variety of reproduction processes at transcriptional, post-transcriptional and epigenetic levels [71; 72]. For example, lncRNAs near the master-switch gene SXL (Sex-lethal) promoter regulates sex determination [73], while OSKAR and TXS involve in gametogenesis in the Drosophila [74; 75]. In sheep, LNC_001056, LNC_000322 and LNC_000207 associated with hormone secretion and pituitary gland development [76]. Moreover, lncRNAs act as ceRNAs by miRNA competition with mRNAs. However, there is no study about the potential role of lncRNAs involved in DSD pig hypothalamus. Here, we found 38 DELns in the sub-ceRNA network, and TCONS_00198346 and TCONS_00157173 were hub lncRNAs. Both of them can target ssc-miR-181a to regulate the PKA6, CAMK4, FOS and GRIA2 gene expression, suggesting these two lncRNAs may be involved in GnRH secretion regulation in XX-DSD pigs. Interestingly, key genes in the GnRH secretion pathway, including AVP, KISS1 and GNRH, were not part of the ceRNA network, and further experiment will be conducted to explore the regulation of these genes.

DSD is a severe threat in the pig breeding. Whole-transcriptome RNAseq identified 1,086 EMs in the hypothalamus of XX-DSD pigs, and AVP, KISS1 and GNRH are considered to be the candidate genes for the XX-DSD phenotype in pig. We also found 61 DEMis and 1,258 DELns in the hypothalamus. 625 lncRNA-miRNA-mRNA co-expression competing triplets were identified and used to construct the ceRNA networks. TCONS_00198346, TCONS_00157173 and ssc-miR-181a may involve in the regulation of the GnRH release in the XX-DSD pigs.

Availability of data and materials

All the raw sequencing data were deposited in the BIG DataCenter (http://bigd.big.ac.cn/) with GSA accession No. CRA002945.

Acknowledgments

We thank the commercial fattening farm for providing theXX-DSD pigs. We also acknowledge Li Julang and Tan Yuansheng for the language editing.

Funding

This work was financially supported by the Guangdong Provincial Key Laboratory of Animal Molecular Design and Precise Breeding (2019B030301010, 2019KSYS011), the Natural Science Foundation of the Hubei Province (2019CFA075),the Innovation Team of Precise Animal Breeding (2019KCXTD004), and the National Natural Science Foundation of China (31730089 and 31722051).

Author information

Affiliations

Department of Ecology, Tibetan Centre for Ecology and Conservation at WHU-TU, Hubei Key Laboratory of Cell Homeostasis, College of Life Science, Wuhan University, Wuhan, Hubei, China

Shuwen Tan & Huabin Zhao

Guangdong Provincial Key Laboratory of Animal Molecular Design and Precise Breeding, Foshan University, Foshan, China.

Haiquan Zhao, Jinhua Wu, Hui Yu, Yin Yang, Yalan Yang& Hua Li

College of Basic Medicine, Zunyi Medical University, Zunyi, China

Yi Zhou

Department of Animal Science and Biotechnology, Chungnam National University, Daejeon, South Korea

College of Science, Tibet University, Lhasa, China

Huabin Zhao

Contributions

The contributions of the authors were as follows: Shuwen Tan, Yi Zhou, Haiquan Zhao, Jinhua Wu and Ying Yang conducted theanimal experimentation and the laboratory work. Shuwen Tan and Yalan Yang performed data analysis. Shuwen Tanwrotethe manuscript. Hui Yu provided experimental animal. Huabin Zhao performed the scientific input and proof reading. Hua Li designed the experiment, oversaw the development of the study.The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Hua Li and Huabin Zhao.

Ethics declarations

Ethics approval and consent to participate

The protocol for ourexperiments was reviewed and approved by the Institutional Animal Care and Use Committeesof Foshan Universityand Wuhan University.

Consent for publication

All authors provide consent to publish this paper.

Competing interests

The authors declare that they have no competing interests.

Terribile M, Stizzo M, Manfredi C, Quattrone C, Bottone F, Giordano DR, Bellastella G, Arcaniolo D, De Sio M. 46,XX yesticular disorder of sex development (DSD): A case report and systematic review. Medicina. 2019;55:371.
Looijenga LHJ, Kao CS, Idrees MT. Predicting gonadal germ cell cancer in people with disorders of sex development; Insights from developmental biology. Int J Mol Sci. 2019;20:5017.
Stachowiak M, Szczerbal I, Nowacka-Woszuk J, Jackowiak H, Sledzinski P, Iskrzak P, Dzimira S, Switonski M. Polymorphisms in the SOX9 region and testicular disorder of sex development (38, XX; SRY-negative) in pigs. Livestock Science. 2017;203:48–53.
Eyarefe OD, Atawalna J, Emikpe BOo, Folitse R, Dei D, Duduyemi B, Okungbowa S, Okai D. Intersex piglet with bilobed urinary bladder in KumasiI, Ghana: A case report. Animal Research International. 2017;14:2720–4.
Thackray VG, Mellon PL, Coss D. Hormones in synergy: Regulation of the pituitary gonadotropin genes. Mol Cell Endocrinol. 2010;314:192–203.
Jin JM, Yang WX. Molecular regulation of hypothalamus-pituitary-gonads axis in males. Gene. 2014;551:15–25.
Dagklis T, Ravanos K, Makedou K, Kourtis A, Rousso D. Common features and differences of the hypothalamic–pituitary–gonadal axis in male and female. J Gynecological Endocrinology. 2015;31:14–7.
Esteller M, Pandolfi PP. The epitranscriptome of noncoding RNAs in cancer. J Cancer Discovery. 2017;7:359–68.
McFarlane L, Wilhelm D. Non-coding RNAs in mammalian sexual development. Sex Dev. 2009;3:302–16.
Cocquet J, Pannetier M, Fellous M, Veitia RA. Sense and antisense Foxl2 transcripts in mouse. Genomics. 2005;85:531–41.
Bao L, Tian C, Liu S, Zhang Y, Elaswad A, Yuan Z, Khalil K, Sun F, Yang Y, Zhou T, Li N, Tan S, Zeng Q, Liu Y, Li Y, Li Y, Gao D, Dunham R, Davis K, Waldbieser G, Liu Z. The Y chromosome sequence of the channel catfish suggests novel sex determination mechanisms in teleost fish. BMC Biol. 2019;17:6.
Pailhoux E, Vigier B, Schibler L, Cribiu EP, Cotinot C, Vaiman D. Positional cloning of the PIS mutation in goats and its impact on understanding mammalian sex-differentiation. Genetics selection evolution: GSE. 2005;37:55–64.
Real FM, Sekido R, Lupianez DG, Lovell-Badge R, Jimenez R, Burgos M. A microRNA (mmu-miR-124) prevents Sox9 expression in developing mouse ovarian cells. Biol Reprod. 2013;89:78.
Ye RS, Xi QY, Qi Q, Cheng X, Chen T, Li H, Kallon S, Shu G, Wang SB, Jiang QY, Zhang YL. Differentially expressed miRNAs after GnRH treatment and their potential roles in FSH regulation in porcine anterior pituitary cell. PLoS One. 2013;8:e57156.
Zhou M, Wang X, Shi H, Cheng L, Wang Z, Zhao H, Yang L, Sun J. Characterization of long non-coding RNA-associated ceRNA network to reveal potential prognostic lncRNA biomarkers in human ovarian cancer. Oncotarget. 2016;7:12598–611.
Chen X, Shi W. Genome-wide characterization of coding and non-coding RNAs in the ovary of honeybee workers and queens. Apidologie. 2020.
Capra E, Lazzari B, Frattini S, Chessa S, Coizet B, Talenti A, Castiglioni B, Marsan PA, Crepaldi P, Pagnacco G, Williams JL, Stella A. Distribution of ncRNAs expression across hypothalamic-pituitary-gonadal axis in Capra hircus. BMC Genom. 2018;19:417.
Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA. language? Cell. 2011;146:353–8.
Tay Y, Rinn J, Pandolfi PP. The multilayered complexity of ceRNA crosstalk and competition. Nature. 2014;505:344–52.
Huang H, Tian Q. Primary amenorrhea after bone marrow transplantation and adjuvant chemotherapy misdiagnosed as disorder of sex development: A case report. Medicine. 2016;95:e5190.
Zhou Y, Zhao HQ, Liu YQ, Yu H, Shuai ISR, Li H. Study of a true hermaphrodite pig. Scientia Agricultura Sinica. 2014;47:2021–9.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome biology. 2009;10:R25.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11:1650–67.
Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, Gao G. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017;45:W12–6.
Finn RD, Bateman A, Clements J, Coggill P, Eberhardt RY, Eddy SR, Heger A, Hetherington K, Holm L, Mistry J, Sonnhammer EL, Tate J, Punta M. Pfam: the protein families database. Nucleic Acids Res. 2014;42:D222–30.
Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41:e166.
Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52.
Wen M, Shen Y, Shi S, Tang T. miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012;13:140.
Planet E, Attolini CS, Reina O, Flores O, Rossell D. htSeqTools: high-throughput sequencing quality control, processing and visualization in R. Bioinformatics. 2012;28:589–90.
Love MI, Huber W, Anders S, Huber W. Beginner's guide to using the DESeq2 package. J Theflippintruth.wordpress.com. 2014.
Zhou L, Chen J, Li Z, Li X, Hu X, Huang Y, Zhao X, Liang C, Wang Y, Sun L, Shi M, Xu X, Shen F, Chen M, Han Z, Peng Z, Zhai Q, Chen J, Zhang Z, Yang R, Ye J, Guan Z, Yang H, Gui Y, Wang J, Cai Z, Zhang X. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS One. 2010;5:e15224.
Wang C, Jiang W, Li W, Lian B, Chen X, Hua L, Lin H, Li D, Li X, Liu Z. Topological properties of the drug targets regulated by microRNA in human protein-protein interaction network. J Drug Target. 2011;19:354–64.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome biology. 2003;5:R1.
Smith JT, Cunningham MJ, Rissman EF, Clifton DK, Steiner RA. Regulation of Kiss1 gene expression in the brain of the female mouse. Endocrinology. 2005;146:3686–92.
Zhang C, Wang C, Jia Z, Tong W, Liu D, He C, Huang X, Xu W. Differentially expressed mRNAs, lncRNAs, and miRNAs with associated co-expression and ceRNA networks in ankylosing spondylitis. Oncotarget. 2017;8:113543–57.
Gargouri M, Park JJ, Holguin FO, Kim MJ, Wang H, Deshpande RR, Shachar-Hill Y, Hicks LM, Gang DR. Identification of regulatory network hubs that control lipid metabolism in Chlamydomonas reinhardtii. J Exp Bot. 2015;66:4551–66.
Lin Y, Sibanda VL, Zhang HM, Hu H, Liu H, Guo AY. MiRNA and TF co-regulatory network analysis for the pathology and recurrence of myocardial infarction. Sci Rep. 2015;5:9653.
Ozdemir M, Kavak RP, Yalcinkaya I, Guresci K. Ovotesticular disorder of sex development: An unusual presentation. J Clin Imaging Sci. 2019;9:34.
Pailhoux E, Vigier B, Vaiman D, Schibler L, Vaiman A, Cribiu E, Nezer C, Georges M, Sundstrom J, Pelliniemi LJ, Fellous M, Cotinot C. Contribution of domestic animals to the identification of new genes involved in sex determination. J Exp Zool. 2001;290:700–8.
Philibert P, Biason-Lauber A, Gueorguieva I, Stuckens C, Pienkowski C, Lebon-Labich B, Paris F, Sultan C. Molecular analysis of WNT4 gene in four adolescent girls with mullerian duct abnormality and hyperandrogenism (atypical Mayer-Rokitansky-Kuster-Hauser syndrome). Fertil Steril. 2011;95:2683–6.
Eozenou C, Gonen N, Touzon MS, Jorgensen A, Yatsenko SA, Fusee L, Kamel AK, Gellen B, Guercio G, Singh P. Testis formation in XX individuals resulting from novel pathogenic variants in Wilms’ tumor 1 (WT1) gene. Proceedings of the National Academy of Sciences. 2020.
Gomes NL, de Paula LCP, Silva JM, Silva TE, Lerario AM, Nishi MY, Batista RL, Faria Junior JAD, Moraes D, Costa EMF, Hemesath TP, Guaragna-Filho G, Leite JCL, Carvalho CG, Domenice S, Costa EC, Mendonca BB. A 46,XX testicular disorder of sex development caused by a Wilms' tumour Factor-1 (WT1) pathogenic variant. Clin Genet. 2019;95:172–6.
Gimelli G, Giorda R, Beri S, Gimelli S, Zuffardi O. A 46,X,inv(Y) young woman with gonadal dysgenesis and gonadoblastoma: cytogenetics, molecular, and methylation studies. Am J Med Genet A. 2006;140:40–5.
Mozhui K, Lu L, Armstrong WE, Williams RW. Sex-specific modulation of gene expression networks in murine hypothalamus. Front Neurosci. 2012;6:63.
Lenz KM, McCarthy MM. Organized for sex - steroid hormones and the developing hypothalamus. Eur J Neurosci. 2010;32:2096–104.
Gao Y, Sun T. Molecular regulation of hypothalamic development and physiological functions. Mol Neurobiol. 2016;53:4275–85.
Hou LP, Chen H, Tian CE, Shi WJ, Liang Y, Wu RR, Fang XW, Zhang CP, Liang YQ, Xie L. The progestin norethindrone affects sex differentiation and alters transcriptional profiles of genes along the hypothalamic-pituitary-gonadal and hypothalamic-pituitary-adrenal axes in juvenile zebrafish Dario renio. Aquat Toxicol. 2018;201:31–9.
Kim KW, Sohn JW, Kohno D, Xu Y, Williams K, Elmquist JK. SF-1 in the ventral medial hypothalamic nucleus: a key regulator of homeostasis. Mol Cell Endocrinol. 2011;336:219–23.
Grgurevic N, Budefeld T, Spanic T, Tobet SA, Majdic G. Evidence that sex chromosome genes affect sexual differentiation of female sexual behavior. Horm Behav. 2012;61:719–24.
Navarro VM. New insights into the control of pulsatile GnRH release: the role of Kiss1/neurokinin B neurons. Front Endocrinol (Lausanne). 2012;3:48.
Thompson IR, Kaiser UB. GnRH pulse frequency-dependent differential regulation of LH and FSH gene expression. Mol Cell Endocrinol. 2014;385:28–35.
Decety J, Yoder KJ, Lahey BB. Sex differences in abnormal white matter development associated with conduct disorder in children. Psychiatry Res. 2015;233:269–77.
Harter CJL, Kavanagh GS, Smith JT. The role of kisspeptin neurons in reproduction and metabolism. J Endocrinol. 2018;238:R173–83.
Piet R, Kalil B, McLennan T, Porteous R, Czieselsky K, Herbison AE. Dominant Neuropeptide Cotransmission in Kisspeptin-GABA Regulation of GnRH Neuron Firing Driving Ovulation. J Neurosci. 2018;38:6310–22.
Piet R, Fraissenon A, Boehm U, Herbison AE. Estrogen permits vasopressin signaling in preoptic kisspeptin neurons in the female mouse. J Neurosci. 2015;35:6881–92.
Ronnekleiv OK, Qiu J, Kelly MJ. Arcuate Kisspeptin Neurons Coordinate Reproductive Activities with Metabolism. Semin Reprod Med. 2019;37:131–40.
Padilla SL, Qiu J, Nestor CC, Zhang C, Smith AW, Whiddon BB, Ronnekleiv OK, Kelly MJ, Palmiter RD. AgRP to Kiss1 neuron signaling links nutritional state and fertility. Proc Natl Acad Sci U S A. 2017;114:2413–8.
Smith MS, True C, Grove KL. The neuroendocrine basis of lactation-induced suppression of GnRH: role of kisspeptin and leptin. Brain Res. 2010;1364:139–52.
Luo Q, Li W, Li M, Zhang X, Zhang H. Leptin/leptinR-kisspeptin/kiss1r-GnRH pathway reacting to regulate puberty onset during negative energy balance. Life Sci. 2016;153:207–12.
Cardoso RC, Alves BR, Prezotto LD, Thorson JF, Tedeschi LO, Keisler DH, Amstalden M, Williams GL. Reciprocal changes in leptin and NPY during nutritional acceleration of puberty in heifers. J Endocrinol. 2014;223:289–98.
Swelum AA, Ayadi M, Alhidary I, Alowaimer A, Abouheif M. The relationships between body fatness, leptin, testosterone, and reproductive performance in ram lambs as affected by level and frequency of feeding. Theriogenology. 2017;89:79–85.
Guzman A, Hernandez-Coronado CG, Rosales-Torres AM, Hernandez-Medrano JH. Leptin regulates neuropeptides associated with food intake and GnRH secretion. Ann Endocrinol. 2019;80:38–46.
Talbi R, Navarro VM. Novel insights into the metabolic action of Kiss1 neurons. Endocr Connect. 2020;9:R124–33.
Xu J, Kirigiti MA, Cowley MA, Grove KL, Smith MS. Suppression of basal spontaneous gonadotropin-releasing hormone neuronal activity during lactation: role of inhibitory effects of neuropeptide Y. Endocrinology. 2009;150:333–40.
Roy D, Angelini NL, Belsham DD. Estrogen directly respresses gonadotropin-releasing hormone (GnRH) gene expression in estrogen receptor-alpha (ERalpha)- and ERbeta-expressing GT1-7 GnRH neurons. Endocrinology. 1999;140:5045–53.
Frattarelli JL, Krsmanovic LZ, Catt KJ. The relationship between pulsatile GnRH secretion and cAMP production in immortalized GnRH neurons. Am J Physiol Endocrinol Metab. 2011;300:E1022–30.
Zayed Y, Qi X, Peng C. Identification of Novel MicroRNAs and Characterization of MicroRNA Expression Profiles in Zebrafish Ovarian Follicular Cells. Front Endocrinol (Lausanne). 2019;10:518.
Goyette SR, Schott E, Uwimana A, Nelson DW, Boganski J. Detection of the steroid receptor interacting protein, PAK6, in a neuronal cell line. Heliyon. 2019;5:e01294.
Mahfouz A, Lelieveldt BP, Grefhorst A, van Weert LT, Mol IM, Sips HC, van den Heuvel JK, Datson NA, Visser JA, Reinders MJ, Meijer OC. Genome-wide coexpression of steroid receptors in the mouse brain: Identifying signaling pathways and functionally coordinated regions. Proc Natl Acad Sci U S A. 2016;113:2738–43.
Schrantz N, da Silva Correia J, Fowler B, Ge Q, Sun Z, Bokoch GM. Mechanism of p21-activated kinase 6-mediated inhibition of androgen receptor signaling. J Biol Chem. 2004;279:1922–31.
Mercer TR, Mattick JS. Structure and function of long noncoding RNAs in epigenetic regulation. Nat Struct Mol Biol. 2013;20:300–7.
Yoon JH, Abdelmohsen K, Gorospe M. Posttranscriptional gene regulation by long noncoding RNA. J Mol Biol. 2013;425:3723–30.
Mulvey BB, Olcese U, Cabrera JR, Horabin JI. An interactive network of long non-coding RNAs facilitates the Drosophila sex determination decision. Biochim Biophys Acta. 2014;1839:773–84.
Ryu YH, Macdonald PM. RNA sequences required for the noncoding function of oskar RNA also mediate regulation of Oskar protein expression by Bicoid Stability Factor. Dev Biol. 2015;407:211–23.
Anguera MC, Ma W, Clift D, Namekawa S, Kelleher RJ 3rd, Lee JT. Tsx produces a long noncoding RNA and has general functions in the germline, stem cells, and brain. PLoS Genet. 2011;7:e1002248.
Li X, Li C, Wureli H, Ni W, Zhang M, Li H, Xu Y, Rizabek K, Bolatkhan M, Askar D, Gulzhan K, Hou X, Hu S. Screening and evaluating of long non-coding RNAs in prenatal and postnatal pituitary gland of sheep. Genomics. 2020;112:934–42.

SupplementaryFig.2.pdf
Supplementary Fig.2.The view of the CeRNA network. The rhombus represents lncRNAs, the squarerepresents mRNAs, and the circlerepresents miRNAs.
SupplementaryFigure1.pdf
Supplementary Fig. 1.The PPI network of candidate genes according to the results of function analysis
Additionalfile1.docx
Additional file 1: Supplementary Table 1.Primers for QRT-PCR;Supplementary Table 2. Supplementary Table 2. Significantly enriched GO terms of all the DEMs (p<0.05);Supplementary Table 3. The nodes topological characters in the PPI network; Significantly enriched GO terms of DEMs in the ceRNA network (p<0.05).

Download PDF

Journal Publication

published 01 Jun, 2021

Read the published version in The Journal of Steroid Biochemistry and Molecular Biology →

Version 1

posted

You are reading this latest preprint version

Comprehensive Transcriptome Analysis of Hypothalamus Reveals Genes Associated With Disorders of Sex Development in Pigs

Status:

Journal Publication

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

Background

Methods

Experimental Design And Sample Collection

RNA isolation

Library construction and sequencing

LncRNA and miRNA Identification

Identification of significantly differentially expressed mRNA, lncRNAs and miRNAs

Functional analysis of DEMs

Protein-protein interaction (PPI) network construction

CeRNA network construction

Confirmation and quantification of miRNAs, lncRNAs and mRNAs by qRT-PCR

Results

Identification of mRNA, lncRNA and miRNA in the DSD hypothalamus

Expression profiling of mRNAs, lncRNAs and miRNAs

PPT network construction

CeRNA network construction and Function analysis

Validation of differentially expressed miRNAs, mRNAs and lncRNAs

Discussion

Conclusions

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1