LncRNAs regulate the cyclic growth and development of hair follicles in Dorper sheep

Introduction Hair follicles in Dorper sheep are characterized by seasonal cyclic growth and development, consequently resulting in hair shedding during spring. The cyclic growth and development of hair follicles are regulated by several influencing factors such as photoperiods, hormones, age of the animal, genes, long non-coding RNAs (lncRNAs), and signaling pathways. Methods In the present study, skin samples of five shedding sheep (S), used as experimental animals, and three non-shedding sheep (N), used as controls, were collected at three time points (September 27, 2019; January 3, 2020; and March 17, 2020) for RNA sequencing (RNA-seq) technology. Nine different groups (S1-vs-S2, S1-vs-S3, S2-vs-S3, N1- vs-N2, N1-vs-N3, N2-vs-N3, S1-vs-N1, S2-vs-N2, and S3-vs-N3) were compared using FDR < 0.05 and log 21 FC >as thresholds to assess the differences in the expression of lncRNAs. Results and discussion In total, 395 differentially expressed (DE) lncRNAs were screened. Cluster heatmap analysis identified two types of expression patterns, namely, high expression during the anagen phase (A pattern) and high expression during the telogen phase (T pattern). Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that the target genes were largely enriched in the Estrogen signaling pathway, PI3K-Akt signaling pathway, Fc gamma R-mediated phagocytosis, and cell adhesion molecules (CAMs), which are associated with hair follicle cyclic growth and development-related pathways. In addition, 17 pairs of lncRNAs-target genes related to hair follicle cyclic growth and development were screened, and a regulatory network was constructed. Altogether, candidate lncRNAs and their regulated target genes were screened that contributed to sheep hair follicle cyclic growth and development. We believe these findings will provide useful insights into the underlying regulatory mechanisms.


Introduction
The hair follicles of certain animal species, such as Inner Mongolian cashmere goats and Dorper sheep, are characterized by seasonal cyclic growth and development. This process results in hair shedding during spring. The cyclic growth and development of hair follicles occur in three phases, namely, anagen, catagen, and telogen (1). Anagen is the most active phase of hair follicle growth, during which hairs grow rapidly and form a complete hair shaft (HS) (2), and the duration of anagen determines the length of the hair. During the catagen phase, HS stops growing; cell proliferation and differentiation capacity begin to decline; cells start undergoing apoptosis; and hair follicles rapidly degenerate (3). Subsequently, the telogen phase begins, during which the biological activity of hair follicles is the weakest and HS decreases (3). During the cyclic growth and development of hair follicles, a large and complex regulatory network comprising genes, microRNAs (miRNAs), long non-coding RNAs (lncRNAs), and other regulators is formed through mutual inhibition or synergistic effects (3).
Several pathways, such as WNT (4), MAPK (5), BMP (6), and TGF-β (7) signaling pathways, have been implicated in the regulation of hair follicle growth and development. Any abnormality in the ligands, receptors, and signal transduction molecules of these signaling pathways can affect the development of hair follicles in animals and changes in hair growth (3). In these pathways, BMP [BMP4 (8), BMP7 (9)], FGFs [FGF20 (10), FGF10 (11), FGF4 (12)], Sox2 (13), Sox18 (14) are intricately associated with hair follicle development and wool bending. NF-κB (15), WNT10b (16), LEF1, IGF1R (17), and Noggin (18) have been reported to promote early hair follicle morphogenesis. Similarly, DKK1 (16), PF4 (19), FGF5, and FGF18 (20) have been found to inhibit the growth of hair follicles. Numerous studies have demonstrated the role of lncRNAs in the complex regulatory network. For instance, Yue et al. determined the functions of lncRNAs and mRNAs in sheep skin by strand-specific RNA sequencing (ssRNA-seq) during the initiation phase of secondary hair follicle (SHF) development. They concluded that crucial differentially expressed genes (DEGs) and lncRNAs influence hair follicle initiation (21). Yang et al. screened lambda cyclic growth-specific lncRNAs and found that 6,127 lncRNAs were expressed during anagen and telogen phases, of which 54 were significantly differentially expressed (32 up-regulated and 22 down-regulated) (22). LncRNA15479 was identified as a crucial molecule that was intricately associated with the anagen phase of hair follicles and involved in the formation of keratin by target knockdown experiments (22). In addition, PlncRNA-1 was found to be responsible for the biological changes in the Wnt/β-catenin signaling pathway by regulating TGF-β1, thereby participating in the biological regulation of hair follicle stem cell proliferation and differentiation (23). In the present study, we screened and identified lncRNAs related to hair follicle cyclic growth and development. However, the regulatory mechanisms warrant further investigation.
Shedding sheep (S) with cyclic hair follicle growth and development were selected as experimental animals and non-shedding (N) sheep as the control. The expression of lncRNAs was determined at three different time points by RNA sequencing technology (RNAseq), followed by the screening of lncRNAs by differential expression analysis and expression pattern analysis. Finally, candidate lncRNAs were screened by lncRNAs-target genes synergy analysis. The specific analysis flow of lncRNA is shown in Supplementary Figure 1. We investigated the lncRNAs affecting anagen and telogen phases of sheep hair follicles and explored the specific regulatory mechanisms to provide a theoretical basis for breeding automatic shedding meat sheep.

