孫 爽，阿布來(lái)提·蘇來(lái)曼，李 彬，李志強(qiáng)，侯晨曦，決肯·阿尼瓦什，依明·蘇來(lái)曼*

(1.新疆農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)學(xué)院，烏魯木齊830052；2.新疆畜牧科學(xué)院畜牧研究所，烏魯木齊830011)

和田羊?qū)儆诙讨伯愘|(zhì)半粗毛羊，具有耐干旱炎熱、耐粗飼、抗病性強(qiáng)等優(yōu)點(diǎn)[1，2]，其毛色潔白，富有光澤，紋長(zhǎng)，纖維彈性強(qiáng)，是織制地毯和提花毛毯的優(yōu)質(zhì)原料[3，4]?，F(xiàn)有的和田羊經(jīng)過(guò)多年的保種選育仍存在體型較小、產(chǎn)毛量低、底絨細(xì)少、干死毛多等缺陷[5]。研究表明，羊毛的品質(zhì)與產(chǎn)量主要由毛囊決定[6]，而后者存在周期性生長(zhǎng)的特點(diǎn)，一個(gè)生長(zhǎng)周期包括生長(zhǎng)期、退行期和休止期[7]。對(duì)于牛、羊哺乳動(dòng)物來(lái)說(shuō)，生長(zhǎng)期一般在每年的4-11月；退行期一般是在每年的12月到次年的1月[8]，休止期一般在每年的2-3月，哺乳動(dòng)物皮膚毛囊的生長(zhǎng)周期主要是由外界溫度變化誘導(dǎo)，一般生長(zhǎng)周期與每年的季節(jié)交替保持一致。王麗[9]研究發(fā)現(xiàn)和田羊毛囊形態(tài)結(jié)構(gòu)在2月以后迅速變化，2-8月，和田羊總毛囊密度整體呈增長(zhǎng)趨勢(shì)，處于生長(zhǎng)期，2月時(shí)毛囊活性較弱，4-8月毛囊群豐富，生長(zhǎng)密集，活性較好，為和田羊毛纖維的生長(zhǎng)發(fā)育提供良好的生理基礎(chǔ)。山區(qū)型和田羊毛囊密度在6月到8月增長(zhǎng)趨勢(shì)尤為顯著，因山區(qū)型和田羊生活在海拔1 900 m左右的山腰上，氣溫較低且溫差大，結(jié)合本課題組的長(zhǎng)期觀察研究，山區(qū)型和田羊的毛囊生長(zhǎng)期一般在5月，退行期一般在10月，休止期一般在1月。

毛囊周期性生長(zhǎng)發(fā)育與許多代謝過(guò)程有關(guān)，毛囊在特定的細(xì)胞中，可以選擇性的表達(dá)或抑制某些基因，在這過(guò)程中有大量調(diào)控因子參與其中，如miRNA、LncRNA等，涉及多種復(fù)雜的調(diào)控通路，如Wnt、Notch、TGF-β、Hedgehog、MAPK、TNF等通路，還有大量細(xì)胞參與其中，相互協(xié)調(diào)，共同完成了毛囊生長(zhǎng)發(fā)育過(guò)程[10]。miRNAs是一類機(jī)體自身基因轉(zhuǎn)錄產(chǎn)生的、長(zhǎng)度約為18～25個(gè)核苷酸的內(nèi)源性小分子非編碼RNA，能通過(guò)反向互補(bǔ)配對(duì)的方式與mRNA的3'UTR靶向結(jié)合，進(jìn)而在轉(zhuǎn)錄后水平影響mRNA的穩(wěn)定性或蛋白翻譯過(guò)程[11]。研究表明，miRNAs在毛囊發(fā)育過(guò)程中發(fā)揮著重要的調(diào)控作用。Zhang[12]等利用芯片和生物信息學(xué)技術(shù)比較內(nèi)蒙古絨山羊和美利奴綿羊的耳尖和體側(cè)區(qū)皮膚中的差異表達(dá)miRNAs，發(fā)現(xiàn)有19個(gè)miRNAs在體側(cè)皮膚區(qū)特異表達(dá)，15個(gè)miRNAs可能在毛囊發(fā)育中上調(diào)。曲海娥[13]發(fā)現(xiàn)miRNA-125b在細(xì)毛羊和絨山羊中作用于FGFR2和MXD4這兩個(gè)靶基因然后調(diào)控毛囊的生長(zhǎng)發(fā)育。李勤群[14]發(fā)現(xiàn)miRNA-148b和miRNA-320在藏綿羊毛囊發(fā)育生長(zhǎng)期的表達(dá)量顯著增加。唐曉慧[15]通過(guò)對(duì)毛囊周期性變化中3個(gè)時(shí)期轉(zhuǎn)錄組測(cè)序的分析，在毛囊生長(zhǎng)期、退行期、休止期發(fā)現(xiàn)差異表達(dá)miRNAs分別有10、27、7個(gè)，并同時(shí)發(fā)現(xiàn)miRNA-184、miRNA-205的靶基因富集在Wnt通路中。

