绒山羊胚胎期次级毛囊发生发育相关miRNA的筛选及靶基因验证

尚方正 马 荣 戎友俊 潘剑锋 王 敏 牛舒冉 齐云鹏 李彦伯 李金泉 张燕军

(内蒙古农业大学 动物科学学院/农业农村部肉羊遗传育种重点实验室/内蒙古自治区动物遗传育种与繁殖重点实验室/ 内蒙古自治区山羊遗传育种工程技术研究中心，呼和浩特 010018)

MicroRNAs (miRNAs)是一类短链非编码RNA，是动物、植物和某些病毒中发现的一类主要的非编码小RNA，在mRNA(信使RNA)水平上负调控基因表达[1]，并在细胞生长、分化、发育和凋亡在内的各种细胞活动中发挥着重要的调控作用。目前，已有研究筛选鉴定了马、绵羊、鸭、小鼠和兔子等动物皮肤毛囊发生发育相关的miRNA[2-6]。当然在绒山羊皮肤毛囊上也取得了较多的研究进展，Liu等[7]鉴定了内蒙古绒山羊成年期22个新的microRNA与316个保守的microRNA，推测可能对皮肤和毛囊的生长有重要作用；Yuan等[8]鉴定了陕北白绒山羊399个保守的microRNA，其中326个microRNA在生长期、退行期和休止期都有表达，而在生长期、退行期和休止期分别有3、12和11个microRNA特异性表达；郝玉等[9]开展了绒山羊胚胎期和出生后羔羊皮肤差异表达microRNA筛选；袁超[10]对陕北白绒山羊次级毛囊周期性变化相关的microRNA进行鉴定和差异性分析，继而对表达差异性显著的miR-125b进行了初步的功能研究；张桂山[11]开展了microRNA在辽宁绒山羊毛囊生长周期中的差异表达及其靶基因功能鉴定研究。

绒山羊毛囊分为初级毛囊和次级毛囊，初级毛囊产生粗毛，次级毛囊产生羊绒，研究表明绒山羊初级毛囊与次级毛囊形态发生的启动处于胎儿发育的不同时期，且初级毛囊的启动早于次级毛囊[12]。本课题组研究发现，内蒙古绒山羊(阿尔巴斯型)胚胎期45～55 d，上皮细胞开始聚集形成基板，并在此过程中，上皮细胞诱导成纤维细胞定向聚集；在55～65 d，胎儿头顶、肩部和颈部等部位的成纤维细胞在基板诱导下已经形成真皮聚集；到65 d时，各部位的真皮聚集向下延伸，并且体侧部形成了能够明显观察的初级毛囊原始体；在65～75 d，初级毛囊进入高速分化期，而次级毛囊在胎儿各部位开始发生，一般从初级毛囊附近的表皮部位生长出来，与初级毛囊类似，体侧部的次级毛囊形成也晚于其他部位，直到75 d体侧部才可观察到明显的次级毛囊原始体；95 d初级毛囊中出现毛纤维，且初级毛囊在95～120 d发育速度最快，115 d时部分初级毛囊发育成熟，次级毛囊在90～120 d发育迅速，135 d部分次级毛囊发育基本成熟[13-14]。

miRNA通常与基因的3′UTR进行靶向结合来充当基因表达的转录后负调控元件，从而影响若干个主要的生物学途径[15-16]。miRNA-1-3p与FGF14 mRNA的3′UTR特异性结合调控辽宁绒山羊的周期性生长[17]，miRNA-203显著下调DDOT和NAE1 mRNA和蛋白质的表达进而调节绒山羊周期发育[18]。上述研究主要集中在miRNA在绒山羊毛囊周期性生长发育的调控，在绒山羊胎儿期次级毛囊发生发育中的调控作用还相对匮乏。因此，本研究基于课题组前期对内蒙古绒山羊胎儿期皮肤(45、55、65、75、95、115和135 d)进行高通量测序，得到miRNA表达谱的基础上，运用生物信息学分析得到差异表达的miRNA[19]；结合绒山羊毛囊形态发生发育规律，筛选绒山羊次级毛囊发生发育相关miRNA,并利用生物学软件预测其靶基因，结合前期次级毛囊发生发育相关基因的筛选结果，构建次级毛囊发生发育关键miRNA-mRNA调控网络，最后运用双荧光素酶报告系统验证其靶向关系，为解析miRNA在绒山羊次级毛囊中发生发育中的分子机制奠定基础，为绒山羊的良种选育提供理论依据。

