Comprehensive transcriptome analysis of hypothalamus reveals genes associated with disorders of sex development in pigs

lncRNA-miRNA-mRNA co-expression triplets and constructed three competing endogenous RNA (ceRNA) potentially associated with XX-DSD. Functional enrichment studies suggested that TCONS_00340886 , TCONS_00000204 and miR-181a related to GnRH secretion. Further, miR-181a inhibitor up-regulated GNRH1 , PAK6 , and CAMK4 in the GT1-7 cells. Conversely, transfection of miR-181a mimics obtained the opposite trends. The expression levels of FSHR , LHR , ESR1 and ESR2 were significantly higher in XX-DSD gondas than those in normal sows. Taken together, we proposed that the balance of endocrine had broken in XX-DSD pigs. The current study is the first to examine the transcriptomic profile in the hypothalamus of XX-DSD pigs. It provides new insight into coding and non-coding RNAs that may be associated with DSD in pigs.


Introduction
In domestic animals, sex reversal is a rare condition characterized by inconsistencies in genetic/chromosomal sex, gonadal, and genital sex development [1].Sex reversal animals usually have normal karyotype, and are categorized into XX and XY-DSD [2].XX-DSDs are relatively common in pigs, which were also called intersexes in the past [3].Typical features of this condition, include 38, XX karyotype, the intersex phenotype of genitalia, and the testis-like gonads (male pseudohermaphrodites) [4] or ovotestis (true hermaphrodites) [5].Many reports have claimed that XX-DSDs relate to the SRY-box transcription factor 9 (SOX9) mutation, leading to abnormal up-regulation during the pigs' ovarian differentiation process [6][7][8].However, the sex reversal region of SOX9 [9] or other gene in pigs needs further study.
Hormone level is a criterion for DSD diagnosis since it usually accompanies by hormone secretion disorder.Progesterone in XX-DSD pigs is 5-90 times higher than normal sows [10].Androstenone produced in testes or ovotestes of XX-DSD pigs can aff ;ect meat quality [11].Research on the molecular mechanism of XX-DSD hormone regulation can thus better our understanding of the pathogenesis of XX-DSD pig.
The hypothalamus is the integration center of HPGA mediating neural, hormonal, and environmental stimuli to sex development.Gonadotropin-releasing hormone (GnRH), an essential modulator for HPGA, regulates gonad development and sex hormone secretion [27].Our previous study found that the number of organelles in hypothalamus GnRH neurons increased and fused in XX-DSD pigs [5].However, the molecular mechanism was not clear.Here, whole-transcriptome RNA sequencing (RNA-seq) of the hypothalamus in XX-DSD pigs were performed to uncover the gene regulation mechanism at the mRNAs, lncRNAs and miRNAs levels.Moreover, ceRNA networks were constructed to identify critical genes associated with XX-DSD.

Experimental design and sample collection
Yorkshire pigs were housed indoor at an intensive farm in Guangdong, China, and fed with compound premix.After inspection of external genitalia, SRY detection and karyotype analysis [5], unrelated XX-DSD (SRY-negative) pigs and normal female pigs (NF) (five months old, 50 kg ± 5 kg) were selected for hypothalamus and gonad 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, RNA was extracted from the frozen hypothalamus and gonad using TRIzol reagent (Invitrogen, USA).RNA concentration was measured using the Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer (Life Technologies, USA).The quality of extracted RNA was determined using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA).

Library construction and sequencing
A total of 3 μg of RNA from three XX-DSD and three NF hypothalamus was used for RNAs library preparation.Ribosomal RNA was removed from each sample using the epicenter Ribo-zero TM rRNA Removal Kit (Epicenter, USA).The cDNA libraries were prepared using the NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina® (NEB, USA) per manufacturer's recommendations.After the purity and quality check, the libraries were sequenced on an Illumina HiSeq 2500 platform, then 150 bp paired-end reads were generated.Clean data of RNA-seq were obtained after removing reads containing adapter, poly-N, and of low quality from raw data.STAR (v2.5.1b) was used to map the cleaned reads to the reference genome (Sscrofa 11.1).
The 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-160 nt (the length of small non-coding RNA plus the 3' and 5' adaptors).The quality was assessed with the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips.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, and ploy A, T, G or C. Reads without a 3' adapter, or the insert tag were also removed.Furthermore, the cleaned data were matched to the reference sequence (Sscrofa 11.1) using Bowtie [28].

