基于RNA-Seq 技術鑒定撒壩豬與杜洛克豬肝臟組織差異表達基因

吳正明，邱丙姍，李志勛，程 娜，翁亞煩，李再磊，趙桂英*

（云南農業大學動物科學技術學院，云南昆明 650201）

我國很多地區都分布有獨具特色的地方豬種資源，它為我國養豬生產可持續發展奠定了基礎，也是滿足人類需求的一份重要基因庫。近年來，隨著組學技術的不斷發展，地方豬種資源的保護及開發利用越來越受到重視。有相關研究表明，大多數中國地方豬具有肉質好、適應性強和繁殖力高等特性[1]。研究調控中國地方豬種優質性狀的分子機制，發掘和鑒定中國地方豬種優質性狀的主效基因對中國地方豬種資源的保護及開發利用具有重要意義。轉錄組（Transcriptome）指某一物種在某一特定生理狀態下的單個細胞或細胞群體內幾乎全部轉錄本，即轉錄后所有mRNA 的總稱[2]，它不僅是研究基因結構、功能和不同基因間表達差異的基礎，也是探索未知功能基因的方法之一[3]。轉錄組測序即RNASeq（Ribonucleic Acid Sequencing），是一種與高通量測序技術相結合從而對轉錄組進行分析的技術[4]，其具有高靈敏度、成本低和不受物種限制等特點，是一種被運用在轉錄組學上比較精確的測定分析方法[5]。現今RNA-Seq 普遍運用在各種畜禽重要分子機制和鑒定相關功能基因的研究領域。基于RNA-Seq 技術，研究者采用了不同的策略開展了不同品種豬各組織間存在的差異表達基因（Differentially Expressed Genes，DEGs）的鑒定和篩選，并對不同豬種間的種質差異進行相關分析。Liu 等[6]利用RNA-Seq 技術對民豬和大白豬背最長肌建立了轉錄組文庫，通過基因表達圖譜鑒定了1 371 個DEGs，為研究不同品種豬肉質差異的形成機制提供了參考。Ropka 等[7]通過對5 個品種944 個樣本的RNA-Seq 數據分析，找到了2 個與肌纖維直徑和密度極顯著相關的單核苷酸多態性位點（Single Nucleotide Polymorphisms，SNPs），說明基于RNA-Seq數據的SNPs 分析可用于尋找和定位基因。Liu 等[8]采用RNA-Seq 技術分別對妊娠期不同日齡（40、55、63、70、90 d）通城豬與約克夏豬的骨骼肌轉錄圖譜進行了綜合比較，分別檢測到了1 317 個和691 個DEGs，結果顯示55 d 是導致通城豬和約克夏豬發育差異的關鍵階段，并闡明了2 個品種間發育差異的6 個關鍵候選基因。Huang 等[9]比較了金華豬和長白豬肝臟中環狀RNA 的差異表達情況，挖掘到336 個DEGs，結合GO（Gene Ontology）分析發現其宿主基因大部分與代謝途徑相關。上述研究為轉錄組層次闡明豬種間的差異表達奠定了分子基礎。盡管通過RNA-Seq 技術挖掘了大批與生產和肉質性狀相關的候選基因，但只在國外豬種和我國少數地方豬種的研究中有相關報道，而云南地方豬種的相關研究更是鮮有報道，大量優良地方豬種的遺傳基礎仍很薄弱，仍需深入研究。本研究擬選用生產性能和肉質性能存在較大差異的云南地方豬——撒壩豬（肉脂兼用型品質）和國外引進的杜洛克豬（瘦肉型品種）為材料，采用RNA-Seq 技術篩選撒壩豬和杜洛克豬肝臟組織的DEGs，并對篩選的DEGs 進行GO 功能富集分析、KEGG 通路富集分析和可變剪接分析，對調控撒壩豬生長性狀和肉質性狀的相關候選基因進行挖掘，為進一步闡明中國地方豬種生產及肉質性能遺傳機制及探索肝臟組織的轉錄譜提供理論依據。

1 材料和方法

1.1 實驗動物和樣品制備 選擇胎次、產仔日齡相近的撒壩豬和杜洛克豬肥豬各4 頭，飼養于云南農業大學動物科學技術學院實習豬場，在相同飼養條件下飼養至260日齡（撒壩豬體重近90 kg，杜洛克豬體重近110 kg）后屠宰，屠宰后用無RNA 酶的凍存管采集肝臟樣品，樣品在低溫條件下運回實驗室，并存于-80℃的冰箱中。

