miR-29a-5p Inhibits Prenatal Hair Placode Formation Through Targeting EDAR by ceRNA Regulatory Network

Hair placode formation is an important stage of hair follicle morphogenesis and it is a complex process facilitated by non-coding RNAs. In this study, we conducted whole transcriptome sequencing analysis of skin, heart, liver, lung, and kidney tissues of day 41 (E41) normal and hairless pig embryos, and respectively detected 15, 8, and 515 skin-specific differentially expressed (DE) lncRNAs, miRNAs, and mRNAs. Furthermore, 18 competing endogenous RNA (ceRNA) networks were constructed. Following weighted gene co-expression network analysis (WGCNA) of stages E39, E41, E45, E52, and E60, between normal and hairless pig embryos, only two ceRNAs (lncRNA2162.1/miR-29a-5p/BMPR1b and lncRNA627.1/miR-29a-5p/EDAR) that showed period-specific differential expression in E41 skin were retained. Dual-luciferase reporter assays further indicated that EDAR was a direct, functioning target of miR-29a-5p and that no binding site was found in BMPR1b. Moreover, miR-29a-5p overexpression inhibited the mRNA and protein expression of EDAR while no significant differential expression of BMPR1b was detected. In addition, over-expressed lncRNA627.1 reduces the expression of miR-29a-5p and increase EDAR expression while inhibits lncRNA627.1 resulted in a opposite expression trend. Cell proliferation result demonstrated that lower expression of EDAR and lncRNA627.1 inhibited hair placode precursor cells (HPPCs) proliferation in a manner similar to that shown by over-expressed miR-29a-5p. This study identified that miR-29a-5p inhibited HPPCs proliferation via the suppression of EDAR expression in the EDA/EDAR signaling pathway, while lncRNA627.1 rescues EDAR expression. Our study provides a basis for a better understanding of the mechanisms underlying the ceRNA complex, miR29a-5p/EDAR/lncRNA627.1, that could regulate hair placode formation, which may help decipher diseases affecting human hair.


