Transcriptome Dynamics Reveal Stage-Specic and Melatonin-Triggered Gene Expression Patterns During The Cashmere Growth Cycle in Capra Hircus

Cashmere goat is famous for its high-quality bers. The growth of cashmere in secondary hair follicles exhibits a seasonal pattern arising from circannual changes in the natural photoperiod. Although several studies have compared and analyzed the differences in gene expression between different cashmere growth stages, the selection of samples in these studies relies on research experience or morphological evidence. Distinguishing cashmere growth cycle according to gene expression patterns may help to explore the regulation mechanisms related to cashmere growth and the effect of melatonin from a molecular level more accurately.


Abstract Background
Cashmere goat is famous for its high-quality bers. The growth of cashmere in secondary hair follicles exhibits a seasonal pattern arising from circannual changes in the natural photoperiod. Although several studies have compared and analyzed the differences in gene expression between different cashmere growth stages, the selection of samples in these studies relies on research experience or morphological evidence. Distinguishing cashmere growth cycle according to gene expression patterns may help to explore the regulation mechanisms related to cashmere growth and the effect of melatonin from a molecular level more accurately.

Results
In this study, we applied RNA-sequencing to the hair follicles of three normal and three melatonin-treated Inner Mongolian cashmere goats sampled every month during a whole cashmere growth cycle. A total of 3559 and 988 genes were subjected as seasonal changing genes (SCGs) in the control and treated groups, respectively. The SCGs in the normal group are divided into three clusters, and their speci c expression patterns help to group the cashmere growth cycle into anagen, catagen and telogen stages.
Some canonical pathways such as Wnt, TGF-beta and Hippo signaling pathways are detected as promoting the cashmere growth, while Cell adhesion molecules (CAMs), Cytokine-cytokine receptor interaction, Jak-STAT, Fc epsilon RI, NOD-like receptor, Rap1, PI3K-Akt, cAMP, NF-kappa B and many immune-related pathways are detected in the catagen and telogen stages. The PI3K-Akt signaling, ECMreceptor interaction and Focal adhesion are found in the transition stage between telogen to anagen, which may serve as candidate biomarkers for telogen-anagen regeneration. Pairwise comparisons between the control and melatonin-treated groups also indicate 941 monthly differentially expressed genes (monthly DEGs). These monthly DEGs are mainly distributed from April and September, which reveal a potential signal pathway map regulating the anagen stage triggered by melatonin. Enrichment analysis shows that Wnt, Hedgehog, ECM, Chemokines and NF-kappa B signaling pathways may be involved in the regulation of non-quiescence and secondary shedding under the in uence of melatonin.

Conclusions
Our study decodes the key regulators of the whole cashmere growth cycle, laying the foundation for the control of cashmere growth and improvement of cashmere yield.