LncRNA and miRNA identification
StringTie 1.3.2d[29] was used to assemble transcripts based on the clean reads mapped to the reference genome.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 with 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 [30], Pfam [31] and CNCI [32], and (4) The novel lncRNAs were named based on HGNC (The HUGO Gene Nomenclature Committee) standards and the location of its coding gene.
MiRBase 20.0 (http://www.mirbase.org/)was used to identify known miRNA, and Mirdeep2 [33] was used to obtain the potential miRNA and draw the secondary structures.The miREvo [34] and miR-Deep2 [33] were integrated to predict the novel miRNAs by exploring the secondary structure, dicer cleavage site, and minimum free energy of the small RNA tags not annotated in the step as mentioned above.

Identification of significantly differentially expressed mRNAs (DEMs), lncRNAs (DELns) and miRNAs (DEMis)
HTSeq v0.6.0 [35] was used to count the numbers of reads that were mapped to each gene.And then, the fragments per kilobase millions (FPKMs) of each gene were calculated based on the length of the gene.After that, the reads count was mapped to this gene.Differential expression analysis between DSD and NF group was performed using the DESeq2 R package (1.10.1)[36].The P.adj < 0.05 and |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) [37].Differential expression analysis was examined 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 value less than 0.05 was 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.The network topological features, including degree and betweenness centrality (BC) were also analyzed using NetworkAnalyzer plugin [38].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 [39].To identify the significant modules in the network, Cytoscape plugin MCODE (Molecular Complex Detection) was conducted with a score > 2.

CeRNA network construction
The ceRNA network was built by the DEMis, DELns and mRNAs DEMs.The target lncRNAs and mRNA of DEMis were predicted by mi-Randa [40].Then the pearson correlation coefficient (PCC) was used to filter the paired 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 criteria mentioned above, lncRNA and mRNA in a co-expressed pair are the target genes of a common miRNA.The qualified miRNA-mRNA-lncRNA was identified as one co-expression competing triplet.A ceRNA network was produced 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 quantitative real-time PCR (qRT-PCR)
We performed qRT-PCR for DEMis, DELns and DEMs in 10 XX-DSD and 10 NF pigs to validate the sequencing datasets.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.The M-MLV RTase cDNA Synthesis kit (Takara, China) was used to synthesize the cDNAs from lncRNA or mRNA.The cDNAs of miRNAs were reverse transcribed using One Step PrimeScript miRNA cDNA Synthesis Kit (Takara, China).According to the manufacturer's protocol, the qRT-PCR was carried using Quant-Studio 3 (Applied Biosystems, USA) with SYBR Premix Ex TaqTM II (Takara, China).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 and evaluated by one-way analysis of variance (ANOVA) followed by Duncan's post-hoc test.

Cell culture and transfection
GT1-7 cells were purchased from Cell Bank of the Chinese Academy of Sciences, China.GT1-7 cells were maintained in DMEM supplemented with 10 % FBS, 1 % penicillin-streptomycin (Gibco, China) at 37 • C in 5 % CO 2 humidified atmosphere.The culture medium was changed twice a week, and cultures were passaged at about 80 % confluence.Changes in cell morphology and growing conditions were carefully monitored under an inverted microscope.GT1-7 cells (2*10 5 viable cells/mL) were plated in the six-well plate, treating with 100 nmol concentrations of 17 β-estradiol and testosterone (T) (Aladdin, Shanghai, China).After 12 h of culture, the medium was removed, and total cellular RNAs were extracted using TRIzol reagent (Invitrogen, USA).The expression of miR-181a and GNRH1 was detected by qRT-PCR analysis.

Endogenous reproduction-related genes in XX-DSD gonad
To investigate the effect of inappropriate GnRH secretion on XX-DSD pig gonad, the transcription levels of the endogenous reproductionrelated genes were also measured using qRT-PCR, including FSHR, luteinizing hormone receptor (LHR), estrogen receptor 1 (ER1), ER2 and androgen receptor (AR).Primer pairs for qRT-PCR are listed in Supplementary Table 2 [41].

