野生盤羊×巴什拜羊雜交二代尾部脂肪差異表達(dá)lncRNA的鑒定與分析

蘇曉慧，何海迎，吐爾遜阿依·牙力坤，王 瓊,2，馬海玉，劉玲玲，劉武軍

(1.新疆農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)學(xué)院，烏魯木齊 830052；2.吐魯番職業(yè)技術(shù)學(xué)院，吐魯番 838000)

目前，市場上大部分脂尾型綿羊都是由瘦尾型綿羊經(jīng)過長期選擇進(jìn)化而來[1-2]，但是，由于當(dāng)代人的飲食和社會(huì)結(jié)構(gòu)的變化，綿羊尾部脂肪被視為一種能量浪費(fèi)，現(xiàn)階段小尾甚至無尾成為當(dāng)下綿羊尾部育種的目標(biāo)性狀[3]，因此，了解綿羊尾部脂肪沉積的分子機(jī)制對綿羊育種至關(guān)重要。近年來常通過研究長鏈非編碼RNA(long non-coding RNA，lncRNA)的作用來探討動(dòng)物脂肪沉積的過程，lncRNA是一類長度>200 nt的非編碼RNA分子，在多種層面上調(diào)控基因的表達(dá)水平[4-5]。Huang等[6]研究發(fā)現(xiàn)，萊蕪豬和長白豬背最長肌內(nèi)脂肪中存在55個(gè)差異lncRNAs，其中XLOC_046142、XLOC_004398和XLOC_015408可通過靶向MAP活化蛋白激酶2(MAPK2)等調(diào)節(jié)豬肌內(nèi)脂肪沉積；李富原等[7]研究了不同喂養(yǎng)期朗德鵝肝臟、腹部脂肪和胸肌中l(wèi)ncRNA的差異表達(dá)，結(jié)果發(fā)現(xiàn)，LOC106047490在3種組織中的表達(dá)量均顯著下降，推測可能是高糖、胰島素和不飽和脂肪酸誘導(dǎo)的，表明LOC106047490與脂代謝、糖代謝有關(guān)；Zhang等[8]對雞腹前脂肪細(xì)胞lncRNA研究發(fā)現(xiàn)，在前脂肪細(xì)胞中有1 336個(gè)lncRNAs在不同分化階段差異表達(dá)，且lncRNA是通過順式作用于鄰近的編碼基因控制腹部前脂肪細(xì)胞的分化，說明lncRNA是參與調(diào)控動(dòng)物脂肪沉積過程的重要因子，但報(bào)道主要集中在豬、雞等家養(yǎng)動(dòng)物。Ma等[9]通過對蘭州肥尾羊、小尾寒羊和小尾藏羊進(jìn)行高通量RNA測序發(fā)現(xiàn)，在尾部脂肪沉積過程中發(fā)揮重要作用的核心lncRNAs有：TCONS_00372767、TCONS_00171926、TCONS_00054953和TCONS_00373007；王鑫悅等[10]對野生盤羊與巴什拜羊雜交二代群體(F2代羊)的脂尾形狀進(jìn)行了K-means聚類分析，結(jié)果分成4個(gè)類型：小尾型、中尾型、大尾型和肥尾型，是理想的尾部脂肪變異和進(jìn)化的研究模型。鑒于此，本試驗(yàn)利用全轉(zhuǎn)錄組測序技術(shù)對新疆2個(gè)不同脂尾型綿羊群體巴什拜羊(肥尾型)和F2代羊(小尾型)的尾部脂肪進(jìn)行測序，篩選差異表達(dá)lncRNA，對其進(jìn)行生物信息學(xué)分析并表達(dá)驗(yàn)證，以期為解析綿羊脂肪沉積調(diào)控機(jī)制提供一定的理論依據(jù)。

1 材料與方法

1.1 材料