1 材料与方法

1.1 样品采集

课题组前期于内蒙古金莱牧业科技有限责任公司选取21只饲喂条件相同的内蒙古绒山羊(阿尔巴斯型)3岁母羊分为7组，每组3只，进行同期发情处理，记录配种时间。采取胎儿期45、55、65、75、95、115和135 d各3只羊皮肤样品。剖腹产后用焦碳酸二乙酯水清洗胎儿。使用无菌、无酶的一次性手术刀和镊子快速采集体侧皮肤，放入冷冻储存管中，编号，在液氮中快速冷冻，并储存在-80 ℃冰箱中，用于转录组测序。

1.2 miRNA的鉴定

共鉴定了三类miRNA：miRbase中已收录的该物种miRNA为已存在miRNA；与miRbase中已知miRNA进行比对鉴定得到已知miRNA;结合参考序列，进行发卡结构预测鉴定到的miRNA为新miRNA。

1.3 差异表达miRNA分析

根据miRNA的表达情况，对miRNA进行差异分析，差异的筛选标准为表达量变化2倍以上并且P值小于0.05。

1.4 次级毛囊发生发育相关miRNA的筛选

为了筛选绒山羊次级毛囊发生发育相关的miRNA和mRNA，将(55 d vs 45 d、65 d vs 55 d)差异miRNA或mRNA取并集(stage A)并将(75 d vs 65 d、95 d vs 75 d、115 d vs 95 d和135 d vs 115 d)差异miRNA或mRNA取并集(Stage B)然后对Stage A和Stage B做维恩图。将Stage B与Stage A共有的差异miRNA和mRNA筛除，Stage B剩余的差异miRNA和mRNA作为与次级毛囊发生发育相关的miRNA和mRNA。

1.5 次级毛囊发生发育相关miRNA靶基因的预测及筛选

运用RNAhybrid、miRanda和TargetScan软件预测筛选到的miRNA靶基因,将每款软件预测得到的靶基因取交集，得到靶基因预测最终结果。进一步将预测的靶基因与课题组前期筛选的次级毛囊发生发育相关的基因取交集，得到关键靶基因。

1.6 关键靶基因的富集分析

对上述筛选到的关键基因进行GO功能分析和KEGG Pathway分析，一方面对它们进行功能注释，另一方面为了能够找出这些候选基因主要相关的功能。将相关基因映射到GO数据库(http:∥www.geneontology.org/)的各term，并计算每个term中的基因数，从而得到具有某个GO功能的基因列表及基因数目统计。GO总共有3个ontology(本体)，分别描述基因的分子功能(Molecular function)、细胞组分(Cellular component)、参与的生物过程(Biological process)。同时基于KEGG(http:∥www.genome.j p/kegg/)的分析有助于更进一步了解基因的生物学功能。筛选富集到毛囊发生发育的关键基因，构建次级毛囊发生发育相关的miRNA-mRNA调控网络。

1.7 双荧光素酶报告基因检测