Identification of mRNA, lncRNA and miRNA in the DSD hypothalamus
The external genitalia of XX-DSD pigs showed typical female sexual characteristics.The vulva located below the anus, and the clitoris swelling protruded outside the vulva.After external genitalia inspection and molecular identification (Supplementary Fig. 1), six hypothalamus from XX-DSD and NF groups were sampled for RNA-seq and small RNAseq analysis.A total of 651,894,346 raw paired-end reads and 630,405,324 clean reads were generated from the sequencing of six transcriptome libraries (Supplementary Table 3).After mapping the clean reads to the reference genome (Ssrofa 11.1), we detected a total of 22,409 mRNAs, 474 annotated lncRNAs and 13,248 novel lncRNAs.A comparison of the genomic characterizations of the lncRNAs with mRNAs showed that lncRNAs had shorter transcripts and open reading frame (ORF) than mRNAs.A higher percentage (86.73%) of lncRNAs had 2-4 exons (Fig. 1A-C).A total of 118,510,461 clean reads were obtained from 116,170,500 raw reads generated in all the miRNA libraries (Supplementary Table 3).The lengths of most clean reads were 20-24 nt (Fig. 1D).After mapping and annotation, about 70 % of clean reads were classified as miRNAs, including 363 known miRNAs and novel miRNAs (Fig. 1E).

Expression profiling of mRNAs, lncRNAs and miRNAs in the DSD hypothalamus
Expression levels of mRNA and lncRNA were quantified using FPKM, and TPM was used to determine the expression levels of miRNA.A total of 1,086 DEMs were identified between the XX-DSD and NF pigs (Fig. 2A).Most genes from the hormone biosynthesis and secretion pathway, such as oxytocin/neurophysin I prepropeptide (OXT), GNRH1, agouti related neuropeptide (AGRP), and FB significantly up-regulated in DSD pigs.For ncRNAs, 61 miRNAs (37 up-regulated and 24 downregulated; Fig. 2B) and 1,258 lncRNAs (755 up-regulated and down-regulated; Fig. 2C) were differentially expressed in XX-DSD pigs, respectively.
GO and KEGG method was performed to determine the biological significance of the DEMs.Our results showed that 595 significantly annotated GO functions (P < 0.05) (Supplementary Table 4).Interestingly, the DEMs enriched in sex development and hormone signaling were almost up-regulated, such as doublesex and mab-3 related transcription factor 1 (DMRT1), ESR1, FB, Wnt family member 5A (WNT5A), steroidogenic factor 1 (NR5A1; SF-1) and GNRH1.The expression profiles of DEMs enriched in both sex development and hormone signaling pathways were listed in Fig. 3A.Fig. 3B exhibited the top 20 significant pathways.According to the pathological study before [5], we focused on the GnRH secretion pathway (P = 0.0004).KEGG-GSEA showed that genes involved in the GnRH secretion pathway were up-regulated (Fig. 3C).The Kiss1 metastasis-suppressor (KISS1) gene, the critical factor of mediating the negative feedback of sex steroids upon GnRH release [42], was also up-regulated (Fig. 3D).

PPI network construction
We established the PPI network with 62 DEMs basing on the S. Tan et al.
candidate GO and KEGG analysis (Supplementary Fig. 2).The genes with the highest BC ≥ 0.01 [43] and the top five percent of highest degree [44,45] were considered as the hub components.Four hub genes were identified and listed in Table 1, including pro-neuropeptide Y (NPY), GNRH1, arginine vasopressin (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 5).NPY, GNRH1, AVP, and CGA occurred in the same cluster, suggesting these genes had similar functions in the DSD hypothalamus.

CeRNA network construction and function analysis
A total of 117 DEMs and 101 DELns were targeted by one or multiple DEMis.A total of 23,300 DELns-DEMs pairs were discovered with 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 lncRNA nodes (Supplementary Fig. 3A).CeNET2 consisted with miRNA, 45 mRNA and 55 lncRNA (Supplementary Fig. 3B).In the CeNET3, there was only one miRNA, two mRNAs and two lncRNAs (Supplementary Fig. 3B).
To further validate the potential functional implication of CeNET in XX-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 6).The DEMs involved in the cAMP signaling pathway, GnRH secretion and ErB signaling pathway (Fig. 4A) were selected for the subnetwork construction (Fig. 4B).In the subnetwork, TCONS_00198346 and TCONS_00157173 were the hub lncRNA, while miR-181a was the hub miRNA.Moreover, TCONS_00198346 and TCONS_00157173 shared miR-181a with DEMs in the cAMP signaling pathway and ErbB signaling pathway.They were   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.

Validation of DEMs, DELns and DEMis
The qRT-PCR and the stem-loop qRT-PCR were performed to measure the expression levels of seven DEMs, four DELns and two DEMis to validate the RNA-seq results.The expression profiles of these RNAs were consistent with those obtained by sequencing (R 2 = 0.867), which confirmed our sequencing results (Supplementary Fig. 4A).Four coexpression competing triplets were selected from the sub-CeNET for qRT-PCR to ensure the expression correlation of miRNA and their targets.The expression profiles of co-expression competing triplets were also related, as shown in (Supplementary Fig. 4B).