INTRODUCTION
Hair follicle (HF) development is a complex morphogenetic process which occurs in two stages: prenatal hair morphogenesis development, and postnatal hair cycle development (Choi, 2018). Embryonic hair follicle morphogenesis development, which involves mesenchymal and epithelial interactions, begins with the formation of a hair placode in the early epidermis (Duverger and Morasso, 2009;Rishikaysh et al., 2014;Saxena et al., 2019). Therefore, placode formation is key to successful hair follicle morphogenesis. Morphological studies have indicated that hair placodes are formed on day 14.5 (E14.5) in mouse embryos, and during the 12th week [E55-E65 (~E60)] in humans and Cashmere goats (Duverger and Morasso, 2014;Gao et al., 2016). Previous findings of ours have indicated that induction stage of pig hair follicle morphogenesis took place at E41 which the characteristic structure of this stage is hair placode formation (Jiang et al., 2019;Jiang et al., 2021). Thus, in pigs, E41 is the critical point which pig hair follicle morphogenesis begin to develop and hair placode start to form in this point and formation. Placode formation in developing embryonic skin requires stepwise signaling between the epithelial epidermis and the mesenchymal dermis. In recent years, many functional studies have uncovered the essential roles played by major signaling pathways, such as the Wnt/β-catenin, EDA/EDAR/NF-κB, and BMP pathways, in placode formation (Rishikaysh et al., 2014;Saxena et al., 2019). Moreover, numerous key genes that play a crucial role in hair placode formation, including Wnt10b (Andl et al., 2002), CTNNB1 (Gay et al., 2015), LHX2 (Tomann et al., 2016), and EDA/EDAR (Zhang et al., 2009) which induce and maintain placode formation, as well as BMP4 (St-Jacques et al., 1998) and DKK4 (Fliniaux et al., 2008) which exert inhibitory effects on placode formation, have been identified.
Previous studies were primarily focused on gene function associated with the genetic regulation of placode formation. However, in gene regulation process, whether noncoding RNAs [long noncoding RNA (lncRNA) and microRNAs (miRNA)] are involved by influence gene function at various stages of HF morphogenesis and HF cycling development remains unclear (Andl and Botchkareva, 2015;Ge et al., 2019). Several specific miRNAs that play multiple roles in the regulation of HF development include miR-22-5p (Yan et al., 2019a) and miR-203a-3p (Luo et al., 2021) in hair follicle stem cells (HFSCs), miRNA-203 (Ma et al., 2021) and miR-199a-5p  in hair follicle development, as well as miR-195-5p in dermal papilla (Zhu et al., 2018), and miR-218-5p in hair shaft growth (Zhao et al., 2019), among others. Furthermore, biogenesis of canonical miRNAs depends on the cytoplasmic processing of pre-miRNAs to mature miRNAs by the Dicer endoribonuclease (Vishlaghi and Lisse, 2020) while Dicer deletion in epidermal result in hair follicle stem cell markers and degenerated (Andl et al., 2006), hair shaft depigmented (Liu et al., 2018a), failure of dermal papilla or hair-follicle maintenance (Teta et al., 2012), a decline in cell proliferation and increased apoptosis in hair follicle (Yi et al., 2006;Yi et al., 2009). In addition, lncRNAs, such as MSTRG.223165/miR-21/SOX6 in wool HF development (Zhao et al., 2020), as well as lncRNA5322/miR-19b-3p/MAPK1 (Cai et al., 2019) and XR_310320.3/chi-miR-144-5p/HOXC8  in hair follicle stem cells (HFSCs), also act as competitive endogenous RNAs (ceRNAs) to sponge miRNA and relieve the inhibitory effects exerted by these miRNAs on target genes during the regulation of HF development. Although the role of ncRNAs in HF development has been gradually deciphered, molecular mechanisms underlying the role played by ncRNAs, particularly the ceRNA regulatory networks, in the formation and maintenance of placodes remain unclear. Although a few studies have revealed that lncRNA2437 (Yue et al., 2016), miRNA195 (Liu et al., 2018b) and XLOC297809 (Nie et al., 2018) are involved in hair follicle initiation and placode induction, the genetic mechanisms underlying these processes have not yet been elucidated. A previous study of ours showed that hairless pigs were mainly caused by the blockage of hair placode formation at E41 in induction stage of pig hair follicle morphogenesis (Jiang et al., 2019). However, the precise molecular basis of blocked hair placode formation in pigs remains elusive, while the role of miRNA and lncRNA in hair placode formation remains unexplored as well.
Therefore, the objective of the current study was to screen out key candidate ceRNAs which may exert significant effects on hair placode formation via whole-transcriptome sequencing. Different tissues (skin, heart, liver, lung, and kidney) at E41 and different periods (E39, E41, E45, E52, and E60) in skin were sequenced to screened tissue-and periods-specific expression RNAs at E41 skin. Furthermore, the ceRNA based regulatory mechanism underlying hair placode formation was verified at the cellular level. Thus, our study of pig models should enhance future exploration of the hitherto unrecognized role played by the ceRNA network in hair placode formation in humans.

Ethics Approval
All animal procedures were evaluated and authorized by Institutional Animal Care and Use Committee (IACUC). The whole procedure for samples collected was carried out in strict accordance with the protocol approved by the IACUC at the China Agricultural University. The IACUC of the China Agricultural University specifically approved this study (permit number DK996).

Animals and Phenotype
E41 is a crucial stage for prenatal placode formation during hair follicle morphogenesis. In this study, tissues, including skin, heart, liver, lung, spleen, and kidney, from hairless (n = 3) and normal (n = 3) pig embryos at day 41 of gestation (E41). One embryo each from a litter, were sampled. In addition, mRNA sequencing of skin from hairless (n = 3) and normal (n = 3) embryos, one each from a litter, were also sampled at days 39 (E39), 45 (E45), 52 (E52), and 60 (E60) of gestation. The phenotype of every sample was identified following Hematoxylin and Eosin (HE) staining (Supplementary Figure  S1). Hair follicle density was measured and those embryos with less than one hair follicle per cm 2 and larger than four hair follicles per cm 2 were determined as hairless and normal, respectively (Jiang et al., 2019).

RNAs (mRNA and ncRNA) Sequencing and Analysis
Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, United States), according to the manufacturer's instructions. Quantity and purity were analyzed using a Bioanalyzer 2,100 and an RNA 6,000 Nano Labchip Kit (Agilent, Palo Alto, CA). Only RNA samples with suitable RNA electrophoresis results (28S/18S ≥ 1.0) and RNA integrity number (RIN) ≥ 7.5 could be analyzed further. The skin, heart, liver, lung and kidney tissues of normal (n = 3) and hairless (n = 3) embryos from the same litter were selected for mRNA, lncRNA and miRNA sequencing at E41. Embryo skins (normal = 3, hairless = 3) from the other four periods, namely E39, E45, E52, and E60, were selected for mRNA sequencing. mRNA and lncRNA were sequenced on Illumina HiSeq 4,000 and miRNA on HiSeq 2,000 platforms, which generated 150-bp paired-end and 50-bp single-end reads, respectively. Differentially expressed (DE) ncRNAs and mRNAs were identified using HISAT2 (Kim et al., 2019) and DESeq2 (Love et al., 2014) with an adjusted p value (false discovery rate, FDR) < 0.05 and a log 2 |FoldChange| > 1.