研究和田羊毛囊周期性的差異表達(dá)miRNA，發(fā)掘綿羊miRNA圖譜，探索miRNA在周期表達(dá)的變化，再通過(guò)靶基因預(yù)測(cè)及功能富集分析，為進(jìn)一步研究其在綿羊周期調(diào)控過(guò)程中的功能提供幫助，為篩選和田羊毛性狀關(guān)鍵基因提供前期基礎(chǔ)，最終為在實(shí)際育種生產(chǎn)中通過(guò)調(diào)控miRNA及其靶基因的表達(dá)，人工干預(yù)綿羊羊毛周期生長(zhǎng)發(fā)育實(shí)現(xiàn)增產(chǎn)提供科學(xué)依據(jù)。目前和田羊毛囊周期性變化中差異表達(dá)miRNAs的篩選及相關(guān)分析對(duì)和田羊毛生長(zhǎng)發(fā)育和高品質(zhì)地毯毛形成的分子機(jī)制方面的研究還未見(jiàn)報(bào)道。

本研究利用RNA測(cè)序(RNA-sequencing)技術(shù)鑒定了和田羊毛囊周期性變化過(guò)程中不同時(shí)期毛囊的差異表達(dá)miRNAs，并通過(guò)生物信息學(xué)分析鑒定了差異表達(dá)miRNAs的潛在生物學(xué)功能，結(jié)果表明miRNAs可能參與調(diào)控和田羊毛囊周期性變化，為進(jìn)一步探究miRNAs對(duì)和田羊毛囊周期性變化中的調(diào)控機(jī)制提供理論基礎(chǔ)，同時(shí)為和田羊品種的選育提供了一定的科學(xué)依據(jù)。

1 材料與方法

1.1 試驗(yàn)材料

1.1.1 試驗(yàn)動(dòng)物

在新疆和田地區(qū)于田縣昆侖種羊場(chǎng)選擇3只周歲、健康無(wú)病、體況均一的山區(qū)型和田母羊作為試驗(yàn)對(duì)象。經(jīng)過(guò)統(tǒng)一飼養(yǎng)管理，分別于毛囊生長(zhǎng)的生長(zhǎng)期(5月份)、退行期(10月份)和休止期(次年1月份)采集數(shù)塊同一側(cè)肩胛部位大小約1 cm2的皮膚組織，分三份置于液氮中保存?zhèn)溆谩?/p>

1.1.2 試劑及儀器

Infinite M200酶標(biāo)儀，DYY-6C型電泳儀，7900熒光定量PCR儀，Gel Doc 2000型凝膠成像儀，反轉(zhuǎn)錄試劑盒(諾唯贊公司)，Trizol試劑(Invitrogen公司)，熒光定量PCR試劑盒(Tiangen公司)。

1.2 試驗(yàn)方法

1.2.1 提取皮膚組織RNA將三個(gè)時(shí)期的皮膚組織樣分別在加入液氮的滅菌研缽中研磨至粉末狀，用Trizol法提取RNA，檢測(cè)提取RNA的濃度與質(zhì)量。合格的RNA置于-80℃保存，備用。本試驗(yàn)所有操作均在低溫下進(jìn)行.

1.2.2 cDNA的制備和文庫(kù)的構(gòu)建

cDNA制備：于滅菌的eppendorf管內(nèi)依次加入2X miRNA RT Reaction、10 μL Buffer、2 μL miRNA RT Mix、500 ng RNA，最后加RNase-free water配成20 μL的體系，按42℃、60 min；95℃、3 min程序進(jìn)行反轉(zhuǎn)錄，產(chǎn)物置于冰箱-20℃保存。

