绒山羊胎儿期毛囊发生发育关键circRNA的鉴定及功能预测分析

尚方正 马 荣 戎友俊 王 敏 潘剑锋 梁丽丽 张燕军* 李金泉

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

circRNA是一类比较特殊的RNA，它没有游离的5′帽子结构和3′poly(A)结构，并且对核酸酶不敏感[1]。circRNA的作用机制主要包括：调控亲本基因的表达[2-3],与 RNA 结合蛋白相互作用[4],翻译蛋白质[5],作为竞争性内源RNA调控基因的表达[6-8]。绒山羊毛囊发育受多种信号分子的调控，特别是近年来较多的circRNA在绒山羊毛囊发生发育有了部分研究进展。研究发现circRNA-1926 可与miR-148a/b-3p 竞争性结合，促进CDK19基因的表达和绒山羊次级毛囊干细胞的分化[9-10]；通过对成年的辽宁绒山羊和内蒙古绒山羊生长期皮肤样本转录组测序数据分析，在辽宁绒山羊中鉴定出17个上调的circRNAs和15个下调的circRNAs，进一步研究发现ciRNA128、circRNA6854、circRNA4154和circRNA3620在不同毛被辽宁绒山羊中存在显著差异表达，揭示circRNA可能对羊绒细度有调控作用[11]。其他方面研究发现在湖羊皮肤毛囊中共鉴定出5 527个circRNA，其中在卷毛和直毛比对组中有114个差异表达的circRNA。circRNA在绒山羊胎儿期毛囊发生发育中的研究还相对匮乏。

绒山羊毛囊可分为初级毛囊和次级毛囊。初级毛囊与次级毛囊形态发生的启动处于胎儿发育的不同时期，张燕军等[12-13]研究发现初级毛囊的启动早于次级毛囊。在胎儿期45～55 d，皮肤形成完整的表皮结构，毛囊还未发生；在55～65 d，胎儿各部位初级毛囊开始发生；在65 d时，体侧部可观察到明显的初级毛囊毛芽；在65～75 d，在胎儿各部位可观察到次级毛囊原始体，次级毛囊开始发生；在75 d时，体侧部可观察到明显的次级毛囊毛芽。

研究发现，定位于细胞质中的外显子circRNA主要和mRNA竞争miRNA的靶标结合位点，从而调控mRNA的表达[14]；而内含子circRNA主要通过参与其亲本基因的表达，如circ-ankrd52是来自于ANKRD52的第二内含子，ci-ankrd52的敲低会导致ANKRD52 mRNA的表达显著降低，表明ci-ankrd52 参与了其亲本基因的的顺式表达[15]。因此，本研究在课题组前期得到绒山羊胎儿期(45、55、65和75 d)4个时期皮肤组织circRNA表达谱的基础上[16]，筛选鉴定绒山羊胎儿期毛囊发生发育相关的外显子circRNA和内含子circRNA，并初步探究其功能，为后续解析circRNA在绒山羊胎儿期毛囊发生发育中的分子机制提供理论基础。

1 材料与方法

1.1 样品采集

于内蒙古金莱牧业科技有限责任公司选取12只饲喂条件相同的内蒙古绒山羊阿尔巴斯型三岁母羊分为四组，每组3只，进行同期发情处理，记录配种时间。采取胎儿期45、55、65和75 d各3只羊皮肤样品。剖腹产后用焦碳酸二乙酯水清洗胎儿。使用无菌、无酶的一次性手术刀和镊子快速采集体侧部皮肤，放入冷冻储存管中，编号，在液氮中快速冷冻，并储存在-80 ℃冰箱中，用于qRT-PCR 试验。所有的胎儿皮肤标本均按照国际动物研究指导原则采集，并获得内蒙古农业大学科学研究和学术道德委员会和内蒙古农业大学生物医学研究伦理的批准(批准文号[2020]056)。

1.2 绒山羊胎儿期毛囊发生发育相关circRNA的筛选

基于全转录组测序得到的circRNA在绒山羊胎儿期(45、55、65和75 d)4个时期表达谱RNA-Seq数据被提交到SRA数据库，注册号(SRR13306949、SRR13306948、SRR13306947、SRR13306946、SRR13306945、SRR13306944、SRR13306943、SRR13306942、SRR13306941、SRR13306940、SRR13306939、SRR13306938)及各比对组中差异表达的circRNA(log2(Fold Change)绝对值≥1且P<0.05)。因此结合毛囊在不同时期的发生发育研究结果，筛选绒山羊毛囊发生发育关键circRNA。

1.3 总RNA的提取及反转录DNA制备