試驗(yàn)動(dòng)物來自于新疆塔城裕民縣，隨機(jī)挑選相同飼養(yǎng)條件下6月齡巴什拜羊及野生盤羊×巴什拜羊的雜交二代羊各6只，屠宰后立即采集尾部脂肪組織，并將新鮮樣品放入2 mL無菌無酶的低溫保存管中暫存于液氮中，樣品帶回實(shí)驗(yàn)室后-80 ℃保存?zhèn)溆谩ｋs交二代是以野生盤羊?yàn)楦副尽褪舶菅驗(yàn)槟副具M(jìn)行雜交，獲得F1代，再以巴什拜羊?yàn)楦副尽1代為母本進(jìn)行雜交，獲得F2代。

1.2 方法

1.2.1 總RNA提取及檢測 采用Trizol法提取尾部脂肪組織總RNA，利用1.0%瓊脂糖凝膠電泳檢測RNA完整性，使用NanoDrop分光光度儀檢測RNA濃度，巴什拜羊樣品編號為CG1～CG6；F2代羊樣品編號為TG1～TG6。

1.2.2 文庫構(gòu)建和測序 總RNA去除rRNA后純化回收，被隨機(jī)打斷成短片段，以片段化去除rRNA的RNA為模板，合成第一鏈cDNA，形成雙鏈結(jié)構(gòu)，利用T4 DNA聚合酶和Klenow DNA聚合酶將DNA的黏性末端修復(fù)為平末端，3′-末端加堿基A并加接頭，經(jīng)過PCR擴(kuò)增獲得最終測序文庫，文庫質(zhì)檢合格后采用Illumina Hiseq4000進(jìn)行測序。

1.2.3 轉(zhuǎn)錄組測序數(shù)據(jù)分析 篩選得到原始數(shù)據(jù)，與NCBI中參考基因組進(jìn)行序列比對，獲取比對數(shù)據(jù)，使用FastQC軟件對序列質(zhì)量(QC)進(jìn)行評估，過濾含有適配器、高N含量或大量低質(zhì)量堿基的reads，得到clean reads。利用Hisat2工具將reads映射到參考基因組[11]，使用StringTie[12]對Hisat2比對結(jié)果進(jìn)行拼接，并通過ASprofile[13]軟件獲得備選拼接的類型和每個(gè)樣本的相應(yīng)表達(dá)量。轉(zhuǎn)錄水平的計(jì)算采用FPKM(Fragments Per Kilobase of transcript per Million fragments mapped)片段法，采用DESeq軟件[14]分析樣本組間差異表達(dá)，將log2|FoldChange|≥2且FDR<0.05作為篩選標(biāo)準(zhǔn)。分析巴什拜羊和F2代羊尾部脂肪差異表達(dá)lncRNA，以及這些差異表達(dá)lncRNA功能的注釋分析。

1.2.4 lncRNA靶基因預(yù)測 利用順式(cis)或者反式(trans)方法對差異表達(dá)的lncRNA進(jìn)行靶基因預(yù)測。選取與lncRNA距離10 kb的上游或下游基因作為順式靶基因，采用Miranda方法對反式靶基因進(jìn)行預(yù)測。

1.2.5 GO功能和KEGG通路富集分析 GO數(shù)據(jù)庫用于預(yù)測和闡明基因產(chǎn)物在分子功能、生物過程及細(xì)胞組分方面的功能[15]。這些基因首先被映射到數(shù)據(jù)庫中的GO項(xiàng)，然后以校正后的P<0.05為閾值，計(jì)算每個(gè)GO項(xiàng)的基因數(shù)，以確定顯著富集的GO項(xiàng)。利用KEGG進(jìn)行通路富集分析[16]，確定基因參與的主要生化和信號通路，以校正后的P<0.05為閾值，確定顯著富集的KEGG通路。

1.2.6 引物設(shè)計(jì) 挑選顯著差異的7個(gè)與脂肪沉積相關(guān)的lncRNAs，根據(jù)測序所得序列設(shè)計(jì)引物；從NCBI數(shù)據(jù)庫中查找β-actin內(nèi)參基因序列(登錄號：HM067830.1)設(shè)計(jì)引物，引物序列見表1。