文庫(kù)的構(gòu)建：由上海鯨舟基因科技有限公司完成。

1.2.3 轉(zhuǎn)錄組的高通量測(cè)序及分析

(1)由上海鯨舟基因科技有限公司采用illumina Xten測(cè)序平臺(tái)的雙端150 bp測(cè)序模式對(duì)樣本進(jìn)行高通量測(cè)序，共構(gòu)建3個(gè)時(shí)期9個(gè)轉(zhuǎn)錄組測(cè)序文庫(kù)。

(2)對(duì)測(cè)序原始數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估，并利用Seqtk軟件過(guò)濾原始數(shù)據(jù)(raw reads)得到高質(zhì)量數(shù)據(jù)(clean reads)。

(3)將clean reads分別比對(duì)到miRBase、Genome數(shù)據(jù)庫(kù)羊的基因組(Oar-v3.1.91)上，確定已知的miRNAs，將未知的smallRNA與Rfam非編碼數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，對(duì)小RNA片段進(jìn)行分類注釋和數(shù)量統(tǒng)計(jì)，再利用miRDeep軟件對(duì)未知的小分子rRNA根據(jù)miRNAs前體的標(biāo)志性發(fā)夾結(jié)構(gòu)來(lái)預(yù)測(cè)新的miRNAs。

(4)以下將以P1、P2、P3，Ⅰ、Ⅱ、Ⅲ分別代表生長(zhǎng)期、退行期、休止期，P1、P2、P3為同一時(shí)期3(如01180635Ⅰ、00617149Ⅰ、01180656Ⅰ)個(gè)樣本的平均值，分三組對(duì)比組對(duì)比，P2/P1對(duì)比組，P2為處理組，P1為參考組；P3/P1對(duì)比組，P3為處理組，P1為參考組；P3/P2對(duì)比組，P3處理組，P2參考組。將每個(gè)樣本的reads比對(duì)到miRBase上，計(jì)算miRNAs表達(dá)量(count數(shù))。miRNAs表達(dá)量的篩選采用CPM計(jì)算度量指標(biāo)(counts per million)，再利用DESeq軟件對(duì)樣本處理組與樣本參考組進(jìn)行差異表達(dá)分析，在兩兩對(duì)比比較組選擇Pvalue≤0.05并且log2(foldchange)≥1的miRNAs作為差異表達(dá)miRNAs。

1.2.4 差異表達(dá)miRNAs靶基因的預(yù)測(cè)

利用miRanda、miRDB、TargetScan等工具對(duì)差異表達(dá)miRNAs進(jìn)行靶基因預(yù)測(cè)，取預(yù)測(cè)結(jié)果的交集作為預(yù)測(cè)的最終結(jié)果。

1.2.5 靶基因的GO和KEGG富集分析

將預(yù)測(cè)的靶基因分別利用GO seqR和KOBAS進(jìn)行GO和KEGG功能富集來(lái)分析靶基因的功能及信號(hào)通路以研究差異表達(dá)miRNAs潛在的功能。

1.2.6 qRT-PCR

為了驗(yàn)證轉(zhuǎn)錄組測(cè)序的準(zhǔn)確性，利用qRT-PCR技術(shù)檢測(cè)了P3與P1兩個(gè)時(shí)期差異表達(dá)較大、CPM值較高的7個(gè)miRNAs進(jìn)行相對(duì)表達(dá)量分析。(1)根據(jù)7個(gè)miRNAs的成熟序列設(shè)計(jì)qRT-PCR的上游引物，下游引物由試劑盒(miRcute Plus miRNA qPCR Kit(Tiangen FP411)提供，引物序列見(jiàn)表1，由上海生工生物工程技術(shù)服務(wù)有限公司合成。(2)以1.2.2獲得的cDNA為模板，以U6為內(nèi)參基因，以SYBR Green染料法進(jìn)行7個(gè)miRNAs相對(duì)表達(dá)水平的實(shí)時(shí)定量分析。qRT-PCR反應(yīng)總體系為10 μL：Mix(2X)5 μL、Amplification primer 1(2 μM)1 μL、Reverse Primer(2 μM)1 μL、cDNA 1 μL、ddH2O 2 μL，擴(kuò)增程序?yàn)?0℃2min；95℃，10 min；95℃15 s，60℃1 min，40個(gè)循環(huán)；95℃，15 s，60℃，15 s，95℃，15 s，每個(gè)樣品設(shè)置3個(gè)重復(fù)。結(jié)果以2-ΔΔCt方法進(jìn)行目的基因相對(duì)表達(dá)量的計(jì)算，其中ΔCt=Ct目的基因-Ct內(nèi)參基因，以01180635Ⅰ作為對(duì)照組，ΔΔCt=ΔCt實(shí)驗(yàn)組-ΔCt對(duì)照組。與測(cè)序結(jié)果進(jìn)行對(duì)比，用excel作圖表示。