总RNA的提取按照Trizol Reagent试剂操作说明进行。将提取的总RNA用NanoDrop 2000紫外分光光度计测定总RNA的质量浓度。利用1%的琼脂糖凝胶电泳检测RNA的完整性。将总RNA等量分成两份，一份不做处理(供内参定量使用)，另一份经RNase R(RNA外切酶)37 ℃ 25 min，70 ℃ 10 min处理去除线性RNA。按照Takara的 PrimeScriptRT reagent Kit with gDNA Eraser反转录试剂盒置于PCR仪中37 ℃ 15 min；85 ℃ 5 s；4 ℃∞获得 cDNA模板。反转录后的cDNA放置于-20 ℃冰箱保存备用。

1.4 引物设计

以circRNA的全长序列为模板设计绒山羊毛囊发育相关circRNA的引物，β-actin作为内参基因，相关引物信息详见表1，引物由广州锐博生物有限公司合成。

表1 circRNA qRT-PCR引物序列信息Table 1 circRNA qRT-PCR primer sequence information

1.5 circRNA 的实时荧光定量PCR验证

qRT-PCR反应体系为20 μL，其中TB Green premix EXTapⅡ10 μL，cDNA 2 μL，上下游引物各0.8 μL，ROX 0.4 μL，ddH2O 6 μL。使用LightCycler®96 PCR扩增仪进行扩增，扩增程序为：预变性：95 ℃ 30 s；PCP反应：95 ℃ 10 s，62 ℃ 30 s，72 ℃ 10 s，40个循环；熔解：95 ℃ 10 s，60 ℃ 60 s，97 ℃ 1 s；冷却：37 ℃ 30 s。利用2-ΔΔCt进行计算。

1.6 外显子circRNA调控网络的构建

运用TargetScan和miRanda(分析TargetScan和miRanda值的大小，TargetScan阈值为50，该值越大说明两者存在互作的可能性越大；miRanda阈值为-10，该值越小说明两者存在互作的可能性越大)，初步得到circRNA-miRNA靶向关系，进一步对靶向miRNA的靶基因进行预测，将富集到毛囊发育相关重要通路WNT、TGF-β和Notch等信号通路中的靶基因构成的调控网络作为毛囊发生发育相关circRNA-miRNA-mRNA重要的调控网络。并用Cytoscape将其可视化。

1.7 内含子circRNA宿主基因的功能分析

综合分析22个差异表达内含子circRNA宿主基因在KEGG的富集情况，探究其富集的通路是否有与绒山羊胎儿期毛囊发育相关的通路。

1.8 生物信息学分析及相关数据分析

运用Hisat2和Tophat-fusion软件进行参考基因组比对分析，运用CIRCExplorer和CIRI鉴定识别circRNA，得到circRNA其序列全长及来源基因，进一步在NCBI查找相应circRNA的结构组成(Ensembl)。最后使用Cytoscape将circRNA-miRNA-mRNA调控网络可视化。

运用SPSS 18.0进对相关数据进行Kruskal-Wallis秩和检验及Nemenyi检验，其中P<0.05表示差异显著，P<0.01表示差异极显著，P>0.05表示差异不显著。

2 结果与分析

2.1 绒山羊胎儿期毛囊发生发育相关circRNA的筛选

在绒山羊胎儿期45、55、65和75 d皮肤中分别鉴定出2 910、3 710、3 228和3 059个circRNA。其中在绒山羊胎儿D75组与D45组中上调59个，下调33个circRNA；胎儿D65组与D45组中上调96个，下调63个circRNA；胎儿D55组与D45组中，上调76个，下调42个circRNA。将上述3个比对组中均差异表达的circRNA作为毛囊发生发育相关的circRNA(图1)。共筛选到4个绒山羊胎儿期毛囊发生发育相关的circRNA，包括circRNA2049、circRNA2225、circRNA3411和circRNA5681。

图1 绒山羊胎儿期毛囊发生发育相关circRNA的筛选Fig.1 Screening of circRNAs related to hair follicle development in cashmere goats during fetal period

2.2 绒山羊胎儿期毛囊发生发育相关circRNA的鉴定

circRNA2049长度为536 bp，由山羊ATRNL1基因第2～5外显子的全部序列组成，为外显子circRNA；circRNA2225长度为791 bp，由山羊SERINC5基因第3～8外显子的全部序列组成，为外显子circRNA；circRNA5681序列长度为265 bp，由山羊基因102 184 877第3～5外显子的全部序列组成，为外显子circRNA；circRNA3411长度为333 bp，由山羊102 173 234基因第175～178外显子的全部序列组成，为外显子circRNA(表2)。

表2 绒山羊胎儿期毛囊发生发育相关circRNA的结构组成Table 2 The structural composition of circRNA related to hair follicle development in the fetal period of cashmere goats

表2(续)

2.3 绒山羊胎儿期毛囊发生发育相关circRNA差异表达分析