由中国上海汉生物科技有限公司合成了Chi-miR-26b-5p mimic和Chi-miR-495-3p mimic。psiCHECK2-SMAD1-WT结构是通过将含有miRNA结合序列的SMAD1片段插入到Renilla荧光素酶基因3′末端的psiCHECK-2载体(Promega)中而产生的，psiCHECK2-SMAD1-MUT结构是通过将miRNA结合序列突变为互补序列产生的。对于SMAD1荧光素酶检测，将miRNA mimic和psiCHECK2-SMAD1-WT或突变的psiCHECK2-SMAD1-Mut报告质粒转染HEK 293T细胞。在转染后48 h，使用双荧光素酶报告检测系统(Promega)测量荧光素酶活性。通过比较萤火虫/肾虫荧光素酶的活性比，计算荧光素酶的相对活性。

2 结果与分析

2.1 不同时期皮肤毛囊miRNA长度分布

测序样品miRNA长度分析表明7个时期检测样品中小RNA长度主要集中在21～23 bp(图1)。该长度与之前研究中报道的miRNA长度一致，表明测序结果可靠，可用于后续分析。

图1 不同时期miRNA长度分布图Fig.1 miRNA length distribution in different periods

2.2 不同时期皮肤毛囊miRNA的鉴定

通过与miRbase数据库和参考序列比对，在每个时期的样本中得到约390个已存在miRNA，530个已知miRNA和90个新发现miRNA(表1)。已知miRNA(数量同mirbase中已知的植物或动物miRNA进行比对，鉴定到的miRNA)多于已存在miRNA(mirbase中已收录的该物种miRNA)和新发现miRNA(结合参考序列，进行发卡结构预测，鉴定到的miRNA)，符合miRNA的鉴定规律。后续研究中主要集中在已存在的miRNA。

2.3 不同时期皮肤毛囊miRNA中差异表达的miRNA

筛选了绒山羊胎儿期毛囊发生发育的7个阶段形成的6个比对组中差异表达的已存在的miRNA。每个比对组自左向右依次为表达量差异倍数大于2且(P<0.01,P<0.05和P<0.1)3个标准下上调和下调的miRNA数目(图2)。结果发现(55 d vs 45 d，65 d vs 55 d和75 d vs 65 d)比对组中差异的miRNA多于(95 d vs 75 d，115 d vs 95 d和135 d vs 115 d)比对组。而在毛囊形态研究中发现45～75 d为毛囊的发生阶段，95～135 d为毛囊的发育过程，推测绒山羊毛囊的发生过程受到更多miRNA的调控。

2.4 不同时期皮肤毛囊miRNA中差异表达miRNA靶基因富集分析

为更好地探究皮肤毛囊形态发生的生物学过程，将各比对组中差异表达miRNA靶基因进行GO和KEGG富集分析。结果显示，每个组中富集前20的信号通路均有Wnt和Mapk信号通路(图3)，而先前众多研究表明Wnt和Mapk信号通路是绒山羊毛囊发生发育的重要通路，表明筛选到的miRNA可能是毛囊发生发育的关键miRNA。

表1 各样本中鉴定到的各类miRNA数量统计Table 1 Statistics on the number of various miRNAs identified in each sample

图2 6个比对组中差异表达的miRNA数目统计图Fig.2 Statistical diagram of the number of miRNAs differentially expressed in six comparison group

图3 差异miRNA靶基因 KEGG 富集图Fig.3 KEGG enrichment map of differential miRNA target gene

2.5 绒山羊次级毛囊发生发育相关miRNA的筛选

为探究次级毛囊发生发育相关miRNA,根据毛囊形态发生规律对差异表达的miRNA进行进一步筛选。结果发现StageA(55 d vs 45 d、65 d vs 55 d)包含166个差异miRNA， Stage B(75 d vs 65 d、95 d vs 75 d、115 d vs 95 d和135 d vs 115 d)包含265个差异miRNA,进而取交集得到153个次级毛囊发生发育关键miRNA(图4)。

图4 绒山羊次级毛囊发生发育相关miRNA upset图Fig.4 MiRNA upset map related to secondary hair follicle development in Cashmere Goats

2.6 绒山羊次级毛囊发生发育相关miRNA靶基因预测