Experimental animals and sample collection
The study was approved by the Experimental Animal Welfare and Ethics Committee of Ningxia University (Animal Ethics approval no. NXU-19-018). The experimental sheep were obtained from the sheep farm of Ningxia China Animal Husbandry Yilin Livestock Co., Ltd. The skin samples were collected from five 2-yearold shedding Dorper ewes and three 2-year-old non-shedding Dorper ewes on September 27, 2019; January 3, 2020; and March 17, 2020, from the side of the sheep's body (the junction of the posterior edge of the last rib and the midline of the body). The sample size was approximately 2 cm 2 of skin tissue. Samples were split and stored in liquid nitrogen. The hair follicles were in the anagen phase on September 27, 2019, for both shedding and non-shedding sheep. The hair follicles were in the telogen phase on January 3, 2020, for shedding sheep, whereas the hair follicles of non-shedding sheep were still in the anagen phase. The hair follicles were in the early anagen phase on March 17, 2020, for shedding sheep, whereas those of non-shedding sheep were in the slow anagen phase. The samples were divided and stored in liquid nitrogen.

RNA isolation, library preparation, and sequencing
First, 50-100 mg of sheep skin samples were obtained from eight Dorper sheep separately and placed in mortars. They were ground into powder using a pestle under liquid nitrogen. Next, the powder was put in a 1.5 mL centrifuge tube and 1 mL of the TRIzol reagent (Invitrogen, USA) was added. It was kept at room temperature for 5 min. Next, 0.2 mL chloroform (Tianjin Regent Chemical Co., Ltd) was added, shaken for 15 s, and allowed to stand for 10 min. Next, the supernatant was centrifuged at 4°C using a benchtop high-speed refrigerated centrifuge (Kate's Laboratory Instruments Co., Ltd., TGL 18 M) at 12,000 g for 15 min, following which the supernatant was placed in another centrifuge tube. Third, 0.5 mL isopropanol (Tianjin reagent Chemical Co., Ltd.) was added to the supernatant, the liquid was mixed in the tube gently, and it was made to stand for 10 min at room temperature and centrifuged (using the method above). After centrifugation, the supernatant was discarded and the precipitate was retained. Finally, 1 mL of 75% ethanol (Yantai Shuang Chemical Co., Ltd.) was added to the precipitate, which was gently obtained. It was centrifuged at 4°C at 7500 g for 10 min and the supernatant was discarded. It was air dried (15 min) and an appropriate volume of RNase-free water (Tiangen Biochemical Technology; Beijing) was added to dissolve the precipitate (65°C for 10-15 min). The purity and concentration of the extracted RNA were measured using the NanoDrop 2000 spectrophotometer (NanoDrop 2000, UV-Vis spectrophotometer, USA) and Agilent 2100 Bioanalyzer RNA Nano 6000 kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The extracted total RNA was reverse transcribed into cDNA using a reverse transcription kit (Nanjing Novozymes Biotechnology Co., Ltd.). The library was constructed using the illumina Ribo-Zero Gold (Human/Mouse/Rat) Kit (MRZG12324, illumine) and NEB#7490 Kit (NEB E7490L, New England Biolabs) according to the manufacturer's instructions. Subsequently, 24 cDNA libraries from eight sheep at three different time points were generated. The quality of the constructed library was evaluated using Agilent 2100, and after passing library inspection, 24 cDNA libraries were sequenced together using the Illumina HiSeqTM 4000 (Illumina; San Diego, CA, USA).
Frontiers in Veterinary Science 03 frontiersin.org

Quality control
The quality of the raw data was assessed using fastp (24). The procedure included the following steps: (1) removing reads containing adapters, (2) removing reads with a percentage of unknown bases (N) greater than 10%, (3) removing reads consisting entirely of A bases, and (4) removing low-quality reads (reads with Q phred ≤ 20 bases accounting for more than 50% of the entire read length). The remaining clean reads were mapped to the reference genome Oar_ rambouillet_v1.0 using HISAT2 (25).

Comprehensive analysis of lncRNAs 2.4.1. Identification of lncRNAs
The transcripts of each biological sample were reconstructed using the merge function of StringTie (26), and the obtained transcripts were compared with known transcripts from the genome annotation to exclude known transcripts. The coding ability of the new transcripts was then performed using both CPC2 (27) and CNCI (28) software, and the intersection of transcripts without coding potential was taken as a reliable prediction, and the existing lncRNA annotations of the reference genome were excluded. Rfam is a database containing information on various lncRNA families, including conserved regions of RNA secondary structure, mRNA cis-acting elements, and other RNA elements. It classifies ncRNAs into different families based on their common ancestry at the evolutionary level, and each family consists of secondary structures and covariates predicted from multiple sequence comparisons. To better annotate lncRNA at the level of evolution, we used Infernal (v1.1.2) (29) to classify all predicted lncRNAs according to their conserved sequence and secondary structure through multiple sequence alignment, secondary structure, and a covariance model. In addition, we used Cuffcompare (30) to compare the positions of lncRNAs with the positions of protein-coding RNAs. Based on this lncRNA were then classified into intergenic lncRNA (cuffcompare class code u,), intronic lncRNA (cuffcompare class code i), sense lncRNA (cuffcompare class code j, o), antisense lncRNA (cuffcompare class code x) and other lncRNA (remaining cuffcompare class codes).

Analysis of differentially expressed lncRNAs
The expression of lncRNAs was displayed as raw read count and FPKM (Fragments Per Kilobase of transcript per Million mapped reads) value. Many factors influence raw reads such as transcript length and the total number of reads. Raw reads were not conducive to the comparison of differential lncRNAs between the samples. To ensure the accuracy of the subsequent analysis, we first corrected the sequencing depth and subsequently, the length of the transcript to obtain the FPKM value of lncRNA before further analyses. The transcripts were quantified with RSEM (31). Based on the lncRNA expression information, we conducted principal component analysis (PCA) using R 1 , and calculate the Pearson correlation coefficient (cor) between samples. The more similar the samples were, the closer they were reflected in the PCA plot, and the samples from different effective 1 http://www.r-project.org/ treatments tended to show their own aggregated distribution. Differential expression analysis of lncRNAs was performed using DESeq2 (32), which provided statistical analyses using a model based on the negative binomial distribution. Nine different groups (S1-vs-S2, S1-vs-S3, S2-vs-S3, N1-vs-N2, N1-vs-N3, N2-vs-N3, S1-vs-N1, S2-vs-N2, S3-vs-N3) were compared using FDR < 0.05 and log 2 1 FC > as thresholds to assess the statistical significant differences in the expression of lncRNAs. The differentially expressed (DE) lncRNAs between the final different comparison groups were considered a concurrent set.

Cluster heatmap analysis and expression pattern analysis of differentially expressed lncRNAs
A cluster heatmap analysis was conducted using the FPKM values of DE lncRNAs for three phases of two groups to explore the expression pattern of lncRNAs using the OmicShare tools 2 , a free online platform for data analysis. Whether the lncRNAs have similar expression patterns were determined using the high and low expression (FPKM) and trend directions of lncRNAs at three different time points in the S and N groups. For example, both lncRNA 1 and lncRNA 2 displayed a high-low-high expression trend in the S group, and both exhibited an all-time elevated expression trend in the N group, which we considered as similar expression patterns. LncRNAs with similar expression patterns often have similar functions or are involved in the same metabolic processes or pathways.

Target gene prediction and correlation analysis of differentially expressed lncRNAs 2.5.1. Target genes prediction of differentially expressed lncRNAs
To explore the function of lncRNAs, we predicted cis-target genes (neighboring genes), trans-target genes (co-expressed genes), and antisense lncRNA genes of identified lncRNAs using mRNAs sequence and expression information from mRNA sequencing data (33) of the same experimental sheep as in this study. Trans-acting lncRNAs influenced the gene expression in diverse biological processes at the transcriptional or post-transcriptional level, whereas cis-acting lncRNAs epigenetically regulated the nearby genes only at the transcriptional level (34). LncRNAs can trans-regulate distant target genes because they have the same expression pattern (35). Therefore, genes with |cor| ≥ 0.95 with this lncRNA were selected as trans-target genes. A cis-acting lncRNA is defined as the one regulating proteincoding genes within the proximity of its genomic locus (10 kb in our case) (36). Subsequently, the cor between lncRNA expression and cis-target genes expression was calculated, and |cor| ≥ 0.8 was retained. Antisense lncRNAs are transcription products from the opposite strand of the protein-coding or sense strand, which regulate the corresponding gene by gene silencing or by degrading sense transcripts (37,38). The software RNA-plex (39) was used to predict the complementary correlation of antisense lncRNAs and target genes Frontiers in Veterinary Science 04 frontiersin.org (antisense-role). Subsequently, the cor between lncRNA and gene expression was calculated, and |cor| ≥ 0.8 was retained.

GO and KEGG enrichment analysis of target genes and construction of the regulatory network
Gene Ontology (GO) enrichment analysis of target genes after target gene prediction and screening were performed using the GOseq R software package. Signaling pathway enrichment analysis was performed using the Kyoto Encyclopedia of Gene and Genomes (KEGG) database. GO terms and pathways with q-value (adjusted p-value) < 0.05 were considered significantly enriched. Candidate lncRNAs-target genes-pathways were screened and the regulatory network was visualized using Cytoscape (v3.9.1) (40).

qRT-PCR
To validate the accuracy of the RNA-seq results, six DE lncRNAs were selected for quantitative real-time polymerase chain reaction (qRT-PCR). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the reference gene for analysis. The qRT-PCR primers were designed using the Primer Premier 5.0 software (41), and primer information is shown in Table 1. The 2 −ΔΔCt method was used to process the qRT-PCR results (42).

Overview of RNA-seq
After extracting total RNA from eight Dorper sheep during three periods, we constructed 24 transcriptome libraries from skin samples and sequenced them individually. A total of 2,184,682,876 raw reads were produced using the Illumina HiSeqTM 4000. After trimming to remove adaptor sequences, and discarding low-quality sequences, 2,178,092,054 clean reads were retained (accounting for 99.69% of raw reads) ( Table 2). The percentage of Q20 of each sample was not less than 97.28% and the percentage of the Q30 base of each sample was not less than 92.75%. The GC content ranged between 48.46 and 53.21%. The average GC content value with a standard deviation was 1.16%. Subsequently, we mapped the clean reads to the reference genome Oar_rambouillet_v1.0 using HISAT2. The matching rate was above 93.65% (Table 2). Quality control results confirmed the reliability of sequencing results and adequacy for further data processing. We reconstructed 9,939 transcripts and 6,970 genes using StringTie. After coding potential analysis using the software CNCI and CPC2, 1,699 novel lncRNA transcripts were identified ( Figure 1A). Five different categories of lncRNAs were classified based on their positions relative to the protein-coding genes, including intergenic lncRNAs, bidirectional lncRNAs, intronic lncRNAs, antisense lncRNAs, and sense lncRNAs ( Figure 1B).

Principal component analysis
Principal component analysis (PCA) was performed using FPKM values of all lncRNAs to identify the sample clusters and distribution patterns (Figure 2A). For the PCA, the closeness of samples was associated with a stronger positive correlation. For the S group, the samples of S1 were closer to those of S3, whereas those of S2 were away from both because both S1 and S3 were in the Frontiers in Veterinary Science 05 frontiersin.org For example, S1-1, S represent shedding group, the first 1 is the first stage of S, the second 1 is the first sheep, and so on.
Frontiers in Veterinary Science 06 frontiersin.org anagen phase and S2 was in the telogen phase. Therefore, they were distributed differently in Figure 2A. The three phases in the N group were all anagen, and the difference was that N1 and N2 were in the anagen phase and N3 was in the slow anagen phase, resulting in the samples being close together. There existed outlier samples, such as S1-1 and N3-2, probably because the cluster was affected by the gene's expression from other tissues in skin samples instead of hair follicles.

Analysis of differential expression lncRNAs
A total of 395 DE lncRNAs were identified. Among them, 313 DE lncRNAs were in the S group (S1-vs-S2, S1-vs-S3, S2-vs-S3), 49 DE lncRNAs in the N group (N1-vs-N2, N1-vs-N3, N2-vs-N3), and 155 DE lncRNAs in the S-vs-N group (S1-vs-N1, S2-vs-N2, S3-vs-N3). The number of DE lncRNAs with the up−/down-regulation in each group is shown in Figure 2B. The S2 group had the highest number of DE lncRNAs compared with others due to the larger difference between S2 in the telogen phase and other phases. However, the number of DE lncRNAs between the other groups was low. These findings are consistent with the results of PCA clustering.

Cluster heatmap and expression pattern analysis of differentially expressed lncRNAs
Cluster heatmap analysis revealed clustering of similar development stages of hair follicles ( Figure 3A). S1, N1, and N2 groups were clustered because they were present in the same active anagen phase in both groups. S3 and N3 groups were clustered together because early anagen and slow anagen phases had similar hair follicle growth states. S2 was a single cluster because the telogen phase was considerably different from other stages.
Cluster heatmap analysis revealed two types of important expression patterns. Pattern one, in which the expression of DE lncRNAs was the highest in the anagen phase (S1, N1, and N2), lowest in S2 (telogen), and middle in the early anagen phase and the slow anagen phase (S3 and N3), was defined as the A (Anagen) pattern. Pattern two, in which the expression of DE lncRNAs was the highest in

Target gene prediction results of lncRNAs
To investigate the functions of lncRNAs, potential target genes of 129 screened lncRNAs were predicted according to trans-role (|cor| ≥ 0.95), cis-role (upstream and downstream 10 kb), and antisenserole. The 118 (trans-role), 29 (cis-role), and 29 (antisense-role) lncRNAs-target genes pairs were identified by three prediction methods, respectively (Supplementary Table 1). The cor of lncRNAs-target genes pairs in cis-role and antisense-role were further calculated, and 19 pairs with |cor| ≥ 0.8 were left in each. All lncRNAs-target gene pairs screened by three prediction methods were positively correlated, suggesting that the lncRNAs probably positively regulated the corresponding target genes to affect the cyclic growth and development of hair follicles.

Expression pattern analysis of target genes
A prediction of the target genes of lncRNAs resulted in 99 targets (Supplementary Table 2), all of which were differentially expressed genes (DEGs) (FDR < 0.05 and | log 2 | 1 FC > ). Their expression patterns were analyzed. Among them, 63 fitted the A pattern ( Figure 4A), 28 fitted the T pattern ( Figure 4B), and 8 DEGs did not fit the A and T patterns (Figures 4C,D). The 91 DEGs of A and T patterns were left for subsequent analysis.

GO and KEGG enrichment analysis of target genes
In terms of biological processes, genes were concentrated in cellular process, single-organism process, regulation of biological process, and biological regulation. In terms of molecular functions, genes preferentially corresponded to binding, including protein binding, nucleic acid binding, heterocyclic compound binding, and organic cyclic compound binding-related genes. In terms of cellular components, genes were enriched in organelle, cell and cell part, including genes related to cell structure, connectivity, and communication. The GO enrichment results displayed specific patterns in the growth and development of Dorper sheep hair follicles, thereby contributing to an in-depth understanding of the biological mechanism of hair follicles.
Next, the KEGG enrichment analysis was performed on 91 DEGs that fitted the A and T patterns. A total of 40 pathways were enriched, and the target genes were significantly enriched in Systemic lupus erythematosus (ko05322), Alcoholism (ko05034), Shigellosis (ko05131), Transcriptional misregulation in cancers (ko05202), Staphylococcus aureus infection (ko05150), Estrogen signaling pathway (ko04915), Viral carcinogenesis (ko05203), and Necroptosis (ko04217). The top 20 enriched canonical pathways are shown in Figure 6. Among these 40 pathways, 9 lncRNAs that fitted the A pattern were identified to regulate 12 pathways by targeting 16 corresponding target genes, whereas 8 lncRNAs that fitted the T pattern were identified to regulate 31 pathways by targeting 8 corresponding target genes. The lncRNAs and their corresponding target genes identified in these pathways are summarized Frontiers in Veterinary Science 08 frontiersin.org Results of GO enrichment analyses of target genes. Blue, orange, and purple words represent the GO terms in biological processes, molecular functions, and cellular components, respectively. The x-axis represents the different GO terms, and the y-axis represents the number of genes. Frontiers in Veterinary Science 09 frontiersin.org in Supplementary Table 3. In addition, the Estrogen signaling pathway (ko04915) and PI3K-Akt signaling pathways (ko04151) were found to be associated with the anagen phase of hair follicles. Similarly, Fc gamma R-mediated phagocytosis (ko04666), Chemokine signaling pathway (ko04062), Bacterial invasion of epithelial cells (ko05100), Cell adhesion molecules (CAMs) (ko04514), and Regulation of actin cytoskeleton (ko04810) are associated with the telogen phase of hair follicles.

Construction of the regulatory network
As shown in Table 3

qRT-PCR
Six DE lncRNAs were randomly selected for quantitative real-time polymerase chain reaction (qRT-PCR) verification analysis. The verification results were generally consistent with the transcriptome sequencing results, indicating the reliability of our sequencing results (Figure 8).

Discussion
The hair follicle undergoes three phases of cyclic growth and development, namely, anagen, catagen, and telogen phases (43). LncRNAs have been implicated in the regulatory role of the cyclic growth and development of hair follicles (44). In this study, shedding sheep and non-shedding sheep were used as experimental animals, and lncRNAs affecting the anagen and telogen phases were effectively screened. In addition, their target genes and corresponding pathways enriched to the regulatory network were analyzed. Results of KEGG enrichment analyses of target genes. The x-axis and y-axis indicate the rich factor and the KEGG pathway names, respectively. The size of the circle represents the number of genes enriched in the pathway. The purple color represents a high enrichment (−log10 of Q value).
Frontiers in Veterinary Science 10 frontiersin.org

Validity of lncRNA screening
To study hair follicle cyclic growth and development, the majority of the researchers have used different developmental stages of the same population. For example, Wu et al. selected three 2-yearold female Jiangnan Cashmere goats as experimental animals to screen the candidate lncRNAs, mRNAs, and related pathways affecting the three developmental stages of secondary hair follicles by RNA-seq (45). Similarly, Wang et al. selected five 1-year-old female Capra hircus in their comprehensive analysis of genes encoding and non-coding RNAs of the hair follicle cycle, of which three were used to construct and analyze the lncRNA library and two for constructing miRNA library (46). Li et al. screened three Inner Mongolian Cashmere goats for seasonal genes during the cyclic growth and development of hair follicles (47). These studies selected only one population as the experimental animals and screened regulatory molecules such as mRNAs, lncRNAs, and miRNAs related to hair follicle cyclic growth and development. These studies were prone to adulteration of false-positive genes caused by seasonality and other factors. In contrast, we used the shedding population (S1: anagen, S2: telogen, S3: early anagen) of Dorper sheep as experimental animals and the non-shedding population (N1, N2: anagen, N3: slow anagen) as the control. We revealed two lncRNAs expression patterns with high expression in anagen or telogen phases using the cluster heatmap analysis. This experimental design effectively screened the candidate lncRNAs affecting the anagen and telogen phases of hair follicles. The use of a non-shedding population as a control allowed the removal of certain false positives caused by seasonal differential genes, such as MSTRG. 11915 Figure 2A) and other 76 lncRNAs identified in this experiment (Supplementary Table 6). These lncRNAs displayed a consistent expression pattern in both populations. As shown in Figure S1B, 31 lncRNAs exhibited a highlow-high expression pattern at all three time points in both S and N groups. Thus, these DE lncRNAs were affected by the season instead of the cyclic growth and development of hair follicles. We removed these DE lncRNAs with consistent expression patterns. Finally, 129  LncRNAs-target genes-pathways association network. A total of 10 lncRNAs, 10 target genes, and 7 signaling pathways are involved in this network. The circle represents pathways, the diamond represents target genes, and V represents lncRNAs.
Frontiers in Veterinary Science 11 frontiersin.org DE lncRNAs significantly associated with hair follicle anagen and telogen phases were screened for subsequent analysis, which effectively excluded the interference of seasonal factors.