表1 引物信息

1.2.7 實(shí)時(shí)熒光定量PCR 使用濃度和純度合格的總RNA，根據(jù)PrimeScript RT Reagent Kit with gDNA Eraser說明書反轉(zhuǎn)錄為cDNA。以cDNA為模板進(jìn)行目的lncRNA的PCR擴(kuò)增。PCR反應(yīng)體系20 μL：2×PCR MasterMix 10 μL，上、下游引物各0.4 μL，cDNA 1.5 μL，RNase-free ddH2O 7.7 μL。取4 μL PCR產(chǎn)物與2 μL 6×Loading Buffer混勻，用1.5%瓊脂糖凝膠電泳檢測，180 V 25 min，使用凝膠成像系統(tǒng)采集圖像。

以β-actin為內(nèi)參基因，采用CFX96型熒光定量PCR儀進(jìn)行實(shí)時(shí)熒光定量PCR檢測。PCR反應(yīng)體系20 μL：TB Green Premix ExTaqⅡ 10 μL，上、下游引物各0.5 μL，cDNA 1 μL，RNase-free ddH2O 8 μL。PCR反應(yīng)程序：95 ℃ 30 s；95 ℃ 5 s，54 ℃ 30 s，共39個(gè)循環(huán)。每個(gè)樣品重復(fù)3次，采用2-ΔΔCt法計(jì)算目的lncRNA的相對表達(dá)量。

1.2.8 統(tǒng)計(jì)分析 運(yùn)用SPSS 25.0軟件的單因素方差分析(One-way ANOVA)，組間采用獨(dú)立樣本t檢驗(yàn)。結(jié)果以平均值±標(biāo)準(zhǔn)誤表示，P<0.05表示差異顯著；P<0.01表示差異極顯著。

2 結(jié) 果

2.1 總RNA檢測及測序數(shù)據(jù)統(tǒng)計(jì)

提取的總RNA具有較好的完整性，28S和18S條帶明顯，5S條帶暗淡，未發(fā)生降解，可用后續(xù)試驗(yàn)。

對巴什拜羊和F2代羊的尾部脂肪分別進(jìn)行RNA-Seq高通量測序后，共獲得1 438 500 450條原始序列數(shù)據(jù)(raw reads)；對原始數(shù)據(jù)進(jìn)行過濾后，共獲得1 338 966 230條高質(zhì)量序列數(shù)據(jù)(clean reads)；Q20在所有樣品中均是99%以上，說明單個(gè)堿基的錯(cuò)誤率非常低；GC含量在47%～49%之間(表2)，更加驗(yàn)證了測序結(jié)果的可靠性。

表2 測序數(shù)據(jù)統(tǒng)計(jì)表

2.2 序列比對

由表3可知，各樣品的原始序列數(shù)據(jù)與參考基因組的比對率在85.28%～91.06%之間，單一映讀取率在64.11%～71.57%之間，多個(gè)映讀取率在19.12%～21.43%之間，說明巴什拜羊和F2代羊相似性較高，測序結(jié)果質(zhì)量可靠，可用于后續(xù)分析。

表3 基因組比對情況

2.3 差異表達(dá)lncRNA的篩選

使用DESeq2軟件篩選巴什拜羊和F2代羊中尾部脂肪的差異表達(dá)lncRNA，并繪制火山圖，結(jié)果見圖1。由圖1可知，共篩選到728個(gè)差異lncRNAs，與巴什拜羊相比，F(xiàn)2代羊中270個(gè)lncRNAs表達(dá)下調(diào)，458個(gè)lncRNAs表達(dá)上調(diào)。

2.4 lncRNA靶基因預(yù)測

