Twelve 6-month-old male Wanxi Angora rabbits were used to establish the HF synchronization model. They were all housed under the same conditions, including temperature, and were fed the same diet (feed pellet and grass). Animals were reared in a controlled environment and had the same length of the hair coat phenotypes. The experimental procedures in this study were approved by the Animal Care and Use Committee of Yangzhou University.

To estimate the wool growth rate and to determine the onset of the anagen phase, the dorsal area of experimental animals was shaved with electronic clippers and entry into anagen was determined by the appearance of light pink skin and by hair regrowth. The length of the hair coat was measured, skin samples were collected after shaving, samples were fixed in 4% formaldehyde, and paraffin sections were stained with hematoxylin–eosin (HE) for histological observations. Longitudinal sections of the HFs showed the skin status and the phase of the HF cycle.

Rabbits were anesthetized via ear vein injections of 0.7% pentobarbital sodium (6 mL/kg), dorsal skin samples (1 cm²) were collected, and placed immediately in liquid nitrogen for RNA extraction. Iodine solution was applied on the wound to prevent bacterial infection. Samples were harvested at different phases of the HF cycle for gene expression profiling: growth (anagen), cessation (catagen), and rest (telogen). Three sample replicates were collected at days 90, 130, and 150 of the HF cycle for ncRNA and mRNA sequencing analysis.

Total RNA from nine samples was extracted from skin tissue using Trizol reagent (Invitrogen, Carlsbad, CA, United States), according to the manufacturer’s instructions. RNA degradation and contamination were monitored by running samples on 1% agarose gels. RNA purity was analyzed via a NanoPhotometer^® spectrophotometer (IMPLEN, CA, United States). RNA concentration was measured using the Qubit^® RNA Assay Kit and a Qubit^® 2.0 Fluorometer (Life Technologies, CA, United States). RNA integrity was assessed via the RNA Nano 6000 Assay Kit and a Bioanalyzer 2100 system (Agilent Technologies, CA, United States). lncRNAs and miRNAs were quantified following the same procedure used for conventional mRNAs. Quantification of circRNAs was performed adding an exonuclease to degrade non-circRNAs. Briefly, two samples containing the same amount of RNA were collected. In one sample, linear RNA was digested with RNase R (Cat. No. RNR07250, Epicentre Company, United States), leaving only the circRNAs, while the other sample was not treated with RNase R. The two RNA samples were reverse transcribed. The samples subjected to RNase treatment were used to detect circRNAs, whereas the untreated samples were used to detect β-actin.

A total amount of 3 μg of RNA per sample was used for lncRNA sequencing and of 5 μg for circRNA sequencing. First, ribosomal RNAs were removed with the Epicentre Ribo-zero^TM rRNA Removal Kit (Epicentre, United States) and the rRNA-depleted samples were purified by ethanol precipitation. Subsequently, sequencing libraries were generated using the rRNA-depleted RNA and the NEBNext^® Ultra^TM Directional RNA Library Prep Kit for Illumina^® (NEB, United States), following the manufacturer’s recommendations. First strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNaseH). Second strand cDNA synthesis was performed using DNA Polymerase I and RNase H. After adenylation of the 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. To select cDNA fragments with a preferential length of 150∼200 bp, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, MA, United States). Then, 3 μl of USER Enzyme (NEB, United States) was used with size-selected, adaptor-ligated cDNA before the PCR. Finally, the PCR products were purified (AMPure XP system) and the library quality was assessed with the Agilent Bioanalyzer 2100 system.

