中國(guó)美利奴羊羊毛彎曲相關(guān)microRNA的篩選

于麗娟, 張艷花，拉扎特·艾尼瓦爾，徐新明，瑪爾孜婭·亞森，狄 江

(新疆畜牧科學(xué)院畜牧研究所，烏魯木齊 830011)

0 引 言

【研究意義】彎曲是羊毛最主要的物理性質(zhì)之一，與毛紡工藝有著密切的關(guān)系。研究表明，降低羊毛纖維彎曲度能減小纖維抗彎曲剛性和紗線(xiàn)厚度并進(jìn)而增加織品柔軟度，其對(duì)紡織品柔軟度的影響與纖維直徑有同等重要的作用[1]。了解羊毛彎曲形成的分子機(jī)理，將有助于利用分子育種培育符合市場(chǎng)需求羊毛彎曲類(lèi)型細(xì)毛羊新群體?！厩叭搜芯窟M(jìn)展】羊毛彎曲與毛囊中細(xì)胞的不對(duì)稱(chēng)分裂、毛纖維中正副皮質(zhì)的雙邊排列、不同角蛋白的不對(duì)稱(chēng)組成及角質(zhì)化過(guò)程等相關(guān)[2]。對(duì)小鼠、人等哺乳動(dòng)物毛發(fā)纖維彎曲的形成機(jī)制方面的研究取得了很大進(jìn)展，Wnt、EDA/EDAR、EGFR、BMP、NF-κB、Shh等經(jīng)典信號(hào)通路與毛發(fā)彎曲有著直接或間接的聯(lián)系[3]。【本研究切入點(diǎn)】MicroRNA(miRNA)為一類(lèi)大小在18～24個(gè)核苷酸的非編碼單鏈RNA，通過(guò)與靶mRNA結(jié)合位點(diǎn)的特異性結(jié)合，抑制靶mRNA翻譯或降解靶mRNA來(lái)發(fā)揮其重要的生理調(diào)節(jié)作用[4]。研究表明，miRNA幾乎參與了機(jī)體所有的代謝調(diào)節(jié)活動(dòng)，包括肌肉生長(zhǎng)發(fā)育、細(xì)胞周期變化、細(xì)胞增殖與凋亡、脂肪代謝、器官形成以及機(jī)體的免疫應(yīng)答等反應(yīng)[5]。隨著miRNA對(duì)皮膚和毛囊發(fā)育形成的研究的不斷深入，miRNAs是毛囊發(fā)育過(guò)程中重要的轉(zhuǎn)錄后調(diào)控因子，參與毛囊的發(fā)生和周期性生長(zhǎng)[6-9]。目前，miRNA對(duì)羊毛彎曲性狀形成的相關(guān)調(diào)控機(jī)制仍不清楚，需要研究miRNAs在羊毛彎曲形成過(guò)程中發(fā)揮的作用?！緮M解決的關(guān)鍵問(wèn)題】通過(guò)采用高通量測(cè)序技術(shù)對(duì)羊毛彎曲頻率有極端差異表型值細(xì)毛羊個(gè)體皮膚的miRNA組進(jìn)行測(cè)序并比較分析，篩選出與細(xì)毛羊羊毛彎曲性狀相關(guān)的差異表達(dá)miRNA，利用miRNAda與TargetScan軟件對(duì)差異表達(dá)的miRNA進(jìn)行靶基因預(yù)測(cè)，并對(duì)靶基因進(jìn)行簡(jiǎn)單的Gene Ontology與KEGG Pathway富集分析，分析差異 miRNA的預(yù)測(cè)靶基因參與的生物學(xué)功能及調(diào)控途徑，為進(jìn)一步解析羊毛彎曲形成分子機(jī)理、以及借助分子育種手段培育符合市場(chǎng)需求羊毛彎曲類(lèi)型細(xì)毛羊新群體提供理論指導(dǎo)。

1 材料與方法

1.1 材 料

綿羊?yàn)樾陆莩欠N羊場(chǎng)的中國(guó)美利奴羊，羊群在基本一致的條件下飼養(yǎng)。采集樣品主要以羊毛彎曲頻率性狀為依據(jù)。在羊毛品質(zhì)鑒定中，羊毛彎曲頻率一般分為大、中、小3個(gè)彎曲類(lèi)型，其中大彎為1 cm 羊毛纖維長(zhǎng)度內(nèi)有小于3 個(gè)彎曲、中彎為1 cm 羊毛纖維長(zhǎng)度內(nèi)有3 ～ 5 個(gè)彎曲、小彎為1 cm 羊毛纖維長(zhǎng)度內(nèi)大于5 個(gè)彎曲。個(gè)體鑒定時(shí)，選取大彎周歲母羊個(gè)體與小彎周歲母羊個(gè)體各10只。體側(cè)肩胛部皮膚剪毛消毒后，用皮膚采樣器分別采集其肩胛部皮膚組織1 cm2左右，投入裝有RNA保存液(RNAlater)的凍存管中。冷藏條件下帶回實(shí)驗(yàn)室，-80℃冰箱保存?zhèn)溆?。?duì)其毛樣進(jìn)行細(xì)度與羊毛彎曲頻率測(cè)定，選取羊毛彎曲頻率有極端差異的大彎、小彎個(gè)體各4只的皮膚樣品用于后續(xù)實(shí)驗(yàn)。

