The experimental goats were all from the Jining Grey Goats Breeding Farm (Jiaxiang, Shandong, China). Under the same feeding management conditions, 20 healthy and disease-free female Jining grey goats were selected. The selected goats were divided into four groups according to age: 1 day old (D1, n = 5; body weight (BW): 2.08 ± 0.11 kg), 2 months old (M2, n = 5; BW: 4.42 ± 0.24 kg), 4 months old (M4, n = 5; BW: 7.62 ± 0.50 kg), and 6 months of age (M6, n = 5; BW: 8.82 ± 0.53 kg). The body condition of Jining grey goats was similar in each group. The experimental goats were slaughtered on the same day, and after the electric shock, the hypothalamic tissue was quickly slaughtered and collected, and stored at −80°C.

Total RNA was extracted from 20 hypothalamic tissues using TRIzol^® reagent (Thermo Fisher Scientific, Waltham, MA, United States). Screening of qualified RNA samples for RNA strand-specific library construction. The rRNA was removed from total RNA samples using the Ribo-Zero rRNA Removal Kit (Illumina, Inc., San Diego, United States), and then a sequencing library was generated using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB E7420) for Illumina to generate sequencing libraries. The enriched RNA was fragmented using a fragmentation buffer to yield small fragments. Then, the fragmented RNA served as a template for reverse transcription with the addition of 6 bp random primers (random hexamers) to synthesize the first cDNA strand. This was followed by the addition of buffer, dNTPs (with dTTP replaced by dUTP), DNA polymerase I, and RNase H to synthesize the second cDNA strand. The synthesized double-stranded cDNA was purified and enriched via PCR. The PCR product was then purified to obtain the final strand-specific library. After reverse transcription and PCR amplification, 150 bp paired-end reads were sequenced using the Illumina Novaseq6000 platform (Illumina, Inc., San Diego, United States).

To obtain high-quality sequencing data, we utilized Fastp (v0.23.1) to eliminate sequences containing poly-N, low-quality reads, and adapters from the obtained sequencing data. The high-quality reads obtained are used for downstream data analysis. We generated an index of the reference genome by employing HISAT2 (v2.0.5.) Subsequently, we aligned the clean reads with the goat reference genome (GCF_001704415.2_ARS1.2) using HISAT2 (20). Transcript assembly is performed using Stringtie (v1.3.3b), and gene expression levels are calculated (21). Gene expression levels were normalized using fragments per kilobase of exon model per million mapped reads (FPKM).

The novel lncRNAs were identified in hypothalamic tissue following the steps shown in Figure 1: (1) removing transcripts with an exon number of 1, (2) removing transcripts less than 200 nt in length, and (3) screening out transcripts that overlapped the exon region annotated in the database by gffcompare software (v0.10.6) (22); (4) CPC2 (v3.2.0) (23), Pfam (v1.6) (24), and CNCI (v2.0) (25) were used to predict the encoding potential of lncRNAs, and the transcripts that were predicted in the three software without coding potential were intersected, and (5) the low-expression lncRNAs (FPKM <0.5) were filtered out.

The DEseq2 (v1.20.0) package was utilized to examine the differential expression of lncRNAs (DELs) across various developmental stages. The readcounts from the sequencing data were used as the input matrix, and the p-value was adjusted utilizing the Benjamini & Hochberg method. The DELs were screened according to the threshold |log₂ (Fold change)| ≥1 and False Discovery Rate (FDR) <0.05. Cluster analysis of FPKM values of lncRNAs was performed using the ggplot2 package (v3.4.4). Data for lncRNAs were normalized [log₂(X + 1)] and then standardization (z-score).

lncRNAs can regulate the expression of potential target genes through cis-or trans-regulatory methods. Based on the location information of lncRNAs, mRNAs within 100 kb upstream and downstream of lncRNAs are defined as cis-target genes of lncRNAs (26). There will be a significant correlation with lncRNA expression (|R| > 0.95 and p < 0.05) is defined as a potential trans-target gene for lncRNA.

