The experimental procedures have been reviewed and approved by the Experimental Animal Ethics and Welfare Committee of Shenyang Agricultural University with the ethical code: 2023030208. Also, all experiments were carried out following the approved protocol guidelines. Nine clinically healthy cashmere goats (adult, female, and no traceable genetic relationships) were selected from the breeding center of Liaoning cashmere goat (Liaoyang, Liaoning, China). Using sterile scalpel blades, we collected skin tissues from the body side part of the selected goats at three stages of SHF cycles: anagen, catagen, and telogen as described in our previous study [12]. The skin tissue sample from each goat was subjected to clearing with alcohol of 75%. Subsequently, the skin tissue was cut into pieces (5 mm²) with phosphate-buffered saline washing three times. The samples were further subjected to digesting with dispase II of 0.25% (Roche, Mannheim, Germany) at 4°C overnight. Ultimately, the SHFs were isolated from digested skin tissue using disposable syringe needle and sterilized scalpel under a stereo-microscope (Leica, Wetzlar, Germany). The total RNA was isolated from SHFs of each sample using the RNAiso reagent kit (TaKaRa, Dalian, China), and the genomic DNA were isolated using the TIANamp Genomic DNA Kit (TIANGEN, Beijing, China).

The circHECA has been identified from skin tissues of cashmere goats in a previous investigation [8]. The BioEdit software (Version 7.0.5) [13] was used for the managing and displaying of circHECA sequence. The transcription source of circHECA was defined based on an aligning of its linear sequence against the goat genome datasets (Genome assembly ARS1.2, https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_001704415.2, last access: 08 November 2023). We predicated the potential m⁶A sites on circHECA sequence with the SRAMP online program ( https://www.cuilab.cn/sramp, last access: 09 November 2023). The potential binding miRNAs within circHECA sequence were predicted using the custom sub-procedure of miRDB online service program ( https://www.mirdb.org, last access: 09 November 2023), along with the online service database: miRNAsong ( https://www2.med.muni.cz/histology/miRNAsong/index.php, last access: 09 November 2023). In the prediction analysis, here, due to the unavailability of goat data in miRDB database, alternatively, we used the human data for performing the prediction with removing the potential target miRNAs unidentified in goats. The coding potential of circHECA sequence was analyzed using two programs: Coding Potential Assessment Tool (CPAT, https://lilab.research.bcm.edu/cpat/index.php, last access: 09 November 2023) and Coding Potential Calculator (CPC, https://cpc.cbi.pku.edu.cn, last access: 09 November 2023). The open reading frame (ORF) of circHECA sequence was analyzed using ORF Finder program at NCBI website ( https://www.ncbi.nlm.nih.gov/orffinder, last access: 09 November 2023) along with the BioEdit software (Version 7.0.5) [13].

The m⁶A site validation of circHECA was performed using methylation immunoprecipitation (Me-RIP) along with reverse transcription-quantitative polymerase chain reaction (RT-qPCR) technique. We carried out the Me-RIP analysis of m⁶A-circHECA following the described methods in previous publication [14]. In short, the total RNA sample (100 μg) was incubated with RNase R (Geneseed, Guangzhou, China) at 37°C for 3 min, and was further subjected to concentration using the Monarch RNA Cleanup Kit (NEB, Ipswich, MA, USA). And then, the obtained RNA sample was fragmented by the NEBNext Magnesium RNA Fragmentation Module (NEB, USA), followed by a concentration using the Monarch RNA Cleanup Kit (NEB, USA). The fragmented RNA sample of 2 μg was used as input control. At 4°C, half fragmented RNA product was incubated with anti-m⁶A antibody (2 μg; Synaptic Systems, Gottingen, Germany) or immunoglobulin G (2 μg; Cell Signaling Technology, Danvers, MA, USA) for 4 h. Subsequently, The Dynabeads Protein A (Thermo Scientific, Rockford, IL, USA) was incubated with the RNA-antibody complex at 4°C for 2 h. Finally, the RNA was isolated with the RNAiso reagent kit (TaKaRa, China). Using random primers, the first strand cDNAs were synthesized by the M MuLV cDNA synthesis kit (Sangon, Shanghai, China). Based on the site-specific primers (Table 1), the enrichment abundance of each m⁶A site within circHECA molecule was measured by RT-qPCR technique where the relative enrichment of each m⁶A site of circHECA was normalized to the input control [14].

The nucleus and cytoplasm localizations of m⁶A-circHECA were performed using SHF stem cells of cashmere goats stored in our laboratory as described in our previous publication [1]. The cells were cultured in Dulbecco’s modified Eagle medium (DMEM)/F12 medium with 10% fetal bovine serum (Hyclone, Logan, UT, USA). The incubator was set as 37°C with the concentration of CO₂ being 5%. The medium was changed every 2 days. The cytoplasmic and nuclear RNA was isolated from SHF stem cells by Cytoplasmic and Nuclear RNA Purification Kits (AmyJet, Wuhan, China). Subsequently, the relative expression of circHECA, GAPDH and snRNA-U6 was measured in cytoplasm and nucleus of SHF stem cells using RT-qPCR technique. Here, the expression level of GAPDH and snRNA-U6 was measured to serve as the cytoplasm and nuclear control, respectively. In RT-qPCR analysis, the relative expression of the analyzed genes was calculated with the 2^−ΔΔCt method.

As described above, six potential target miRNAs of m⁶A-circHECA had been predicted using the custom sub-procedure of miRDB along with the miRNAsong procedure including chi-miR-20b, chi-miR-27a-5p, chi-miR-129-3p, chi-miR-187, chi-miR-449a-5p, and chi-miR-449b-3p. For the competitive endogenous RNA (ceRNA) regulatory network of m⁶A-circHECA, the potential target genes of the six miRNAs (chi-miR-20b, chi-miR-27a-5p, chi-miR-129-3p, chi-miR-187, chi-miR-449a-5p, and chi-miR-449b-3p) were further predicted using the custom prediction sub-procedure of miRDB along with the miRNAsong procedure as recommended in previous publication [11]. Finally, the ceRNA regulatory network of m⁶A-circHECA was generated and visualized using the Cyotoscape (Version 2.8.3) procedure [15]. The CluePedia built-in plugin in Cyotoscape procedure was used for enriching the m⁶A-circHECA regulatory genes mediated by the six predicted miRNAs into signaling pathways with the default settings ( https://www.ici.upmc.fr/cluepedia/, last access: 20 November 2023). CluePedia generated a regulatory network connecting the analyzed genes and implicated signaling pathways which indicated novel potential regulatory relationships among the analyzed genes and involved pathways [16].

For the interactional network of m⁶A-circHECA with potential target proteins, we performed a prediction within m⁶A-circHECA sequence to screen the potential interactional proteins using RNA-protein Interaction Prediction (RPISeq). The RPISeq ( https://pridb.gdcb.iastate.edu/RPISeq/batch-rna.html, last access: 24 November 2023) can computationally screen potential interactional proteins within given RNA sequence based on sequence-derived pattern [17]. The obtained interactional relationships of m⁶A-circHECA with potential target proteins were presented as a network that was further deeply explored using an online service procedure: FunRich ( www.funrich.org, last access: 24 November 2023). The potential target proteins directly or indirectly regulated by m⁶A-circHECA were further enriched into signaling pathways using the CluePedia procedure with default settings [16]. Ultimately, the significantly enriched pathways were presented as a chordmap that was generated by an online platform for data analysis and visualization ( https://www.bioinformatics.com.cn, last access: 24 Nov 2023) [18].

Using the M-MuLV cDNA synthesis kit with random primers (Sangon, Shanghai, China), the reverse-transcriptions were performed on total RNA extracted from cashmere goat SHFs during different stages (telogen, anagen and catagen). For an expression test of the m⁶A-circHECA in SHFs of cashmere goats, we performed the RT-qPCR analysis with the divergent primers (Table 1) in a light Cycler 480 real-time PCR system (Roche Diagnostics, Germany). In a final volume of 25 μL, the qPCR reactions were performed where the reaction system were comprised of Green Premix Ex Taq II of 12.5 μL TB (Tli RNaseH Plus; TaKaRa, Chnia), the first-strand cDNA solution of 2.0 μL, each primer of 1.0 μL (10 μM), and ddH₂O water of 8.5 μL. In qPCR reactions, the thermal cycling parameters were set as a single cycle of 95°C for 3 min, followed by 40 cycles (95°C for 5 s, 57°C for 30 s, and 72°C for 30 s). The relative express of m⁶A-circHECA was normalized to the GAPDH and calculated with the 2^−ΔΔCt method [19].

For assessing the methylation level of promoter region of m⁶A-circHECA host gene (HECA), we hunted potential a CpG island within the 800-bp region immediately upstream of transcription start sites (TSS) in HECA gene using Methyl Primer Express program (Version 1.0; Applied Biosystems, Foster City, CA, USA). We analyzed the potential binding sites of transcription factors using JASPAR program ( https://jaspar.elixir.no, last access: 25 Nov 2023). The extracted genomic DNA from each goat was treated by the MethylCode Bisulfite Conversion Kit (Invitrogen, Shanghai, China), and pooled into three groups: telogen, anagen and catagen of SHF cycle of cashmere goats. Subsequently, under the described above assay, bisulfite sequencing PCR (BSP-PCR) reactions were performed with BSP-primers (Table 1). After purifying with DNA purification kit (TaKaRa, China), the obtained BSP-PCR products were ligated to pMD18-T Vector (TaKaRa, China), followed by a propagation in competent Escherichia coli DH5α cells. For each analyzed stage (telogen, anagen, and catagen) of SHF cycle of cashmere goats, we sequenced ten positive clones. The resultant results were displayed by the QUMA procedure [20].

The SPSS 17.0 software (SPSS Inc., Chicago, IL, USA) was used to conduct the statistical analysis for the obtained data. We presented the analyzed results using the GraphPad Prism version 8.3.0 for Windows, GraphPad Software, San Diego, CA, USA ( www.graphpad.com, last access: November 27, 2023). The difference between analyzed groups was compared with Student’s t-test where p-value less than 0.05 was considered as significant difference statistically.

The circHECA (also known as chi_circ_2092) was previously identified from skin tissue of cashmere goat with a spliced length of 1,049-nt [8]. To define the transcription source of m⁶A-circHECA, we performed an alignment of its linear sequence against the goat genome datasets. As a result, the HECA on chromosome 9 was revealed to be its host gene (Figure 1a). As annotated well in National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov ), the goat HECA gene is composed of four exons, including exons 1, 2, 3 and 4. Whereas, m⁶A-circHECA is formed by reverse splicing of its entire exon 2 with position nos. 63,980,390–63,981,438 of the NC_030816.1 sequence (Figure 1a).

Based on in-silico analysis, four potential m⁶A modification sites were revealed in m⁶A-circHECA including m⁶A-213, m⁶A-297, m⁶A-780, and m⁶A-927 with the motif of GAACA, GGACU, GGACC, and GAACU, respectively (Figure 1b). These motif structures of m⁶A-circHECA are fully in line with the m⁶A motif: RRACH (R: A/G and H: A/C/U) previously found in linear RNAs [21]. On the other hand, bioinformatively, we also found that m⁶A-circHECA harbored potential binding target sites of several miRNAs, including chi-miR-449a-5p, chi-miR-129-3p, chi-miR-187, chi-miR-449b-3p, chi-miR-20b, and chi-miR-27a-5p (Figure 1b), which implies that these miRNAs may be implicated in the functional performance of m⁶A-circHECA in SHF physiological process of cashmere goats as shown in previous investigation [1,12].

Based on the above predicted m⁶A sites within circHECA sequence (Figure 1b; Figure 2a), we further performed a validation for each potential m⁶A modification site by Me-RIP along with qPCR technique (Me-RIP-qPCR). As shown in Figure 2a, the four potential m⁶A modification sites were confirmed in circHECA molecule including m⁶A-213, m⁶A-297, m⁶A-780, and m⁶A-927 (Figure 2a). Interestingly, it is thought that the m⁶A modifications sites within m⁶A-circRNA molecules are essential to their functional performance in biological cells [9,22]. Thus, it can be suggested that the verified m⁶A modifications of circHECA may mean further biological significance in regulating its functional roles SHFs of cashmere goats as reported in a recent publication by Yin et al [12]. On the other hand, as well known, circRNAs were originally regarded as non-coding RNAs that had no ability of encoding protein/peptide, however, increasing lines of evidence showed that some of circRNAs have capability of encoding a protein/peptide with biological significance [23]. Based on bioinformatics analysis, therefore, we assessed the encoding potential of m⁶A-circHECA by both CPAT and CPC programs with a screening of potential ORF within its sequence. As shown in Figure 2b, m⁶A-circHECA harbored one ORF of 960-nt with encoding a potential novel HECA-319aa protein (Figure 2b). Also, there is evidence that the m⁶A motifs are heavily implicated in driving an effective translation initiation of m⁶A modification circRNAs [9]. Taken together with our results, we strongly suggest that the potential novel HECA-319aa protein should be further identified with its functional validation, which might mean further biological significance in SHF physiology of cashmere goats.

In previous investigations, it has been demonstrated that the subcellular localization of a given circRNA can provide significant clues on its functional roles in biological cells [24]. In this study, we also evaluated the subcellular localization of m⁶A-circHECA in SHF stem cell of cashmere goats. As shown in Figure 2c, although the expression of m⁶A-circHECA was detected in both nucleus and cytoplasm of the analyzed SHF stem cells, it was detected mainly in cytoplasm of the analyzed cells (Figure 2c). On the whole, it is widely accepted that the majority of circRNAs located in cytoplasm regulate the availability of miRNAs binding to target mRNA molecules [24]. Based on in-silico analysis, therefore, we further investigated the binding features of m⁶A-circHECA with the above predicted six potential target miRNAs (chi-miR-449a-5p, chi-miR-129-3p, chi-miR-187, chi-miR-449b-3p, chi-miR-20b, and chi-miR-27a-5p), including binding structure, site types, and free energy (ΔG). As shown in Figure 2d, a good sequence matching structure can be formed between m⁶A-circHECA and its target miRNA with binding site type of 5 to 7 mer and free energy (ΔG) of −30.7 to −19.3 kcal/mol (Figure 2d). Thus, we speculate that the m⁶A-circHECA may exert its functional roles via miRNA mediated pathways in SHFs of cashmere goats.

It is well known that circRNAs can sequester target miRNAs to inhibit the binding of the miRNAs with corresponding target mRNAs thereby ultimately regulating the expression of protein-coding genes via a ceRNA network mechanism [10]. To explore the potential functional mechanisms of m⁶A-circHECA in SHFs of cashmere goats, here, we generated a ceRNA network of m⁶A-circHECA with potentially interactive miRNAs and their corresponding target genes based on the use of bioinformatics tools. As presented in Figure 3, six miRNAs were revealed to have potential interactive relationships with m⁶A-circHECA including chi-miR-449a-5p, chi-miR-129-3p, chi-miR-187, chi-miR-449b-3p, chi-miR-20b, and chi-miR-27a-5p. Whereas, these miRNAs may be further implicated in potentially regulating the expression of corresponding protein-coding genes (Figure 3). Interestingly, the chi-miR-449a-5p, chi-miR-129-3p, chi-miR-187, chi-miR-20b, and chi-miR-27a-5p were previously identified in hair follicles of cashmere goats [25]. Moreover, the chi-miR-129-3p and chi-miR-20b were recorded to be significantly upregulated at telogen hair follicles of cashmere goats compared with the counterpart of anagen with a p-value of 0.00066 and 0.00026, respectively [25]. As well known, the telogen is a relative quiescent stage of SHFs in cashmere goats with low proliferation rate and small volume [26]. Thus, it can be suggested that the chi-miR-129-3p and chi-miR-20b may play significant roles in establishing the optimal expression balance of protein-coding genes at telogen SHFs of cashmere goats, which may be essentially important for cashmere goat SHFs to maintain the quiescent status at telogen stage. Whereas this biological process may be ultimately regulated by m⁶A-circHECA via the miRNAs (chi-miR-129-3p and chi-miR-20b) meditated pathways.

On the other hand, some of the protein-coding genes in the generated ceRNA network have been verified to play essential roles in physiological process of hair follicle. As an example, the LEF1 was found to direct the fate of hair follicle stem cells through promoting the β-catenin translocation [27]. Also, it was reported that the inhibition of CDK6 along with CDK4 in skin tissue by a membrane-transducible TAT-p16INK4a protein impeded the growth and differentiation of hair follicle [28]. Here, both LEF1 and CDK6 were revealed to be potential target genes of chi-miR-449a-5p that might interact with m⁶A-circHECA of cashmere goat (Figure 3). Similarly, the FOXC1, AAK1, and ZNF704 were predicted as potential target genes of chi-miR-27a-5p, chi-miR-20b, and chi-miR-449b-3p, respectively. Previously, it was demonstrated that FOXC1 could not only maintain the niche of hair follicle stem cells, but also reinforce their quiescence to preserve the long-term potential of hair follicle regeneration [29]. Whereas the AAK1 and ZNF704 were also implicated in the development and growth of hair follicle [22]. Taken together with our results, it can be inferred that m⁶A-circHECA might play multiple roles in maintaining the normal SHF physiology of cashmere goat through precisely regulating the expression of its miRNA mediated target genes.

To gain a further insight into the potential functional mechanisms of m⁶A-circHECA, we performed a pathway enrichment upon its miRNA mediated target genes using the CluePedia procedure. As observed from Figure 4, the potential miRNA mediated target genes of m⁶A-circHECA were significantly enriched into multiple signaling pathways, including MAPK signaling pathway, autophagy, systemic lupus erythematosus, olfactory transduction, ribosome, Axon guidance, and Parathyroid hormone synthesis, secretion, and action (Figure 4). Of them, interestingly, several signaling pathways were previously verified to play important roles in the development and growth of hair follicles. For example, it was demonstrated that the activation of MAPK signaling pathway was heavily implicated in hair regeneration and growth with heightening expression of involved growth factors [30,31]. Also, there is evidence that the activation of axon guidance signal is required for the morphogenesis of hair follicles where it exerts function through driving the rearrangement of localized cells [32]. More recently, it was reported that the activation of autophagy pathways was also involved in promoting the hair growth of 5-Bromo-3,4-dihydroxybenzaldehyde along with activating Wnt/β-catenin pathway and inhibiting of the transforming growth factor-β (TGF-β) pathway in dermal papilla cells [33]. Thus, we inferred that the miRNA mediated regulatory genes of m⁶A-circHECA might be involved in the SHF physiological processes of cashmere goat, where their biological function might be ultimately subjected to the regulation of the m⁶A-circHECA molecules.

Over past few years, increasing evidence indicated that the interactions of non-coding RNAs and RNA binding proteins are widely implicated in a variety of cellular processes, such as, transcription, post-transcriptional regulation of gene expression and host defense against pathogens [34]. Accordingly, to further explore the potential interactional proteins of m⁶A-circHECA in SHFs of cashmere goat, we generated an interactional regulatory network of m⁶A-circHECA of cashmere goat with its potential RNA-binding proteins. As a result, there are six proteins having a direct binding relationship with m⁶A-circHECA of cashmere goats, including RBMY1A1, RBMX, SNRPA, ZRANB2, EIF4B, and FUS (Figure 5a). As shown in Figure 5a, moreover, each of them has been revealed to have multiple further regulatory relationships with other proteins. Although, it is not yet known whether these six proteins play directly functional roles in SHFs of cashmere goats, several of their proteins were implicated in hair follicle-related key signaling pathways, such as RBMX in Wnt/β-catenin signaling pathway [35], EIF4B in MAPK signaling pathway [36] and FUS in TGF-β signaling pathway [37]. Additionally, the SNRPA is previously known as an important paralog of SNRPB2 because of its structure closely related with SNRPB2 [38]. Although it is not yet known whether SNRPA plays similar functional roles with SNRPB2 in SHFs of cashmere goats, it was demonstrated that SNRPB2 was involved in Wnt signaling pathway [39]. Whereas it is well known that the Wnt signal play an important role in the growth and development of hair follicles [40].

Also, we performed a signaling pathway enrichment on the potential target proteins directly or indirectly regulated by m⁶A-circHECA molecules. As a result, we found that the analyzed regulatory proteins of m⁶A-circHECA were significantly enriched into multiple signaling pathways including T cell receptor signaling pathway, mRNA surveillance pathway, Spliceosome, MAPK signaling pathway, pathways in cancer, and viral carcinogenesis (Figure 4b). Of them, the T cell receptor signaling pathway was implicated in T cell-mediated developmental defect of hair follicles [41]. Also, the MAPK signaling pathway was demonstrated to play important roles in regeneration and growth of hair follicles [30,31]. Taken together, we speculate that the significantly enriched proteins may be essentially implicated in the SHF growth and development of cashmere goats which may be directly or indirectly regulated by m⁶A-circHECA molecules.

It is generally thought that the biogenesis of circRNAs is highly spatio-temporal specific across different cell types and developmental stages of tissues and organs, and their generation in cells is precisely controlled at multiple layers [42]. In previous investigation, it was recorded that the expression of m⁶A-circHECA along with its host gene HECA in SHFs of cashmere goats was significantly higher at both anagen and catagen than the counterpart of telogen [8]. Here, this was further validated in SHFs of nine individuals from Liaoning cashmere goats (Figure 6a; 6b). In addition, as shown in Figure 6a and b, m⁶A-circHECA and its host gene HECA linear mRNA exhibited a highly similar expression pattern in SHFs of cashmere goat during SHF cycles: telogen, anagen and catagen. Moreover, we found that a significantly positive correlation relationship existed in expression pattern between m⁶A-circHECA and its host gene HECA linear mRNA molecules (Figure 6c). Although the precise biological functions of HECA linear mRNAs in SHF physiology need to be further clarified in cashmere goats, it is thought that HECA is deeply implicated in the main aspects of cellular behavior (such as growth and differentiation) through Wnt/β-catenin signaling pathway [8], while the Wnt/β-catenin signaling pathway plays key role in the growth and development of hair follicle [11]. Thus, we inferred that the HECA along with its m⁶A-circHECA might form a novel regulatory layer in the SHF physiological process of cashmere goats.

It is well recognized that the promoter methylation of protein-coding gene is essentially implicated in its expression regulation without the alteration in promoter DNA sequence [43]. Moreover, emerging evidence has indicated that biosynthesis regulation of circRNAs is similar to the counterpart of protein-coding genes with their expression regulated by corresponding promoter methylation status [44]. These findings drove us to ask whether the promoter methylation may be implicated in the revealed expression changes of m⁶A-circHECA in SHFs of cashmere goats during SHF cycle. Therefore, we further investigated the promoter methylation status of m⁶A-circHECA host gene (HECA) at three SHF stages of cashmere goats including telogen, anagen and catagen. As observed from Figure 6d, a CpG island of 649-bp was revealed in promoter region directly upstream to the TSS of HECA gene. Within the CpG island region, we amplified a fragment with a length of 276-bp that contained 26 CpG sites and spanned multiple potential binding sites of transcription factors, such as AP-2α, Sp1, c-Fos, NF-1, c-Jun, and SRF (Figure 6e).

We performed the detection on methylation status of analyzed HECA gene promoter region in nine analyzed DNA samples pooled into three stage groups: telogen, anagen and catagen. As a result, the ratios of methylated CpG sites are 12.31%, 1.92%, and 2.69% at telogen, anagen and catagen, respectively (Figure 6f). Taken together with the express change of m⁶A-circHECA in SHFs of cashmere goats among telogen, anagen and catagen (Figure 6a), the telogen exhibited lower expression of m⁶A-circHECA, but had the higher methylation level in analyzed promoter region of HECA gene. In contrast, both anagen and catagen exhibited higher expression of m⁶A-circHECA, but had lower methylation level in analyzed promoter region of HECA gene. Thus, it appears that a higher methylation level of promoter region of HECA gene in SHFs of cashmere goats may cause the lower expression of m⁶A-circHECA (Figure 6a; 6f). Recently, a highly similar case was also reported in multiple myeloma cells where the decreasing expression of circRNA-ciRS-7 was found to be associated with promoter hypermethylation of its host gene LINC00632 [45]. On the other hand, however, it was also reported that some circRNAs could regulate the promoter methylation of their host genes via other epigentica factors [46]. As an example, the circRNA-FECR1 transcribed from FLI1 gene was found to recruit TET1 (a demethylase) binding to the promoter of its host gene further to cause the hypomethylation in cis [47]. Therefore, we strongly recommend further investigate whether the m⁶A-circHECA in turn regulates the promoter methylation level of its host gene (HECA) in SHFs of cashmere goats through forming a feedback loop mechanism.

Meanwhile, although the precise mechanism on how epigenetic machinery regulates the biogenesis of circRNAs in cells is largely unknown, taken together with our results, it can be inferred that the lower promoter methylation of HECA gene most likely is implicated in higher expression of m⁶A-circHECA at both anagen and catagen SHFs of cashmere goats, which might explain, at least in part, the changed expression of m⁶A-circHECA in SHFs of cashmere goat among the analyzed three stages: telogen, anagen and catagen (Figure 6a; 6d). These results provided further information from epigenetic layer for understanding the regulatory mechanism on m⁶A-circHECA expression in SHFs of cashmere goats.