Skip to main content
Download PDF

Combining ATAC-seq and RNA-seq reveals key genes for gonadal abnormalities in one-month-old XX-DSD pigs

BMC Genomics volume 26, Article number: 447 (2025) Cite this article

Abstract

Background

Disorders of Sex Development (DSD) are caused by congenital abnormalities in the chromosomes, and subsequent development of gonads or sexual anatomy. XX-DSD pigs exhibit a series of adverse symptoms such as sterility, genital infections, and decline in meat quality, leading to significant economic losses in the breeding industry. However, the understanding of the etiology and pathogenesis of XX-DSD in pigs remains limited. To investigate the molecular mechanisms underlying abnormal gonadal development in XX-DSD pigs, we analyzed the gonads of 1-month-old XX-DSD pigs, normal females, and normal males using RNA-seq and ATAC-seq techniques.

Results

From RNA-seq, we identified potential genes involved in gonadal development in XX-DSD pigs, including SOX9, HSD3B1, CYP19A1, CCNB2, CYP11A1, DMRT1, and MGP. Following this, we analyzed ATAC-seq data and identified 14,820 differential accessible chromatin peaks. Then, by integrating the ATAC-seq and RNA-seq analysis results, we identified several candidate genes (SOX9, COL1A1, COL1A2, FDX1, COL6A1, HSD3B1, FSHR, and CYP17A1) that might be associated with sex development. Through PPI (Protein-Protein Interaction Networks) analysis, we found that SOX9 gene was the top hub gene. Furthermore, we confirmed the effect of the open chromatin region on SOX9 gene expression by a dual-luciferase reporter assay, thus further validating the critical role of this open region in regulating SOX9 expression.

Conclusions

This study elucidates the critical regulatory role of specific open chromatin structures in the SOX9 gene promoter region (8647563–8648475) in gonadal development of XX-DSD pigs. Additionally, we identify that genes such as SOX9, HSD3B1, and CYP19A1 act in concert to participate in gonadal development. These findings provide molecular evidence for the dynamic chromatin regulatory network underlying gonadal dysgenesis in XX-DSD and lay the foundation for subsequent mechanistic studies.

Peer Review reports

Background

Disorders of Sex Development (DSD) refer to atypical features in genetic, gonadal, and phenotypic sex, resulting in inconsistencies in the development of internal and external genitalia [1]. In swine, karyotypes designated as 38, XX-DSD, reflect a genetic disorder resulting from the deficiency of the SRY gene. It is prevalent in some breeds, such as Large White pigs and Landrace pigs, or in their hybrids. An XX-DSD pig often has diminished reproductive performance and increased susceptibility to genitourinary infections, which frequently hinder the advancement of high-quality pig production. Although, the etiology of gonadal development in XX-DSD pigs has achived some progress [2], but the molecular mechanisms still remain unclear.

Chromatin accessibility refers to the extent of physical contact between nucleosomes and chromatin DNA. This property is influenced by nucleosome composition and post-translational modifications, thereby affecting gene expression [3,4,5]. Recent genomic studies in porcine XX-DSD have identified high-penetrance variants in IFITM1, NOBOX, and WWOX through whole-genome sequencing, suggesting that these genetic alterations may be associated with DSD [6]. Open chromatin regions (OCRs) were specific areas in the genome typically associated with the binding of transcription factors (TFs) and other regulatory proteins, which could influence gene expression. In the process of sexual differentiation and sex reversal, changes in chromatin accessibility were key factors [7]. The accessible regions where TFs bound were believed to drive differential gene expression and to be involved in determining gonadal fate [8,9,10]. For example, using ATAC-seq and ChIP-seq techniques, researchers identified changes in chromatin accessibility in support cells of both female and male mouse gonads before and after sex determination, and discovered a large number of cis-regulatory elements (CREs) involved in this process [11,12,13,14]. Therefore, based on the relationship between chromatin accessibility and gene expression, it is important to understand the dynamic patterns and regulatory functions of open chromatin regions in pig sex differentiation and sex reversal.

The chromatin accessibility characteristics of gonads in DSD pigs have not been studied yet. In this study, we conducted RNA sequencing coupled with functional and pathway enrichment analyses on gonadal tissues from one-month-old normal males (M), normal females (FM), and XX-DSD pigs to elucidate the molecular regulatory networks associated with abnormal sexual development. Moreover, we performed ATAC sequencing to analyze chromatin accessibility characteristics in gonadal tissues from both normal females and XX-DSD pigs. This study provides essential data and theoretical insights into the molecular mechanisms underlying XX-DSD pig pathogenesis by systematic analysis of chromatin accessibility and gene expression profiles in these pigs.

Methods

Animal collection and assessment

One-month-old XX-DSD pigs, normal females, and normal males (Landrace×Yorkshire crossbred hybrids) were provided from a pig farm in Guangdong province. XX-DSD pigs were preliminarily identified through sex organ appearance examination and SRY gene testing (Table S1). Normal females (FM, n = 3), normal males (M, n = 3), and XX-DSD pigs (n = 5) were collected as test subjects. Their gonadal samples were collected and frozen in liquid nitrogen, then transferred to an ultra-low temperature refrigerator at -80 °C for long-term storage. In addition, serum samples were collected for hormone analysis.

The levels of three hormones, testosterone (T), estradiol (E2), and progesterone (PROG), in pig serum were determined using an enzyme-linked immunosorbent assay kit provided by Shanghai Enzyme-linked Biotechnology Co., Ltd.

RNA extraction, library construction and sequencing