lncRNA可調(diào)控臨近mRNA表達(dá)，利用lncRNA與蛋白編碼基因的共定位預(yù)測lncRNA的靶基因，以便進(jìn)一步對獲得的lncRNA進(jìn)行功能研究，其中上、下游100 kb作為共定位的閾值，結(jié)果顯示，728個(gè)差異表達(dá)的lncRNAs靶向17 850個(gè)靶基因，其中MSTRG.37980、MSTRG.38164、MSTRG.36912、MSTRG.36913、MSTRG.8169、MSTRG.24995和MSTRG.31389可以靶向與脂肪沉積和脂肪代謝相關(guān)的基因，如SCD、GPAM、THRSP和FASN等。lncRNA靶向與脂肪沉積和脂肪代謝相關(guān)基因的部分結(jié)果見表4。

Y軸>1.3且X軸≥1區(qū)域代表上調(diào)；Y軸>1.3且X軸≤-1區(qū)域代表下調(diào)The area of Y>1.3 and X≥1 means up-regulated；The area of Y>1.3 and X ≤-1 means down-regulated圖1 差異表達(dá)lncRNA火山圖Fig.1 Volcano diagram of differentially expressed lncRNA

表4 與脂肪沉積相關(guān)的部分lncRNAs及其靶基因

2.5 GO功能和KEGG通路富集分析

通過GO功能分析，差異顯著lncRNA的靶基因參與到634個(gè)功能分類，注釋到生物過程、細(xì)胞組分和分子功能三大類，主要富集于T細(xì)胞受體信號通路(T cell receptor signaling pathway)、T細(xì)胞受體復(fù)合體(T cell receptor complex)、SH3/SH2適配器的活動(dòng)(SH3/SH2 adaptor activity)、正向調(diào)節(jié)T細(xì)胞增殖(positive regulation of T cell proliferation)、GTPase活性的正調(diào)控(positive regulation of GTPase activity)、陰性胸腺T細(xì)胞選擇(negative thymic T cell selection)、肽-酪氨酸磷酸化的負(fù)調(diào)控(negative regulation of peptidyl-tyrosine phosphorylation)、MHCⅡ類蛋白復(fù)合物(MHC class Ⅱ protein complex)、膜的積分分量(integral component of membrane)、先天性免疫應(yīng)答(innate immune response)、炎癥應(yīng)答(inflammatory response)、免疫突觸(immunological synapse)、免疫反應(yīng)(immune response)、質(zhì)膜的外側(cè)(external side of plasma membrane)、對病毒的防御反應(yīng)(defense response to virus)、染色質(zhì)沉默(chromatin silencing)、趨化性(chemotaxis)、chemokine-mediated信號通路(chemokine-mediated signaling pathway)、細(xì)胞表面(cell surface)、B細(xì)胞受體信號通路(B cell receptor signaling pathway)。其中顯著富集到膜的積分分量、GTPase活性的正調(diào)控、炎癥應(yīng)答、免疫反應(yīng)、質(zhì)膜外側(cè)、對病毒的防御反應(yīng)等，最顯著富集到了膜的積分分量(圖2)。

KEGG通路富集結(jié)果顯示，巴什拜羊和F2代羊尾脂組織差異表達(dá)lncRNA的靶基因共富集到76條通路，主要涉及病毒致癌作用(viral carcinogenesis)、Ⅰ型糖尿病(type Ⅰ diabetes mellitus)、T細(xì)胞受體信號通路(T cell receptor signaling pathway)、全身性紅斑狼瘡(systemic lupus erythematosus)、金黃色葡萄球菌感染(Staphylococcusaureusinfection)、原發(fā)性免疫缺陷(primary immunodeficiency)、破骨細(xì)胞分化(osteoclast differentiation)、自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性(natural killer cell mediated cytotoxicty)、麻疹(measles)、利什曼病(Leishmaniasis)、IgA生產(chǎn)的腸道免疫網(wǎng)絡(luò)(intestinal immune network for IgA production)、造血玻璃紙血統(tǒng)(hematopoietic cel lineage)、移植物抗宿主病(graft-versus-host disease)、Fcγ-R介導(dǎo)的吞噬作用(Fc gamma R-mediated phagocytosis)、細(xì)胞因子受體相互作用(cytokine-cytokine receptor interaction)、趨化因子信號通路(chemokine signaling pathway)、B細(xì)胞受體信號通路(B cell receptor signaling pathway)、自身免疫性甲狀腺疾病(autoimmune thyroid disease)、同種異體移植物排斥(allograft rejection)、醇中毒(alcoholism)，主要通路見圖3。