1.2 總RNA 的提取和文庫的構建 用Trizol 試劑分別提取每頭撒壩豬和杜洛克豬肝臟組織中的總RNA；Nanodrop、瓊脂糖凝膠電泳、Qubit、Aglient 2100 等軟件檢測總RNA 的濃度、純度和完整性，提取完整性和純度檢測符合要求后的總RNA 樣本。向Oligo（dT）的磁珠富集過的mRNA 中加入fragmentationbuffer，把mRNA 打斷片段。以mRNA 為模板、六堿基為隨機引物合成第1 鏈cDNA；加入dNTPs、緩沖液、DNA polymerase I 和RNase H 合成第2 鏈cDNA。先 對cDNA 進行純化并洗脫、做末端修復和連接測序接頭，然后用瓊脂糖凝膠電泳進行片段大小選擇，最后進行PCR 擴增富集cDNA 即完成文庫制備。

1.3 數據質控 為了知曉測序數據質量的優劣，對在Illumina 測序平臺上完成測序的轉錄組數據下機后進行質控，質控包含3 個方面：堿基質量的分布、堿基平衡性和重復序列水平。對原始序列進行質量評估的目的是通過評估數據質量的情況來引導是否有必要進行下游分析。堿基平衡性分析可以清晰檢測到數據中是否具有AT、GC 分散情況，正常情況下4 種堿基出現頻率是一致的，并且位置上也沒有差異，當個別位置堿基的比例發生偏斜現象時，往往提示有超出比例的核酸污染。測序深度越高，越容易產生一定程度的重復，這屬于正常現象；但如果重復的程度很高，可能有大量核酸受到了污染。

1.4 基因差異表達與富集分析 以每百萬序列中來自于某基因每千堿基長度的序列數（Reads Per Kilobase per Million reads，RPKM）作為基因表達量的衡量水平，計算公式如下：

為了對假陽性率（False Positive rate）的控制，需結合假發現率（False Discovery Rate，FDR）即P-value和RPKM 值倍數變化（Fold Change，FC）對DEGs 進行篩選。以|logFC|>1 且FDR<0.05 為條件篩選出撒壩豬與杜洛克豬肝臟組織表達量顯著差異的基因。為了比較2 個豬種間基因表達量的差異，以lgRPKM 值進行聚類分析。使用edgeR[7]對撒壩豬與杜洛克豬肝臟組織轉錄組數據進行差異分析，并利用DESeq[8]軟件對DEGs進行GO 富集分析和KEGG 通路顯著富集分析。

1.5 差異剪接事件檢測 rMATS 軟件具有對轉錄組數據進行自動檢測和差異可變剪接分析的優點，因此可用該軟件對成對樣品進行可變剪接事件分析，主要通過rMATS 統計模型的方式檢測不同樣品的可變剪切類別及相應的表達量（ψ），并設定|Δψ|=|ψ1 ? ψ2|<0.05 為差異表達閾值。

2 結果

2.1 轉錄組測序數據 利用RNA-Seq 進行撒壩豬和杜洛克肝臟組織轉錄組測序，每個文庫產生1 億多個原始序列，過濾掉低質量、帶接頭和重復的原始序列后，8 個RNA-Seq 文庫中均獲得8 900 萬個以上的過濾序列，經過濾處理后，約有97.13%的序列被定位到杜洛克豬的參考基因組（Sus scrofa 11.1），低于撒壩豬的97.34%。杜洛克豬的外顯子、內含子和基因組間區域分別占85.83%、3.62%和10.28%。而約克夏豬區域的外顯子、內含子和基因組區間分別占89.56%、4.23% 和0.588%。2 種豬所比對到參考基因組上的序列比例錯誤率小于0.001（Q30）的均大于99%。結果表明：RNASeq 得到的數據利用率高，所選擇的參考基因組能夠滿足本實驗后續分析的需求。

2.2 篩選DEGs 從差異基因火山圖（圖1）中可以非常直觀并且合理地篩選出撒壩豬和杜洛克豬肝臟組織的DEGs。其中x 軸為lgFC 值，y 軸為lgFDR 值。結果顯示從撒壩豬和杜洛克豬肝臟組織間共鑒定出2 326 個差異表達基因，在撒壩豬上調顯著表達基因有1 049 個，在杜洛克豬顯著上調表達的基因有1 277 個。