Background
The Inner Mongolian cashmere goat (Capra hircus), an excellent cashmere goat breed in China, is famous for producing cashmere with superior quality and high yield. The cashmere goat has two different brous hair structures: thick and coarse hairs forming the guard layer and soft hairs forming the ground cashmere. The cashmere comes from secondary hair follicles (SHFs) in the skin [1], and the coarse hair comes from primary hair follicles [2,3]. The cashmere obtained from goats is specially used for production of expansive textile products [4]. Under the periodic changes of the natural photoperiod, the growth of cashmere in the Inner Mongolian cashmere goat presents a seasonal pattern. Normally, a typical cashmere growth cycle starts in July and stops in March of the following year with shedding of the cashmere at April [5].
Melatonin is an important mediator between photoperiod and the cashmere growth. The cyclical uctuation of melatonin levels directly affects the cashmere growth [6]. Previous studies have shown that the use of exogenous melatonin could stimulate cashmere growth during the resting period. However, a number of studies have shown that the implantation and duration of melatonin may lead to earlier shedding of cashmere [6-9], or increase production by promoting the growth of cashmere [10,11].
Although the positive role on cashmere growth of exogenous melatonin has been con rmed, previous experiments have not been able to dynamically show the gene expression pro les related to the whole cycle of cashmere growth and the potential effect of exogenous melatonin on the cashmere growth cycle.
Long noncoding RNAs (lncRNAs) are a type of RNAs that are longer than 200 nucleotides but do not encode proteins. However, lncRNAs can regulate the expression of protein-coding genes at various levels to in uence biological processes, including epigenetic regulation, transcriptional regulation and posttranscriptional regulation. The regulation mechanisms of lncRNAs in hair growth have been reported by some recent studies. For example, several important hair follicle development signals (lncRNAs and mRNAs) are involved in primary wool follicle induction in carpet sheep [22]. Yin et al. indicated that lncRNA-599554 contributes the inductive property of dermal papilla cells in cashmere goat, which might be achieved through sponging chi-miR-15b-5p to promote the WNT3A expression [23]. Wang et al. integrated analysis of lncRNA, miRNA and mRNA in cashmere goat skin during anagen and telogen stages and revealed potential ceRNA regulatory networks [24,25]. Sulayman et al. performed a comprehensive analysis of lncRNA and mRNA expression pro les during sheep fetal and postnatal hair follicle development and demonstrated that the interaction between lncRNA and their target genes may regulate the development of hair follicles [26]. However, the roles of lncRNAs in controlling the whole cashmere growth cycle have not been well described.
To clarify the regulatory mechanism of cashmere growth cycles and gain insight into the gene regulatory network perturbed by exogenous melatonin, we use RNA-seq analysis to investigate the expression patterns of seasonal changing genes (SCGs) among different gene clusters. The interactions between lncRNAs and mRNAs were also explored using co-expression network analysis. The monthly differentially expressed genes obtained from pairwise comparisons between the control and melatonin-treated groups were detected to identify the key regulators associated with the growth of secondary hair follicles and reveal potential signaling pathways which may be involved in melatonin-affected growth patterns.

