絨山羊胎兒期毛囊發生發育關鍵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在絨山羊毛囊發生發育中的分子機制奠定了基礎，也為絨山羊分子輔助育種提供了理論依據。