本研究同時還通過熱圖顯示了撒壩豬和杜洛克豬差異表達基因的表達水平（圖2），可以清楚地觀察到這2 個品種豬肝臟組織的基因表達量存在差異。

圖2 差異表達基因熱圖

2.3 DEGs 的功能分析 將2 豬種篩選到的2 326 個表達顯著差異的基因與GO 數據庫和KEGG 數據庫進行對比。在2 個數據庫中共有926 個DEGs 被注釋，其中在GO 數據庫中共檢測到521 個差異表達基因，在KEGG數據庫中共檢測到651 個差異表達基因。

2.3.1 GO 富集分析 將521 個DEGs 注釋到GO 數據庫進行功能分類，分類內容包括：分子功能（Molecular Function，MF）、生物過程（Biological Process，BP）及細胞組分（Cellular Component，CC）3 個方面。其中MF 分成7 個小類，BP 分成33 個小類，CC 分成12 個小類。在BP 分類中參與代謝過程的DEGs 最多；在CC分類中參與調控細胞內環境和細胞外環境的DEGs 最多；在MF 分類中調控催化活性的DEGs 最多（圖3）。

圖3 差異表達基因GO 注釋

首先，通過與整個基因組背景相比并應用超幾何檢驗，挑選出撒壩豬與杜洛克豬肝臟組織的DEGs 顯著富集的GO 條目。圖4 為按P-value 值排序選擇的20 個最顯著GO 條目作出的GO 富集分析散點圖。用-log10（P-value）表示富集的顯著性（橫軸），該值越大富集得越顯著，縱軸表示富集到的GO 條目數。圓點大小及顏色的深淺表示該GO 條目包含的差異基因數目和富集程度。

圖4 差異表達基因GO 富集分析散點圖

2.3.2 KEGG 富集分析 首先，通過與整個基因組背景相比并應用超幾何檢驗，篩選出撒壩豬和杜洛克豬肝臟組織DEGs 中顯著性富集的通路。共有47 條顯著富集的通路，其中有24 條通路發生上調，有23 條通路發生下調。對DEGs 做KEGG 通路富集后，統計出它們參與的KEGG 代謝通路，如圖5 所示：橫軸表示富集的顯著性，P值越高表示富集得越顯著，縱軸表示富集的KEGG 通路數。圓點大小表示KEGG 通路包含的差異基因數目，顏色深淺表示KEGG 富集因子富集的程度。這里按P值進行排序，挑選顯著性最高的前20 個KEGG 條目作出相應的KEGG 通路富集散點圖；20 條有關生產和肉質性能的代謝通路具體情況如表1 所示。

表1 撒壩豬與杜洛克豬之間表達差異顯著的有關生產性能及肉質性能的相關通路

圖5 差異表達基因KEGG 富集分析散點圖

2.4 可變剪接分析 本研究通過對撒壩豬和杜洛克豬肝臟組織DEGs 可變剪接事件的檢測發現，共有5 種可變剪接事件類型：外顯子互斥（Mutually Exclusive Exons，MEX）、外顯子跳躍（Skipped exon，SE）、內含子保留（Retained intron，RI）、5'端可變剪接（Alter native 5’splice site，A5SS）和3'端可變剪接（Alternative 3'splice site，A3SS）。交叉點讀數（Junction reads）中檢測到總可變剪接事件數分別為21 487、3 489、438、746 和1 319；交叉點讀數檢測到的撒壩豬和杜洛克豬差異可變剪接事件比例分比為369:266、820:498、21:20、7:24 和14:22；目標讀數（Reads on target）和交叉點讀數檢測到的總可變剪接事件數分別為21 491、3 489、439、747 和1 329；目標讀數和交叉點讀數檢測到的撒壩豬和杜洛克豬存在差異可變剪接事件比例為389:282、817:495、20:20、9:29 和13:20（表2）。

表2 可變剪接事件分析

3 討 論

本研究通過對GO 數據庫和KEGG 數據庫顯著富集通路的分析，初步篩選出PPAR 信號通路等20 條有關生產和肉質性能的代謝通路，撒壩豬和杜洛克豬肝臟組織間共鑒定出2 326 個差異表達基因，其中大部分通路內富集的DEGs 與脂類代謝調控有關，包含脂蛋白酯酶（Lipoprotein Lipase,LPL）、長鏈脂酰輔酶A 合成酶3（Long-chain Acyl-CoA Synthetase3，ACSL3）、6-磷酸果糖激酶（6-phesphofructokinase,PFKL）、脂肪酸結合蛋白4（Fatty Acid Binding Protein 4,FABP4）、瘦素（Leptin，LEP）、脂肪轉運蛋白4（Latty Acid Trans port Proteins 4，FATPs4）、脂聯素（Adiponectin，ADIPOQ）、脂肪酸結合蛋白1（Fatty Acid Binding Protein1，FABP1）、長鏈輔酶A 脫氫酶（Long-chain Acyl-CoA dehydrogenase，ACADL）等。