1.2 方 法

1.2.1 RNA提取、小RNA文庫(kù)構(gòu)建及Solexa測(cè)序

按照Trizol法提取細(xì)毛羊皮膚組織總RNA，采用Qubit 2.0 (美國(guó))試劑盒進(jìn)行RNA的精確定量；利用NanoPhotometer?分光光度計(jì)進(jìn)行純度檢測(cè)，同時(shí)利用Bioanalyzer 2100 系統(tǒng)的 RNA Nano 6000 Assay Kit進(jìn)行質(zhì)量完整性的檢測(cè)。樣品檢測(cè)合格后，使用Small RNA Sample Pre Kit構(gòu)建文庫(kù)，利用Small RNA的3’及5’端特殊結(jié)構(gòu)(5’端有完整的磷酸基團(tuán)，3’端有羥基)，以總RNA為起始樣品，直接將Small RNA兩端加上接頭，然后反轉(zhuǎn)錄合成cDNA。隨后經(jīng)過(guò)PCR擴(kuò)增，PAGE膠電泳分離目標(biāo)DNA片段，切膠回收得到cDNA文庫(kù)。利用安捷倫2100高靈敏度DNA芯片檢測(cè)文庫(kù)質(zhì)量。庫(kù)檢合格后，使用 Illumina Hiseq 2500 上機(jī)測(cè)序。

1.2.2 測(cè)序數(shù)據(jù)

對(duì)測(cè)序得到的原始數(shù)據(jù)(raw data)進(jìn)行質(zhì)控和過(guò)濾，包括去除低質(zhì)量、有5’接頭污染、沒(méi)有3’接頭序列、沒(méi)有插入片段及包含poly A/T/G/C的reads。進(jìn)行質(zhì)控后得到的數(shù)據(jù)為clean reads，用于后續(xù)分析。采用bowtie軟件與羊的參考基因組(Oar_v4.0)進(jìn)行比對(duì)，分析small RNA在參考序列上的分布情況。且利用miRNA所在的位置對(duì)已知(保守)miRNA進(jìn)行鑒定，而新miRNA則通過(guò)miREvo和mirdeep2這2個(gè)軟件進(jìn)行預(yù)測(cè)。利用DESeq軟件篩選差異表達(dá)的miRNA，利用TPM分別計(jì)算樣本huACF，huBDE表達(dá)量取log2，滿(mǎn)足Pvalue≤0.05篩選2組之間的差異miRNA。利用miRNAda、TargetScan靶基因預(yù)測(cè)軟件對(duì)miRNA進(jìn)行靶基因預(yù)測(cè)。根據(jù)miRNA與其靶基因的對(duì)應(yīng)關(guān)系，對(duì)差異miRNA的預(yù)測(cè)靶基因進(jìn)行簡(jiǎn)單的Gene Ontology與KEGG Pathway富集分析。

1.2.3 實(shí)時(shí)熒光定量PCR驗(yàn)證

隨機(jī)選擇4個(gè)差異表達(dá)的miRNA，檢測(cè)其在2組樣品的表達(dá)量。以U6作為內(nèi)參基因，分別根據(jù)4個(gè)miRNAs的成熟序列設(shè)計(jì)熒光定量PCR的上游引物，下游引物由試劑盒提供，引物由上海生工生物科技有限公司合成。表1

表1 qRT-PCR引物Table 1 The list of primers of qRT-PCR

采用poly(A)加尾法進(jìn)行cDNA的反轉(zhuǎn)錄，具體操作步驟參照康為世紀(jì)miRNA cDNA Synthesis Kit進(jìn)行。反應(yīng)體系為20 L：2×miRcute增強(qiáng)型miRNA定量預(yù)混試劑(含SYBR和ROX)10 μL，反向引物(10 μmol/L )0．4 μL，正向引物(10 μL )0.4 μL ，cDNA 2 μL，無(wú)核酶ddH2O。補(bǔ)至20 μL。PCR反應(yīng)條件：95℃ 預(yù)變性15 min；94℃ 變性20 s，60℃ 退火延伸34 s，共40個(gè)循環(huán)。每個(gè)樣品每個(gè)miRNA均進(jìn)行3次技術(shù)重復(fù)，結(jié)果采用2-ΔΔCt法計(jì)算miRNA的相對(duì)表達(dá)量，與測(cè)序結(jié)果對(duì)比，檢驗(yàn)測(cè)序結(jié)果的可靠性。

