基于轉(zhuǎn)錄組測序篩選和田羊毛囊周期性變化過程中差異表達miRNAs以及初步分析

孫 爽，阿布來提·蘇來曼，李 彬，李志強，侯晨曦，決肯·阿尼瓦什，依明·蘇來曼*

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

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

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

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

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

1 材料與方法

1.1 試驗材料

1.1.1 試驗動物

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

1.1.2 試劑及儀器

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

1.2 試驗方法

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

1.2.2 cDNA的制備和文庫的構(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程序進行反轉(zhuǎn)錄，產(chǎn)物置于冰箱-20℃保存。

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

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

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

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

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

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

1.2.4 差異表達miRNAs靶基因的預測

利用miRanda、miRDB、TargetScan等工具對差異表達miRNAs進行靶基因預測，取預測結(jié)果的交集作為預測的最終結(jié)果。

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

將預測的靶基因分別利用GO seqR和KOBAS進行GO和KEGG功能富集來分析靶基因的功能及信號通路以研究差異表達miRNAs潛在的功能。

1.2.6 qRT-PCR

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

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

2 結(jié)果與分析

2.1 提取RNA檢測

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

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

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

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

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

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

圖1 S0617149I樣本的reads長度占比情況

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

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

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

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

2.2.3 miRNAs鑒定結(jié)果

共鑒定了143個已知的miRNAs，并發(fā)現(xiàn)這些miRNAs屬于75個家族，其中oar-let-7家族的最多，有7個：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個：oar-miR-200a、oar-miR-200b、oar-miR-200c；并預測了759新miRNAs。

2.2.4 差異表達miRNAs的篩選

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

圖2 新預測miRNAs中P2/P1比較組差異表達miRNAs的火山圖

為了進一步鑒定差異表達miRNA在和田羊3個不同時期中的分布情況，對已知的差異表達miRNAs進行了維恩圖分析(圖3)，在維恩圖分析中可直觀展現(xiàn)出各差異比較組合之間共有與特有的差異表達miRNAs數(shù)量，在P2/P1時期檢測到了6個時期特異性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時期檢測到了2個時期特異性miRNAs：oar-miR-134-3p、oar-miR-323a-5p。提示這些時期特異性表達miRNAs在毛囊周期性變化過程有重要作用。

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

2.3 靶基因預測

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

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

2.4 靶基因GO和KEGG富集分析

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

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

2.5 熒光定量PCR

結(jié)果如圖6所示，所有選擇的7個miRNAs表達趨勢與轉(zhuǎn)錄組測序結(jié)果一致，說明轉(zhuǎn)錄組測序結(jié)果準確可靠。

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

3 討論

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

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

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

4 結(jié)論

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