表1 引物序列參數(shù)

2 結(jié)果與分析

2.1 提取RNA檢測(cè)

提取的質(zhì)量檢測(cè)結(jié)果如表2，OD260/OD280在1.8～2.1之間，濃度均在50 ng/μL以上，RIN在7.0以上，表示所提RNA濃度高且無(wú)蛋白質(zhì)、DNA等污染，符合后期反轉(zhuǎn)錄的要求，可用于后續(xù)的實(shí)驗(yàn)。

表2 提取RNA質(zhì)量檢測(cè)結(jié)果

2.2 轉(zhuǎn)錄組測(cè)序及分析

2.2.1 測(cè)序質(zhì)量與數(shù)據(jù)處理情況

平均每個(gè)樣品測(cè)序獲得約24 M的數(shù)據(jù)量，平均每個(gè)樣本獲得24 768 122個(gè)原始數(shù)據(jù)，處理后平均每個(gè)樣本獲得21 499 205個(gè)clean reads，過(guò)濾各項(xiàng)有影響的reads得到的結(jié)果見(jiàn)表3，高質(zhì)量數(shù)據(jù)占比較高為86.80%，30%在所有樣本中均為90%以上，說(shuō)明單個(gè)堿基的錯(cuò)誤率非常低，表明產(chǎn)出的reads質(zhì)量可靠；9個(gè)文庫(kù)所得序列長(zhǎng)度大部分位于19～25 nt之間，其中長(zhǎng)度為22 nt的序列所占比例最大，符合miRNAs長(zhǎng)度分布特征。圖1為樣本S0617149I的reads長(zhǎng)度占比情況。以上均表明測(cè)序數(shù)據(jù)質(zhì)量較好，符合試驗(yàn)要求，可用于后續(xù)結(jié)果分析。

表3 轉(zhuǎn)錄組的數(shù)據(jù)統(tǒng)計(jì)

圖1 S0617149I樣本的reads長(zhǎng)度占比情況

2.2.2 基因組比對(duì)結(jié)果

如表4，平均每個(gè)樣本有8 173 966條reads比對(duì)到miRBase mature數(shù)據(jù)庫(kù)上，平均每個(gè)樣本可比對(duì)reads占比為37.90%；平均每個(gè)樣本有18 544 011條reads比對(duì)到Genome數(shù)據(jù)庫(kù)上，平均每個(gè)樣本可比對(duì)reads占比為86.35%，73.09%在參考序列上只有唯一比對(duì)位置，占比越高說(shuō)明轉(zhuǎn)錄組的測(cè)序數(shù)據(jù)利用率越高，可用于后續(xù)試驗(yàn)；平均每個(gè)樣本有17 160 325條reads比對(duì)到Rfam數(shù)據(jù)庫(kù)，未被注釋的reads有2 116 948條，各類小RNA的分類注釋和數(shù)量統(tǒng)計(jì)結(jié)果如表5，可用于后續(xù)新miRNAs的預(yù)測(cè)。

表4 比對(duì)到各基因組數(shù)據(jù)庫(kù)的情況

表5 不同長(zhǎng)度的數(shù)據(jù)庫(kù)讀取映射

2.2.3 miRNAs鑒定結(jié)果

共鑒定了143個(gè)已知的miRNAs，并發(fā)現(xiàn)這些miRNAs屬于75個(gè)家族，其中oar-let-7家族的最多，有7個(gè)：oar-let-7a、oar-let-7b、oar-let-7c、oar-let-7d、oar-let-7f、oar-let-7g、oar-let-7i，oar-miR-200家族的有3個(gè)：oar-miR-200a、oar-miR-200b、oar-miR-200c；并預(yù)測(cè)了759新miRNAs。