为构建绒山羊次级毛囊发生发育相关miRNA-mRNA调控网络，运用RNAhybrid、miRanda和TargetScan三款软件预测筛选上述筛选到的153个miRNA的靶基因。结果表明，153个miRNA共预测到32 978个靶基因，与课题组先前研究中筛选到的次级毛囊发生发育相关的1 548个基因取交集，得到1 352个关键靶基因(图5)。

2.7 绒山羊次级毛囊发生发育相关miRNA关键靶基因富集分析

为进一步探究关键靶基因作用途径，将1 352个关键靶基因进行GO和KEGG富集分析。结果显示在426条 GO term中显著富集，包括301条生物学过程(Biological process)，46条细胞组分(Cellcular component)和79条分子功能(Molecular function)(图6)；同时显著富集在29条KEGG通路中，其中包括毛囊发生发育重要通路TGF-β信号通路(图7)，并有22个基因富集在此通路中(表2)。

2.8 绒山羊次级毛囊发生发育相关miRNA-mRNA调控网络的构建

基于上述富集分析结果，预测富集在毛囊发生发育相关的TGF-B信号通路中的21个基因靶向的miRNA,构建了232对次级毛囊发生发育可能相关的miRNA-mRNA靶向关系对(图8)，为进一步探究绒山羊胚胎期次级毛囊发生发育相关miRNA和mRNA的分子机制奠定了理论基础。

2.9 关键mRNA-miRNA序列分析及双荧光素酶报告基因载体的构建

通过分析miRNA和SMAD1-3′UTR的序列信息，预测到了chi-miR-26b-5p和chi-miR-495-3p与SMAD1-3′UTR的结合位点，这表明SMAD1是这2个miRNA的潜在靶基因。再将合成好的SMAD1-3′UTR序列克隆到pSI-checK2双荧光素酶报告基因载体中的特定位置(图9(a))，然后再对miRNA与SMAD1-3′UTR的结合位点利用PCR突变的方法进行突变(图9(b)和图9(c))，构建突变载体。

图5 次级毛囊发生发育相关miRNA关键靶基因upset图Fig.5 Upset map of key target genes of miRNA related to secondary hair follicle development

图6 关键靶基因GO富集柱状图Fig.6 Histogram of GO enrichment of key target genes

图7 关键靶基因KEGG富集图Fig.7 Enrichment map of key target gene KEGG

表2 显著富集的前11条KEGG信号通路及相关基因Table 2 The significantly enriched first 11 KEGG signaling pathways and related genes

表2(续)

2.10 双荧光素酶报告基因检测chi-miR-26b-5p和chi-miR-495-3p与SMAD1的靶向关系

在chi-miR-26b-5p试验组中(图10(a))，结果显示chi-miR-26b-5p未能下调SMAD1-3′UTR的Luciferase的表达，说明本研究中两者间不存在结合作用。在chi-miR-495-3p试验组中(图10(b))，实验结果显示chi-miR-495-3p未能下调SMAD1-3′UTR的Luciferase的表达，说明本研究中两者间不存在结合作用。

3 讨 论

近年来，miRNA在绒山羊皮肤毛囊中调控作用受到了研究人员的广泛关注，但主要集中在毛囊的周期性生长[18]，而关于miRNA在绒山羊胚胎次级毛囊发生发育中的研究极其匮乏。本研究基于课题组前期建立的绒山羊胚胎期(45、55、65、75、95、115和135 d)皮肤组织的miRNA表达谱的基础上，得到390个已存在miRNA，530个已知miRNA和90个新发现miRNA；6个比对组中共鉴定出525个差异表达的miRNA;根据绒山羊毛囊形态在不同时期的发生发育状况，筛选到次级毛囊发生发育相关miRNA 153个。前期通过对 mRNA和miRNA联合分析，筛选出与内蒙古绒山羊毛囊发生发育相关的 miRNA，并进行了功能验证[20-21]。同时，研究人员在绵羊胎儿发育的11个时期共检测到了179个新的miRNA和152个已知的miRNA，其中37个差异表达miRNAs与次级毛囊形态发生相关[22]。表明miRNA可能参与调控毛囊的发生发育。