A total amount of 3 μg of RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using the NEBNext^® Multiplex Small RNA Library Prep Set for Illumina^® (NEB, United States), following the manufacturer’s recommendations. Briefly, NEB 3′ SR Adaptor was directly and specifically ligated to the 3′ end of miRNAs, siRNAs, and piRNAs. After the 3′ ligation reaction, the SR RT Primer was hybridized to the excess of 3′ SR Adaptor, transforming the single-stranded DNA adaptor into a double-stranded DNA molecule. Then, the 5′ ends adapter was ligated to the 5′ ends of the miRNAs, siRNAs, and piRNAs. First strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H–). DNA fragments of 140–160 bp length (the length of small non-coding RNAs plus the 3′ and 5′ adaptors) were recovered and dissolved in 8 μL of elution buffer. Finally, the library quality was assessed using the Agilent Bioanalyzer 2100 system and DNA High Sensitivity Chips.

Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina), according to the manufacturer’s instructions. After cluster generation, the lncRNA and circRNA libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired end reads were generated. The miRNA library was sequenced on an Illumina Hiseq 2500 platform and 50 bp single-end reads were generated.

For lncRNA and circRNA sequencing, raw data (raw reads) in fastq format were first processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N, and low-quality reads from the raw data. For miRNA sequencing, raw data (raw reads) in fastq format were first processed through custom perl and python scripts. In this step, clean data (clean reads) were obtained by removing reads containing ploy-N, with 5′ adapter contaminants, without 3′ adapter or the insert tag, containing ploy A or T or G or C, and low-quality reads from raw data. At the same time, the Q20, Q30 scores, and GC-content of the raw data were calculated. A specific length range from the clean reads was selected to conduct all the downstream analyses, based on clean data of high quality.

For lncRNA and circRNA sequences, the reference genome (Oryctolagus cuniculus genome obtained from Ensembl OryCun2.0) and annotation files were directly downloaded from the genome website. An index of the reference genome was built using bowtie2 (Langmead and Salzberg, 2012), and paired-end clean reads were aligned to the reference genome using HISAT2 v2.0.4 (Pertea et al., 2016). Also, the small RNA tags were mapped to the reference sequence with bowtie2 (Langmead and Salzberg, 2012) without mismatch to analyze the expression and distribution of miRNA sequences in the reference genome.

The mapped lncRNA and mRNA reads from each sample were assembled by means of StringTie (v1.3.1) (Pertea et al., 2016), following a reference-based approach. The circRNAs were detected and identified using find_circ (Memczak et al., 2014). Alignment of the small RNA tags to miRBase20.0 identified known Oryctolagus cuniculus and Mus musculus (near-source species) miRNAs. Mirdeep2 software (Friedländer et al., 2011) was used to identify potentially novel miRNAs and to draw the secondary structures and the characteristics of the hairpin structures of miRNA precursors.

Cuffdiff (v2.1.1) was used to calculate fragments per kilo-base millions of exon per million fragments mapped (FPKM) of both lncRNAs and mRNA in each sample (Trapnell et al., 2010). FPKMs of genes were computed by summing the FPKMs of transcripts in each gene group. Also, the raw counts were first normalized using transcripts per million (TPM) (Zhou et al., 2010) and normalized expression levels = (read count^∗1,000,000)/lib size (lib size is the sum of circRNA read counts). This was used to determine the circRNA expression levels. On the other hand, miRNA expression levels were estimated by TPM based on the following criteria: Normalization formula: Normalized expression = mapped read count/total reads^∗1,000,000. The differential expression of ncRNAs was determined using the DESeq R package (1.10.1) (Wang et al., 2010).

In cis regulation, lncRNAs can act on neighboring target genes. Coding genes 10 k/100 k upstream or downstream of the lncRNA gene were searched for and their function was analyzed. For trans regulation, lncRNAs and their target genes were analyzed based on their expression levels. The correlation between lncRNAs and coding gene expression levels were calculated with custom scripts; then, the genes from different samples were clustered using WGCNA (Langfelder and Horvath, 2008) to search for common expression modules and to analyze the function via functional enrichment analysis. The target genes of miRNAs and miRNA target sites in exons of circRNA loci were identified using miRanda (version 3.3a, main parameter: -sc 140; -en -10; -scale 4; -strict) (Enright et al., 2004). Differentially expressed (DE) ncRNAs were annotated by gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses to investigate their biological functions. Briefly, GO analysis was applied to elucidate genetic regulatory networks of interest by forming hierarchical categories according to the molecular function (MF), cellular component (CC), and biological process (BP) aspects of the differentially expressed genes¹. KEGG pathway analysis was performed to explore the significantly enriched pathways of DE genes².