Regulatory role of lncRNAs and corresponding target genes in the cyclic growth and development of hair follicles
The normalized FPKM values of 395 DE lncRNAs were analyzed by clustering heatmap analysis, and 129 DE lncRNAs were analyzed to fit two kinds of expression patterns (A and T patterns), among which 57 lncRNAs were found to play regulatory roles in the anagen (A pattern) and 72 lncRNAs in telogen (T pattern) phases. The target genes of these lncRNAs were predicted using three methods (co-localization, co-expression, and antisense effects), combined with GO and KEGG functional analyses, to explore the regulatory functions of candidate lncRNAs in anagen and telogen phases of hair follicles.
Finally   (53). Among them, MSTRG.13836.1-KRT28 was screened by both trans-role and antisense-role analyses and was used for further study. Wang et al. identified crucial pathways and genes associated with hair follicle cyclic growth and development in cashmere goats and reported that genes associated with the anagen phase were enriched in the PI3K-Akt signaling pathway (ko04151) (49). The PI3K-Akt signaling pathway (ko04151) is essential for maintaining and restoring hair induction in dermal papilla cells (54) and promoting proliferation and inhibiting apoptosis in dermal papilla cells (55). Lpar6 was detected to be enriched in the PI3K-Akt signaling pathway (ko04151) in this study. Therefore, it was hypothesized that MSTRG.12818.1 plays a regulatory role in the PI3K-Akt signaling pathway (ko04151) by targeting Lpar6 to affect the anagen phase of hair follicles.
Frontiers in Veterinary Science 12 frontiersin.org phase (48). Nagao et al. demonstrated that chemokines can induce immune cell migration following the entry of hair follicles into the catagen phase, thus further promoting the apoptosis of hair follicle cells (56). Interestingly, lncRNA MSTRG.15931.1 was predicted to target PTPRM, ELMO1, and Pip5k1c. Among them, PTPRM has been reported to regulate intercellular communication in keratinforming cells (57), and ELMO1 has been demonstrated to play an unexpected role in clearing apoptotic cells from hair follicles and maintaining homeostasis (58). Both ELMO1 (59) and Pip5k1c (60) were implicated in the remodeling of the actin cytoskeleton via the fifth pathway, i.e., the Regulation of the actin cytoskeleton (ko04810). Therefore, the specific regulatory function of lncRNA MSTRG.15931.1 in the telogen phase of hair follicles will be the focus of our further research.

Conclusion
In this study, 129 lncRNAs were screened for their expression patterns to elucidate their regulatory roles in hair follicle growth and development during the anagen and telogen phases in Dorper sheep. In total, 10

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/bioproject; PRJNA963059.

Ethics statement
The animal study was reviewed and approved by Experimental Animal Welfare and Ethics Committee of Ningxia University.

Author contributions
XL and HS contributed to the conception and design of the study. HS organized the database, performed statistical analysis, and wrote the first draft of the manuscript. KM, YFW, YYW, and XY provided the usage of some of the software. XL completed the revision of the manuscript. All authors contributed to the manuscript revision, read, and approved the submitted version.

Funding
This study was financially supported by the National Natural Science Foundation of China "Screening of key molecules and construction of regulatory network about sheep shedding trait" (No. 31960650).