2.6 差異表達(dá)lncRNA的實(shí)時(shí)熒光定量PCR驗(yàn)證

為了驗(yàn)證RNA-Seq數(shù)據(jù)的準(zhǔn)確性，選擇巴什拜羊和F2代羊尾部脂肪之間差異表達(dá)的7個(gè)lncRNAs進(jìn)行了實(shí)時(shí)熒光定量PCR驗(yàn)證，結(jié)果見圖4。 由圖4可知，MSTRG.37980、MSTRG.38164、MSTRG.36912、MSTRG.36913、MSTRG.24995、MSTRG.8169和MSTRG.31389的實(shí)時(shí)熒光定量PCR與RNA-Seq測序結(jié)果趨勢相同，說明RNA-Seq結(jié)果可靠。實(shí)時(shí)熒光定量PCR檢測結(jié)果顯示，7個(gè)差異表達(dá)lncRNAs中有5個(gè)(MSTRG.37980、MSTRG.38164、MSTRG.36912、MSTRG.36913和MSTRG.24995)在F2代羊中表達(dá)下調(diào)，其中MSTRG.37980、MSTRG.38164和MSTRG.36912表達(dá)量極顯著低于巴什拜羊(P<0.01)，MSTRG.36913和MSTRG.24995表達(dá)量顯著低于巴什拜羊(P<0.05)；MSTRG.8169和MSTRG.31389在F2代羊中表達(dá)上調(diào)，且顯著高于巴什拜羊(P<0.05)。

圖2 差異表達(dá)lncRNA靶基因的GO功能富集散點(diǎn)圖Fig.2 Rich distribution map of GO function with differentially expressed lncRNA target genes

圖3 差異表達(dá)lncRNA靶基因的KEGG通路富集散點(diǎn)圖Fig.3 Rich hub diagram of KEGG pathway with differentially expressed lncRNA target genes

3 討 論

脂肪不僅是一種重要的儲能物質(zhì)，而且在調(diào)控動(dòng)物體能量平衡過程中發(fā)揮著重要作用，并與動(dòng)物產(chǎn)肉量、肉品質(zhì)、肉風(fēng)味等產(chǎn)肉性狀密切相關(guān)。隨著人們生活水平的不斷提高，飲食結(jié)構(gòu)發(fā)生改變，對食物質(zhì)量的要求也越來越高，也逐漸發(fā)現(xiàn)食用過多的脂肪對人體健康有害，對現(xiàn)代消費(fèi)者而言尾部脂肪并不理想，并降低了綿羊肉的價(jià)值，因此小尾型綿羊備受青睞[17-18]。但是，脂肪沉積是一個(gè)復(fù)雜的過程，涉及脂肪細(xì)胞的增殖與分化、發(fā)生與凋亡、甘油三酯的合成與水解等生化過程，受到多種分子和途徑的調(diào)控。此外，各種生物大分子之間也存在相互作用，如RNA-RNA、DNA-蛋白質(zhì)、RNA-蛋白質(zhì)、蛋白質(zhì)-蛋白質(zhì)等，因此關(guān)于脂肪沉積相關(guān)機(jī)制的探索十分重要。