2.2.4 差異表達(dá)miRNAs的篩選

在已知的miRNAs中，相較于P1組，P1與P2兩個(gè)時(shí)期共鑒定了62個(gè)差異表達(dá)miRNAs，全部下調(diào)；而在P1與P3兩個(gè)時(shí)期鑒定了差異表達(dá)miRNAs 71個(gè)，包含P3時(shí)期3個(gè)顯著上調(diào)的miRNAs與68個(gè)顯著下調(diào)的miRNAs；在P2與P3時(shí)期未發(fā)現(xiàn)差異表達(dá)miRNAs，其中oar-miR-370-3p、oar-miR-541-3p、oar-miR-134-3p、miR-125b等miRNA在三個(gè)時(shí)期差異較明顯；在新預(yù)測(cè)的759個(gè)miRNAs中，相較于P1組，P1與P2兩個(gè)時(shí)期共鑒定了111個(gè)顯著差異表達(dá)miRNAs，包含P2期7個(gè)顯著上調(diào)的miRNAs與104個(gè)顯著下調(diào)的miRNAs；而在P1與P3兩個(gè)時(shí)期比較中鑒定了172個(gè)差異表達(dá)miRNAs，包含P3期25個(gè)顯著上調(diào)的miRNAs與147個(gè)顯著下調(diào)的miRNAs。相較于P2組，P3組共鑒定了71個(gè)顯著差異表達(dá)miRNAs，其中包含23個(gè)顯著上調(diào)的miRNAs與48個(gè)顯著下調(diào)的miRNAs。以上結(jié)果表明休止期和生長(zhǎng)期比較組獲得的差異miRNAs明顯多于其他兩個(gè)比較組，說(shuō)明由休止期向生長(zhǎng)期轉(zhuǎn)變過(guò)程中，有更多的miRNAs參與調(diào)控毛囊周期性變化過(guò)程；同時(shí)在三個(gè)比較組中，差異表達(dá)miRNA都是下調(diào)表達(dá)的比上調(diào)表達(dá)的多，說(shuō)明在和田羊毛囊周期性變化過(guò)程中miRNA主要發(fā)揮負(fù)性調(diào)控作用。將上述數(shù)據(jù)繪制成火山圖，如圖2，為新預(yù)測(cè)的miRNAs中P2/P1對(duì)比組差異表達(dá)miRNAs的火山圖。

圖2 新預(yù)測(cè)miRNAs中P2/P1比較組差異表達(dá)miRNAs的火山圖

為了進(jìn)一步鑒定差異表達(dá)miRNA在和田羊3個(gè)不同時(shí)期中的分布情況，對(duì)已知的差異表達(dá)miRNAs進(jìn)行了維恩圖分析(圖3)，在維恩圖分析中可直觀展現(xiàn)出各差異比較組合之間共有與特有的差異表達(dá)miRNAs數(shù)量，在P2/P1時(shí)期檢測(cè)到了6個(gè)時(shí)期特異性miRNAs：oar-miR-329a-5p、oar-miR-411b-3p、oar-miR-433-5p、oar-miR-496-3p、oar-miR-539-5p、oar-miR-665-5p，在P3/P2時(shí)期檢測(cè)到了2個(gè)時(shí)期特異性miRNAs：oar-miR-134-3p、oar-miR-323a-5p。提示這些時(shí)期特異性表達(dá)miRNAs在毛囊周期性變化過(guò)程有重要作用。

圖3 已知的miRNA中差異表達(dá)miRNAs維恩圖分析

2.3 靶基因預(yù)測(cè)