2 結(jié)果與分析

2.1 總RNA質(zhì)量檢測(cè)

研究表明，提取的總RNA進(jìn)行質(zhì)量檢測(cè)，其OD260/280值均在1.9～2.0，RIN值均高于8.0，總RNA純度及完整性均達(dá)到構(gòu)建文庫(kù)要求。表2

表2 總RNA質(zhì)量檢測(cè)Table 2 The quality test results of total RNA

2.2 中國(guó)美利奴羊皮膚組織測(cè)序

研究表明，經(jīng)Solexa測(cè)序和生物信息學(xué)分析,去除接頭序列以及低質(zhì)量序列，8個(gè)文庫(kù)共獲得81 820 492條總序列，經(jīng)過(guò)濾后共獲得79 455 352條純凈序列。細(xì)毛羊皮膚組織小RNA序列長(zhǎng)度大部分位于21～24 nt，其中長(zhǎng)度為22 nt的序列所占比列最大。表3，圖1

表3 數(shù)據(jù)過(guò)濾信息Table 3 The information of raw data filtering

圖1 小RNA的長(zhǎng)度分布Fig.1 Length distribution of sRNAs

2.3 小RNA的分類(lèi)注釋

利用Bowtie 軟件將得到的純凈序列在數(shù)據(jù)庫(kù)中比對(duì)后進(jìn)行分類(lèi)注釋?zhuān)玫礁鞣N類(lèi)型sRNA,包含有已知miRNA、新miRNA、rRN A、snRNA、snoRNA、tRNA和其他RNA。其中，各樣品中已知的miRNA所占的比例最大。表4

表4 各樣本sRNA分類(lèi)注釋Table 4 sRNA classification annotation results for each sample

2.4 小RNA序列基因組比對(duì)

針對(duì)每個(gè)樣本預(yù)處理后所有序列和去重復(fù)后的唯一序列，分別采用Bowtie軟件與羊的參考基因組(Oar_v4.0)進(jìn)行比對(duì)分析。獲取其在參考基因組上的位置信息。比對(duì)結(jié)果表明81%以上的片段能比對(duì)到綿羊基因組上。61%以上的種類(lèi)比對(duì)到綿羊基因組上。各樣品間的比對(duì)率存在一定的差異說(shuō)明樣品間小RNA存在一定差異。表5

表5 與參考序列比對(duì)Table 5 Alignment results of the reference sequence

2.5 已知miRNA分析和新miRNA預(yù)測(cè)

將mapped到參考序列上的reads，與miRBase中指定范圍序列進(jìn)行比對(duì)，得到各樣品匹配上的sRNA的詳細(xì)情況，包括匹配上的已知miRNA的二級(jí)結(jié)構(gòu)，各樣本中miRNA的序列、長(zhǎng)度、出現(xiàn)的次數(shù)等信息。比對(duì)上的已知成熟miRNA143條，已知 miRNA前體有106條。

整合miREvo(Wen et al.，2012)和mirdeep2(Friedlander et al.，2011)這些miRNA預(yù)測(cè)軟件來(lái)進(jìn)行新miRNA的分析。在小RNA文庫(kù)中總共預(yù)測(cè)得到306條miRNA前體，成熟miRNA也有282條。

2.6 miRNA差異表達(dá)

利用DESeq軟件篩選差異表達(dá)的miRNA，利用TPM分別計(jì)算樣本huACF，huBDE表達(dá)量取log2，滿(mǎn)足Pvalue≤0.05篩選2組之間的差異miRNA。共鑒定出差異表達(dá)miRNA25個(gè)，其中上調(diào)表達(dá)8個(gè)，下調(diào)表達(dá)17個(gè)。

2.7 差異miRNA靶基因預(yù)測(cè)與功能富集

利用miRNAda與TargetScan靶基因預(yù)測(cè)軟件對(duì)差異表達(dá)的miRNA進(jìn)行靶基因預(yù)測(cè)。預(yù)測(cè)結(jié)果顯示，大彎曲與小彎曲組的25個(gè)差異表達(dá)的miRNA共預(yù)測(cè)到了15 761個(gè)無(wú)重復(fù)的靶基因。根據(jù)miRNA與其靶基因的對(duì)應(yīng)關(guān)系，對(duì)差異miRNA的預(yù)測(cè)靶基因進(jìn)行簡(jiǎn)單的Gene Ontology與KEGG Pathway富集分析。GO分析結(jié)果表明，其主要參與跨膜信號(hào)受體活動(dòng)、信號(hào)轉(zhuǎn)導(dǎo)活動(dòng)、分子轉(zhuǎn)導(dǎo)活動(dòng)、膜的組成、細(xì)胞交流、生物過(guò)程、細(xì)胞過(guò)程的調(diào)節(jié)等過(guò)程。15 761個(gè)靶基因注釋到252個(gè)信號(hào)通路。其中，107個(gè)靶基因參與了Wnt信號(hào)通路，93個(gè)靶基因參與了AMPK信號(hào)通路，34個(gè)靶基因參與了Notch信號(hào)通路，161個(gè)靶基因參與了MAPK信號(hào)通路，68個(gè)靶基因參與了NF-kappa B信號(hào)通路。圖2