Eight mRNAs, four lncRNAs, and five circRNAs associated with skin and the HF cycle were selected for validation by qRT-PCR analysis. Approximately 1μg of total RNA was used to synthesize cDNA using HiScript II Q Select RT SuperMix for qPCR (Vazyme). qRT-PCR was performed using the AceQ qPCR SYBR^® Green Master Mix (Vazyme), according to the manufacturer’s instructions, and data were analyzed via QuantStudio^® 5 (Applied Biosystems). The specific primer sequences are listed in Supplementary Table S1. The expression levels were calculated using the 2^-ΔΔCt method (Schmittgen and Livak, 2008), with glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as reference gene.

To confirm the miRNA transcriptome data, three miRNAs were selected for qRT-PCR analysis. Approximately 2μg of total RNA was used to synthesize cDNA after adding a poly (A) tail to the 3′ end of the miRNAs using the miRcute Plus miRNA First-Strand cDNA Synthesis Kit (Tiangen). qRT-PCR was performed using the miRcute miRNA qPCR Detection Kit (SYBR Green), according to the manufacturer’s instructions. The specific primers were designed by Beijing Tiangen Co., Ltd. and the product code sets are listed in Supplementary Table S1. The U6 small nuclear RNA gene was chosen as internal control. The expression levels were calculated using the 2^-ΔΔCt method (Schmittgen and Livak, 2008), and the results of the experiments were normalized to the expression levels of the constitutively expressed U6 gene.

To investigate the role and interactions between ncRNAs and mRNAs during the HF cycle, ncRNAs regulatory networks were constructed. For the interaction network of lncRNA–miRNA, DE lncRNAs were filtered out according to the homology between lncRNA and miRNA precursor; then, the targeted relationships between lncRNA and miRNA were predicted by miRanda. Then, the regulatory networks of lncRNA–miRNA–mRNA pairs and circRNA–miRNA–mRNA pairs were constructed according to the following steps: (i) the ncRNAs and mRNAs that were upregulated or downregulated were retained; (ii) the interactions of lncRNA–miRNA, miRNA–mRNA, and miRNA–circRNA were predicted by miRanda, which predicts miRNA binding seed sequence sites, as well as overlapping the same miRNA binding site in lncRNAs, circRNAs, and mRNAs; (iii) The lncRNA–miRNA–mRNA pairs network covered two cases: one was the upregulated lncRNA-downregulated miRNA-upregulated mRNA, the other was the downregulated lncRNA-upregulated miRNA-downregulated mRNA. The circRNA–miRNA–mRNA pairs network covered two cases: one was the upregulated circRNA-downregulated miRNA-upregulated mRNA, the other was the downregulated circRNA-upregulated miRNA-downregulated mRNA. Cytoscape software was used to build and visually display the networks.

The dual-luciferase reporter system E1910 (Promega, Madison, WI, United States) was used to perform luciferase activity assays. The miR-320-3p mimic and miR-320-3p negative control mimics were purchased from Shanghai GenePharma Co., Ltd. Wild-type luciferase reporter vectors (pMir-HTATIP2-3’UTR-WT, pMir-LNC_002919-WT, and pMir-novel_circ_0026326-WT) were constructed using the primers shown in Supplementary Table S2. Their substitution mutants (pMir-HTATIP2-3’-UTR-MUT, pMir-LNC_002919-MUT, and pMir-novel_circ_0026326-MUT) were synthesized by Beijing Tsingke Co., Ltd. Briefly, the skin fibroblast cells of rabbit (RAB-9, ATCC^® CRL-1414^TM) were cultured in 24-well tissue culture plates. Cells were co-transfected with the pMir-report luciferase reporter, the miRNA (miR-320-3p) mimics and pRL-TK using Lipofectamine^TM 2000 (Invitrogen). After 48 h of culture at 37°C, transfected cells were lysed with 100 μl of passive lysis buffer. Next, 20 μl of lysates were mixed with 100 μl of LAR II, and firefly luciferase activity was measured by using a luminometer. As an internal control, 100 μl of Stop & Glo reagent was added to the sample. Firefly luciferase activity was normalized to the corresponding Renilla luciferase activity.

