Intersex goats show different gene expression levels in the hypothalamus and pituitary compared with non-intersex goats based on RNA-Seq.

The conditions for sex reversal in vertebrate species have been extensively studied, and the results highlighted numerous key factors involved in sex differentiation. However, the transcriptomes in hypothalamic and pituitary tissues from intersex goats have rarely been studied. The aim of this study was to screen candidate genes and signalling pathways related to sex reversal in Huai goats by analyzing gene expression in hypothalamic and pituitary tissues via transcriptome sequencing and bioinformatics analyses. In total, 612 and 139 differentially expressed genes (DEGs) were identified between the intersex and non-intersex groups in the hypothalamus and pituitary, respectively. The DEGs in the hypothalamus and pituitary were significantly enriched in 41 and 16 signalling pathways, respectively, including the calcium signalling pathway, neuroactive ligand-receptor interaction signalling pathway, and oestrogen signalling pathway, which might be related to intersex sex development disorders. A candidate gene from the tachykinin family (TACR1) was significantly enriched in the calcium signalling pathway. Thirty-one DEGs were shared between these two comparisons and were enriched in several acetyl-CoA-related processes and the oestrogen signalling pathway. The results of the real-time PCR analysis show that the transcriptome sequencing results were reliable. The transcriptome data indicate that the regulation of various physiological systems is involved in intersex goat development. Therefore, these results provide helpful data enhancing our understanding of the molecular mechanisms underlying intersex syndrome in goats.

INTRODUCTION

Cases of intersex development have been described in several domestic species, including pigs, dogs and goats (Cribiu & Chaffaux 1990). In goats, intersexuality refers to a group of developmental abnormalities of the sexual system ranging from only subtle changes to complete sex reversal, consequently preventing breeding improvement and reproduction. Overall, the prevalence of intersexuality in goats is reportedly 5%–6% (McEntee 1990). Intersex goats are anatomically classified into five types based on their whole genital organs.

Polledness intersex syndrome (PIS) is a phenomenon associated with polledness and intersexuality. The genetic basis of PIS syndrome has been identified as a deletion of 11.7 kb of DNA on goat chromosome 1q43 (Pailhoux et al., 2001; Schibler et al., 2000; Zhang et al., 2018). An in‐depth analysis revealed that the PIS interval contains no known coding sequences; this interval is known to be involved in the transcriptional regulation of two adjacent genes (PISRT1 and FOXL2) in goats (Crisponi et al., 2001), suggesting that a complex regulatory mechanism underlies intersexual development in goats. PISRT1 is known to enhance FOXL2 expression in females, and PISRT1 insufficiency in PIS(–/–) female goats is associated with significantly reduced FOXL2 expression (Boulanger et al., 2008). The transcriptional activity of the SOX9 gene is clearly associated with the functional state of Sertoli cells. These cells are completely degraded, and the SOX9 expression level is strongly reduced in 12‐month‐old females with PIS (Szatkowska et al., 2014). A transcriptome study investigating differentially expressed genes (DEGs) in pituitary tissues from intersexual and non‐intersex goats revealed a total of 10,063 DEGs. The DEGs were clustered into 56 gene ontology (GO) categories and determined to be significantly enriched in 53 Kyoto Encyclopaedia of Genes and Genomes (KEGG) signalling pathways (E et al., 2019).

The phenotypic diversity of goats with PIS is most likely the result of the complexity of the sexual differentiation process in which numerous molecular factors may potentially play a crucial role in regulating development (Gao et al., 2018). Thus, a genetic analysis of intersex goats from multiple perspectives is necessary to better understand the molecular genetics underlying sex development and reversal. In this study, we comparatively analyzed DEGs and their functions in pituitary and hypothalamic tissues from intersexual and non‐intersex goats using transcriptome sequencing (RNA‐Seq).

MATERIALS AND METHODS

Tissue collection

The Huai goats used in this experiment were obtained from the Jierui Huai goat conservation farm, Shenqiu County, Henan Province, China. The hircine hypothalamic and pituitary tissue samples were obtained from five intersexual (hypothalamus, XJ1‐5; pituitary, CJ1‐5) and five non‐intersex (hypothalamus, XN1‐5; pituitary, CN1‐5) Huai goats aged > 1 year under the approval of the University Institutional Animal Care and Use Committee. All fresh tissue samples were collected, divided into 1.5 ml plastic centrifuge tubes (each sample weighed approximately 100 mg), rapidly frozen in liquid nitrogen and stored at −80°C.