本研究运用qRT-PCR检测了circRNA2049、circRNA2225、circRNA3411和circRNA5681在绒山羊胎儿期(45、55、65和75 d)皮肤组织中的表达量。Kruskal-Wallis秩和检验表明上述circRNA在4个时期中均差异表达(P<0.05)，Nemenyi组间多重比较发现circRNA2049、circRNA2225和circRNA5681在D65vsD55差异显著(P<0.05)(图2)。表明circRNA2049、circRNA2225和circRNA5681可能调控初级毛囊的发生；circRNA3411在D75vsD45差异显著(P<0.05)，推测circRNA3411可能调控次级毛囊的发生。

*P<0.05表示差异显著。* P<0.05 indicates significant difference.图2 绒山羊胎儿期毛囊发生发育相关circRNA qRT-PCR验证分析Fig.2 Validation analysis of circRNA qRT-PCR related to hair follicle development in the fetal period of cashmere goats

2.4 绒山羊胎儿期毛囊发生发育相关circRNA靶向miRNA的预测

本研究运用TargetScan和miRanda预测了筛选到的4个外显子circRNA靶向的miRNA(TargetScan score percentile>50且miRanda energy<-10)。其中circRNA2049有19个靶向的miRNA，circRNA2225有8个靶向miRNA，circRNA3411有4个靶向miRNA，circRNA5681有1个靶向miRNA(表3)。

表3 TargetScan和miRanda 预测circRNA靶向的miRNATable 3 Prediction of miRNAs targeted circRNA by Targetscan and miRanda

表3(续)

2.5 绒山羊胎儿期毛囊发生发育外显子circRNA调控网路构建

本研究基于上述circRNA靶向miRNA的预测结果，同样利用TargetScan和miRanda预测靶向miRNA的靶基因。并结合mRNA的富集情况，筛选富集在Notch signaling pathway 和TGF-β signaling pathway与毛囊发生发育相关的重要信号通路中的mRNA构建了调控绒山羊毛囊发生发育的调控网络。其中DLL4位于NOTCH信号通路的起始段，同时研究表明DLL4与羊毛囊发育相关。因此circRNA2049-chi-miR-27b-3p-DLL4，circRNA2225-chi-miR-145-5p-DLL4和circRNA9181-chi-miR-145-5p-DLL4成为后续绒山羊胎儿期毛囊发生发育的重要调控网络(图3)。

绿色三角形为circRNA;粉色四边形为miRNA;蓝色圆形为mRNA。The green triangle was circRNA; The pink quadrilateral was miRNA; Blue circle was mRNA.图3 绒山羊毛囊发生发育重要circRNA-miRNA-mRNA调控网络Fig.3 The important circRNA-miRNA-mRNA regulatory network of hair follicle development in cashmere goats

2.6 绒山羊胎儿期毛囊发生发育相关内含子circRNA宿主基因功能分析

本研究为探究内含子circRNA宿主基因的功能，分析了6个对照组中22个差异表达的ciRNA(表4)。其中ciRNA450的宿主基因100 860 901富集在Wnt信号通路中；ciRNA587的宿主基因102 189 934富集在Mapk信号通路中(表5)，而Wnt 信号通路和Mapk信号通路是毛囊发生发育相关的重要通路，推测ciRNA可能通过调控其宿主基因的表达来影响毛囊的发生发育。

表4 绒山羊胎儿期毛囊发生发育差异表达的内含子circRNA统计Table 4 Statistics of intron circRNA differentially expressed during hair follicle development in the fetal period of cashmere goats

表5 绒山羊胎儿期毛囊发生发育内含子circRNA宿主基因富集的KEGG通路Table 5 KEGG pathway enriched in circRNA host genes in the development of fetal hair follicles in cashmere goats

3 讨 论

3.1 绒山羊胎儿期毛囊发生发育相关circRNA结构

毛囊的发生发育受到多种信号通路分子的调控，除了众多编码RNA外，非编码也是调控毛囊发生发育的重要因子，包括miRNA、lncRNA和circRNA，但circRNA在毛囊的发生发育中调控作用的研究还相对匮乏。因此，本研究根据课题组前期研究结果，筛选出4个绒山羊毛囊发生发育相关的circRNA，进一步运用生物信息学分析出这4个长度分别为536、791、265和333 bp，这与Li等[17]在绵羊circRNA测序中鉴定到circRNA长度在200～400 bp的结果基本一致。同时，筛选的4个circRNA均由不同个数的外显子组成，均为外显子circRNA。与Li等[18]在circRNA的鉴定中发现外显子circRNA占80%以上的研究结果相符。qRT-PCR表明筛选的4个circRNA均为差异表达的circRNA，组间差异分析发现部分结果与测序结果不一致，但整体表达趋势一致，推测可能与实验误差有关。

3.2 绒山羊胎儿期毛囊发生发育相关内含子circRNA功能预测