目前已對動(dòng)物脂肪沉積進(jìn)行過轉(zhuǎn)錄組測序研究，Miao等[19]研究了lncRNA在金華豬和長白豬肌內(nèi)脂肪組織中的表達(dá)水平，發(fā)現(xiàn)119個(gè)差異表達(dá)lncRNAs，其中有6個(gè)共表達(dá)lncRNAs參與脂肪沉積和脂質(zhì)代謝通路；鄭竹清等[20]通過CPC(coding potential calculator)和CPAT(coding potential assessment tool)預(yù)測發(fā)現(xiàn)1 472個(gè)lncRNAs在山羊脂肪細(xì)胞中表達(dá)，其中29個(gè)在內(nèi)蒙古絨山羊成熟前后肌內(nèi)脂肪細(xì)胞中顯著表達(dá)，推測這29個(gè)lncRNAs影響細(xì)胞的分泌和生長調(diào)節(jié)，從而在脂肪細(xì)胞的發(fā)育中發(fā)揮重要作用；晁天樂[21]以杜泊羊和小尾寒羊?yàn)椴牧希瑱z測到47個(gè)lncRNAs參與調(diào)節(jié)十二指腸和大腸的轉(zhuǎn)錄后表達(dá)，發(fā)現(xiàn)這些lncRNAs可能在碳水化合物代謝和脂質(zhì)代謝中發(fā)揮主要的作用；Ma等[9]對蘭州肥尾羊、小尾寒羊和藏羊的尾部脂肪進(jìn)行高通量測序，在差異表達(dá)的6個(gè)lncRNAs中發(fā)現(xiàn)4個(gè)lncRNAs(TCONS-00372767、TCONS-00171926、TCONS-00054953和TCONS-00373007)可能作為核心lncRNA在尾部脂肪沉積過程中發(fā)揮關(guān)鍵作用。本研究通過轉(zhuǎn)錄組測序共篩選得到728個(gè)差異lncRNAs，共鑒定到17 850個(gè)靶基因，并篩選到至少7個(gè)與綿羊尾脂沉積相關(guān)的lncRNAs，說明lncRNA在動(dòng)物脂肪沉積中發(fā)揮著至關(guān)重要的調(diào)控作用。本研究采集的是不同尾型的巴什拜羊和F2代羊的尾部脂肪組織，這2種羊與前人研究所用的試驗(yàn)材料有較近的親緣關(guān)系，相對降低了種間差異性，保證了結(jié)果的準(zhǔn)確性和可靠性。

本研究鑒定到的靶基因可富集到76條信號通路，多條與脂代謝相關(guān)(如趨化因子信號通路)，且一些差異靶基因?qū)d羊尾部脂肪沉積的調(diào)控作用得到了部分驗(yàn)證，如，MSTRG.37980的靶基因SCD編碼硬脂酰輔酶A去飽和酶，此酶在PPAR信號通路中發(fā)揮關(guān)鍵作用，以激活脂肪細(xì)胞的成脂分化[22-23]，且在動(dòng)物肝臟和脂肪組織中高度表達(dá)[24]；MSTRG.38164的靶基因GPAM可以催化磷脂生物合成過程中甘油三酯合成的第一反應(yīng)，調(diào)節(jié)磷脂和甘油三酯的水平[25-26]，且GPAM基因在所有脂肪組織中均有表達(dá)，在肝臟和脂肪組織中表達(dá)最高[27]。本試驗(yàn)結(jié)果為從lncRNA的角度分析綿羊尾部脂肪在轉(zhuǎn)錄組中的差異及綿羊脂肪沉積的進(jìn)一步研究提供理論依據(jù)。

4 結(jié) 論

本研究共篩選出728個(gè)差異表達(dá)lncRNAs，并對差異顯著的lncRNA靶基因進(jìn)行分析，得到參與疾病、免疫、蛋白修飾、細(xì)胞代謝等相關(guān)功能分類634個(gè)，共涉及76條信號通路，靶向與脂肪沉積相關(guān)的基因SCD、GPAM、THRSP、FASN等，表明lncRNA對綿羊脂肪沉積具有重要的調(diào)節(jié)作用。