For the HF cycle synchronization model, Angora rabbits were used. The obtained observations showed that the length of the hair coat increased steadily until day 110. Between days 120 and 150, the growth rate of wool declined rapidly. Then, between days 160 and 180, the wool recovered and once again showed an increased growth rate (Figure 1A). Histological analysis showed rapid growth of the hair shaft and increasing depth of the HF between days 0 and 110. Then, the growth of the hair shaft and the depth of the HF decreased between days 120 and 130. Finally, the hair shaft started to fall off and the hair bulbs atrophied between days 140 and 150. After the HF cycle ended, a new HF appeared, the growth of the hair shaft recovered and the HFs moved into a new cycle (Figure 1B). In conclusion, the hair cycle of Angora rabbits is characterized by an anagen phase between days 0 and 110, a catagen phase between days 120 and 130, and a telogen phase between days 140 and 150.

A summary of the lncRNA-seq, miRNA-seq, and circRNA-seq data from the three HF cycle phases is shown in Supplementary Table S3, indicating the relatively high quality of the transcriptome data. The lncRNA-seq, miRNA-seq, and circRNA-seq data were deposited in the Short Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under the bioproject numbers PRJNA479733, PRJNA495446, and PRJNA495449. DE ncRNAs and mRNAs were analyzed using Cuffdiff software with a criterion of p < 0.05. Volcano plots, clustering maps, and Venn diagrams were used to illustrate the distribution of the DE ncRNAs and mRNAs between the three groups (Figure 2–5). Table 1 summarizes the number of DE ncRNAs and mRNAs. Differential expressions of 111 lncRNAs (60 upregulated and 51 downregulated), 247 circRNAs (128 upregulated and 119 downregulated), 97 miRNAs (38 upregulated and 59 downregulated), and 1,168 mRNAs (750 upregulated and 418 downregulated) were found between the three HF cycle stages. Complete information on all DE lncRNAs, circRNAs, miRNAs, and mRNAs is listed in Supplementary Tables S4–S7. Several lncRNAs were found to be associated with the HF cycle, such as LNC_002694, LNC_002919, LNC_003354, LNC_003790, LNC_008354, LNC_008931, and LNC_005484, which could regulate gene expression by recognizing their target mRNAs. Based on analysis of their biological function, the candidate lncRNAs associated with the HF cycle are listed in Supplementary Table S8. Moreover, analysis of the relationships between circRNAs and genes allowed identification of novel_circ_0004876, novel_circ_0005177, novel_circ_0026326, novel_circ_0034968, and novel_circ_0036671, which may play a role during the HF cycle. In addition, several miRNAs, including miR-128-3p, miR-200a-3p, miR-27a-3p, miR-30e-5p, and miR-320-3p; mRNAs, such as BMP2, CSNK2B, KRT17, LAMB1, FZD4, SMAD2, HTATIP2, and SIAH1 were identified to play pivotal roles during the HF cycle and during skin development.

To validate the lncRNAs, mRNAs, miRNAs, and circRNAs differential expression results, the relative expression of four DE lncRNAs (LNC_002694, LNC_002919, LNC_003354, and LNC_005484), five DE circRNAs (novel_circ_0004876, novel_circ_0005177, novel_circ_0026326, novel_circ_0034968, and novel_circ_0036671), four DE miRNA (miR-128-3p, miR-200a-3p, miR-27a-3p, and miR-320-3p), and eight DE mRNAs (BMP2, CSNK2B, FAM45A, FUOM, HTATIP2, KRT17, ME1, and SIAH1) were measured by qRT-PCR (Figure 6–9). The qRT-PCR results were consistent with the transcriptome sequencing data.