Results
Transcriptome analysis and differential gene expression overview A total of 72 samples spanning 12 months were analyzed in the control (D) and melatonin (M) groups ( Figure 1A). The experiments produced 2,959,842,480 clean reads in total (~926G). A total of 22,404 genes were detected from reads counting through the RNA-seq data analysis (Additional le 1). A total of 2365 novel lncRNA genes were identi ed after the protein-coding-potential test (Additional le 2). Using DESeq2 time-series data analysis, 3559 and 988 genes were subjected as SCGs by fold change ≥2 and adjusted p-value ≤0.05 in the D and M groups, respectively ( Figure 1B) (Additional le 3 and Additional le 4). Among these SCGs, a total of 345 novel lncRNAs and 211 annotated lncRNAs were detected in the D group. We found that SCGs in the M group were much less than those in the D group, and 76.8% (759) of the SCGs in the M group also existed in the D group ( Figure 1C). According to the PCA diagrams in Figures 1D and 1E, the samples of the D group showed a certain periodic distribution, but the sample distribution was disturbed after the melatonin treatment. It can be seen in Figure 1F and 1G that the goats in two groups shed the cashmere in May, while goats in the M group had another shedding in September-November. An interesting nding is that in the M group, new cashmere has grown during the shedding in May, that is, the SHFs have not entered the resting period and the cashmere growth cycle has restarted in advance ( Figure 1G).
Stage-speci c gene expression dynamics in the cashmere growth cycle The correct division of the cashmere growth cycle is the basis for subsequent identi cation and analysis of regulatory factors. The existing studies are mostly based on morphological evidence or prior knowledge. A recent study divided the cashmere growth cycle into growth period, regression period and resting period based on skin tissue sections and transcriptome data [5]. However, their analysis is based on static differential expression analysis, the dynamic change of gene expression has not been well studied. In our study, we analyzed the time-series transcriptome data covering the entire cashmere growth cycle and obtained 3559 SCGs. Through WGCNA analysis, three key gene clusters (DC1, DC2 and DC3) that may be involved in the growth cycle of villi were identi ed in the D group (Figure 2A). The detailed genes of three clusters are listed in Additional le 5, which encompass 48.2% (1717) of the SCGs detected in the D group, including 91 annotated lncRNAs and 131 novel lncRNAs ( Figure 2B). Combining the expression patterns and functional enrichment results of three clusters (Figure 2A and 2C), we inferred the effects of these clustered genes on the cashmere growth cycle. The genes in DC1 may promote cashmere growth and the genes in DC2 may be related to the regression of the SHFs. There may be an antagonistic relationship between DC1 and DC2, which can also be inferred from the gene expression trends in Figure 2D. The gene expression in DC1 began to be up-regulated in April and reached a peak in October and started to be down-regulated, while the gene expression pattern in DC2 was exactly the opposite. In addition, it can also be found that the high expression of DC3 occurs when the expression of DC1 genes begins to be up-regulated and the expression of DC2 genes begins to down-regulate. Based on the above analysis, we divided the cashmere growth cycle into three stages associated with three gene clusters: (1) anagen (April-October); (2) catagen and telogen (October-December and January-April); (3) telogenanagen regeneration (February-May). The detailed functional analysis of these three clusters and their relationships with the cashmere growth cycle are discussed below.
The DC1 cluster shows high expression in anagen progression stage A total of 394 genes in DC1 (including 22 annotated lncRNAs and 52 novel lncRNAs) are positively correlated with the anagen progression stage (April-October), which may be involved in promoting the cashmere growth (Figure 2A and 2D). Some canonical pathways such as the Wnt, TGF-beta and Hippo signaling pathways, which have been proven to be closely related to hair growth [14,27,28] are also enriched in DC1 ( Figure 2C) (Additional le 6). Three genes of the Wnt family WNT2/WNT2B/WNT11, together with their receptors FZD3 and FZD10 are involved in the regulation of the SHF development in cashmere goats. Lymphoid enhancer binding factor-1 (LEF1) is an essential transcription factor in the Wnt signaling, and it was strongly expressed during the anagen progression stage in this study. Its function in hair cell differentiation and follicle morphogenesis has already been discussed [29]. Genes in the TGF-beta signaling (BMP2, BMP8A, BAMBI and SMAD6) also showed similar patterns with the Wnt family. It is found that the genes enriched in the Hippo signaling overlap with those in the Wnt and TGFbeta signaling, indicating the crosstalk between Wnt and TGF-beta signaling formed by the joint regulation of downstream pathways [14,30]. In addition, some downstream regulatory mechanisms of hair growth and cycling were also con rmed in this study. For example, the homeobox transcription factor DLX3 and one of its regulating genes, HOXC13, were also highly expressed in the anagen stage. The regulatory cascade positions DLX3 downstream of Wnt signaling and regulates other transcription factors related to hair follicle (HF) differentiation (such as HOXC13) [31]. KRTs (KRT26, KRT35, KRT36, KRT39, KRT6A, KRT74 and KRT84) and KRTAPs (KRTAP3-1 and KRTAP11-1) in DC1 were also associated with HF development [32]. SHH and its receptor PTCH2 in the sonic hedgehog (Shh) signaling pathway were also found in DC1. The function of Shh signaling is indicated as an essential regulator for controlling ingrowth and morphogenesis of the HFs [33][34][35], but it is not necessary for initiating the HF development [33].
The DC2 cluster prefers to be highly expressed in catagen and telogen A total of 919 genes in DC2 (including 45 annotated lncRNAs and 53 novel lncRNAs) were up-regulated from October to April of the following year ( Figure 2A and 2D), which corresponded to the degenerative and resting periods of the cashmere growth cycle. It should be pointed out that due to the limitation of sampling interval, the catagen and telogen stages were not further distinguished in this study. The expression pattern of DC2 is negatively correlated with DC1, suggesting that there may be an antagonistic relationship between the genes in DC2 and DC1. DC2 genes are mainly enriched in pathways such as Cell adhesion molecules (CAMs), Cytokine-cytokine receptor interaction, Jak-STAT, Fc epsilon RI, NOD-like receptor, Rap1, PI3K-Akt, cAMP, NF-kappa B and many immune-related pathways ( Figure 2C) (Additional le 6). A previous study also found that differentially expressed genes between anagen and telogen SHF-derived dermal papilla cells of the Cashmere goat were also enriched in CAMs and Cytokinecytokine receptor interaction pathways [36]. A number of interleukin (IL) superfamily genes were involved in the enriched pathways of DC2, such as the IL1 family (IL7 and IL18), IL10, IL15, and IL receptors (IL2RG, IL3RA, IL6R, IL7R, IL11RA and IL20RA). The CC chemokine subfamily (CCL5, CCL21, CCL22 and CCL26) and receptors (CCR4 and CCR6) were also found in the DC2 cluster. These chemokines are mainly involved in cell migration, immunity and in ammation [37]. During the HF regression, these chemokines may guide the migration of immune cells such as dendritic cells [38] and Regulatory T cells [39], thereby regulating the immune response to apoptotic cells. The JAK3 and STAT4 genes in the JAK-STAT signaling were highly expressed in catagen and telogen, which have been found to maintain HF stem cell quiescence and inhibit hair growth [40][41][42]. In addition, Dickkopf1 (DKK1), a Wnt signaling inhibitor, was also found in DC2. DKK1 has been strongly suggested to promote regression of HFs by suppressing Wnt/β-catenin signaling and inducing apoptosis in follicular keratinocytes [43,44].
The DC3 cluster is speci cally expressed in telogen-anagen regeneration stage DC3 cluster contains a total of 404 genes (including 24 annotated lncRNAs and 26 novel lncRNAs) and shows speci c high expression during the transition period from February to May (Figure 2A and 2D). Due to the overlap of the early anagen and telogen-anagen regeneration stages, the genes of DC3 and DC2 were partially enriched in several same pathways like CAMs, Focal adhesion, extracellular matrix (ECM)-receptor interaction, PI3K-Akt and NF-kappa B signaling ( Figure 2C) (Additional le 6). The PI3K-Akt signaling has been proved to be essential for HF regeneration [45,46]. A previous study found that the Toll-Like Receptor 3 (TLR3) activated by a dsRNA was able to promote HF regeneration [47]. Several collagen genes (COL1A1, COL1A2, COL3A1, COL5A2, COL6A3, COL6A5 and COL6A6) in the PI3K-Akt signaling, ECM-receptor interaction and Focal adhesion were found in DC3. These collagen genes may serve as candidate biomarkers for telogen-anagen regeneration. For example, a kind of self-assembling peptide hydrogel scaffold was used to build the ECM environment in vitro to promote HF regeneration Pathway crosstalk through mRNA-lncRNA co-expression network in cashmere growth cycle A total of 16 signaling pathways, 145 pathway genes, and 93 co-expressed lncRNAs (Pearson correlation ≥0.8) (Additional le 7) are enrolled in this pathway-mRNA-lncRNA network ( Figure 3). It can be seen that the DC1 sub-network has no positive correlation with DC2 and DC3, which can be explained by the possible antagonistic relationship between DC1 and DC2. The telogen-anagen transition phase of DC3 overlaps with the telogen stage, so the DC3 subnet is closely connected to DC2. Through this network, the function of lncRNAs can be inferred by their co-expressed coding genes. Expression differences triggered by melatonin reveal a potential signal pathway map regulating cashmere growth A total of 908 genes were detected to have different expression patterns between the D group and the M group ( Figure 4A) (Additional le 8). When treated with melatonin, 80 (24.6%) genes in DC3 maintained the same pattern as group D, but the expression patterns of most genes in DC1 and DC2 have changed ( Figure 4B). Among the 908 differential changing genes (DCGs), 369 genes belonging to the three clusters DC1, DC2 and DC3 were divided into three clusters MC1 (159), MC2 (144) and MC3 (66), respectively ( Figure 4C) (Additional le 9). After melatonin treatment, MC1 and MC2 in the M group lost the periodic expression uctuations as in the D group ( Figure 4D). In Apr-May, gene expression of MC1 dropped to the lowest while expression of MC2 rose to the highest in the D group and shedding appeared. Pathway enrichment showed that MC1 and MC2 were involved in the promotion and regression of HFs, respectively ( Figure 4E) (Additional le 10). In the M group, the expression of MC1 was still rising and MC2 was falling from January to April, which may lead to a non-resting period after shedding in the M group. Therefore, genes in MC3, which may be responsible for restarting HF growth, did not show a signi cant increase in expression in February-May. In September-November, the expression of MC1 in the M group was relatively lower than that in the D group, while the expression of MC2 in the M group was higher than that in the D group, which may disrupt the growth and maintenance of HFs and cause a local shedding.
To reveal the differential expression pattern every month, we use DESeq2 to obtain monthly differentially expressed genes (monthly DEGs) between melatonin and control groups with adjusted p-value (padj) ≤0.05 and |log2FoldChange| ≥1. A total of 941 monthly DEGs were identi ed from monthly pairwise comparisons. The monthly DEGs are mainly distributed from April and September ( Figure 5A), which exactly cover the whole anagen stage, and 96% of them are protein-coding genes. The KEGG analysis results ( Figure 5B) (Additional le 11) showed that Hedgehog related genes (SHH, PTCH2, PTCH1) and Wnt related genes (FZD10, WIF1, LEF1, WNT11) were up-regulated in April, while the GO results ( Figure 5C)  [54]. KRT71 is an inner root sheath keratin, and the mutant of KRT71 can disrupt keratin intermediate lament formation [55]. Therefore, the Hedgehog related genes (SHH, PTCH2, PTCH1, FOXE1), Wnt related genes (FZD10, WIF1, LEF1, WNT11), and other hair development related genes FOXN1, HOXC13, KRT25 and KRT71 may be responsible for the initiation of a fast anagen progressing stage from April to July. Meanwhile, the KEGG results ( Figure 5B) showed that the expression of ECM receptor interaction genes (COL6A3, THBS3 Table 1. seasons and explored the regulation of seasonal variation genes on the cashmere growth cycle of the cashmere goat and milk goat [59]. However, these studies only selected three or more stages determined by experiments or experience at the cellular level. There is still a lack of research on the in-depth exploration of the dynamic pattern of gene expression during different cashmere growth stages on the scale of whole cycles. Therefore, this study performed transcriptome sequencing on the skin samples covering the entire cashmere cycles for 12 months, which aims to explore the dynamics of gene expression in the cashmere growth cycle in more detail. The gene expression pattern for 12 months can provide useful information for distinguishing different cashmere growth stages from the genetic and molecular levels. According to the cluster-month correlations in Figure 2C, we grouped the cashmere growth cycle into three main stages: (1) anagen (April-October); (2) catagen and telogen (October-December and January-April); (3) telogen-anagen regeneration (February-May). The corresponding gene clusters are DC1, DC2 and DC3, respectively. Some canonical pathways such as the Wnt, TGF-beta and Hippo signaling pathways are enriched in DC1. DC2 genes are mainly enriched in pathways such as Cell adhesion molecules (CAMs), Cytokine-cytokine receptor interaction, Jak-STAT, Fc epsilon RI, NOD-like receptor, Rap1, PI3K-Akt, cAMP, NF-kappa B and many immune-related pathways. Interestingly, due to the overlap of the early anagen and telogen-anagen regeneration stages, the genes of DC3 and DC2 were partially enriched in several same pathways like CAMs, Focal adhesion, extracellular matrix (ECM)-receptor interaction, PI3K-Akt and NF-kappa B signaling. Besides, by constructing a co-expression network of genes (that are enriched in key pathways) and lncRNAs in three clusters, we reveal the possible regulators for crosstalk between different signaling pathways, and unearthed novel lncRNAs that may participate in these pathways.
In addition to unraveling the gene expression regulation mechanisms of the transition between different stages of the hair follicle cycle, this study also helps to gure out the role of exogenous melatonin in the speci c stages of the cashmere growth cycle. By identifying genes that exhibit different expression patterns during the cashmere growth cycle under the stimulation of melatonin, we also obtained three gene clusters (MC1, MC2 and MC3) that may affect the cashmere growth cycle. Among them, MC1 genes (BAMBI, BMP2, BMP8A, FZD10, LEF1, PPP2R1B, SMAD6 and WNT11) and MC2 genes (IL6R, IL7R, IL11RA, IL15, IL18, PDE1A, PDE1B and PDE3B) showed opposite periodicity in group D. However, after the melatonin treatment, this regular uctuation has been disordered. MC3 genes (COL1A1, COL1A2, COL3A1, CHAD, CREB3L1 and THBS3) were expressed speci cally in the anagen restart phase (Apr-May) in group D, but there was no similarly signi cant expression pattern in group M.
The relative expression levels of monthly DEGs (Additional le 13) show that the hair development related genes HOXC13, KRT25, KRT71, FOXN1 were generally expressed at higher levels at the beginning of fast anagen progressing period from April to May, implying that they may function to promote the initiation of anagen. Wnt genes (Wif-1, WNT11, FZD10, LEF1, NOTUM, SFRP2, WNT6) together with Hedgehog genes (SHH, PTCH1, PTCH2, FOXE1) showed higher expression levels between April and May, but decreased in August, which implied that Wnt-related genes may promote the rapid transition into anagen phase of hair follicles between April and May, and repress the growth of hair follicles on the eve of the second cashmere shedding period in August. Chemokines (CCL17, CCL22, CCL2, LYN, RAC2, LOC102170772, PIK3CG and VAV1) and NF-κB genes (ZAP70, LYN, BTK, CD40LG, LTB) were highly expressed in September. The NF-κB pathway may facilitate the progress of the subsequent cashmere growth phase.
Meanwhile, chemokines such as LTN may promote the second cashmere shedding.
KEGG pathway could be used as a reference to demonstrate the regulatory relationships of differentially expressed genes. Taking the above results together and collating the relevant KEGG pathway visualization results (Additional le 14), here we proposed a signaling pathway diagram of melatonin in uenced cashmere growth cycle (Figure 6), which covered the main differentially expressed genes related to cashmere growth in anagen phase from April to September. The anagen phase of melatonintreated groups was composed of a fast anagen progressing stage and a second cashmere shedding stage. The fast anagen progressing stage was from April to July, and this period was characterized by the occurrence of the rst massive cashmere shedding at the end of April, and the rapid transition into anagen phase of hair follicles from May to July, when the quick resumption of cashmere growth appeared instead of residing in a resting non-growth period. The rapid resumption of the anagen phase of hair follicles may be due to the high expression of KRT25, HOXC13 and HOXC13's regulatory target FOXN1, high expression of FZD10, WIF1, LEF1, WNT11 in Wnt signaling pathway, and SHH, PTCH1, PTCH2, FOXE1 in sonic hedgehog signaling pathway. The second cashmere shedding period was from August to September. The appearance of the second cashmere shedding may not only be associated with the low expression of ECM signaling molecules such as FREM1, FREM2, FRAS1, COL1A1, COL6A3, THBS3 in June and July, sonic hedgehog signaling pathway genes such as SHH, PTCH2 and WNT signaling pathway genes such as NOTUM, SFRP2, WNT6 in August, but also with the high expression of chemokines such as CCL, LYN, PIK3CG, VAV1, RAC2 in August. In addition, the highly expressed genes in NF-kappa B signaling pathway such as CD40LG, LTB, ZAP, LYN, BTK in September may promote the subsequent growth of cashmere after the second cashmere shedding period.

Conclusion
In summary, this study systematically analyzed RNA-seq data from skin samples of cashmere goats covering the entire cashmere growth cycle and identi ed a series of key regulators (including genes and lncRNAs) that may be involved in the cashmere growth processes. Based on gene expression patterns, we elucidated a more precise division of the cashmere growth cycle from the molecular level. However, due to the lack of sampling points, some key stages (especially the transition state) are still not well identi ed. A possible way is to increase the sampling time density before and after the stage of interest. In addition, differences in individual development, such as different growth rates, may also cause bias in the monthly differential expression analysis.

Skin sample and cashmere collection
We enrolled six cashmere goats of the same gender and age. All goats were randomly divided into two groups: control group (D) and melatonin group (M). Three cashmere goats in the experimental group were subcutaneously implanted with melatonin for one year since December 2014. The control group was untreated. The skin samples from the six cashmere goats were collected monthly (Jan 19th, Feb 13th, Mar 21st, Apr 18th, May 20th, Jun 20th, Jul 18th, Aug 23rd, Sept 19th, Oct 20th, Nov 20th, and Dec 20th ) from both the control and melatonin groups. Immediately after collection, all samples were frozen in liquid nitrogen and stored at −80°C until use. At each sampling time point, 50 cashmeres were taken from each goat to measure the cashmere length.

RNA extraction, library preparation and sequencing
Total RNA from 72 collected skin samples were isolated using the TRIzol TM reagent (Invitrogen, USA) following the manufacturer's instructions (https://assets.thermo sher.com/TFS-Assets/LSG/manuals/trizol_reagent.pdf). RNA purity was veri ed using the NanoPhotometer® spectrophotometer (Implen, USA) (https://www.implen.de/wp-content/uploads/docs/Implen-NanoPhotometer-User-Manual-N120-NP80-N60-N50-C40.pdf Where G ijt is the normalized expression value of gene i in goat j at sampling time t. The thresholds for adjusted p-value and fold change were set to 0.05 and 2, respectively. Besides, the differential changing genes (DCGs) between control and melatonin group were detected using LRT with model "~ time + condition + time: condition" versus "~ time + condition".

Gene cluster detection
The SCGs detected in control group were included for weighted gene co-expression network analysis The annotations were all achieved with clusterPro ler [70] package, with q-value ≤ 0.05 for GO and pvalue ≤ 0.05 for KEGG. Annotation information was retrieved from the Ensembl database using AnnotationHub [71] to generate an OrgDb annotation le.
Pathway-mRNA-lncRNA network construction We selected genes enriched in key signaling pathways as sources, and co-expressed lncRNA genes as targets to construct the pathway-mRNA-lncRNA network. The Pearson correlations between pathway genes and lncRNA genes were calculated, and those mRNA-lncRNA pairs with correlation higher than 0.8 were selected to construct the network. CytoScape [72] was used to visualize the network.

Calculation of relative gene expression for monthly DEGs
To observe the effects of melatonin on the expression of representative differentially expressed gene, we plotted relative expression boxplots to visually compare the expression trends of these genes. The horizontal axis of the boxplot represents months, and the vertical axis represents relative expression levels, which were calculated as follows: The gene counts were normalized with total counts in different groups. Three relative expression values per month were used for boxplot visualization.    Pathway-mRNA-lncRNA co-expression network. A total of 16 signaling pathways, 143 pathway genes, and 93 co-expressed lncRNAs (Pearson correlation ≥0.8) are involved in this network. The size of the pathway node is positively related to its degree.