對(duì)上述篩選的和田羊毛囊周期性變化的3個(gè)時(shí)期的所有差異表達(dá)miRNAs分別使用miRanda、miRDB、TargetScan工具進(jìn)行靶基因預(yù)測(cè)，取三者交集為差異表達(dá)miRNAs調(diào)控的靶基因。在生長(zhǎng)期與退行期比較組173個(gè)差異表達(dá)miRNAs預(yù)測(cè)得到1 883個(gè)靶基因；休止期與生長(zhǎng)期比較組243個(gè)差異表達(dá)miRNAs預(yù)測(cè)得到1881個(gè)靶基因；休止期與退行期比較組71個(gè)差異表達(dá)miRNAs預(yù)測(cè)得到140個(gè)靶基因。我們發(fā)現(xiàn)在這些預(yù)測(cè)得到靶基因中存在眾多與毛囊生長(zhǎng)、發(fā)育密切相關(guān)的功能基因，例如差異表達(dá)miRNA oar-miR-134-5p和oar-miR-432共同調(diào)控的KRT15已知在毛囊發(fā)育期、生長(zhǎng)期啟動(dòng)和靜止期的毛囊隆突區(qū)細(xì)胞表達(dá)，但在生長(zhǎng)期和退化期表達(dá)不一致[16]；差異表達(dá)miRNA oar-miR-495-5p調(diào)控的BMP4屬于骨形態(tài)發(fā)生蛋白質(zhì)(BMP)，而B(niǎo)MP在控制皮膚生長(zhǎng)、出生后的組織重塑和腫瘤發(fā)生上起著重要的作用，是毛囊發(fā)育所涉及的重要信號(hào)分子之一[17]；oar-miR-433-3p和oar-miR-432的靶基因KRT25被發(fā)現(xiàn)可能與灘羊皮膚不同生長(zhǎng)階段的毛囊周期性變化有關(guān)[18]。以上預(yù)測(cè)的靶基因均與毛囊發(fā)育密切相關(guān)，提示調(diào)控這些基因的miRNAs可能通過(guò)靶向調(diào)控靶基因參與毛囊周期性變化過(guò)程，而其中確切的靶向關(guān)系也需要進(jìn)一步的試驗(yàn)來(lái)進(jìn)行驗(yàn)證。同時(shí)，我們還發(fā)現(xiàn)miRNAs與靶基因的調(diào)控關(guān)系可能是一一對(duì)應(yīng)的，也可能是多個(gè)miRNAs共同調(diào)控單個(gè)靶基因，還可能是單個(gè)miRNA同時(shí)調(diào)節(jié)多個(gè)靶基因，如圖4，oar-miR-485-3p既能調(diào)控STEAP4、INA、SMC4、FSD2，同時(shí)還和oar-miR-154a-3p共同調(diào)控GHITM。這說(shuō)明和田羊中miRNAs具有多重靶向性。

圖4 差異表達(dá)miRNAs與靶基因網(wǎng)絡(luò)

2.4 靶基因GO和KEGG富集分析

本研究通過(guò)GO富集分析三個(gè)比較組的靶基因被注釋毛囊發(fā)育相關(guān)的激素水平調(diào)節(jié)、細(xì)胞分解代謝、周期蛋白結(jié)合、骨細(xì)胞發(fā)育等GO條目中；通過(guò)和KEGG富集分析發(fā)現(xiàn)靶基因顯著富集在與毛囊發(fā)育相關(guān)的AMPK、Hedgehog、TNF信號(hào)通路中，圖5為退行期和生長(zhǎng)期比較組靶基因KEGG富集結(jié)果中選取的顯著富集的20個(gè)KEGG通路繪制成的柱形圖。大量研究證實(shí)，出生后毛囊周期性變化激素調(diào)控，如雄激素以及雌激素等，并通過(guò)各相關(guān)調(diào)控信號(hào)分子/通路包括Wnt/β-catenin、骨形成蛋白(BMP)、Sonichedgehog(Shh)、notch、Hox、成纖維細(xì)胞生長(zhǎng)因子(FGF)、胰島素樣生長(zhǎng)因子(IGF)、血小板衍生生長(zhǎng)因子(PDGF)、轉(zhuǎn)化生長(zhǎng)因子(TGF)、AMPK參與毛囊周期性發(fā)育過(guò)程。本研究與前人研究相符，提示在和田羊中皮膚毛囊差異表達(dá)miRNAs通過(guò)各信號(hào)通路調(diào)控靶基因參與毛囊周期進(jìn)程。

圖5 基于KEGG數(shù)據(jù)庫(kù)的靶基因相關(guān)通路功能分布圖

2.5 熒光定量PCR

結(jié)果如圖6所示，所有選擇的7個(gè)miRNAs表達(dá)趨勢(shì)與轉(zhuǎn)錄組測(cè)序結(jié)果一致，說(shuō)明轉(zhuǎn)錄組測(cè)序結(jié)果準(zhǔn)確可靠。