mRNA Sequence Analysis
An Illumina high-throughput platform was used for mRNA sequencing and raw data was obtained using FastQC v0.11.9 . Quality control of the raw reads was conducted using NGSQCtoolkit_v2.3.3 (Patel and Jain, 2012 (Pertea et al., 2016) were performed to obtain clean reads aligned to the reference genome and mRNA reads were assembled for each sample. DESeq2 (Love et al., 2014) package was used for differential gene expression analysis.

lncRNA-Seq Analysis
Quality control and mapping methods of lncRNA-seq data were similar to those used for mRNA-seq data. Candidate lncRNAs were selected using the following conditions: 1) Transcript length ≥200 and exon number ≥2, and Fragments Per Kilobase per Million reads (FPKM) >0.5; 2) Minimal read coverage ≥3 in at least one sample; 3) Filter known non-lncRNA annotations; and 4) Classify selected candidate lncRNAs, using the following four tools to predict coding potential: CPC2 (Kang et al., 2017), CNCI (Sun et al., 2013), PhyloCSF (v20121028) (Lin et al., 2011), and Pfam Scan databases (Finn et al., 2016). DESeq2 package was used to identify differentially expressed lncRNA between different groups at p adj < 0.05. The lncRNA function depends on protein coding genes via cis and trans acting elements. Cisacting (in which lncRNAs act on adjacent genes within a 100 kb distance) and trans-acting (in which the Pearson correlation of mutual expression levels is ≥0.95 or ≤−0.95) are widely adopted to forecast lncRNA gene interactive pairs. Bedtools (Quinlan and Hall, 2010) was used to identify neighboring genes approximately 100 kb upstream and downstream of differentially upregulated and downregulated lncRNAs, respectively.

Weighted Gene Co-Expression Network Analysis
The mRNA sequencing of 30 embryo skins from five different stages (WGCNA_skins_periods; E39, E41, E45, E52, and E60) of normal and hairless pigs were integrated into the gene expression matrix for WGCNA. Then mRNA-seq data, which contained the expression data of 14,555 genes and 15,160 genes (Sum counts >10) for WGCNA_E41_tissues and WGCNA_skins_periods, was used to construct the gene expression matrix. Variance stabilizing transformation (VST) (Anders and Huber, 2010)  gene expression matrices respectively. The construction of two gene co-expression networks was based on the WGCNA package (Pei et al., 2017). Gene co-expression networks must conform to scale-free characteristics and obey power law distribution. Following sensitivity analysis of scale-free topology, the soft threshold power parameters of the two networks were set at 10 and 12, respectively (Ravasz et al., 2002).