Total RNA was isolated from 50 to 100 mg of left-side gonadal tissues using Trizol Reagent kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The analyzed samples included intact ovaries from normal female pigs, whole testes from normal male pigs, and ovotestis tissues from XX-DSD pigs. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase free agarose gel electrophoresis. After total RNA was extracted, eukaryotic mRNA was enriched by Oligo (dT) beads. Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse-transcribed into cDNA by using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #7530, New England Biolabs, Ipswich, MA, USA). The purified double-stranded cDNA fragments were end repaired, A base added, and ligated to Illumina sequencing adapters. The ligation reaction was purified with the AMPure XP Beads (1.0X). Ligated fragments were subjected to size selection by agarose gel electrophoresis and polymerase chain reaction (PCR) amplified. The resulting cDNA library was sequenced using Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

ATAC library Preparation and sequencing

ATAC-seq Libraries were constructed using three samples of normal females and five samples of XX-DSD. Nuclei suspensions were incubated in a transposition mix that includes a transposase. The transposase enters the nuclei and preferentially fragments the DNA in open regions of the chromatin. Simultaneously, adapter sequences were added to the ends of the DNA fragments. Incubate reaction was performed at 37℃ for 30 min. Immediately following transposition, the products were purified using a QIAGEN minielute kit. Following purification, we amplified library fragments using 1x NEBnext PCR master mix and 1.25 µM of custom Nextera PCR primers 1 and 2, using the following PCR conditions: 72 °C for 5 min, 98 °C for 30 s, followed by thermocycling at 98 °C for 10 s, 63 °C for 30 s and 72 °C for 1 min. The libraries were purified using a Qiagen PCR cleanup kit yielding a final library concentration of ~ 30 nM in 20 µL. Clustering of the index-coded samples was carried out on a cBot cluster generation system using HiSeq PE Cluster Kit v4-cBot-HS (Illumina). The libraries were sequenced on an Illumina platform after the cluster generation and 150 bp paired-end reads were produced [15]. The original Illumina high-throughput sequencing image data files were transformed into original sequences after base group recognition using the bcl2fastq software v1.8.4, and were stored in the FASTQ file format.

RNA‑seq data analysis