圖6 差異表達(dá)miRNAs qRT-PCR與轉(zhuǎn)錄組測(cè)序比較

3 討論

和田羊尤其是山區(qū)型和田羊，其羊毛是其經(jīng)濟(jì)價(jià)值的所在，羊毛的品質(zhì)決定著經(jīng)濟(jì)價(jià)值，羊毛的品質(zhì)受毛囊的影響，而毛囊的發(fā)育和周期性變化受到多種因子的調(diào)控。很久之前就有研究證明miRNAs能夠通過(guò)與mRNA的3'UTR序列配對(duì)的形式與mRNA結(jié)合，來(lái)影響mRNA的翻譯過(guò)程，起到轉(zhuǎn)錄后的調(diào)控作用[19]。近年來(lái)有研究表明，miRNAs在多種哺乳動(dòng)物皮膚毛囊中均有表達(dá)，在毛囊發(fā)育中起著重要作用[20]。為了探討miRNAs在和田羊毛囊周期性變化中皮膚組織的轉(zhuǎn)錄組表達(dá)差異及可能發(fā)揮的作用，本研究利用高通量測(cè)序分析發(fā)現(xiàn)miRNAs在三個(gè)時(shí)期比較組有顯著差異表達(dá)，表示miRNAs可能參與調(diào)控和田羊毛囊的周期性變化進(jìn)程。

Liu[21]利用深度測(cè)序技術(shù)鑒定了絨山羊皮膚中的miRNAs，發(fā)現(xiàn)了22個(gè)全新的山羊miRNAs。張桂山[22]在遼寧絨山羊和乾華肉用美利奴羊皮膚毛囊發(fā)現(xiàn)毛囊發(fā)育的三個(gè)時(shí)期—生長(zhǎng)期、退行期和休止期特有的miRNAs數(shù)分別為21個(gè)、23個(gè)和26個(gè)，15個(gè)、30個(gè)和27個(gè)。Yuan等[23]鑒定了絨山羊皮膚毛囊的399個(gè)保守的miRNAs，其中有36個(gè)miRNAs在毛囊發(fā)育的3個(gè)時(shí)期都有表達(dá)，有23個(gè)在皮膚毛囊發(fā)育的生長(zhǎng)期特異表達(dá)，29個(gè)皮膚毛囊發(fā)育的退行期特異表達(dá)，44個(gè)在皮膚毛囊發(fā)育的休止期特異性表達(dá)。本研究共鑒定了143個(gè)已知miRNAs并預(yù)測(cè)了759個(gè)新miRNAs，在退行期和生長(zhǎng)期對(duì)比組鑒定差異表達(dá)miRNA 173個(gè)，休止期和生長(zhǎng)期對(duì)比組發(fā)現(xiàn)差異表達(dá)miRNA 243個(gè)；休止期和退行期對(duì)比組發(fā)現(xiàn)差異表達(dá)miRNA 71個(gè)；并發(fā)在退行期和生長(zhǎng)期檢測(cè)到了6個(gè)時(shí)期特異性miRNAs：oar-miR-329a-5p、oar-miR-411b-3p、oar-miR-433-5p、oar-miR-496-3p、oar-miR-539-5p、oar-miR-665-5p，在休止期和退行期檢測(cè)到了2個(gè)時(shí)期特異性miRNAs：oar-miR-134-3p、oar-miR-323a-5p；而在三個(gè)時(shí)期顯著差異表達(dá)miRNAs中的oar-let-7家族與山羊皮膚毛囊周期性生長(zhǎng)發(fā)育密切相關(guān)[24]，oar-miR-200家族在羊駝皮膚和毛囊發(fā)育中也可能起重要作用[25]，并且發(fā)現(xiàn)顯著差異表達(dá)miR-125b，可以導(dǎo)致皮膚增厚，對(duì)毛囊周期生長(zhǎng)具有調(diào)節(jié)作用，但是有研究報(bào)道的miR-203對(duì)調(diào)控皮膚生長(zhǎng)發(fā)育有重要調(diào)控作用，而在我們的研究中，在和田羊毛囊的3個(gè)時(shí)期對(duì)比組中均未發(fā)現(xiàn)miR-203的差異表達(dá)，說(shuō)明出生后毛囊周期中的miRNA調(diào)控與胚胎期及出生后皮膚中的miRNA調(diào)控是存在差異的，還需要進(jìn)一步研究。以上差異表達(dá)miRNAs可作為和田羊毛囊周期性變化過(guò)程中的關(guān)鍵候選miRNAs用于更進(jìn)一步的研究。