lncRNAs can regulate neighboring protein-coding genes; therefore, a colocalization threshold of 100 kb upstream or downstream of lncRNAs was set for the GO and KEGG analyses. Several GO terms were found that were significantly enriched in the three experimental groups (Supplementary Table S9), including skin and HF-related GO terms like HF development (GO: 0001942), hair cycle (GO: 0042633), hair cycle process (GO: 0022405), regulation of HF development (GO: 0051797), and skin morphogenesis (GO: 0043589), among others. The top 20 KEGG pathways associated with DE lncRNAs between days 90, 130, and 150 of the HF cycle based on the function of colocalized mRNAs (Supplementary Figure S1) and co-expressed mRNAs (Supplementary Figure S2) included the Wnt signaling pathway, TGF-β signaling pathway, MAPK signaling pathway, and JAK/STAT signaling pathway.

In addition, based on the relationship between circRNAs and genes, GO analysis of genes producing DE circRNAs was performed (Supplementary Table S10). The GO terms identified HF development (GO: 0001942), hair cycle process (GO: 0022405), hair cycle (GO: 0042633), and skin development (GO: 0043588), which were all related to skin and HF development. The top 20 KEGG pathways associated with genes producing DE circRNAs between 90, 130, and 150 days (Supplementary Figure S3) of the HF cycle were likewise related to skin and HF development, such as the Hedgehog signaling pathway, Wnt signaling pathway, and MAPK signaling pathway.

Furthermore, GO enrichment analysis of genes targeted by DE miRNA (Supplementary Table S11) identified GO terms related to HF development, such as HF morphogenesis (GO: 0031069), negative regulation of HF development (GO: 0051799), and regulation of HF development (GO: 0051797), among others. The top 20 KEGG pathways associated with DE miRNAs are shown in Supplementary Figure S4. They include pathways related to HF cycle, such as the Hedgehog signaling pathway, NF-κB signaling pathway, and JAK/STAT signaling pathway.

Finally, GO and KEGG analyses of DE mRNAs are shown in Supplementary Table S12. The GO terms identified include, for example, skin morphogenesis (GO: 0043589) and positive regulation of HF development (GO: 0051798). The top 20 enriched KEGG pathways for DE genes between the different stages of the HF cycle are shown in Supplementary Figure S5. These KEGG pathways include the Wnt signaling pathway, the MAPK signaling pathway, and the TGF-β signaling pathway, which participate in skin development and HF cycle. Differentially expressed genes between days 90, 130, and 150 of the HF cycle, as well as their biological functions, are listed in Supplementary Table S13.

Study of the relationship between ncRNAs and mRNAs may increase our understanding of the molecular mechanisms operating during skin development and HF cycle. According to the competing endogenous RNA (ceRNA) regulatory hypothesis, ncRNAs and mRNAs can compete for the same miRNAs, resulting in additional layers of regulation of gene expression. Based on the analysis of DE lncRNAs, circRNAs, miRNAs, and mRNAs, a network of lncRNAs and miRNAs was first constructed (Figure 10). In lncRNA-miRNA-mRNA regulatory networks, miRNA may act as the center, lncRNA as the decoy, and mRNA as the target, which suggests that lncRNAs could act as miRNA sponges to regulate gene expression (Figure 11). In addition, certain circRNAs can competitively bind miRNAs and act as miRNA sponges; therefore, circRNA-miRNA-mRNA triads were constructed with the circRNA as the docoy, miRNA as the center, and mRNA as the target (Figure 12).

LNC_002919 and novel_circ_0026326 were identified as ceRNAs for miR-320-3p, which targets HTATIP2. A dual-luciferase reporter system was used to verify the binding relationships between the identified lncRNA and miRNA, circRNA and miRNA, and mRNA and miRNA. Luciferase assay showed that miR-320-3p could decrease luciferase activity by binding to sites on LNC_002919, novel_circ_0026326, and the HTATIP2 3′UTR (Figure 13). The interactions between ncRNAs and mRNA suggest the existence of novel regulatory mechanisms during skin development and HF cycle.