Paired-end reads were aligned to the pig reference genome (Sscrofa11.1) by HISAT2 v2.1.0 [16] with default parameters, the mapped reads of each sample were assembled by using StringTie v1.3.1 [17, 18]. For each transcription region, a TPM (Transcripts Per Kilobase of exon model per Million mapped reads) value was calculated to quantify its expression abundance and variations, using RSEM software [19]. RNAs differential expression analysis was performed by DESeq2 software v1.30.1 [20] between two different groups (and by edgeR [21] between two samples). The genes with the parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≤ 2 were considered differentially expressed genes. Principal component analysis (PCA) was performed with R package gmodels (http://www.r-project.org/) in this experiment.

ATAC‑seq data processing

Bowtie v2.2.4 [22] with the parameters – X2000 and–m1 was used to align the clean reads from each sample against the reference genome (Sscrofa11.1), Samtools was used to convert the SAM files to the BAM format. Then, the peaks of the above files were detected by MACS2 v2.1.0 [23]. According to the genomic location information and gene annotation information of peak, peak related genes can be confirmed using ChIPseeker [24]. The DiffBind [25] was used to analyse peak differences across groups.

Integration of ATAC-Seq and RNA-Seq

Genes with peaks in ATAC-seq were selected if they were located within the promoter regions, 2 kb upstream and 1 kb downstream of the transcription start site (TSS). TPM values from RNA-seq were combined with these peaks, averaged, and log2 transformed. Box plots were generated in R, and nonparametric tests (Mann-Whitney U test or Wilcoxon signed-rank test) were performed. Overlap analysis and IGV visualization were employed to assess genes differentially expressed in both ATAC-seq and RNA-seq. Finally, transcription factors predicted by Homer motifs were overlaid with RNA-seq data to identify potential transcription factors, and the gene IDs were converted using gProfiler and UniProtKB to obtain pathway commons identifiers.

Footprint prediction

OCRs of promoter regions from genes showing common differences in RNA-seq and ATAC-seq were selected. These sequences underwent motif analysis, followed by footprint analysis based on the identified motifs. The motif footprint analysis was performed using the R package ATAC seq QC and the JASPAR 2022 CORE vertebrates non-redundant database.

RT-qPCR experimental verification

To validate the relative expression patterns identified through RNA-seq analysis, three mRNAs deferentially expressed in the gonads were selected for RT-qPCR analysis (Table S2). RT-qPCR was conducted using the ABI PRISM 7500 Fast Real-Time PCR System, with GAPDH serving as the internal reference gene. Relative gene expression levels were calculated as described by Livak and Schmittgen (2001; 2−ΔΔCt method), and data were presented as mean ± SEM. Differential expression was statistically analyzed using a t-test and one-way ANOVA, with significance denoted by p < 0.05 (*) or p < 0.01 (**). Each group comprised a minimum of three biological replicates, including three normal females, three normal males, and five XX-DSD pigs. Each biological sample underwent triplicate technical qPCR runs to assess experimental consistency.

Double fluorescent reporter enzyme experiment

To validate the transcriptional activity of the SOX9 promoter region (8647563–8648475), a 912-bp fragment was amplified from genomic DNA using primers designed with NCBI Primer-BLAST (Table S3). The psiCheck2 plasmid was digested with XhoI/BglII to remove the SV40 promoter region, generating a negative control plasmid (psiCheck2-ΔSV40) retaining only the basic reporter gene framework. Simultaneously, the SOX9 promoter fragment was cloned into the same restriction sites of the psiCheck2 plasmid to construct the experimental plasmid (psiCheck2-SOX9), forming the mutant reporter vector. These vectors were co-transfected into 293T cells seeded in 24-well plates at a density of 5 × 10⁴ cells per well. Transfection was performed when cell confluence reached 80%. After 48 h, cells were harvested, and the activities of Firefly luciferase (detection wavelength: 560 nm) and Renilla luciferase (detection wavelength: 480 nm) were sequentially measured using a dual-luciferase assay system (Lux-T020). Data were normalized by the Firefly Luc/Renilla Luc ratio. At least three independent biological replicates (each containing three technical replicates) were set up for the experiment. The non-paired t-test was performed using GraphPad Prism 9.0 to compare the statistical differences, with a significance threshold of P < 0.05.

Results

Detection of XX-DSD pigs SRY gene

The amplification of the SRY and GAPDH genes was successfully achieved from the DNA of normal males. The SRY gene was not detected in the DNA of normal females or XX-DSD pigs (Fig. S1). These results indicated that XX-DSD pigs were negative of SRY genes on the chromosomes.

The phenotypes of XX-DSD pigs

Appearance examination revealed that XX-DSD pigs displayed both a scrotum and a vulva, indicating a hermaphroditic malformation. The vulva was located below the anus and appeared morphologically deformed. Additionally, bilateral testicular-like masses were observed within the scrotum. Postmortem examination revealed mixed Müllerian and Wolffian duct derivatives, encompassing epididymis, vas deferens, and uterine structures (Fig. 1). The gonads were hypoplastic testes with dark-brown coloration, oval morphology, and loose texture. Hormone assays revealed no significant differences in PROG concentration. However, the E2 concentration in XX-DSD pigs was the lowest compared to the other two groups. Additionally, the T concentration in XX-DSD pigs was intermediate between that of normal males and normal females (Table S4).

Fig. 1
figure 1

Phenotypic characteristics of normal females, XX-DSD pigs and normal males. 1 Anus. 2 Vulva. 3 scrotum. 4 Kidney. 5 Ureter. 6 Bladder. 7 Ovary. 8 Uterine horn. 9 Uterine body. 10 Rectum. 11 Tail. 12 Epididymis. 13 hypoplastic testes. 14 testis. 15 Vas deferens/Oviduct. 16 Vas deferens. 17 Penis

Full size image

Results of RNA-seq data from XX-DSD pigs

Principal component analysis (PCA) (Fig. 2a) and Pearson’s correlation coefficients (Fig. 2b) revealed both similarities and differences between the biological replicates and differences. The PCA results indicated clear variations in the transcriptome profiles among the three groups.

DEGs were identified through pairwise differential expression analysis among the three groups (Fig. 2c). In the FM vs. XX-DSD comparison, 348 DEGs were up-regulated, while 1,316 DEGs were down-regulated. GO and KEGG pathway enrichment analyses were performed on DEGs to construct PPI networks using key genes. GO enrichment analysis indicated that these DEGs were mainly involved in developmental and reproductive processes (Fig. 3a). KEGG enrichment analysis revealed significant enrichment in pathways related to EMC-receptor interaction, ovarian steroidogenesis, and Wnt signaling (Fig. 3b). PPI network analysis identified CCNB2, FDX1, SOX9, CYP19A1, and HSD3B1 as among the top ten genes based on degree values (Fig. 3c).

Compared with the FM vs. M group, 358 DEGs were found to be up-regulated and 1,178 DEGs were down-regulated. GO enrichment analysis revealed that these DEGs were predominantly involved in reproductive processes, cell differentiation, and reproduction pathways (Fig. 3a). KEGG analysis demonstrated significant enrichment in pathways related to steroid biosynthesis, signaling pathways regulating the pluripotency of stem cells, and ovarian steroidogenesis (Fig. 3b). Screening of genes for PPI analysis of pathways related to reproductive development identified CYP11A1, ESPL1, CDSA8, SOX9, and CYP19A1 as the top five genes based on degree values (Fig. 3c).

Compared the M vs. XX-DSD group, 30 DEGs were up-regulated, and 27 DEGs were down-regulated. GO enrichment analysis revealed that these DEGs were primarily involved in intracellular anatomical structures, binding and cytoplasm pathways (Fig. 3a). KEGG analysis indicated significant enrichment in the MAPK signaling pathway, hippo signaling pathway and other pathways (Fig. 3b). Additionally, PPI network analysis identified ISG15, MX1 and other genes were among the top ten based on degree values (Fig. 3c).

Meanwhile, the Venn diagram (Fig. 2d) showed a total of 1,292 DEGs between FM vs. XX-DSD and FM vs. M groups. Additionally, six DEGs were found to be common across all three comparison groups.

Fig. 2
figure 2

Analysis of RNA-seq data of normal females, XX-DSD and normal males. a PCA Clustering Plot of RNA-seq Sequencing Date. b Pearson Correlation Coefficients of RNA-seq Sequencing Data. c Bar graph showing the number of differentially expressed genes (DEGs). (|log2FC| > 1, p < 0.05) d Venn diagrams between different groups

Full size image
Fig. 3
figure 3

Enrichment analysis of differentially expressed mRNAs. a GO enrichment analysis of differentially expressed genes. b Bubble chart of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of differentially expressed genes. c Protein-Protein Interaction Network Diagram

Full size image

Analysis results of XX-DSD pigs ATAC-seq data

To explore the OCR differences between XX-DSD pigs and normal females, ovaries of normal females were compared with dysplastic testes of XX-DSD pigs. PCA revealed the similarity between biological replicates and the difference between the two groups. The correlation coefficient was 0.99 for normal females and 0.95 for XX-DSD pigs (Fig. 4a, b). These imply that the sequencing data had an excellent quality.

The results indicated that there were 14,820 different OCRs, there are 9,640 increased and 5,180 decreased ORCs. (Fig. 4c). The GO enrichment analysis of the genes corresponding to the screened OCRs revealed significant enrichment in functions such as binding, anatomical structure development, and protein binding. KEGG enrichment analysis revealed significant enrichment in pathways such as axon guidance and the oxytocin signaling pathway.

Analysis of ATAC-seq read distribution across the genome indicated significant enrichment primarily within 2 Kb upstream and downstream of transcription start sites (TSS). To proximity of TSS was positively correlated with ATAC-seq signal intensity, indicating high-quality library construction (Fig. S2a). To identify the different gene functional elements captured by the ATAC-seq data of each sample, the genome was stratified into promoters, exons, 5’and 3’untranslated regions, introns, and other regions. It was found that the proportions of these regions differed significantly (Fig. S2b), along with differences in the proportions of each genomic region (0-1 kb to > 100 kb) between the two groups (Fig. S2c).

OCRs can modulate gene expression by interacting with transcription factors or other regulatory elements. To explore transcription factors potentially regulating differential OCRs between FM and XX-DSD, motif enrichment analysis using homer was conducted. Top enriched transcription factors included NR5A2, ESRRB, ESRRG, NR5A1, TRPS1, GATA4, DMRT1, GATA2, NR2F1, and FOS. These findings suggest that these transcription factors may be critical for gonadal development in XX-DSD pigs, warranting further investigation (Fig. 4f).

Fig. 4
figure 4

ATAC-seq data analysis of normal females and XX-DSD. a PCA clustering plot of ATAC-seq sequencing data. b Pearson Correlation Coefficients of ATAC-seq sequencing data. c Bar graph of differential expression between FM and DSD group. d Bubble chart of top 20 GO enrichment pathways. e Bubble chart of top 20 KEGG enrichment pathways. f Motif enrichment analysis of differential OCRs

Full size image

Correlation analysis results of XX-DSD pigs RNA-seq and ATAC-seq data

Association analysis between ATAC-seq and RNA-seq was conducted. Initially, we compared the impact of promoter OCRs on gene expression across two datasets. We found that genes associated with promoter OCRs had significantly higher expression than others, indicating their regulatory role (Fig. 5a). Subsequently, the identified OCRs were correlated with RNA-seq data to investigate their influence on gene expression. Notably, promoters were found to exert significant effects on gene expression between normal females and XX-DSD pigs (Fig. 5b).

The intersection of differentially expressed genes from RNA-seq and ATAC-seq sequencing results identified a total of 716 common differentially expressed genes. Subsequently, we mapped the regulatory elements of these intersecting genes. The findings revealed that the promoter regions of 141 of these genes contained open chromatin regions, which facilitated further investigation (Fig. 5d). Using the STRING database, we constructed a PPI network diagram for the 141 selected genes. Notably, the top eight genes by degree centrality in the PPI network included SOX9, COL1A1, COL1A2, FDX1, COL6A1, HSD3B1, FSHR, and CYP17A1. Among them, SOX9, COL1A1, and COL1A2 were centrally positioned (Fig. 5c). Based on the previous motif enrichment analysis, footprint analysis further revealed that the transcription factors of ESRRB, ESRRG, FOS, NR5A1, NR5A2, and SP5 had the potential to bind to target genes, indicating that these factors may play important roles in disease pathogenesis (Fig. 5e).

GO enrichment analysis was conducted on 141 genes (Fig. 6). It was indicated by the results that significant enrichment was observed in pathways related to anatomical structure development, developmental processes, multicellular organism development, and system development. It was noted that all these pathways are associated with growth and development.

To further investigate OCRs associated with differential transcription factors and gene promoters, we visually examined several transcription factors (ESRRB, ESRRG, FOS, NR5A2, NR5A1, and SP5) (Fig. 7a) and the promoter regions of AMH, CYPB1, FSHR, SOX9, HSD3B1, and HSD17B3 genes using IGV (Fig. 7b). The reliability of our data was confirmed by these results.

Fig. 5
figure 5

Correlated RNA-seq and ATAC-seq data. a The box plot of the relationship between OCRs located in the promoter and gene expression. b Boxplot of OCRs located in promoters and gene expression. c PPI (Protein-Protein Interaction) Network Interaction Graph. d Venn diagram of differentially expressed genes from RNA-seq and ATAC-seq between normal females (FM) and XX-DSD pigs, including the intersecting co-differentially expressed genes located in the promoter regions. e Results of Footprint analysis

Full size image
Fig. 6
figure 6

Bubble Chart of Top 20 GO enrichment pathways

Full size image
Fig. 7
figure 7

IGV visualization. a Transcription factor IGV visualization. b IGV visualization of differential gene promoter regions

Full size image

Verification results of RT-qPCR experiment

To validate the transcriptomic findings, we performed RT-qPCR on three key genes (DMRT, MGP, and SOX9) showing differential expression across distinct gonadal tissues: ovaries from female pigs, testes from male pigs, and ovotestis tissues from XX-DSD pigs (Fig. 8). Comparison of RT-qPCR results with RNA sequencing data revealed consistency in the expression patterns of the selected coding RNAs, affirming the accuracy of the sequencing data investigated in this study.

Fig. 8
figure 8

Verifications of DEGs by RT-qPCR

Full size image

Verification results of dual luciferase experiment

Dual-luciferase assays confirmed that the SOX9 promoter region (8647563–8648475) exhibits intrinsic transcriptional activity. The assay confirmed that the SOX9 promoter region exhibits significantly elevated transcriptional activity compared to the psiCheck2-ΔSV40 (Fig. 9). This validates that the chromatin accessibility observed in ATAC-seq data directly drives transcriptional activation of SOX9.

Fig. 9
figure 9

Dual-Luciferase Reporter Assay Result

Full size image

Discussion

Sex determination is the process by which the undifferentiated gonads develop into either ovaries or testes, along with their accessory sex organs, during early embryonic development [1]. The undifferentiated genital ridge initially possesses bipotential gonads, and sex determination and differentiation are regulated by genes [26]. Studies have shown that in normal male development, testicular development was induced by SRY gene on the Y chromosome. Androgens such as testosterone are then secreted by the testes, promoting the formation of male reproductive organs [27]. In contrast, during normal female development, the interaction between paracrine signals and transcription factors in the bipotential gonads activates the classical Wnt signaling pathway in the absence of SRY expression [28, 29]. This leads to the accumulation of β-catenin, which induces genes involved in ovary formation and represses the activity of the SOX9 gene [30]. Dysregulation during sex development can result in Disorders of Sex Development (DSD). XX-DSD, characterized by hermaphroditism and malformed reproductive organs, had been previously reported in six-month-old XX-DSD pigs [31], suggesting potential functional abnormalities in their reproductive systems. Notably, a critical 70 bp deletion in intron of the WWOX was identified in 62.5% of XX-DSD pigs through prior investigations, which was hypothesized to contribute to abnormal gonadal development [6]. To systematically explore the molecular mechanisms underlying gonadal dysplasia, this study implemented an integrated approach combining ATAC-seq and RNA-seq analyses. This strategy enabled comprehensive identification of key regulatory elements, transcription factors, and functionally relevant genes associated with gonadal development. Our findings provided novel mechanistic insights into the molecular dysregulation driving gonadal abnormalities in XX-DSD pigs.

In this study, the testosterone level in XX-DSD pigs was higher than that in normal females but lower than that in normal males. Particularly in the combined analysis, several key enzymes involved in testosterone biosynthesis showed high expression levels. As CYP11A1, a key enzyme in the testosterone synthesis pathway, catalyzes the conversion of cholesterol into pregnenolone. Its high expression in XX-DSD pigs suggests an enhanced conversion process from cholesterol to pregnenolone [32]. CYP17A1 further catalyzes the production of testosterone precursors, and its high expression may promote testosterone synthesis in XX-DSD pigs [33]. Additionally, in the gonads, a crucial role was played by CYP17A1. It promoted the metabolism of cholesterol into androgen and estrogen precursors [34]. Furthermore, cholesterol transformation and synthesis were involved, and the synthesis and balance of adrenal corticosteroid and sex hormone were regulated. Defects or mutations in CYP17A1 may lead to abnormal steroid hormone synthesis, thereby affecting gonadal development and hormone balance [35, 36]. Furthermore, we identified SOX9, COL1A1, COL1A2, FDX1, COL6A1, HSD3B1, and FSHR genes as key genes in XX-DSD pigs sex development, which are associated with steroid hormone production, gonadal development, and hormone pathways. The enzyme 3β-HSD1 converts 17α-hydroxyprogesterone into precursor molecules of androgens during testosterone synthesis, a critical step in the testosterone synthesis pathway in the testes [37]. Defects or mutations in HSD3B1 may reduce 3β-hydroxysteroid dehydrogenase activity, thereby affecting hormone synthesis and metabolism [38]. Studies have shown that mutations in the HSD3B1 gene lead to reduced 3β-hydroxysteroid dehydrogenase. This reduction further affected hormone synthesis and metabolism. Consequently, these changes resulted in abnormal gonadal development [39].

The proteins encoded by COL1A1 and COL1A2 are the main components of collagen. A crucial role was play in connective tissue, bone, cartilage, skin, blood vessels, and other structures. Additionally, they were considered essential for maintaining the structure and function of the human body [40]. The expression of the SOX9 gene plays a crucial regulatory role in gonadal development [2, 41, 42]. Once activated by SRY, SOX9 expression in the gonads is significantly enhanced, and it also mutually regulates with other key genes in gonadal development, such as NR5A1 and WNT4. Proper functioning of this reciprocal regulatory network is essential for normal gonadal development [43, 44]. In this study, the transcriptional activity of the SOX9 gene promoter region was confirmed through dual luciferase reporter assays. Combined with previous studies on humans and Amami spiny rats [45,46,47], it was found that Enh13 and Enh14, enhancers upstream of the SOX9 gene, regulate SOX9 expression across different species, playing a crucial role in sex determination and gonadal development. These findings provide new insights into the molecular mechanisms of sex determination and may offer valuable information for the study and treatment of related genetic diseases. According to the interaction relationships depicted in the PPI network diagram, SOX9, COL1A1, and COL1A2 can be explored as novel angles for studying gonadal dysplasia.

ATAC-seq footprint analysis provides detailed information on the actual binding sites of transcription factors by directly identifying where DNA-binding proteins attach. It reflects the binding state of transcription factors under specific conditions and reveals how these factors influence gene regulation. It is associated with open chromatin, and the binding situation of transcription factors can be directly observed in in-situ sequencing. Compared with motif enrichment analysis, it has more direct and comprehensive advantages, which is helpful to further understand the transcriptional regulation mechanism of the genome. In this study, motif prediction was performed on differentially expressed OCRs. To further screen the target transcription factors, a high-confidence footprint analysis was performed and combined with transcriptome data, resulting in the identification of six candidate factors: ESRRB, ESRRG, FOS, NR5A1, NR5A2, and SP5. These factors may play roles in regulating XX-DSD gonadal abnormalities. Among them, NR5A1 has been implicated in regulating the transcription of key genes involved in sex differentiation and gonadal development [48], crucial for ovarian development and maintenance of reproductive function [49]. Deletion of NR5A1 results in incomplete maturation of sertoli cells, thereby affecting the development and function of male reproductive organs [50]. Additionally, previous studies have found differences the FOS gene in the hypothalamic tissue of XX-DSD pigs and normal females [51]. During mouse testis development, FOS gene expression levels were observed to change [52], FOS is involved in regulating the proliferation and function of Leydig cells [53]. Its overexpression may negatively impact ovarian cells, leading to abnormal follicle development [54]. In conclusion, these genes and pathways warrant further attention, with their interactions with TFs and regulation in abnormal gonadal development being a focal point for future studies.

Conclusions

In this study, combined ATAC-seq and RNA-seq analyses revealed that SOX9, COL1A1, COL1A2, FDX1, COL6A1, HSD3B1, FSHR, and CYP17A1 genes are pivotal in the development of XX-DSD pigs sex characteristics, associated with pathways involved in steroid hormone biosynthesis, gonadal development, and hormone responses. Additionally, the footprint analysis identified six transcription factors: NR5A2, ESRRB, ESRRG, NR5A1, FOS, and SP5, which play important roles in regulating gene expression during gonadal dysgenesis. The findings of this study lay the foundation for further elucidating the molecular mechanisms underlying gonadal abnormalities in XX-DSD pigs. Future research will focus on identifying the factors that activate SOX9 expression during the prenatal period, aiming to comprehensively understand the pathogenesis.

Data availability

The original contributions presented in the study are publicly available and The data that support the findings in this study have been deposited into CNGB Sequence Archive (CNSA) of China National GeneBank DataBase (CNGBdb) with accession number CNP0006005.

References

  1. Reyes AP, León NY, Frost ER, Harley VR. Genetic control of typical and atypical sex development. Nat Rev Urol. 2023;20(7):434–51.

    Article  PubMed  Google Scholar 

  2. Wu J, Yu H, Zhang Y, Zhao H, Zhong B, Yu C, Feng Z, Yu H, Li H. Pathological characteristics of SRY-negative 38,XX-DSD pigs: A family case report. Anim Reprod Sci. 2024;270:107579.

    Article  CAS  PubMed  Google Scholar 

  3. Kornberg RD, Lorch Y. Primary role of the nucleosome. Mol Cell. 2020;79(3):371–5.

    Article  CAS  PubMed  Google Scholar 

  4. Liu Y, Zhou K, Zhang N, Wei H, Tan YZ, Zhang Z, Carragher B, Potter CS, D’Arcy S, Luger K. FACT caught in the act of manipulating the nucleosome. Nature. 2020;577(7790):426–31.

    Article  CAS  PubMed  Google Scholar 

  5. Klemm SL, Shipony Z, Greenleaf WJ. Chromatin accessibility and the regulatory epigenome. Nat Rev Genet. 2019;20(4):207–20.

    Article  CAS  PubMed  Google Scholar 

  6. Wu J, Tan S, Feng Z, Zhao H, Yu C, Yang Y, Zhong B, Zheng W, Yu H, Li H. Whole-genome de Novo sequencing reveals genomic variants associated with differences of sex development in SRY negative pigs. Biol Sex Differ. 2024;15(1):68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Zhang X, Li J, Wang X, Jie Y, Sun C, Zheng J, Li J, Yang N, Chen S. ATAC-seq and RNA-seq analysis unravel the mechanism of sex differentiation and infertility in sex reversal chicken. Epigenet Chromat. 2023;16(1):2.

    Article  CAS  Google Scholar 

  8. Raurell-Vila H, Ramos-Rodríguez M, Pasquali L. Assay for transposase accessible chromatin (ATAC-Seq) to chart the open chromatin landscape of human pancreatic Islets. Methods Mol Biol. 2018;1766:197–208.

    Article  CAS  PubMed  Google Scholar 

  9. Cinkornpumin JK, Hossain I, Pastor WA. Mapping chromatin accessibility in human Naïve pluripotent stem cells using ATAC-Seq. Methods Mol Biol. 2022;2416:201–11.

    Article  CAS  PubMed  Google Scholar 

  10. Wu J, Huang B, Chen H, Yin Q, Liu Y, Xiang Y, Zhang B, Liu B, Wang Q, Xia W, et al. The landscape of accessible chromatin in mammalian preimplantation embryos. Nature. 2016;534(7609):652–7.

    Article  CAS  PubMed  Google Scholar 

  11. Wu J, Xu J, Liu B, Yao G, Wang P, Lin Z, Huang B, Wang X, Li T, Shi S, et al. Chromatin analysis in human early development reveals epigenetic transition during ZGA. Nature. 2018;557(7704):256–60.

    Article  CAS  PubMed  Google Scholar 

  12. Lu F, Liu Y, Inoue A, Suzuki T, Zhao K, Zhang Y. Establishing chromatin regulatory landscape during mouse preimplantation development. Cell. 2016;165(6):1375–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Xu R, Li C, Liu X, Gao S. Insights into epigenetic patterns in mammalian early embryos. Protein Cell. 2021;12(1):7–28.

    Article  PubMed  Google Scholar 

  14. Garcia-Moreno SA, Futtner CR, Salamone IM, Gonen N, Lovell-Badge R, Maatouk DM. Gonadal supporting cells acquire sex-specific chromatin landscapes during mammalian sex determination. Dev Biol. 2019;446(2):168–79.

    Article  CAS  PubMed  Google Scholar 

  15. Zhanpeng S, Jingjing L, Li L, Yifei G, Bin W, Tong H: Integration of ATAC-seq and RNA-seq identifies active G-protein coupled receptors functioning in molting process in muscle of Eriocheir sinensis. Frontiers in Marine Science 2022, 9.

  16. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. 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(9):1650–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Love MI, Huber W, Anders S. Moderated Estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Robinson MD, McCarthy DJ, Smyth GK. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  22. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for chip peak annotation, comparison and visualization. Bioinformatics. 2015;31(14):2382–3.

    Article  CAS  PubMed  Google Scholar 

  25. Stark R, Brown G. DiffBind: differential binding analysis of ChIP-Seq peak data. R Package Version 2011, 100(4.3).

  26. Xie Y, Wu C, Li Z, Wu Z, Hong L. Early gonadal development and sex determination in mammal. Int J Mol Sci. 2022;23(14):7500.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Hiramatsu R, Matoba S, Kanai-Azuma M, Tsunekawa N, Katoh-Fukui Y, Kurohmaru M, Morohashi K, Wilhelm D, Koopman P, Kanai Y. A critical time window of Sry action in gonadal sex determination in mice. Development. 2009;136(1):129–38.

    Article  CAS  PubMed  Google Scholar 

  28. Vainio S, Heikkilä M, Kispert A, Chin N, McMahon AP. Female development in mammals is regulated by Wnt-4 signalling. Nature. 1999;397(6718):405–9.

    Article  CAS  PubMed  Google Scholar 

  29. Okashita N, Maeda R, Tachibana M. CDYL reinforces male gonadal sex determination through epigenetically repressing Wnt4 transcription in mice. Proc Natl Acad Sci U S A. 2023;120(20):e2221499120.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Maatouk DM, DiNapoli L, Alvers A, Parker KL, Taketo MM, Capel B. Stabilization of beta-catenin in XY gonads causes male-to-female sex-reversal. Hum Mol Genet. 2008;17(19):2949–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Wu J, Tan S, Zhou Y, Zhao H, Yu H, Zhong B, Yu C, Wang H, Yang Y, Li H, et al. Clinical and gonadal transcriptome analysis of 38,XX disorder of sex development pigs†. Biol Reprod. 2024;111(1):212–26.

    Article  PubMed  Google Scholar 

  32. Miller WL, Auchus RJ. The molecular biology, biochemistry, and physiology of human steroidogenesis and its disorders. Endocr Rev. 2011;32(1):81–151.

    Article  CAS  PubMed  Google Scholar 

  33. Kalra S, Baruah M, Unnikrishnan AG, Sahay R. Publication trends in the Indian journal of endocrinology and metabolism. Indian J Endocrinol Metab. 2011;15(1):27–30.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Unal E, Yıldırım R, Taş FF, Tekin S, Ceylaner S, Haspolat YK. A rare cause of delayed puberty in two cases with 46,XX and 46,XY karyotype: 17 α-hydroxylase deficiency due to a novel variant in CYP17A1 gene. Gynecol Endocrinol. 2020;36(8):739–42.

    Article  CAS  PubMed  Google Scholar 

  35. Zhang D, Yao F, Luo M, Wang Y, Tian T, Deng S, Tian Q. Clinical characteristics and molecular etiology of partial 17α-hydroxylase deficiency diagnosed in 46,XX patients. Front Endocrinol (Lausanne). 2022;13:978026.

    Article  PubMed  Google Scholar 

  36. Yin M, Yang J, Tian Q, Zhang X. Ovarian gonadoblastoma with dysgerminoma in a Girl with 46,XX karyotype 17a-hydroxylase/17, 20-lyase deficiency: A case report and literature review. Front Endocrinol (Lausanne). 2022;13:989695.

    Article  PubMed  Google Scholar 

  37. Kobayashi S, Sata F, Ikeda-Araki A, Miyashita C, Itoh S, Goudarzi H, Iwasaki Y, Mitsui T, Moriya K, Shinohara N, et al. Associations among maternal perfluoroalkyl substance levels, fetal sex-hormone enzymatic gene polymorphisms, and fetal sex hormone levels in the Hokkaido study. Reprod Toxicol. 2021;105:221–31.

    Article  CAS  PubMed  Google Scholar 

  38. Dickson A, Yutuc E, Thornton CA, Dunford JE, Oppermann U, Wang Y, Griffiths WJ. HSD3B1 is an oxysterol 3β-hydroxysteroid dehydrogenase in human placenta. Open Biol. 2023;13(5):220313.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Tseng PW, Wu GC, Kuo WL, Tseng YC, Chang CF. The ovarian transcriptome at the early stage of testis Removal-Induced Male-To-Female sex change in the protandrous black porgy Acanthopagrus schlegelii. Front Genet. 2022;13:816955.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Lu Y, Zhang S, Wang Y, Ren X, Han J. Molecular mechanisms and clinical manifestations of rare genetic disorders associated with type I collagen. Intractable Rare Dis Res. 2019;8(2):98–107.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Vining B, Ming Z, Bagheri-Fam S, Harley V. Diverse regulation but conserved function: SOX9 in vertebrate sex determination. Genes (Basel). 2021;12(4):486.

    Article  CAS  PubMed  Google Scholar 

  42. Gao J, Wang Y, Liu J, Chen F, Guo Y, Ke H, Wang X, Luo M, Fu S. Genome-wide association study reveals genomic loci of sex differentiation and gonadal development in Plectropomus leopardus. Front Genet. 2023;14:1229242.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Jakob S, Lovell-Badge R. Sex determination and the control of Sox9 expression in mammals. Febs J. 2011;278(7):1002–9.

    Article  CAS  PubMed  Google Scholar 

  44. Gonen N, Lovell-Badge R. The regulation of Sox9 expression in the gonad. Curr Top Dev Biol. 2019;134:223–52.

    Article  CAS  PubMed  Google Scholar 

  45. López-Hernández B, Méndez JP, Coral-Vázquez RM, Benítez-Granados J, Zenteno JC, Villegas-Ruiz V, Calzada-León R, Soderlund D, Canto P. Duplication of SOX9 associated with 46,XX ovotesticular disorder of sex development. Reprod Biomed Online. 2018;37(1):107–12.

    Article  PubMed  Google Scholar 

  46. Terao M, Ogawa Y, Takada S, Kajitani R, Okuno M, Mochimaru Y, Matsuoka K, Itoh T, Toyoda A, Kono T, et al. Turnover of mammal sex chromosomes in the Sry-deficient Amami spiny rat is due to male-specific upregulation of Sox9. Proc Natl Acad Sci U S A. 2022;119(49):e2211574119.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Hirata Y, Mizushima S, Mitsukawa S, Kon M, Kuroki Y, Jogahara T, Shinohara N, Kuroiwa A. Identification of a new enhancer that promotes Sox9 expression by a comparative analysis of mouse and Sry-Deficient Amami spiny rat. Cytogenet Genome Res. 2023;163(5–6):307–16.

    Article  CAS  PubMed  Google Scholar 

  48. Yang Y, Workman S, Wilson M. The molecular pathways underlying early gonadal development. J Mol Endocrinol 2018 Jul 24:JME–17.

  49. Luppino G, Wasniewska M, Coco R, Pepe G, Morabito LA, Li Pomi A, Corica D, Aversa T. Role of NR5A1 gene mutations in disorders of sex development: molecular and clinical features. Curr Issues Mol Biol. 2024;46(5):4519–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Kato T, Esaki M, Matsuzawa A, Ikeda Y. NR5A1 is required for functional maturation of Sertoli cells during postnatal development. Reproduction. 2012;143(5):663–72.

    Article  CAS  PubMed  Google Scholar 

  51. Tan S, Zhou Y, Zhao H, Wu J, Yu H, Yang Y, Yang Y, Zhao H, Li H. Comprehensive transcriptome analysis of hypothalamus reveals genes associated with disorders of sex development in pigs. J Steroid Biochem Mol Biol. 2021;210:105875.

    Article  CAS  PubMed  Google Scholar 

  52. Shalini S, Bansal MP. Role of selenium in spermatogenesis: differential expression of Cjun and Cfos in tubular cells of mice testis. Mol Cell Biochem. 2006;292(1–2):27–38.

    Article  CAS  PubMed  Google Scholar 

  53. Sekido R, Lovell-Badge R. Genetic control of testis development. Sex Dev. 2013;7(1–3):21–32.

    Article  CAS  PubMed  Google Scholar 

  54. Richards JS. Hormonal control of gene expression in the ovary. Endocr Rev. 1994;15(6):725–51.

    Article  CAS  PubMed  Google Scholar 

Download references

Funding

This work was supported by the Guangdong Provincial Natural Science Foundation of China (2020B1515120057), the National Natural Science Foundation of China (32072699), and the Key Technologies R&D Program of Guangdong Province (2022B0202090001).

Author information

Author notes

  1. Congying Yu, Bingzhou Zhong these authors contributed equally to this work.

Authors and Affiliations

  1. Guangdong Provincial Key Laboratory of Animal Molecular Design and Precise Breeding, Animal Science and Technology College, Foshan University, Foshan, Guangdong, 528225, China

    Congying Yu, Bingzhou Zhong, Haiquan Zhao, Jinhua Wu, Hui Yu & Hua Li

  2. School of Biological Sciences, Crawley (Perth), The University of Western Australia, 35 Stirling Highway, Perth, WA, 6009, Australia

    Haiyi Yu

  3. Zhongshan Baishi Pig Farm Co., Ltd.r, Zhongshan, 528463, China

    Yuqiao Zhang

Authors
  1. Congying Yu
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Bingzhou Zhong
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Yuqiao Zhang
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Haiquan Zhao
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Jinhua Wu
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Haiyi Yu
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Hui Yu
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Hua Li
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

HL and HY conceived and designed the project. CYY, BZZ, YQZ, HQZ, JHW and HYY performed the experiments and collected samples. CYY and BZZ analyzed the data and wrote the manuscript. All authors read and approved the final draft.

Corresponding authors

Correspondence to Hui Yu or Hua Li.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the Animal Care Committee of Foshan University (Foshan, China; No.682167). Euthanasia was performed using intravenous pentobarbital sodium (> 150 mg/kg). All animal experiments were performed following Chinese National Guidelines for Experimental Animal Welfare.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yu, C., Zhong, B., Zhang, Y. et al. Combining ATAC-seq and RNA-seq reveals key genes for gonadal abnormalities in one-month-old XX-DSD pigs. BMC Genomics 26, 447 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11613-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11613-x

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us