Subsequently, we used clusterProfiler software (v3.8.1) to perform Gene Ontology (GO) functional analysis of these differentially lncRNAs predicted target genes (27). The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis is performed by KOBAS (http://bioinfo.org/kobas; Accessed: 12.19, 2023) (28). p < 0.05 was considered to be significantly enriched.

Six lncRNAs were randomly selected and the accuracy of the lncRNA sequencing results was verified by qRT-PCR. First, we used PrimeScript^™II First strand cDNA synthesis kit (Takara Bio Inc., Dalian) to reverse transcrib total RNA from goat hypothalamus tissue into cDNA. Then qRT-PCR was performed in a Roche LightCycler 96 using the SYBR PrimeScript^™ RT-PCR Kit (Takara Bio Inc., Dalian). GAPDH was used as an internal reference gene to correct gene expression levels and normalize the data. The primers designed using Primer 5.0 (Supplementary Table S1). The relative expression levels of lncRNAs were calculated by 2^−ΔΔCT (29). One-way ANOVA was performed with SPSS 17.0, and the results were expressed as mean ± standard error. Three repetitions are performed for each set. p < 0.05 was considered statistically significant.

RNA was isolated from the hypothalamic tissues of female Jining grey goats at four developmental stages (D1, M2, M4, and M6), and 20 lncRNA libraries were constructed. Sequencing of the libraries was conducted on the Noveseq 6000 platform, resulting in a total of 1,803,310,692 raw reads. Following quality control procedures, we obtained 1,764,056,912 clean reads (Supplementary Table S2). The alignment of these clean reads to the goat reference genome was performed using HiSAT2, achieving alignment rates ranging from 86.05 to 95.99%, with unique mapping reads alignment rates between 75.1 and 91.98% (Supplementary Table S2). Subsequent analysis only considered the uniquely mapped reads.

A total of 10,630 lncRNAs were identified according to the steps shown in Figure 1, of which 1,676 were annotated and 8,954 lncRNAs were newly identified. Cluster analysis showed that most of the lncRNAs were expressed at low levels at the D1 stage (Supplementary Figure S1). Further analysis of the identified lncRNA signatures showed that showed that about 62% of the lncRNAs had 2 exons, and a few lncRNAs (3%) had more than 6 exons (Figure 2A). In addition, the length distribution of the identified lncRNAs ranged from 122 to 73,429 bp. More than 50% of the lncRNAs were less than 1,000 bp in length, about 83% were in the 0–3 kb range, and a few (17%) were greater than 3 kb (Figure 2B). About 35.2% of the lncRNAs are located in the intergenic region. Only 12.4% of the lncRNAs are from the antisense region, and about 32.9% of the lncRNAs are from the intron region (Figure 2C).

According to the screening criteria of |log₂ (Fold change)| ≥1 and p_adj < 0.05, the lncRNAs of four different developmental stages (D1, M2, M4 and M6) in hypothalamic tissue were compared and analysed. This analysis resulted in the identification of 237 differentially expressed lncRNAs (Figure 3A).

A total of 135 DELs (31 up-regulated and 104 down-regulated) were identified in M2 vs. D1, 68 DELs (38 up-regulated and 30 down-regulated) were identified in M4 vs. M2, and 16 DELs (11 up-regulated and 5 down-regulated) were identified in M6 vs. M4. A total of 53 DELs (16 up-regulated and 37 down-regulated) were identified in M4 vs. D1, 33 DELs (15 up-regulated and 18 down-regulated) were identified in M6 vs. D1, and 69 DELs (48 up-regulated and 21 down-regulated) were identified in M6 vs. M2. Interestingly, Figure 3B shows different expression patterns for these DELs, and Figure 3C shows that the M2 vs. D1 group has the most unique DELs, with 72. The M6 vs. M4 group had the lowest number of DELs with 3. Remarkably, 9 DELs were the same in the three comparison groups of M2 vs. D1, M4 vs. D1, and M6 vs. D1. This suggests that TCONS_00076225, TCONS_00148767, TCONS_00191611, TCONS_00102342, TCONS_00192735, TCONS_00032355, TCONS_00032357, TCONS_00053700, and TCONS_00070238 may be significant in the postnatal sexual development of goats.

lncRNAs have been implicated in influencing gene expression through cis-or trans-interactions. To investigate the role of lncRNAs in the goat hypothalamus, we predicted the potential regulatory roles of the identified DELs on both cis and trans target genes. Specifically, based on a distance threshold of 100 kb between lncRNAs and their target genes, 221 DELs were predicted to regulate 693 target genes in a cis manner (Supplementary Table S3).

The Gene Ontology (GO) analysis results showed that the target genes were significantly enriched in 50 categories (Supplementary Table S4) These cis-target genes are involved in many biological processes, such as regulation of catalytic activity, regulation of molecular function, regulation of hydrolase activity, and dephosphorylation (Figure 4A). In addition, KEGG analysis showed that these target genes were significantly enriched in sphingolipid signaling pathway, protein processing in endoplasmic reticulum, progesterone-mediated oocyte maturation, adipocytokine signaling pathway, neurotrophin signaling pathway, oocyte meiosis, sphingolipid metabolism, glutamatergic synapse, cGMP-PKG signaling pathway, inositol phosphate metabolism, phospholipase D signaling pathway and other pathways (p < 0.05) (Figure 4B and Supplementary Table S5).

According to |R| > 0.95 and p < 0.05, 24 lncRNAs were found, which had 63 trans target genes (Supplementary Table S6 and Supplementary Figure S2). The number of trans-target genes identified in lncRNAs was significantly lower than that in cis, suggesting that lncRNAs may function mainly by cis-regulating gene expression. The GO analysis results showed that we found that these cis-target genes were involved in many biological processes, such as signal receptor binding, hormone activity, NADH dehydrogenase activity, transferase activity, transfer aminoacyl, etc. (Figure 4C and Supplementary Table S7). In addition, KEGG analysis showed that these target genes were significantly enriched in neuroactive ligand receptors, multiple apoptosis, adipocytokine signaling pathway, GnRH signaling pathway, JAK-STAT signaling pathway, and p53 signaling pathway (p < 0.05) (Figure 4D and Supplementary Table S8).

The results suggest that these lncRNAs may be involved in hormone secretion, signal transduction processes, and thus regulation of sexual maturation in goats by modulating cis and trans target genes. Interestingly, LOC108634846, LOC108635405, LOC108637396, AGRP, and PIK3C2G were predicted in both cis and trans target genes. Among these, PIK3C2G is the trans-target gene and the cis-target gene of TCONS_00038560.

To better understand the role of goat hypothalamic DELs in the process of sexual maturation, we selected the target gene of DELs involved in reproduction. Construction of lncRNA-mRNA regulatory networks for DELs and their cis-and trans-target genes, respectively (Figure 5). A total of 31 lncRNAs regulated 34 target genes.

To verify the accuracy of the sequencing results, we randomly selected 6 lncRNAs (XR_001297374.2, XR_001917241.1, TCONS_00176496, XR_001919854.1, TCONS_00032357, and TCONS_00074891) for qRT-PCR detection. The results showed that the expression patterns of these lncRNAs and those found in the transcriptome data were consistent with the sequencing results (Figure 6). This further illustrates the high reliability and accuracy of RNA-seq.