圖2 差異表達(dá)miRNA靶基因GO

2.8 miRNA qRT-PCR 驗(yàn)證

隨機(jī)選擇4個(gè)差異表達(dá)miRNA(novel_503、oar-miR-133、oar-let-7g、novel_9)進(jìn)行定量驗(yàn)證，結(jié)果表明：novel_503、oar-miR-133、oar-let-7g、novel_9在大彎曲組中上調(diào)表達(dá)，而在小彎曲組中下調(diào)表達(dá)，與高通量測(cè)序結(jié)果一致，證明高通量測(cè)序結(jié)果的可靠性。圖3，表6

圖3 miRNA 的qRT-PCR 驗(yàn)證Fig.3 qRT-PCR verification of miRNA

表6 用于qRT-PCR 驗(yàn)證的miRNA測(cè)序表達(dá)量數(shù)據(jù)Table 6 Sequencing data of miRNA verified by qRT-PCR

3 討 論

關(guān)于中國(guó)美利奴羊羊毛經(jīng)濟(jì)性狀方面的研究主要集中在對(duì)羊毛細(xì)度、長(zhǎng)度等。在中國(guó)美利奴羊羊毛彎曲性狀方面的研究?jī)H見(jiàn)狄江等[10]通過(guò)全基因組關(guān)聯(lián)分析檢測(cè)與細(xì)羊毛彎曲頻率性狀顯著關(guān)聯(lián)的SNPs 。

高通量測(cè)序技術(shù)已被廣泛應(yīng)用于動(dòng)植物轉(zhuǎn)錄組測(cè)序、小RNA測(cè)序等方面，在中國(guó)美利奴羊羊毛彎曲性狀形成方面的研究還未見(jiàn)報(bào)道。研究運(yùn)用高通量測(cè)序技術(shù)對(duì)羊毛彎曲頻率有極端差異表型值(大彎與小彎)細(xì)毛羊個(gè)體皮膚的miRNA組進(jìn)行測(cè)序并比較，共鑒定出425個(gè)miRNAs，其中包括143個(gè)已知miRNAs和282個(gè)新發(fā)現(xiàn)的miRNAs。鑒定出的miRNAs進(jìn)行表達(dá)量差異分析發(fā)現(xiàn)，在大彎曲組與小彎曲組共篩選到8個(gè)上調(diào)和17個(gè)下調(diào)的miRNAs。其中，差異表達(dá)miRNA中的oar-miR-133、oar-miR-379-5p的表達(dá)量趨勢(shì)與Liu的研究結(jié)果一致[11]。并且，高水平表達(dá)的oar-miR-379-5p在灘羊與藏綿羊皮膚毛囊中也檢測(cè)到[12,13]。差異表達(dá)miRNA中的oar-miR-21可能通過(guò)調(diào)控CNKSR2、KLF3和TNPO1的表達(dá)來(lái)影響毛囊發(fā)育[14]。差異表達(dá)miRNA中的oar-let-7g屬于let-7家族，該miRNA的靶基因與皮膚組織中毛囊的生長(zhǎng)和毛發(fā)的質(zhì)量相關(guān)[15],而且研究表明，let-7a、let-7b、let-7c、let-7i等與灘羊被毛卷曲、湖羊羔皮花紋及生長(zhǎng)相關(guān)[16,17]。推測(cè)let-7家族可能參與羊毛彎曲的形成。對(duì)于其他差異表達(dá)的miRNA,目前還沒(méi)有文獻(xiàn)報(bào)道這些miRNAs在羊毛彎曲形成中的具體作用，其功能尚需進(jìn)一步研究。

4 結(jié) 論

研究利用高通量測(cè)序技術(shù)比較了不同羊毛彎曲頻率中國(guó)美利奴羊皮膚組織中miRNA差異表達(dá)，共鑒定出了25個(gè)差異表達(dá)的miRNAs，其中，8個(gè)上調(diào)，17個(gè)下調(diào)。靶基因預(yù)測(cè)和功能富集分析表明，差異表達(dá)miRNA可能在羊毛毛囊發(fā)育相關(guān)通路中起作用。