Cell Culture
Approximately 2 cm 2 skin samples from the lateral backsides of the pig embryo at days 41 gestation were surgically removed following standard, aseptic procedures. The tissues were cut at the interface of the sub cutaneous adipose layer and dermis. The remaining tissues were sub merged in 0.25% dispase II (Gibco, Grand Island, United States) overnight at 4°C to separate the epidermis from the dermis, following which the epidermis was cut into small blocks of approximately 1 mm 3 in size, which were used as explants to culture hair placode precursor cells (HPPCs

Cell Proliferation Assay
Cell proliferation was assessed using a CCK-8 Cell Counting kit (Beyotime, Shanghai, China). FBCs were seeded in a 96-well plate (2 × 10 3 cells per well) and then cultured for 12, 24, 36, 48, 60, or 72 h before adding 10 μl of CCK-8 (5 mg/ml) to the culture medium in each well. After incubating for 2 h, absorbance at 450 nm was measured using a Thermo-max microplate (ThermoFisher, Massachusetts, United States) reader. All experiments involving each transfection were performed in triplicate. HPPCs were transferred to culture medium with 50 µM 5-ethynyl-20-deoxyuridine (Edu, Beyotime, Beijing, China) for 2 h at 37°C after 36 h transfection. Afterwards, cells were fixed in 4% paraformaldehyde for 15 min at room temperature (RT), and then permeabilized with 0.3% Triton X-100 for 10 min. To block non-specific binding, cells were incubated in blocking buffer (PBS containing 3% bovine serum albumin, 0.3% Triton X-100) for 1 h at room temperature. Next, the cells were incubated with a solution containing 10 mM Edu for 30 min in the dark. Nuclei were stained with 10 μg/ml 4, 6-diamidino-2-phenylindole (DAPI, Solarbio, Beijing, China) solution in the dark for 10 min. A Leica SP8 confocal microscope was used to capture three randomly selected fields to visualize the number of Edustained cells.

Western Blotting
Western blotting (WB) was used to detect protein quantity. In brief, the sample was lysed in RIPA lysis solution, to extract total protein, and denatured. Next, the sample proteins were separated via SDS-PAGE gel electrophoresis under conditions involving a constant voltage of 100 V. In this study, the proteins were transferred from SDS-PAGE gel to PVDF membranes using a semi-dry transfer method. Then the membranes were blocked for 4 h, and incubated with primary antibodies against EDAR (ABclonal, 1:2000) and BMPR1b (ABclonal, 1:1500) overnight at 4°C. Next, the membranes were washed thrice with TBST (Trisbuffered saline) and incubated with secondary antibodies (1: 1200 dilution) at room temperature for 1.5 h. Finally, a BeyoECL plus kit (Beyotime) was used to detect the protein signal. GAPDH (ABclonal, 1:1500) was used as a reference protein and blots were analyzed using IPWIN software.

Quantitative RT-PCR Analysis
Total RNA was isolated from dorsal pig skin with TriZol (Invitrogen, SanDiego, CA) following the manufacturer's instructions. Next, cDNA was reverse transcribed using a PrimeScriptTM RT reagent kit with gDNA Eraser (Takara, Kyoto, Japan). RT-PCR was performed with LightCycler 480 SYBR Green I Master (Roche, Mannheim, Germany) mix on a LightCycler 480 real-time PCR system. GAPDH was used as a normalized control and relative gene expression was calculated based on the 2 −ΔΔCT formula. Measurements were recorded in triplicate. The following PCR conditions used: 95°C "hot start" for 10 min; 35 cycles of 95°C for 10 s, 60°C for 10 s, and 72°C for 10 s; and 72°C for 5 min. Primer sequences are provided (Supplementary Table S3). Differences between the gene expression of normal and hairless pigs were measured via a t test. A summary of the descriptive statistics of lncRNA-seq, miRNAseq, and mRNA-seq data pertaining to different tissues, indicating the relatively high-quality of transcriptome data in this study is shown (Supplementary Table S1). Furthermore, mRNA and ncRNA (lncRNA and miRNA) sequencing of skin, heart, liver, lung, and kidney tissues of three normal and three hairless pig embryos at E41 stage were performed, and the number of differentially expressed (DE) mRNAs and ncRNAs are shown ( Table 1). The DE RNAs (mRNAs, miRNAs, and lncRNAs) are illustrated in detail (Supplementary Figure S2) and sample-to-sample correlation in the five tissues was high (Supplementary Figure S2).

Specific Differential Expression of mRNA and ncRNAs in Pig Embryos Skin at E41
Overlap between DE mRNAs, DE miRNAs and DE lncRNAs are illustrated ( Figures 1A-C). In addition, 515 DE mRNAs ( Figure 1A), 8 DE miRNAs ( Figure 1B), and 15 DE lncRNAs ( Figure 1C) were specifically differentially expressed between E41 skins of hairless and normal embryos. Next, qRT-PCR indicated that expression levels of randomly selected DE mRNA and ncRNA transcripts were highly consistent (Supplementary Figure S3), indicating the reliability of our RNA-seq data. Eight GO terms ( Figure 1D) of the 515 DE mRNAs were significantly enriched in skin development and epidermal development and differentiation, while four hair follicle related signaling pathways, including the Wnt signaling pathway and the BMP signaling pathway, were identified ( Figure 1E). Based on GO and pathway enrichment analysis, 128 of the 515 DE mRNAs (Supplementary Table S2) were selected which play important roles in HF development and hair placode formation. A total of 1,530 cis or trans target genes were predicted for the 15 DE lncRNAs in skin (Table 1), and among that, 68 target genes overlapped with those for the 128 DE mRNAs (Supplementary Figure S4A). Similarly, 1296 overlapped target genes were predicted for the 8 DE miRNAs by miRanda (3,452 target genes) and RNAhybrid (4,507 target genes); (Supplementary Figure S4B). Furthermore, epithelium development, the Wnt signaling pathway, the EDA/EDAR signaling pathway, and the BMP signaling pathway, among others, were enriched in these overlapped target genes (Supplementary Figures S3D,F), of which 60 genes overlapped with those of 128 DE mRNAs ( Figure 1F; Supplementary Figure S4C). Thirty-one genes that overlapped between normal and hairless pig embryos at E41 stage skin tissue were identified based on DE mRNAs and target genes of DE miRNAs and lncRNAs ( Figure 1F). These were enriched in hair follicle development, including skin epidermis development and epithelial cell differentiation ( Figure 1G). Similarly, enriched KEGG pathways included the classical HF morphogenesis pathways such as the Wnt, TGF-β, and EDA/EDAR signaling pathways ( Figure 1H). As well, the key genes EDAR, DKK4, BMPR1B, BMP3, WNT3, and FGFR2, among others, were identified.

Construction and Weighted Gene Co-Expression Network Analysis of ceRNA Regulatory Networks
Based on expression trend models (lncRNA−|miRNA+| mRNA−) pertaining to the expression trends of 31 overlapping target genes as well as the lncRNAs and miRNAs regulating the expression of these genes, 18 ceRNA regulatory networks were constructed. These networks contained 18 overlapped DE genes, five DE lncRNAs, and three DE miRNAs (Figure 2A). Function and pathway enrichment analysis of GO ( Figure 2B) showed that EDAR, FGFR2, BMPR1b, and VANGL1 were enriched in classical hair follicle pathways, such as Wnt, EDAR/NF-κB, TGF-β, BMP, and FGF signaling pathways, and that the GO terms were mainly related to skin morphology and formation, hair and skin development function and fibroblast cell line differentiation.
To further screen ceRNAs for skin-specific expression at E41, WGCNA was performed on normal and hairless pig skin at five stages, namely E39, E41, E45, E52, and E60, using mRNA sequencing. The independence degree approximated 0.8 while the average connectivity degree was higher ( Supplementary Figures S4A,C). In total, 24 and 28 distinct gene co-expression modules were constructed for the normal group and the hairless group, respectively (Supplementary Figures S4B,D). Modules showing p < 0.05 at E41 were selected, and interestingly, four modules, MEbrown (p = 0.03), MEgreen (P = 4E-04), MEpink (p = 0.03) and MEgray (p = 0.0083), were identified in normal pigs ( Figure 2C), while no significant module was found in hairless pigs (Supplementary Figure S6). Among the hub genes of the four significant modules, only EDAR ( Figure 2D) and BMPR1b ( Figure 2E), of the 18 DE target genes in ceRNA, were found to be involved. Additionally, both EDAR and BMPR1b were interacted with many hub genes (such as KRT5, WNT8A, and KRT4) which were related to hair follicle development, in the most significant MEgreen module ( Figure 2D). Moreover, the expression levels of EDAR ( Figures 1A, 2F) and BMPR1b ( Figures 1A, 2G) in skin at different periods, were found to be not only tissue-specific but also period-specific at E41 (p < 0.001). Therefore, two ceRNAs, lncRNA2162.1/miR-29a-5p/BMPR1b and lncRNA627.1/miR-29a-5p/EDAR, were retained. The results of ssc-miR-29a-5p sequence alignment in different species revealed that ssc-miR-29a-5p is highly conserved among species ( Figure 2H), suggesting the importance of the role played by ssc-miR-29a-5p. A previous study of ours indicated that E41 constituted the critical stage of hair placode formation, and the results of our current study indicated that EDAR and BMPR1b may play an important role in hair placode formation via the ceRNA network.

EDAR Was a Target Gene of miR-29a-5p in ceRNA
Among the ceRNAs, only miR-29a-5p was able to regulate EDAR and BMPR1b. In order to validate direct regulation of the expression levels of EDAR and BMPR1b by miR-29a-5p, a construct containing the 5′UTR of the potential target genes, or the sequence with the mutant seed region, was co-transfected together with miR-29a-5p mimics in human 293T cells. As predicted, miR-29a-5p mimics significantly reduced EDAR-WT-5′UTR (p < 0.01) luciferase activity, whereas no Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 902026 6 significant inhibition of EDAR-MUT-5′UTR was detected ( Figure 3A). Furthermore, the mRNA and protein expression levels of EDAR were significantly decreased (p < 0.05) by miR-29a-5p mimics, whereas these levels were significantly increased (p < 0.01) by a miR-29a-5p inhibitor, compared with the effect exerted by miR-29a-5p on EDAR in the control group (NC) (Figures 3A,B). However, in regard to BMPR1b, no luciferase activity was detected in either BMPR1b WT or MUT-5′UTR, suggesting that miR-29a-5p was unable to bind BMPR1b 5′UTR ( Figure 3D). In addition, no significant change was detected in the mRNA ( Figure 3D) and protein ( Figure 3E) expression of BMPR1b treated with miR-29a-5p mimics and an inhibitor, suggesting the absence of a binding site between miR-29a-5p and BMPR1b. Our results confirmed that EDAR is a direct and functionally relevant target of miR-29a-5p, the expression trend of which is illustrated ( Figure 3C) and that both are involved in hair placode formation, possibly via ceRNA activity.

miRNA-29a-5p Inhibits Hair Placode Formation Through Targeting EDAR via the EDA/EDAR Signaling Pathway
To further verify the function of EDAR in hair placode formation, protein expression of EDAR in E41 skin of normal and hairless embryos was determined. EDAR was highly expressed in the hair placodes of normal embryos which consistent with previously reported (Schmidt-Ullrich et al., 2006;Zhang et al., 2009;Wu et al., 2020) while weakly expression in hairless embryos skin, suggesting that EDAR plays a key role in hair placode formation at E41 ( Figure 4A). The competitive mechanism involving lncRNA627.1-miR-29a-5p-EDAR is illustrated ( Figure 4B), lncRNA627.1 mimics significantly upregulated (p < 0.05) EDAR expression while downregulating miR-29a-5p expression (p < 0.01). Meanwhile, the lncRNA627.1 inhibitor group displayed a converse expression trend, where EDAR was downregulated  and miR-29a-5p was upregulated ( Figure 4B). To gain further insights into how lncRNA627.1 and miR-29a-5p regulate EDAR, related genes, EDA, NF-ĸB, Wnt10b, BMP4, EDARADD, LHX2, and FGF20, associated with the EDA/EDAR signaling pathway, were also detected after transfecting HPPCs with miR-29a-5p mimics and a lncRNA627.1 inhibitor. Among the miR-29a-5p mimics group, BMP4 (p < 0.05) was remarkably upregulated while EDA (p < 0.01), NF-κB (p < 0.01), Wnt10b (p < 0.01) and LHX2 (p < 0.01) were significantly downregulated ( Figure 4C). Similar results were also observed in the lncRNA627.1 inhibitor group ( Figure 4D). Our results demonstrated that the similar expression trend in EDAR related genes by lncRNA627.1 expression inhibition and miR-29a-5p overexpression which suggested that a regulation might existed between lncRNA627.1 and miR-29a-5p and need furthermore investigation to explore whether regulated by ceRNA ( Figures 4C,D). In addition, proliferation of HPPCs further confirmed the regulation by miR-29a-5p and lncRNA627.1. Thirty-six h following the transfection of miR-29a-5p mimics and lncRNA627.1 inhibitor, the number of viable cells in the miR-29a-5p mimics group was significantly decreased (p < 0.01), and continued to decrease further for 48, 60, and 72 h, a result consistent with findings for the lncRNA627.1 inhibitor group ( Figure 4E). We further explored the role played by the association between EDAR and miR-29a-5p in HPPCs proliferation, via siRNA-mediated EDAR silencing as well as via miR-29a-5p mimics mediated suppression of EDAR expression. Edu staining assays indicated that following the reduction of EDAR expression, the proportion of Edu + HPPCs cells were reduced when transfected with EDAR siRNA (p < 0.001) ( Figure 4F) or miR-29a-5p mimics (p < 0.05) ( Figure 4G). Meanwhile, compared with that of mimics/ inhibitor NC groups, EDAR + HPPCs cells were upregulated following transfection miR-29a-5p inhibitor (p < 0.01), while it was downregulated in the miR-29a-5p mimics group (p < 0.001) ( Figure 4H). Considered together, our results indicated that miRNA29a-5p inhibits HPPCs proliferation by suppressing the expression of the targeting gene, EDAR, in the EDA/EDAR signaling pathway, while lncRNA627.1 removes such inhibition. Previous studies of ours (Jiang et al., 2019) indicated that impairment of hair placode formation at E41 constitutes the main cause of hairlessness in hairless pigs. In this study, whole transcriptome sequencing (mRNA, lncRNA, and miRNA) of normal and hairless pig embryos at E41 screened out EDAR and an associated regulatory mechanism involving ceRNAs. Of the known key genes involved in hair placode formation, EDAR appears to play an important, early role by regulating signaling molecules involved in the establishment, formation, and regulation of placodes (Saxena et al., 2019). High expression levels of EDAR increased the number and size of hair placodes (Mustonen et al., 2004;Mohri et al., 2008) whereas hair placode formation could neither be induced nor maintained in the absence of EDAR (Schmidt-Ullrich et al., 2006;Fliniaux et al., 2008). Impaired hair placode formation was associated with a reduction in the epidermal expression levels of EDAR, Lef1, and Shh, which are essential for hair follicle morphogenesis (Mohri et al., 2008). In our study, EDAR was highly expression at the hair placodes of normal skin while no expression in hairless ( Figure 4A) which not only confirmed the key role of EDAR in hair placode formation, but also showed consistency with the results of previous studies (Schmidt-Ullrich et al., 2006;Zhang et al., 2009;Wu et al., 2020). Furthermore, this study also demonstrated that high EDAR expression promotes the proliferation of HPPCs, while low EDAR expression inhibits such proliferation ( Figure 4F). Previous studies had also reported that EDAR exerts its effect on epithelial (Jaskoll et al., 2003) and matrix progenitor cell proliferation (Wang et al., 2017a) via the EDA/EDAR/NF-κB signaling pathway. In hair follicles, Gata6 regulates the rapid proliferation of epithelial (matrix) progenitors during adult mouse hair follicle regeneration, by stimulating the activation of EDAR and NF-κB pathways (Wang et al., 2017a). Furthermore, EDAR and NF-κB play a role in keratinocyte proliferation and hair placode down-growth (Schmidt-Ullrich et al., 2006).
Besides that, it was necessary to use HPPCs in this study, due to the absence of a specific cell line that could be utilized for hair placode formation research. Although, it is advantageous to use tissue block culture and organoid culture methods to detect hair follicle related phenotype at the cellular level while still exist many problems. In tissue block culture for hair placode, the culture Proliferation of fibroblasts at 0, 24, 36, 48, 60, and 72 h after transfection of miR-29a-5p mimics and a lncRNA627.1 inhibitor, as well as miR-29a-5p mimics NC and lncRNA627.1 inhibitor NC. (F) Proliferating fibroblasts labeled with Edu after transfecting EDAR siRNA; (G) Proliferating fibroblasts labeled with Edu following transfection with miR-29a-5p mimics. The click-iT reaction revealed Edu staining (red). Cell nuclei were stained with DAPI (blue). (H) Protein expression of EDAR in fibroblasts following transfection of miR-29a-5p mimics, a miR-29a-5p inhibitor, miR-29a-5p mimics NC and lncRNA627.1 inhibitor NC, cell nuclei were stained with DAPI (blue) while EDAR proteins were stained green. **p < 0.05, **p < 0.01, ***p < 0.001.
Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 902026 environment and culture additive such as activator, enzymes, small molecules, growth factors are still unknown in pig [Zhou et al. (2018), Yan et al. (2019b)] while the combination selection of multiple cell lines and precises space culture environment limit the development of organoid culture (Lee et al., 2018;Lee et al., 2020;Lee and Koehler, 2021). HPPCs at the E41 stage were used in this study because another research of us (unpublished) found that the formation of hair placode mainly comes from a kind of hair placode precursor cell (HPPCs) in epidermal (Supplementary Figure S7A), which not only specifically expresses EDAR, but also highly expresses fibroblast cells related marker genes (COL1A1, LUM, Vimentin), just similar FIGURE 5 | The molecular mechanism of lncRNA627.1/miR-29a-5p/EDAR involved in hair placode formation and expression of EDAR related genes detected in the EDA/EDAR signaling pathway. (A) miR-29a-5p inhibited cell proliferation to suppress the formation of hair placodes by targeting EDAR; (B) Schematic illustrating the association between downstream/upstream target gene expression and EDAR in the EDA/EDAR signaling pathway; (C) The expression of EDAR related genes, Wnt5a, CTNNB1, EDA, EDARADD, NFKB1, CCND1, FGF20, Wnt10b, DKK4, BMP2, and BMP4, in normal and hairless pig skins, at stages E39, E41, E45, E52, and E60.
Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 902026 10 to a pseudo fibroblast cell lineage (Supplementary Figure S7B). We speculated that this pseudo fibroblasts (HPPCs) might be the key cells leading to the formation of hair placode which could also screened from epidermal tissue and highly expressed FB marker, Vimentin and LUM (Supplementary Figure S7C). Thus, these results were the main basis for us to select HPPCs to research in this study. In addition, according to Chen et al., fibroblast proliferation plays a functional role in hair placode initiation  while proliferation capacity was the determining factor in the initiation of placode formation (Duverger and Morasso, 2009;Schneider et al., 2009;Rishikaysh et al., 2014). Our study provided a speculation which EDAR could directly promote hair placode development via a pseudo fibroblasts lineage, HPPCs cells proliferation.
Increasingly, evidence indicates that ncRNA participates in the regulation of hair follicle morphogenesis and development via the ceRNA network which imposes mutual regulation among miRNAs, lncRNAs, and mRNAs (Cai et al., 2019;Zhao et al., 2020). According to ceRNA theory, miRNA suppresses gene expression by recognizing specific target mRNAs, while lncRNAs, which act as natural miRNA sponges, inhibit miRNA function and modulate the expression of target mRNAs by interacting with miRNA response elements (Salmena et al., 2011). In our study, significantly high expression levels of miR-29a-5p suppressed the expression of EDAR and lncRNA627.1 in hairless pig skin tissue at E41 ( Figure 2B). Meanwhile, overexpression lncRNA627.1 elevated the expression of EDAR while the expression of miR-29a-5p was reduce which suggested that lncRNA627.1 may disrupt the binding between miRNA and mRNA while competitive binding displacement assay was needed in further investigation ( Figure 4B). Another similar regulatory mechanism based on ceRNAs is that of miR-144-5p, which indirectly affects HFSCs by targeting lncRNA310320.3 and lncRNA311077.2 . lncRNA679 promoted Wnt3 expression by targeting miR-221-5p in hair follicle cycle development (Wang et al., 2017b). Generally, ncRNA function involves regulation of the expression of targeted mRNAs (Salmena et al., 2011;Tay et al., 2014); for example, lncRNA5322 serves as a ceRNA of miR-19b-3p to elevate MAPK1 expression, ultimately promoting proliferation of HFSCs and epidermal regeneration in mouse models (Cai et al., 2019). Previous studies have reported that miR-29a plays a positive role in skin fibroblast protection (Zhou et al., 2016), hair follicle regeneration (Zhu et al., 2020), and HFSC proliferation (Guo et al., 2018;Ge et al., 2019) by regulating target genes. The results of the current study further revealed that miR-29a-5p plays a role in placode formation by regulating EDAR, where miR-29a-5p downregulated EDAR expression, thereby inhibiting HPPCs proliferation and hair placode formation. In contrast, lncRNA627.1 upregulated EDAR expression by suppressing miR-29a-5p expression ( Figure 5A). It's worth noting that current data supported a model where ncRNA627.1 acts upstream of miR-29a expression ( Figures 4B-E), but no direct evidence that lncRNA627.1 could release the EDAR expression by competitively binding miR-29a-5p based on ceRNA regulatory mechanism. Thus, the relation between lncRNA627.1 and miR-29a-5p, such as competitive binding displacement assay, needs to be further verification.
EDAR influences hair placode formation via its participation in the EDA/EDAR signaling pathway (Schmidt-Ullrich et al., 2006). Therefore, we constructed an interaction network of EDAR-related genes based on previous reports (Fliniaux et al., 2008;Zhang et al., 2009;Wu et al., 2020); ( Figure 5B), and detected the expressions of these genes at E39, E41, E45, E52 and E60 stages of normal and hairless pig embryo skins ( Figure 5C). The expression trends of EDA, NF-κB, Wnt10b and BMP4 were consistent with their expression in HPPCs, in which EDAR expression was inhibited by miR-29a-5p and lncRNA627.1 (Figures 4C,D). Reportedly, EDA (Zhang et al., 2009), Wnt10b (Andl et al., 2002), and NF-κB (Schmidt-Ullrich et al., 2006) were required for the initiation of hair placodes by the EDA/EDAR signaling pathway, while BMP4 acted as an inhibitor for EDAR (Rishikaysh et al., 2014). On the other hand, the differences between the expression levels of some EDAR related genes ( Figure 5C) at E41, were not significant ( Figure 5B), and this may be attributed to heterogeneity of skin cells (Takahashi et al., 2020) and spatio-temporal gene expression. Our findings not only revealed the importance of EDAR in hair placode formation, but also the role played by multi-gene co-regulation which is exerted via participation of the EDAR signaling pathway in hair placode formation. However, further research on molecular communication between the multiple genes involved in coregulation is needed.

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 below: https://ngdc.cncb.ac. cn/, CRA004672.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee of the China Agricultural University.