為了探究差異表達(dá)miRNAs可能和田羊周期性變化中進(jìn)程發(fā)揮的作用，篩選出各個(gè)時(shí)期差異表達(dá)的miRNAs之后，對(duì)差異表達(dá)miRNAs進(jìn)行靶基因預(yù)測(cè)并進(jìn)行功能富集分析。閆海龍等人[26]通過(guò)對(duì)毛囊干細(xì)胞與毛乳頭細(xì)胞進(jìn)行轉(zhuǎn)錄組和miRNAs測(cè)序分析，以及毛乳頭細(xì)胞和毛乳頭細(xì)胞來(lái)源外泌體進(jìn)行的miRNAs測(cè)序分析鑒定了4個(gè)與毛囊干細(xì)胞增殖分化相關(guān)的關(guān)鍵候選miRNAs，分別為miR-1、miR-151-3p、miR-182、miR-143-3p，并發(fā)現(xiàn)let-7b-5p能夠靶向調(diào)節(jié)預(yù)測(cè)的FOXC和MSI2的靶位點(diǎn)；江倩等人[27]發(fā)現(xiàn)miR-Let 7a在退行期的表達(dá)量高于生長(zhǎng)期的表達(dá)量，并進(jìn)一步證實(shí)miR-Let 7a可通過(guò)抑制兩個(gè)靶基因Cmyc和FGF5的作用進(jìn)而調(diào)控遼寧絨山羊的毛囊發(fā)育。Liu等[21]發(fā)現(xiàn)部分miRNA的靶基因可能通過(guò)Wnt信號(hào)通路參與調(diào)控絨山羊皮膚毛囊發(fā)育。本研究利用miRanda等軟件對(duì)差異表達(dá)miRNAs的靶基因預(yù)測(cè)[28，29]，發(fā)現(xiàn)KRT15是oar-miR-134-5p和oar-miR-432共同靶基因，oar-miR-495-5p靶基因是BMP4、oar-miR-433-3p和oar-miR-432有共同的靶基因KRT25。并對(duì)所有靶基因進(jìn)行GO和KEGG富集分析，發(fā)現(xiàn)這些靶基因主要富集在毛囊發(fā)育相關(guān)的激素水平調(diào)節(jié)、細(xì)胞分解代謝、周期蛋白結(jié)合、骨細(xì)胞發(fā)育等GO條目中以及與毛囊發(fā)育相關(guān)的AMPK、Hedgehog、TNF信號(hào)通路中，說(shuō)明miRNAs在毛囊周期性變化中可能通過(guò)這些通路來(lái)調(diào)節(jié)下游靶基因的表達(dá)來(lái)發(fā)揮調(diào)控作用，但miRNAs對(duì)預(yù)測(cè)的靶基因的具體靶向關(guān)系還需要進(jìn)行更進(jìn)一步的探究，需要進(jìn)行更深層次的驗(yàn)證。

4 結(jié)論

本研究通過(guò)轉(zhuǎn)錄組測(cè)序獲得了和田羊毛囊周期性變化中生長(zhǎng)期、退行期、休止期的轉(zhuǎn)錄組文庫(kù)，并對(duì)不同時(shí)期的轉(zhuǎn)錄組進(jìn)行比較，鑒定了143個(gè)已知的miRNAs并預(yù)測(cè)了759個(gè)新miRNAs，進(jìn)一步豐富了綿羊miRNA信息庫(kù)；并在退行期和生長(zhǎng)期比較組鑒定了173個(gè)差異表達(dá)miRNA，休止期和生長(zhǎng)期比較組鑒定了243個(gè)差異表達(dá)miRNA，休止期和生長(zhǎng)期比較組鑒定了71個(gè)差異表達(dá)miRNA，還通過(guò)靶基因預(yù)測(cè)和功能富集分析表明上述差異表達(dá)miRNAs可能參與和田羊毛囊周期性變化過(guò)程，但差異表達(dá)miRNAs及靶基因的互作調(diào)控關(guān)系還有待進(jìn)一步驗(yàn)證。