RNA extraction

The total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. The RNA concentration and purity were assessed by measuring the optical density at 260/280 nm using a Thermo Scientific NanoDrop Lite Spectrophotometer. Then, equal amounts of high‐quality RNA samples from each tissue were used for the cDNA library construction and Illumina sequencing.

cDNA library construction and sequencing

RNA libraries from each tissue sample were prepared following the Illumina mRNA‐Seq protocol by Nanjing Personal Gene Technology Co., Ltd. (Nanjing, China). In total, 3 μg RNA per sample were used to construct a cDNA library. The library quality was assessed using an Agilent Bioanalyzer 2100 system. The cDNA library was sequenced using the Illumina HiSeq2000 platform (San Diego, CA, USA). The raw RNA‐Seq data of the 10 pituitary and hypothalamic tissues are available in the NCBI Sequences Read Archive under accession number SRR13302433∼13302452.

Raw read filtering, quality control and mapping to the reference genome

The raw reads were cleaned by removing adaptors, poly‐N‐containing reads and reads of low quality using Cutadapt. A quality control analysis was performed by using FastQC with the default parameters. The clean reads were aligned to the annotated reference genome (Capra_hircus.ARS1.dna.toplevel.fa) from the Ensembl database using HISAT2 (the updated version of TopHat) software (http://ccb.jhu.edu/software/hisat2/index.shtml) with the BWT algorithm (Kim et al., 2013). Bowtie2 (http://bowtie‐bio.sourceforge.net/bowtie2/manual.shtml) was used to align the reads to the transcriptome.

Analysis of DEGs

The read counts were obtained from the mapping results obtained by using HTSEquation (Anders et al., 2015) and performing genome annotation (ftp://ftp.ensembl.org/pub/release‐73/gtf/homo_sapiens/). The gene expression results were normalized by calculating the fragments per kilobase of exon per million fragments mapped (FPKM) value of each gene to eliminate the influence of different gene lengths and sequencing levels. To select expressed genes, a threshold of FPKM ≥ 1 was used. The FPKM values between the biological replicates were analyzed by a Pearson correlation analysis. The Pearson coefficient of gene expression in different replicates was greater than 0.80, indicating consistency between the replicates. A principal component analysis (PCA) was conducted to cluster the replicates in different comparisons using the DESeq2 R package 1.8.1 (Love et al., 2014).

The DESeq R package (1.10.1) was used to identify DEGs between the intersex and non‐intersex tissues. Genes were identified as significantly differentially expressed when they showed a p‐value < 0.05 and a |log2FoldChange| > 1. Pheatmap software was used to cluster the DEGs and samples based on Euclidean and complete linkage methods. A volcano plot of the DEGs identified in each comparison was generated and visualized with the ggplot2 package (Wickham 2009). The common and specific DEGs between two comparisons were identified by a Venn analysis.

Annotation of DEGs

The DEGs were classified according to the cellular component, molecular function and biological process categories through a GO analysis, which was carried out using Blast2go software (http://www.blast2go.com/b2ghome) in the GO database (http://www.geneontology.org/). KEGG pathway analyses were performed using the major public pathway‐related database (http://www.kegg.jp/). The enriched GO and KEGG categories were determined using a p‐value ≤ 0.05 as the cut‐off for significant GO and KEGG categories.

Detection of putative single‐nucleotide polymorphisms

Varscan software was used to perform single‐nucleotide polymorphism (SNP) and indel calling. The filtering thresholds were set as follows: consensus quality of no less than 20 and coverage of no less than 8. To obtain more reliable SNPs, those identified in at least two reads with a p‐value < 0.01 were regarded as SNPs.

Quantitative real‐time PCR validation

The samples used in the qPCR analyses were the same as those used in the RNA‐Seq experiment. To confirm the sequencing results, a real‐time PCR (qRT‐PCR) analysis of 20 randomly selected DEGs was performed. The total RNA was reverse transcribed into cDNA using one‐step gDNA removal and cDNA synthesis SuperMix (TransScript, Beijing, China) according to the manufacturer's instructions. Gene‐specific primers were designed using Primer Premier 5.0 (Premier Biosoft). The mRNA levels were quantified by a qRT‐PCR analysis using an ABI PRISM 7500 fast qRT‐PCR system (Applied Biosystem, CA, USA) and a qPCR detection kit (Top Green qPCR SuperMix, TransScript, Beijing, China). The PCR analyses were performed in 20 μl amplification reactions containing 10 μl of 2 × SuperMix, 0.4 μl of each primer (10 μM), 1 μl cDNA template, 0.4 μl of 50 × Passive Reference Dye and 7.8 μl RNase‐Free ddH2O. The reaction mixtures were subjected to initial denaturation at 94°C for 30 s, followed by 40 cycles of 94°C for 5 s and 60°C for 34 s. The final melting curve analysis was performed at 95°C for 15 s, 60°C for 1 min, 95°C for 30 s and 60°C for 15 s. The specificity of the amplification was confirmed by a melting curve analysis. All reactions were conducted in triplicate, and the cycle threshold (Ct) values of the target genes were normalized to the average value of the reference gene GAPDH. The fold changes in the expression of the target genes were calculated using the 2−ΔΔ formula, where ΔCt is the average Ct of the target gene minus the average Ct of two endogenous controls, and ΔΔCt is the ΔCt of the target sample (intersex tissue) minus the ΔCt of the calibrator sample (normal tissue). All data were subjected to tests of normality and homoscedasticity and a one‐way ANOVA, followed by Duncan's multiple range tests. p < 0.05 was considered significant.

RESULTS

Construction of the mRNA expression profile

We sequenced mRNA from hypothalamic and pituitary tissues collected from five intersex goats. The anatomical structures and histopathological features of ovarian tissues from these five goats can be found in our previous report (Yang et al., 2020). The transcriptome sequencing using the Illumina platform yielded an average of 42,985,856 and 46,921,660 clean reads from the pituitary glands of the intersex and normal goats, respectively. Additionally, the transcriptome sequencing yielded 46,502,082 and 46,734,866 clean reads from the hypothalamic tissues collected from the intersex and non‐intersex individuals, respectively (Table S1). The Q30 values of all samples were higher than 93%. We mapped the sequence reads to the reference Capra hircus genome sequence (Capra_hircus.ARS1.dna.toplevel.fa). The read alignment showed that on average, 97.09% of the reads (44,450,383 reads) were mapped to the reference genome. Among the mapped reads, 92.71% were mapped to unique positions (Table S1). Only the uniquely mapped reads were considered in the subsequent analyses.

Clustering analysis

A Pearson test was applied to the 10 hypothalamic and pituitary samples, and the obtained Pearson's correlation coefficients ranged from 0.9876 to 0.9988 among CJ1, CJ3, CJ4 and CJ5 and from 0.9851 to 0.9969 among CN1‐5, indicating high similarity across these biological repeats. However, the correlation coefficients between CJ2 and the other pituitary samples were lower than 0.45. Furthermore, the correlation coefficients among XJ1, XJ3, XJ4 and XJ5 ranged from 0.8809 to 0.9103, and XN1‐5 were also highly correlated with each other (0.9712–0.9962). XJ2 showed relatively lower correlation coefficients than the other samples. Therefore, the CJ2 and XJ2 samples were removed from the subsequent analysis.

A Pearson test, two‐way clustering analysis and PCA were carried out to verify the biological repeats among the remaining nine hypothalamic samples and nine pituitary samples. Our results showed that the correlation coefficients among the hypothalamic samples were higher than 0.91 (Figure S1A); the XJ samples were clustered together, and the XN samples were clustered together based on the identified DEGs (Figure S1B). According to the PCA plot (Figure S1C), the nine hypothalamic samples were highly correlated and showed different gene expression profiles. The nine pituitary samples showed high correlation coefficients (> 0.99) (Figure S2A). Four samples from the intersex goats and five samples from the normal goats clustered separately and showed different expression profiles (Figure S2B, S2C).

Differential gene expression analysis

The RNA‐Seq technique enabled the analysis of the differential expression profiles via an analysis of the transcript abundance with high sensitivity. In total, 21,341 genes were identified; of these genes, 53.765% (11,474 genes) were defined as expressed genes (FPKM > 1), and 2.01% (428 genes) were defined as highly expressed (FPKM > 100). Five genes, namely, ENSCHIG00000018565, ENSCHIG00000012599, ENSCHIG00000010087 (EEF1A1), ENSCHIG00000022860 (SEPW1) and ENSCHIG00000024457 (CFL1), showed the highest expression levels (FPKM > 1000) in all 20 samples.

When p < 0.05 and a |log2‐fold change|≥1 were used as the cut‐offs, in total, 612 and 139 DEGs were identified in the hypothalamus and pituitary between the intersex and normal groups, respectively. In the comparison of the pituitary tissues, 57 DEGs were found to be downregulated in the intersex tissues relative to the normal tissues, while 82 DEGs were upregulated (Figure 1A). In total, 323 and 289 transcripts were downregulated and upregulated in the hypothalamic tissues from the intersex goats, respectively, relative to those in the non‐intersex goats (Figure 1B). Many genes expressed in the hypothalamus and pituitary exhibited differential expression between the intersex and normal goats as shown in Figure S1B and Figure S2B, indicating that the expression of these genes might be associated with developmental abnormalities of the sexual system. There were 31 shared DEGs identified in both comparisons between the intersex and non‐intersex goats (Figure 1C and Table 1).

FIGURE 1

Volcano plot displaying the differentially expressed genes (DEGs) identified in goat hypothalamic and pituitary tissues and Venn diagram of the DEGs distributed in the two groups. Red dots: upregulated genes. Blue dots: downregulated genes. (a) DEGs in the pituitary between intersex and non‐intersex goats. (b) DEGs in the hypothalamus between intersex and non‐intersex goats. (c) Group 1: DEGs in the hypothalamus between intersex and non‐intersex goats. Group 2: DEGs in the pituitary between intersex and non‐intersex goats

TABLE 1

Differentially expressed genes between intersex and non‐intersex goats shared by the hypothalamus and pituitary

Gene ID	Gene name	Hypothalamus		Pituitary
		log2fc	pval	log2fc	pval
ENSCHIG00000013432	PDK4	−2.427812172	6.63581E‐22	−2.43753152	1.35777E‐15
ENSCHIG00000007606	ARRDC2	−2.165564754	1.12331E‐13	−1.704551636	4.61091E‐22
ENSCHIG00000013289	LOC102175225	2.167072919	9.64985E‐11	2.073546573	2.50607E‐06
ENSCHIG00000012197	RUBCNL	−1.369691591	4.16812E‐08	−1.157475233	1.25094E‐16
ENSCHIG00000010165	–	3.605710593	1.76249E‐05	2.664625694	0.008345249
ENSCHIG00000002934	FOXQ1	1.317826471	0.000125748	1.747153939	8.86787E‐10
ENSCHIG00000014249	DIO2	−2.137857783	0.000209338	−1.267670368	0.005920083
ENSCHIG00000018826	LOC102177708	−2.364483144	0.000213566	−2.203943343	0.009352443
ENSCHIG00000016684	UCP2	1.011722022	0.000254453	1.067308904	0.011238771
ENSCHIG00000017462	–	2.714793328	0.000438798	2.218960657	0.025606302
ENSCHIG00000008160	–	−5.319376219	0.000679637	−4.638708254	0.002026309
ENSCHIG00000014211	HLX	1.35266329	0.001401527	1.588401896	9.03093E‐05
ENSCHIG00000015588	IGFBP1	−3.162415185	0.001533424	−2.847145196	0.001199814
ENSCHIG00000014746	FKBP5	−1.524645718	0.001589511	−1.60966534	0.009906033
ENSCHIG00000025384	FOS	1.105257739	0.00169011	−1.083757283	2.64677E‐17
ENSCHIG00000018361	ID2	1.225944052	0.00188682	1.056150807	5.80544E‐06
ENSCHIG00000012647	ITPRIPL1	1.158646475	0.003034011	1.072796345	0.001956041
ENSCHIG00000015100	LOC102191549	1.936214538	0.003727347	1.2461325	0.022878987
ENSCHIG00000006029	–	4.989606762	0.006256121	5.042526633	0.005298725
ENSCHIG00000006431	–	1.274146263	0.014111449	1.644797567	5.97941E‐05
ENSCHIG00000021088	FAM180A	−1.309410109	0.014276896	−1.008126869	0.032094413
ENSCHIG00000001740	LOC102184922	1.615753031	0.01635515	1.537206406	8.71934E‐05
ENSCHIG00000015403	CXCL10	1.440830258	0.017145139	1.092424552	0.001416087
ENSCHIG00000004128	LOC108636739	4.03596049	0.019940391	2.822482144	0.034620481
ENSCHIG00000004176	ARL11	1.506977442	0.024701215	1.372397871	0.035745909
ENSCHIG00000015444	SCGN	−2.218495693	0.024792096	−1.186203164	0.015807027
ENSCHIG00000022202	LOC102173185	−2.802462903	0.026697706	−2.880216111	0.026985276
ENSCHIG00000006496	FAM124B	1.171587724	0.032447059	1.729987615	0.002257718
ENSCHIG00000024771	LOC102168680	2.220290941	0.03497696	2.362306307	0.029228719
ENSCHIG00000025411	HSPA6	−1.390678161	0.039093708	−1.777478103	0.000175822
ENSCHIG00000013114	FCN1	2.183268543	0.039114686	1.897399214	0.027603116

Note: log2fc, log2‐fold change; pval, p‐value.

Functional enrichment

GO and KEGG analyses were performed to further elucidate the biological functions of the 31 genes. These RNAs were summarized under the three main GO categories and were shown to be significantly enriched in 175 GO terms (p‐value < 0.05) and extremely significantly enriched in 34 terms (p‐value < 0.01). Figure 2A shows the top 20 significant GO terms. Protein folding (GO:0006457), thyroid hormone catabolic process (GO:0042404) and cysteinyl leukotriene receptor activity (GO:0001631) were the most significantly enriched terms. In particular, the GO directed acyclic graph (DAG) of biological processes indicated seven significantly connected GO terms, including thioester biosynthetic process (GO:0035384), acetyl‐CoA metabolic process (GO:0006084), acyl‐CoA biosynthetic process (GO:0071616), acetyl‐CoA biosynthetic process (GO:0006085), regulation of acyl‐CoA biosynthetic process (GO:0050812), acetyl‐CoA biosynthetic process from pyruvate (GO:0006086) and regulation of acetyl‐CoA biosynthetic process from pyruvate (GO:0010510) (Figures 2B, 2C).

FIGURE 2

GO analysis of the 31 common DEGs. (a) Top 20 significantly enriched gene ontology (GO) terms. (b) GO‐directed acyclic graph enrichment of biological processes. (c) An enlarged view of Figure 2B. Note: boxes in Figure 2C indicate the top 10 enriched terms

In addition, 58 signalling pathways were annotated. In total, seven (oestrogen signalling pathway, IL‐17 signalling pathway, Toll‐like receptor signalling pathway, TNF signalling pathway, natural killer cell‐mediated cytotoxicity, influenza A and African trypanosomiasis) of the 58 KEGG pathways were significantly enriched (p‐value < 0.05). However, only three or two enriched DEGs were involved in each pathway. Figure 3A presents the top 20 enriched KEGG pathways. The oestrogen signalling pathway might be associated with intersexuality in goats; in this pathway, the TF gene was upregulated, and the FKBP52 and HSP70 genes were downregulated (Figure 3B).

FIGURE 3

Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways of the 31 common DEGs. (a) Top 20 significantly enriched KEGG pathways. (b) A portion of the estrogen signaling pathway with three DEGs

Real‐time qPCR validation of the RNA‐Seq results

To validate the gene expression results obtained by sequencing, the expression levels of 20 DEGs (FDK4, HSPA6, FKBP5, RUBCNL, DIO2, ARRDC2, CXCL10, FOXQ1, ID2, UCP2, TRIB1, CA14, ZBTB16, GADD45B, BTBD17, SCN4B, DNER, ANXA6, OTOF and COCH) were analyzed through qPCR. The detailed information of the primers employed for the 20 genes is provided in Table S2. The log2‐fold change values determined by RNA‐Seq and the log2 values determined by qPCR through the 2–∆∆ method with normalization to GAPDH in the comparisons of these 20 genes are shown in Table S3. The expression levels of the 20 genes in all test samples and some validation results are displayed in Figure S3. The expression patterns of these 20 genes were consistent with the RNA‐Seq results according to the obtained Pearson correlation coefficients of 0.965712 (XJ vs. XN) and 0.9399 (CJ vs. CN) (Figure 4). The results indicate that the expression trends of these 20 genes according to the two detection methods were almost the same. The qPCR results revealed similar expression profiles, supporting the validity of the RNAseq results.

FIGURE 4

Correlation analysis of the expression patterns of 20 selected DEGs based on RT‐PCR and RNA‐Seq. The x‐axis represents the log2FoldChange value detected by RNA‐Seq, and the y‐axis represents the log2 value measured by the 2–∆∆ method using RT‐PCR with GAPDH as the reference. A represents the XJ vs. XN group, B represents the CJ versus LN group

SNPs associated with intersexuality in goats

In total, 569,740 SNPs were called from all samples, including 89,328 SNPs located in exonic regions. Here, 234 SNPs, including 61 nonsynonymous mutations and 273 synonymous mutations, were detected in ten intersex goat tissues (five hypothalamus and five pituitary) but not in normal goat tissues. Combined with the SNP information from the ovary that we previously reported (Yang et al., 2020), in total, 10 positions were found to be mutated in 15 tissues from intersex goats but not in tissues from nonintersex goats (Table S4). These mutations were located in exonic positions, three of which were nonsynonymous, while seven were synonymous.

DISCUSSION

In this study, the transcriptome datasets of DEGs in pituitary and hypothalamic tissues from intersex and non‐intersex goats were examined using high‐throughput technology. Considering only the pituitary results, in total, 139 DEGs were identified, and 503 significant GO terms and 16 KEGG pathways were annotated. Oestrogen signalling was significant and proven to be associated with gonad development. However, when E et al. (2019) analyzed the digital expression values of coding genes in pituitary tissues from intersex individuals compared to those from healthy goats, many more DEGs (10,063) were detected. Other pathways previously discovered to be related to sex development disorders include the oxytocin signalling, adrenergic signalling, oestrogen signalling and GnRH signalling pathways (Fisher et al. 2016). In addition, pathways, such as the axon guidance, FoxO signalling, MAPK signalling and cAMP pathways, were also identified but were not significantly enriched in our study and previous work (E et al., 2019).

In the comparison of hypothalamic tissues from intersex goats and non‐intersex goats, 612 DEGs were identified. In our study, these DEGs were significantly annotated to 1055 GO categories and were determined to be significantly enriched in 41 signalling pathways. Shi et al. (2019) investigated DEGs in hypothalamic tissue from Henan Huai goats maintained under natural or artificial illumination conditions and identified 448 DEGs enriched in five reproductive‐related signalling pathways (calcium signalling pathway, neuroactive ligand‐receptor interaction signalling pathway, GnRH signalling pathway, oxytocin signalling pathway and oestrogen signalling pathway). Interestingly, the neuroactive ligand‐receptor interaction signalling pathway and calcium signalling pathway were the most significantly enriched pathways identified in our study, involving 39 and 18 DEGs, respectively. The candidate gene TACR1 in the calcium signalling pathway was upregulated and significantly enriched. This finding is consistent with the results reported by Shi et al. (2019), indicating that the calcium signalling pathway and tachykinin family genes may play important roles in goat reproduction.

The hypothalamus, pituitary gland and gonads are very closely interrelated. Some hormones that they produce stimulate the development of the ovaries and testes. Seven pathways (IL‐17 signalling pathway, type II diabetes mellitus, signalling pathways regulating the pluripotency of stem cells, hematopoietic cell lineage, TNF signalling pathway, amphetamine addiction and malaria) were found to be enriched in both the hypothalamus and pituitary in the present study. Therefore, these pathways may play important roles in the developmental regulation of intersexuality. A hormone‐related signalling pathway, that is, the adrenergic signalling pathway in cardiomyocytes, was identified as significant in the hypothalamus in this study and in the pituitary according to the results reported by E et al. (2019). Fisher et al. (2016) previously suggested that adrenergic factors play an important role in sex development disorders. Additionally, the axon guidance and cAMP signalling pathways were significantly enriched in our study and that by E et al. (2019). In addition, some interesting pathways, including the GnRH and PPAR signalling pathways, were not significantly enriched but were identified in both our study and the study by E et al. (2019), supporting the assumption that the abnormal development of sexual organs may show an intrinsic relationship with hormones in intersex animals.

We further investigated the genes that were differentially expressed in both the hypothalamic and pituitary tissues in our study. A Venn diagram analysis was conducted and revealed 31 shared DEGs. As described above, the PIS locus affected the expression of FOXL2 and PISRT1 through transcriptional regulation (Pannetier et al., 2012). Recent studies have also noted that the silencing or uncontrolled expression of SOX9 and FOXL2 may affect gonadal differentiation (Boulanger et al., 2014; Veitia et al., 2001). For example, in mice, the SRY gene specifically upregulates SOX9 expression in the early stage of sex differentiation (Gonen et al., 2018). In goats and other mammals, SRY may inhibit anti‐testicular factors, such as PISRT1 and FOXL2, through the PIS region and indirectly act on SOX9 (Pannetier et al., 2003). The results of the absence of the SRY gene in five intersex samples and the DEGs identified in gonads have been published (Yang et al., 2020). We found that SOX9 was upregulated in both LJA and LJB, indicating that SOX9 can participate in the regulation of gonadal differentiation in goats (Yang et al., 2020). However, the expression of FOXL2 and PISRT1 was not detected in gonad, hypothalamic or pituitary tissues.

The 31 shared DEGs identified in our study were significantly enriched in 175 GO terms. In particular, several acetyl‐CoA‐related processes were connected and significantly enriched according to the GO analysis. Acetyl‐CoA is a key metabolic signalling molecule that influences various important cellular processes, such as gluconeogenesis, glucose oxidation, protein acetylation and steroid and fatty acid biosynthesis (Pietrocola et al., 2015). Sushma et al. (2021) found that silencing the acetyl‐CoA carboxylase beta (ACACB) gene reduced cholesterol synthesis. Cholesterol is a steroid biomolecule in animal cells that performs several important functions in living systems, such as steroid hormone (including androgen, oestrogen and progestin) production (Idoko et al., 2020). PDK4 is a lipid metabolism‐related gene that was downregulated and enriched in acetyl‐CoA‐related processes in our study. However, the functional mechanism of the PDK4 gene in steroid biosynthesis remains unclear. In the KEGG analysis, the oestrogen signalling pathway showed the highest significance. Our results suggest that the biosynthesis of gonadal hormones and their concentrations might be affected in intersex goats.

In our study, in total, 569,740 SNPs were called from all samples; of these SNPs, three mutations were nonsynonymous, but their relationships with intersexuality require further investigation. E et al. (2019) identified 171 SNPs that may contribute to intersexuality, and 53 of these SNPs were determined to be located in coding regions using a general linear model (E et al., 2019). There were no common mutations between the results of these two studies. In a previous report, a deletion of 11.7 kb on chromosome 1 was shown to contribute to PIS in goats (Pailhoux et al., 2001). However, in our research, no deletion fragment containing more than 50 bases was detected. E et al. (2019) also observed no significant correlated candidate site near this chromosomal region using GWAS. Another study demonstrated that the ASR1 goat assembly was based on the genome of a horned San Clemente goat (Bickhart et al., 2017). Thus, it can be assumed that the variant causing polledness was not present in this reference assembly (Simon et al., 2020). In summary, a wide range of genes that participate in determining the intersex phenotype have been identified. Our results could be helpful by providing more information explaining the physiology and development of intersex goats.

CONCLUSION

Our study compared gene expression profiles in hypothalamic and pituitary tissues from intersex and normal goats. Hundreds of DEGs were detected, and 31 genes were shared between these two sets of comparisons. Thirty of the 31 genes showed similar expression patterns and were enriched in the oestrogen signalling pathway and acetyl‐CoA biosynthetic process, suggesting that these genes might perform important functions in the regulation of hormone biosynthesis related to gonad development. Our results contribute to a new understanding of sex reversal and development.

CONFLICT OF INTEREST

The authors declare that they have no competing interests.

ETHICAL APPROVAL

The entire experimental procedure was approved by the Animal Care Commission of the College of Animal Science and Technology, Henan University of Animal Husbandry and Economy, China.

AUTHOR CONTRIBUTIONS

K.Q conceptualized and supervised the study. H.H, S.Y, J.L and X.W completed the whole experiments. H.H. drafted the manuscript. H.H. and J.Z. performed computational analysis. H.W and T.H. prepared experiment material. H.Y. and C.L. collected animals testes, W.L revised the manuscript. All authors read and approved the final manuscript.

SUPPORTING INFORMATION

Click here for additional data file.

SUPPORTING INFORMATION

Click here for additional data file.

SUPPORTING INFORMATION

Click here for additional data file.

SUPPORTING INFORMATION

Click here for additional data file.