全部由内含子环化而成的circRNA称内含子circRNA，这些ciRNAs通常定位于细胞核中，而它们的亲本基因主要定位在细胞质中。研究表明，环状RNA可以通过与线性转录本竞争性剪接，从而影响宿主基因的表达[4]。Zhang等[19]研究发现内含子环状RNA能顺式调控宿主基因的表达。本研究对前期测序结果中22个差异表达的内含子环状RNA的宿主基因富集情况进行整合，发现部分宿主基因富集在Wnt和Mapk信号通路中，其中100 860 901基因富集于Wnt信号通路中。Wu等[20]筛选了Wnt信号通路中Wnt10A，验证发现对绒山羊胎儿期毛囊发生发育有调控作用，马森等[21]研究发现，Wnt/β-catenin 通路与绒山羊绒毛周期性再生相关，并且能够促进绒毛的快速生长，另外发现Wnt蛋白是胎儿毛囊发育启动的最早信号之一[22]；102 189 934基因是本研究中富集到MAPK通路中的基因。研究显示，Tβ4能够影响β-catenin的表达，通过诱导VEGF及MMP-2的表达，使得MAPK中P38的激活，进而影响了毛囊分布和毛干数量[23]；另外研究发现，调控毛囊周期的一些蛋白受MAPK信号激活，肿瘤坏死因子(TNF)作为MAPK信号中的重要成员之一，对毛囊的细胞凋亡有重要作用[24]。因此，推测内含子circRNA可能通过调控宿主基因的表达来影响毛囊的发生发育。

3.3 绒山羊胎儿期毛囊发生发育相关外显子circRNA功能预测

全部由外显子环化而成的circRNA为外显子circRNA。外显子circRNA主要定位于细胞质中，近年来，有大量的研究发现定位于细胞质中的circRNA可以和mRNA竞争miRNA的靶标结合位点，从而调控mRNA的表达[25]。研究表明Notch 信号和毛囊的形态发生有直接的关系[26]。本研究筛选到了绒山羊胎儿期毛囊发生发育关键基因Notch2和DLL4，而先前研究中也发现Notch2对调节羊绒细度具有潜在的重要意义[27]，DLL4与羊毛囊发育相关[28]。TGF-β信号通路同样是毛囊发育相关的重要通路，本研究中筛选到TGF-β信号通路中的FBLN5对毛发生长有影响[29]，本课题组Han等[30]研究表明chi-miR-199a-5p通过调控TGF-β信号通路中的TGF-β2来调控次级毛囊的发生发育通路，另外研究发现TGF-β1对小鼠毛囊退行期有调控作用，敲除TGF-β1小鼠生长期向退行期的转变明显延迟[31]，TGF-β2可以抑制毛干的延长并诱导毛囊退行期的形态学改变[32]，TGF-β2可能通过调控内源性的细胞凋亡引起毛囊角质形成细胞的凋亡，进而抑制毛发的生长[33]。本研究预测关键circRNA靶向miRNA以及靶向miRNA的靶基因，同时部分靶基因富集在Notch和TGF-β信号通路中。当然，circRNA对Notch和TGF-β信号通路调控的研究也相对较多，但主要集中在人类疾病的研究上。研究发现，下调的环状RNA hsa_circ0067301通过miR-141/Notch信号通路调节子宫内膜异位症上皮-间质转化[34]；CircNFIX可通过Notch信号通路介导miR-34a-5p，调控NOTCH1和Notch信号通路，促进胶质瘤的发生[35]；CircAPLP2通过靶向miR-101-3p激活Notch信号通路调控结直肠癌的增殖和转移[36]。circRNA010567通过靶向TGF-β1抑制miR-141促进心肌纤维化[37]；circRNAcESRP1通过分泌miR-93-5p 抑制TGF-β通路在小细胞肺癌化疗敏感性中起着关键作用，提示cESRP1可能是一个有价值的预后生物标志物和潜在的小细胞肺癌治疗靶点[38]；hsa_circ0001085可能通过hsa-miR-196b-5p间接调节TGF-β信号通路从而对前列腺癌细胞起调节作用[39]。因此，认为外显子circRNA主要通过circRNA-miRNA-mRNA途径调控毛囊的发生发育。

4 结 论

本研究筛选了绒山羊毛囊发生发育关键circRNA并发现在山羊胎儿期毛囊发育各阶段存在差异；外显子circRNA主要通过竞争性结合miRNA以circRNA-miRNA-mRNA调控网络来调控毛囊的发生发育；22个内含子circRNA的宿主基因可能通过Wnt和Mapk等信号通路调控毛囊的发生发育。本研究为解析circRNA在绒山羊毛囊发生发育中的分子机制奠定了基础，也为绒山羊分子辅助育种提供了理论依据。