肝臟是研究脂質代謝的重要組織，也是脂肪酸合成最主要的場所，因此肝臟是研究代謝過程、免疫、生長和肉質的首選組織[10]。李建平等[11]通過對育成豬飼喂不同脂肪源飼料飼養屠宰后,采肝臟組織樣品做轉錄組差異分析，結果表明對肝臟組織差異基因的GO 富集分析和KEGG 富集分析中均有與脂類代謝相關基因被注釋；Ge 等[12]選取低肌內脂肪酸和高肌內脂肪酸的巢湖鴨肝臟組織為實驗樣本進行轉錄組測序，篩選出9 個與脂肪沉積的相關基因；本研究結果和以往的研究結果大致相同，表明肝臟確實是研究脂質代謝的重要組織。

撒壩豬LPL基因表達量較高，是杜洛克豬的28.3倍，提示它和脂肪沉積能力較高的性狀有關。脂蛋白脂酶（Lipoprotein Lipase，LPL）位于豬14 號染色體上，屬于水解酶家族，LPL對脂質在機體內的代謝及相關代謝物的轉運發揮重大作用，它能將血乳糜微粒和極低密度脂蛋白中的三酰甘油分解為脂肪酸和甘油，并將脂肪酸和甘油轉運至機體相應組織中儲存[13]。Wang 等[14]對藏豬和滇南小耳豬兩個中國地方豬與大白和長白豬兩個外來豬背最長肌進行同位素標記的相對和絕對定量，找到LPL等10 個與脂肪沉積能力相關的基因，其中藏豬-滇南小耳豬組LPL基因表達量是大白豬-長白豬組的2.84 倍，從而LPL被認為與藏族和滇南小耳豬較高的脂肪沉積能力有關；本研究結果與以往的研究結果具有一致性。

ADIPOQ基因是被人們所發現唯獨與肥胖呈負相關的激素，ADIPOQ基因通過受體1（AdipoR1）、受體2（AdipoR2）以及T-黏鈣蛋白相結合來調控機體內脂肪的含量，當生物體內脂肪含量過高時，就會減少ADIPOQ基因的表達量和ADIPOQ在血液中的含量[15]。Cirera 等[16]分別采集哥本哈根迷你豬和杜洛克豬背最長肌、背部脂肪、肝臟組織，測定ADIPOQ基因在2豬種間3 個組織中的表達量，結果顯示：ADIPOQ基因哥本哈根迷你豬和杜洛克豬背最長肌、背部脂肪、肝臟組織中的表達量并沒有明顯差異，但ADIPOQ基因在杜洛克豬中表達量稍高于哥本哈根迷你豬；而本研究表明ADIPOQ基因在杜洛克豬中的表達量極顯著高于撒壩豬，可能是撒壩豬脂肪沉積能力高于哥本哈根迷你豬所導致的。

PPAR 信號通路可能是影響撒壩豬脂肪沉積能力的關鍵通路。過氧化物酶體增值激活受體是甾體激素受體超家族成員，通過與類色素x 受體相結合而轉變為異二聚體，它是調控脂肪細胞的分化、脂肪酸轉運、炎癥以及調節糖代謝等的相關基因[17-18]。PPARα、PPARβ和PPARγ是PPAR 的3 種不同亞型；其中PPARα和PPARβ主要在肝臟及脂肪酸氧化水平較高的組織中發揮作用，PPARγ主要作用于脂肪細胞，影響其分化和生長，也是I 型糖尿病藥物治療的靶點[19]。

4 結 論

本研究對撒壩豬和杜洛克豬肝臟組織的轉錄組分析，共鑒定出2 326 個DEGs，結合GO 分析和KEGG通路富集分析，本研究初步認為LPL、ACSL3、DIPOQ等9 個基因和PPAR 信號通路等16 條代謝通路在撒壩豬的肉質形成中發揮了重要作用，研究結果為揭示撒壩豬肉質和生長性狀的形成機制提供了重要的參考數據，并提高了對脂質沉積相關分子機制的認識。