Expression profiling of reproductive hormone receptor genes in the XX-DSD gonad
The FSHR, ER1, and ER2 expressions in the DSD gonads were significantly higher (P < 0.01) (Fig. 6), and LHR expression increased (P <0 .05).However, the AR expression slightly decreased in the XX-DSD pigs (P > 0.05).

Discussion
DSD is a congenital genetic disease that is not fully understood.XX-DSD occurs in approximately 0.1− 0.5% of the porcine population [46].Previous studies have focused on the genetic mechanism of XX-DSD.The identified candidate genes include SOX9 [11,[47][48][49][50], WNT4 [51], R-spondin1 (RSPO1), and Wilms tumor-associated gene 1 (WT1) [52,53], Nuclear peceptor subfamily 2 group F member 2 (NR2F2) [54,55] and so on.Many researchers also have shown that the epigenetic participates in various biological processes by regulating gene expression.In XX-DSD dogs, the CpG box located upstream of the SOX3 has higher methylation levels [56].However, the role of ncRNAs during sex development remains elusive.Hormone disorder is one of the main clinical symptoms of DSD.Nevertheless, the regulation mechanism of HPGA, 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 their potential implications for DSD diagnosis.
We used RNA-and miRNA-sequencing to demonstrate the transcriptional profiles of XX-DSD in the hypothalamus.The hypothalamus contains sexually dimorphic structures responsible for driving sex differences in behavior and physiology and controlling many facets and phases of reproduction [57,58].We found that miR-9 and miR-9-1 were the highest expressed miRNAs in the hypothalamus from both XX-DSD and NF groups, which is in accordance with studies in human [59] and goat [24].A total of 1,086 mRNAs, 1,258 lncRNAs, and 61 miRNAs differentially expressed in the hypothalamus between XX-DSD and NF pigs.Function enrichment analysis demonstrated that hormone receptor genes, such as ESR1, prostaglandin F receptor (PTGFR), thyroid stimulating hormone receptor (TSHR), and progesterone receptor (PGR) significantly up-regulated, suggesting that the endocrine system was disrupted in XX-DSD pig.Two male-predominant genes, DMRT1 and SF-1, involving in developing and maintaining of the male sex features also up-regulated [60].SF-1 is essential for ventromedial hypothalamus (VMH) development, energy homeostasis [61], female biased brain circuitry, and behavior [62].In SF-1 KO mice, the number of immunopositive PGR cells was reduced in the VMH [62].We thus hypothesized that the high expression of hormone receptor genes is partially related to the increased expression of SF-1.
GnRH is the main factor that integrates the central and peripheral cues to regulate the synthesis and secretion of gonadotropin during sex development [63].In XX-DSD pigs, the number of organelles in the GnRH neurons cell and nerve fibers surrounding cells was increased, leading to the abnormal increment of GnRH transportation [5].Similarly, the up-regulation of the GnRH secretion pathway and GNRH1 gene was also observed in the present study, which may stimulate CGA and FSHB expression [64].Moreover, LH and FSH contents in XX-DSD pig serum were higher than those in the normal female [5].Generally, the secretion of GnRH is regulated by the sex steroid feedback [65].In the present study, the expressions of FSHR, LHR, ESR1, and ESR2 up-regulated in XX-DSD gonad.It can conclude that the massive secretion of GnRH termed the pituitary LH and FSH secretion [66], in turn, sex hormones promote the secretion by kisspeptin in the rostro periventricular area of the third ventricle (RP3V) in the XX-DSD pigs [67,68].The AVP in suprachiasmatic nucleus (SCN) neurons innervate the RP3V kisspeptin neurons, and orchestrates the downstream GnRH secretion [69].We speculated that the up-regulated of AVP/-KISS1/GNRH1 signaling transduction in XX-DSD pigs might affect the up-regulated of GnRH secretion.Many studies also claimed the coordination between energy states and reproductive functions [70][71][72].Leptin induces GnRH releasing by increasing the synthesis and secretion of kisspeptin [73].It plays a key role in the hypothalamus for reproductive function and energy balance [74].The up-regulated expressions of leptin receptor (LEPR), KISS1, and GNRH1 mRNA may also relate to positive energy balance resulting from the abnormal elevation of testosterone disorder in XX-DSD pig [5,75].In a positive energy balance, leptin has negative effects on the expression of AGRP and NPY [76], suppressing GnRH secretion either directly or indirectly [71,77,78].Controversially, the expressions of AGRP and NPY also increased in our study.We suggested that the up-regulation of GnRH secretion in XX-DSD was mainly due to the endocrine system disruption, leading to the HPGA disorder.Taking together the results of PPI network, we conclude that the up-regulation of AVP, KISS1, and GNRH1 may be the typical feature in XX-DSD pig hypothalamus.
The ceRNA networks were constructed to identify candidate coding and non-coding RNAs associated with pig XX-DSD.By integrating RNAseq and miRNA-seq data, we established three ceRNA networks.DEMs enriched in the cAMP signaling pathway, GnRH secretion, and ErB signaling pathway were selected to construct a sub-CeNET.The ErB signaling pathway [79] and the cAMP signaling pathway [80] positively regulates the GnRH secretion.The GnRH signaling pathway continues to signal through the cAMP signaling pathway.The miR-181a is the hub miRNA in the sub-CeNET, which participates in the ovarian follicle development [81].PAK6 is a target gene to the miR-181a.PKA6 interacts with multiple binding partners, including sex-steroid receptors [82], co-expresses and inhibits AR signaling [83,84].Decreased expression of miR-181a may account for the higher PKA6 protein content in XX-DSD pig, leading to disruption of the sex-steroid feedback of GnRH secretion.The miR-181a also targeted CAMK4, suggesting a possible role of miR-181a in the cAMP signaling pathway.In vitro experiment, the expression level of miR-181a can be induced by the 17 β-estradiol and T stimulation.It suggested that miR-181a may participate in the GnRH secretion regulation.In addition, overexpression of miR-181a decreased the expression of GNRH1, CAMK4 and PAK6, while miR-181a knockdown showed the contrary results.Thus, we suggested that miR-181a involved in GnRH secretion of GT1-7 cells by targeting CAMK4 and PAK6.
A growing number of studies have indicated that lncRNAs regulate gene expression in various reproduction processes at transcriptional, post-transcriptional, and epigenetic levels [85,86].LncRNAs near the master-switch gene Sex-lethal (SXL) promoter regulate sex determination [87], while OSKAR and TXS involve in gametogenesis in the Drosophila [88,89].In sheep, LNC_001056, LNC_000322, and LNC_000207 associate with hormone secretion and pituitary gland development [90].Moreover, lncRNAs may function as ceRNAs and compete with mRNA interacts with miRNAs.However, there are no studies about the potential role of lncRNAs involved in the hypothalamus of DSD pigs.Here, we found 38 DELns in the sub-ceRNA network, and TCONS_00198346 and TCONS_00157173 were hub lncRNAs.Both can target miR-181a to regulate PKA6 and CAMK4 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 GNRH1, were not part of the ceRNA network.Further experiments will be conducted to explore the regulatory role of these genes.It will also be interesting to study the effect of high-level GnRH on gene regulation of the pituitary and gonad in the HPGA.The molecular mechanism of embryonic XX-DSD pigs should be investigated in the future.

Conclusions
DSD poses a severe threat to the pig breeding industry.Wholetranscriptome RNA-seq identified 1,086 DEMs in XX-DSD hypothalamus.Up-regulations of AVP, KISS1 and GNRH1 were detected in XX-DSD hypothalamus, and considered marker genes for XX-DSD diagnosis.Moreover, endogenous reproduction-related genes up-regulated in XX-DSD gonad, indicating the endocrine disorder in XX-DSD pigs.In this study, 61 DEMis and 1,258 DELns were discovered in the hypothalamus.A total of 625 lncRNA-miRNA-mRNA co-expression competing triplets were identified and used to construct the ceRNA networks.TCONS_00198346, TCONS_00157173, and miR-181a may involve in the regulation of the GnRH release in XX-DSD pigs.The biological function of miR-181a was verified in vitro experiment.

Fig. 1 .
Fig. 1.Identification of lncRNA and miRNA.A-C: Comparison of lncRNA and mRNA with respect to the exon number, transcript length and ORF length.D-F: Length and type distribution of all the small RNA in the XX-DSD and NF group.

Fig. 3 .
Fig. 3. Function analysis of DEMs.A: Heatmap of DEMs enriched in the sex development and hormone signaling related GO terms.B: KEGG pathway analyses.C: KEGG-GSEA analysis of GnRH secretion.D: Up-regulated DEMs enriched in the GnRH secretion pathway.

Fig. 4 .
Fig. 4. Function analysis and sub-CeNET construction of candidate RNA.A: KEGG pathway of DEMs in the ceRNA net word.B: Sub-CeNET constructed with 12 related DEMs and all DELns and DEMis in the hypothalamus.

S
. Tan et al.

Table 1
Hub genes in the PPI network.