图8 绒山羊次级毛囊发生发育相关miRNA-mRNA调控网络图Fig.8 miRNA-mRNA regulatory network related to secondary hair follicle development in cashmere goats

(a)pSI-check2载体结构图谱;(b)chi-miR-26b-5p与SMAD1-3′UTR结合位点及突变位点；(c)chi-miR-495-3p 与SMAD1-3′UTR结合位点及突变位点(a) Structural map of pSI-check2 vector； (b) Binding site and mutation site of chi-miR-26b-5p and SMAD1-3′UTR; (c) Binding site of chi-miR-495-3p and SMAD1-3′UTR and mutation sites图9 双荧光素酶报告基因试验载体结构、结合位点及突变位点Fig.9 The structure, binding site and mutation site of the dual-luciferase reporter gene test vector

进一步预测次级毛囊发生发育相关153个miRNA的靶基因，得到32 978个靶基因，取交集后得到1 352个靶向的次级毛囊发生发育相关的基因，且显著富集在29条KEGG通路中，其中21个基因富集在毛囊发生发育相关的重要TGF-β通路中。众多研究表明毛囊的发生发育可能与一些蛋白编码基因相关，目前认为调节毛囊形态发生的信号分子大部分属于NOTCH传导通路[23]，转化生长因子(TGF)家族[24],Wnt传导通路[25],成纤维细胞生长因子(FGF)家族[26]等。其中TGF-β信号通路是一类广泛的存在于各种生物体内的多功能细胞因子，属于转化生长因子超家族。TGF-β共有5种异构体，其结构相关，序列相似。在哺乳动物中只发现3种亚型(TGF-β1，TGF-β2和TGF-β3)[27]，当然在TGF-β信号通路包括许多其他基因，如本研究中富集的TGF-β2和SMAD1等基因。研究发现TGF-β2可以抑制毛干的延长并诱导毛囊退行期的形态学改变；同时发现TGF-β2可能通过调控内源性的细胞凋亡引起毛囊角质形成细胞的凋亡，进而抑制毛发的生长[28-29]；SMAD1在鹅胚胎期毛囊发育分化过程中具有潜在的调控作用.因此，接下来的试验中着重探究了SMAD1的在绒山羊次级毛囊发生发育中的调控作用。

动物体内miRNA通过与mRNA的3′端的非翻译区(Untranslated regions，UTR)相结合，在不减少mRNA数量的情况下，达到抑制mRNA翻译蛋白的功能。具体作用机制为成熟miRNA的种子区域(seed，miRNA 5′端到3′端的2-7位核苷酸)与AGO蛋白形成复合物后结合mRNA3′UTR区域，从而抑制mRNA翻译[30-31]。实际上，只有与Argonaute结合的miRNA的种子序列才可以参与配对[32-34]。chi-miR-26b-5p和 chi-miR-495-3p与SMAD1-3′UTR区域均存在着结合位点，通过双荧光素酶报告基因发现，chi-miR-26b-5p和 chi-miR-495-3p均显示与SMAD1-3′UTR没有结合作用，推测chi-miR-26b-5p和 chi-miR-495-3p可能通过其他途径来影响着SMAD1的表达。

4 结 论

本研究筛选了绒山羊次级毛囊发生发育相关的miRNA，这些miRNA可能通过毛囊发生发育相关的TGF-β等重要信号通路中的基因调控绒山羊次级毛囊发生发育；同时chi-miR-26b-5p和 chi-miR-495-3p与SMAD1不存在靶向结合关系。总之，本研究为解析miRNA在绒山羊次级毛囊中的分子调控机制奠定的了重要的理论依据。