多浪羊與小尾寒羊皮下脂肪組織轉錄組分析

劉天義，馮 卉，Salsabeel Yousuf，解領麗，苗向陽

(中國農業科學院北京畜牧獸醫研究所，北京 100193)

綿羊是全世界主要的肉、奶、皮毛的牲畜資源。羊肉是我國日常膳食結構的重要組成部分。隨著人們生活水平的提高，人們對于羊肉品質的要求也隨之提高，而脂肪含量對肉品質具有重要影響[1-2]。研究脂肪組織的脂肪沉積機制不僅可以為人們提供高品質肉產品，還可以提高養殖業飼料轉化率，避免浪費。同時通過對脂肪組織的研究還能為由肥胖引起的相關疾病的預防和治療提供一定依據[3-5]。

脂肪組織遍布哺乳動物全身，主要分為白色脂肪組織和棕色脂肪組織兩類[6-7]。研究表明，脂肪合成與分解代謝過程受到C/EBPs、PPARγ、BMP4等轉錄因子[8-10]，以及PPAR、MAPK、TGFβ、Insulin等信號通路的調控[11-12]。Xu等[13]研究發現，SCD、FAS、FAP4等基因在脂肪酸代謝過程中發揮作用，其中，SCD1是SCD的一種亞型，可作為脂肪沉積的潛在標志物，SCD1是SREBP1調控的下游基因，當SCD1的啟動子上的固醇調節原件與SREBP1c轉錄因子結合后上調轉錄，脂質合成增加[14]。RNA是遺傳信息的重要載體，mRNA可以編碼蛋白實現遺傳信息的傳遞，lncRNA和miRNA均可以靶向mRNA來調控基因的表達[15-17]。轉錄組測序技術由于具有高靈敏度，已經被廣泛應用。Miao等[18]通過轉錄組測序手段，對多塞特羊與小尾寒羊脂肪組織進行轉錄組測序，以更好地了解兩品種綿羊脂肪代謝的生物學機制。Huang等[19]通過對和牛與荷斯坦牛皮下脂肪組織進行轉錄組測序，篩選得到差異表達基因，進一步對兩品種牛脂肪沉積的分子機制進行解析。為更深入了解mRNA在綿羊脂肪組織中的作用，本研究選用脂肪沉積存在差異的多浪羊與小尾寒羊皮下脂肪組織為試驗材料，其中多浪羊脂臀較大，屬于肉脂兼用型綿羊，小尾寒羊屬于短瘦尾型綿羊，肉用性能較差，皮下脂肪較少，而脂肪主要分布在內臟周圍，以此來適應惡劣的環境條件，二者為研究脂肪的沉積機制提供了良好的試驗材料。目前轉錄組測序技術較為成熟，以往雖有通過轉錄組技術來對脂肪組織進行的研究報道，但對于多浪羊與小尾寒羊皮下脂肪組織沉積的生物學差異還有待于深入研究。

本研究通過對多浪羊與小尾寒羊皮下脂肪組織進行轉錄組測序分析，得到兩品種羊皮下脂肪組織中差異表達的mRNA，并從中篩選出與脂肪沉積相關的候選基因。

1 材料與方法

1.1 試驗動物和樣品制備

本試驗以脂肪沉積存在差異的多浪羊和小尾寒羊為研究對象，選取健康無病、體況良好，種內個體體重相近(約50 kg)的雌性成年多浪羊和小尾寒羊各3只。采集背最長肌皮下脂肪組織為試驗材料，迅速將皮下脂肪組織放入液氮凍存。該過程使用已消毒器械，減少RNA的降解。試驗設置多浪羊皮下脂肪組織和小尾寒羊皮下脂肪組織兩組樣本，每組樣本3個重復，分別為D-PF-1、D-PF-2、D-PF-3、X-PF-1、X-PF-2、X-PF-3。

1.2 方法

1.2.1 脂肪組織RNA的提取及質控 取適量的脂肪組織，采用TRIzol(Invitrogen)法提取脂肪組織中的總RNA。并使用DNase I (RNase-free，TaKaRa)去除DNA。利用1%瓊脂糖凝膠電泳與ND-2000(NanoDrop Technologies)對所提RNA的濃度、純度以及完整性進行初步檢測，再用2100 Bioanalyzer(Agilent)進行檢測，記錄OD260 nm/OD280 nm以及RIN值。

1.2.2 脂肪組織cDNA文庫構建和RNA測序 將質檢合格的總RNA進行文庫構建。采用TruSeqTMstranded total RNA Kit(Illumina，San Diego，CA)試劑盒。首先從2 μg總RNA中除去rRNA；加入Fragmentation Buffer將RNA片段化；采用SuperScript double-stranded cDNA synthesis kit (Invitrogen,CA)試劑盒，加入六堿基隨機引物(Illumina)，利用隨機引物，以mRNA為模板反轉錄合成第一鏈cDNA，隨后合成第二鏈cDNA，形成穩定的雙鏈結構；加入End Repair Mix將粘性末端補成平末端，隨后在3′末端加上一個A堿基，用于Y字形的接頭，擴增前用UNG酶消化cDNA第二鏈；cDNA 經過 PCR 富集(sample preparation Kit(Illumina，San Diego,CA)后，DNA clean beads篩選 200～300 bp 的條 帶。經 TBS380(Picogreen)定量后，文庫使用 Illumina NovaSeq 6000 測序平臺進行高通量測序。

1.2.3 脂肪組織mRNA測序分析 為保證后續分析的準確性，需要使用fastp[20]軟件(https://github.com/OpenGene/fastp)將原始數據(raw data)中包含的測序接頭序列、低質量讀段、不確定堿基信息率較高序列及長度過短序列進行過濾，從而得到高質量的測序數據(clean data)。隨后使用Hisat2[21](https://daehwankimlab.github.io/hisat2/)將篩選后的數據定位到參考基因組，獲取用于后續分析的有效數據，同時將對比的結果用分析軟件RSeQC-2.3.6[22-23]進行質量評估。在已有的參考基因組基礎上使用軟件stringTie[24-25](https://ccb.jhu.edu/software/stringtie/)對每個樣本進行組裝拼接，與已有的轉錄本進行比較，獲得沒有注釋的新轉錄本。用軟件RSEM[26-27](http://deweylab.github.io/RSEM/)進行mRNA的表達量計算，再用軟件DESeq2篩選樣本之間差異表達的mRNA，將條件設定為|Fold change|≥2，Padjust≤0.05進行篩選。

1.2.4 脂肪組織差異表達mRNA功能富集分析 應用Cytoscape中的Clue-GO[28]軟件將差異表達的mRNA進行GO功能注釋和KEGG通路富集分析，GO是基因本體論聯合會建立的數據庫(geneontology，http://geneontology.org/)，分為細胞組分(cellular component，CC)，分子功能(molecular function，MF)和生物學過程(biological process，BP)。KEGG[29]京都基因和基因組百科全書(https://www.genome.jp/kegg/)是系統分析基因功能、聯系基因組信息和功能信息的數據庫，可將基因按照參與的pathway通路或行使的功能分類。將集中的基因顯示在KEGG通路圖上，展示其參與的KEGG注釋通路圖。利用多重檢驗校正方法BH(Benjamini and Hochberg)對P值進行校正，當P值(Padjust)小于等于0.05時則顯著富集。

1.2.5 脂肪組織mRNA蛋白網絡互作分析 將GO和KEGG注釋得到的信息，結合文獻篩選出與脂肪發育相關的條目和信息通路，選擇其中涉及到的基因，利用Cytoscape[30]軟件制作網絡圖，找出網絡中處于節點位置的比較重要的基因，在文獻中找出這些基因具體的功能以及在脂肪組織發育中的作用。

1.2.6 qRT-PCR驗證 采用實時熒光定量PCR驗證基因的表達水平。利用PerfectStartTMGreen qPCR SuperMix試劑盒在LightCycler480 Ⅱ Real-time PCR Instrument(Roche，Swiss)上配置PCR反應體系：2×PerfectStartTMGreen qPCR SuperMix 5 μL，10 μmol·L-1上、下游引物各 0.2 μL，cDNA 1 μL，Nuclease-free H2O 3.6 μL。PCR程序：94 ℃ 30 s；94 ℃ 5 s，60 ℃ 30 s，45個循環。循環結束后利用熔解曲線檢測產物特異性：從60 ℃緩慢升溫至97 ℃，每1 ℃采集5次熒光信號。利用2-ΔΔCt法進行樣本間表達量的計算[31]，用t-檢驗統計分析相對表達量，數據結果表示為“平均數±標準差(Mean±SD)”，P<0.05為差異顯著。對應的引物序列見表1。

表1 基因以及對應的引物Table 1 Gene and the corresponding primer sequences

2 結 果

2.1 多浪羊與小尾寒羊皮下脂肪組織測序數據分析

將原始測序數據總條目數(raw reads)中質檢不合格的條目去除后得到質控后測序數據的總條目數(clean reads)，多浪羊與小尾寒羊皮下脂肪組織clean reads分別為135 393 598、140 587 412、142 895 188、136 746 686、128 433 640和130 999 204。將這些reads與參考基因組對比得到有效數據(mapped reads),以及在參考序列上有多個比對位置的clean reads數目(multiple mapped)和唯一比對位置的clean reads數目(unique mapped)見表2。Mapping率高于94%，說明參考基因組注釋完整，且相關試驗不存在污染情況，為后續數據分析做了很好的鋪墊。

表2 6個樣本比對參考基因組結果Table 2 The result of 6 samples after mapping to the reference genome

2.2 已知與新mRNA分析及表達量計算

將數據進行已知和新mRNA篩選，結果顯示，在多浪羊和小尾寒羊皮下脂肪組織樣本中已知和新mRNA的數量分別為38 672和1 606個。每組設置3個生物學重復計算mRNA的表達量，選取|Fold change|≥2，Padjust≤0.05作為篩選差異表達mRNA的閾值。發現兩組差異表達的mRNA共839個(已知mRNA 828個，新mRNA 11個)，如圖1所示，其中在多浪羊組中上調和下調的差異表達基因數目分別為320和519個。

圖1 差異表達基因火山圖Fig.1 Volcano map of differentially expressed genes

2.3 脂肪組織差異表達mRNA的功能富集分析

對差異表達的mRNA進行GO和KEGG分析。GO注釋結果顯示，在生物學過程中，差異表達的mRNA主要富集于脂質分解代謝過程、脂質生物合成過程、脂質分解代謝負調控過程、MAPK級聯反應調控、低密度脂蛋白顆粒介導的信號傳導、對甘油三酯的反應、脂質代謝過程、MAPK級聯反應負調控等條目；在分子功能方面，主要富集于細胞外基質結構成分、心磷脂結合、甘油酸激酶活性、磷脂酰甘油結合等條目；在細胞組分，主要富集于細胞外基質、脂蛋白顆粒、血漿脂蛋白顆粒等條目(圖2)。在KEGG通路富集中也富集到與脂代謝相關的通路，如PI3K-Akt、MAPK、胰島素以及PPAR等信號通路(圖3)。

圖2 差異表達基因GO富集分析Fig.2 GO enrichment analysis of the differentially expressed genes

圖3 差異表達基因KEGG通路富集Fig.3 KEGG pathway enrichment analysis of the differentially expressed genes

2.4 脂肪組織差異表達基因蛋白互作網絡圖的構建

構建與脂質代謝相關差異表達基因蛋白互作網絡，能夠更直觀的展示多浪羊和小尾寒羊與脂肪沉積相關的基因間的調控關系(圖4)。本研究以網絡中各個節點的聯通性(Degree)為依據，篩選出關鍵節點[32-33]。Degree值越大，節點越大，表明該節點越關鍵。在蛋白互作網絡中發現，COL1A1、AKT3、AKT2、COL1A2、COL3A1等基因編碼的蛋白處于關鍵節點位置，可能在皮下脂肪代謝和成脂分化中具有重要調控作用。

2.5 qRT-PCR驗證

為了驗證測序結果，隨機選擇6個差異表達的mRNAs進行qRT-PCR驗證，結果如圖5所示，MGST3、COL1A1、AKT3在多浪羊皮下脂肪組織中下調，PCK1、PPP2R5A、LOC101113583在多浪羊皮下脂肪組織中顯著上調。以上結果與測序結果一致，表明測序結果可靠。

圖5 差異表達基因qRT-PCR驗證Fig.5 qRT-PCR validation for differentially expressed genes

3 討 論

mRNA是指導蛋白質合成的模板，在生物學過程中發揮著重要作用[34]，可為培育優質畜禽品種以及研究和預防與脂代謝相關的疾病提供依據。本研究對多浪羊與小尾寒羊皮下脂肪組織中的mRNA進行測序和分析，共鑒定到38 672個mRNAs，其中差異表達的mRNAs共839個，在多浪羊中有320個上調表達，519個下調表達。通過qRT-PCR對MGST3、COL1A1、AKT3、PCK1、PPP2R5A、LOC101113583進行了驗證，與測序結果一致，表明測序結果可靠。

有研究表明，細胞外基質可以影響脂肪的形成，而膠原是脂肪細胞外基質的主要成分，對脂肪的形成具有一定的調節作用[35]，推測研究編碼膠原的基因COL1A1對于脂肪的沉積具有重要意義。在秦川牛中的研究表明，COL1A1的表達量與肉質存在一定的相關性[36]。COL1A1在高脂飲食的肥胖人皮下脂肪組織中和高脂飼喂小鼠肌腱內的表達量均降低[37-38]，與本研究趨勢一致，表明COL1A1可能參與了脂肪代謝的調控。本研究中，AKT2顯著富集于PI3K/Akt信號通路中。PI3K/Akt是胰島素信號傳導的主要通路，通過調控胰島素水平參與糖脂代謝，可以抑制脂類分解。AKT2是胰島素信號通路的關鍵激酶，僅在胰島素靶組織中高表達，包括脂肪、肝、肌肉等組織，而胰島素可以調節脂肪細胞的分化，即推測AKT2基因可能調控了脂肪細胞的分化。敲除AKT2基因的小鼠表現出脂肪組織減少、胰島素抵抗和輕度糖尿病等現象[39]，本研究中，AKT2在多浪羊皮下脂肪組織中高表達，表明AKT2基因對多浪羊脂類細胞分解具有促進作用。綜上所述，COL1A1與AKT2不僅在脂肪沉積過程中起到了調控作用，而且還可能成為治療糖尿病等相關代謝疾病的潛在靶基因。

節點的大小代表連通性大小，白色表示下調，灰色表示上調The size of the node represents the level of Degree,white means down-regulation and gray means up-regulation圖4 差異表達基因蛋白互作網絡圖Fig.4 Network diagram of differentially expressed gene proteins interaction

通過差異表達基因GO注釋和KEGG富集分析表明，差異表達基因主要參與了PPAR信號通路、MAPK信號通路、脂肪細胞因子信號通路、脂質代謝過程等與脂類代謝相關的生物學過程。本研究鑒定到SCD、LPL、PCK1、PPP2R5A等基因與脂類代謝、脂類生物合成過程以及胰島素信號通路、PPAR信號通路過程相關。PPAR信號通路是與脂肪酸的生物合成和代謝以及成脂分化相關的通路[40-41]，本研究中，PPAR信號通路顯著富集了10個差異表達基因，其中SCD、LPL和PCK1在多浪羊皮下脂肪組織中上調表達，SCD是一種內質網結合酶，是參與單不飽和脂肪酸生物合成的關鍵酶，并且具有多種亞型。由于PPARγ2對于甘油三酯合成以及脂肪酸吸收的相關基因具有調控作用，所以在脂質代謝過程中具有重要作用，SCD1可以作為其靶基因在脂代謝過程中發揮作用。SCD1的表達水平受類固醇激素、甲狀腺激素、共軛亞油酸等影響。LPL是脂肪分解的關鍵限速酶，有利于脂肪的沉積[42-45]，并且主要在轉錄后蛋白水平起作用。LPL作為脂肪細胞的特征蛋白質之一，在脂肪細胞分化過程中具有重要作用，可以充當前體脂肪細胞向成熟脂肪細胞分化過程的標志。本研究中，LPL在多浪羊皮下脂肪組織中上調表達，與前人研究結果一致，推測LPL同樣對綿羊脂肪的生成具有促進作用。PCK1是調節糖異生和甘油生成的主要酶之一[46]，與肥胖相關，能通過GTP催化草酰乙酸(oxaloacetic acid，OAA)生成磷酸烯醇式丙酮酸(phosphoenolpyruvate，PEP)、GDP以及二氧化碳，在甘油生成以及糖異生過程中均可見到。在白色脂肪組織中，PCK1可以調節游離脂肪酸酯化生成甘油三酯，PPARγ是其啟動子中的結合位點，與野生型小鼠相比，PPAR缺失小鼠患有嚴重胰島素抵抗[47]，還有研究表明，敲除小鼠PCK1基因會導致糖尿病以及肥胖等[48]。因此PCK1對于維持脂質代謝至關重要。

PPP2R5A是蛋白磷酸酶2A的一個亞單位，與細胞增殖分化相關,可以通過調節CDK1/2和CHK2的活性促進有絲分裂，并且可以調節脂質激酶的活性[49-50]。本研究中，PPP2R5A在多浪羊皮下脂肪組織中上調，GO注釋發現其主要參與了細胞過程和脂質代謝過程的調節。PPP2R5A顯著富集于PI3K-Akt信號通路，這是細胞內重要的通路，可以調控細胞生長增殖、分化等過程，因此推測PPP2R5A可能在脂肪細胞增殖分化等過程中發揮作用。MGST3與LOC101113583在脂肪沉積的研究中較為少見，但在本研究中顯著富集到脂質生物合成與代謝生物學過程中，對其功能的研究還需進一步驗證。

本研究對多浪羊與小尾寒羊皮下脂肪組織中差異表達的mRNA進行了分析，結果表明，這些差異表達基因可能參與脂代謝過程，并且從中篩選到影響脂肪沉積的候選基因以及相關通路。但對于篩選到的候選mRNA在生物學功能以及調控的分子網絡還需要進一步探索。

4 結 論

本研究鑒定分析了多浪羊和小尾寒羊皮下脂肪組織中mRNA的表達情況，發現差異表達的mRNAs共839個，在多浪羊皮下脂肪組織中320個上調，519個下調；通過對差異表達基因進行GO和KEGG富集分析，篩選出COL1A1、AKT2、SCD、LPL、PCK1與PPP2R5A等與綿羊脂質代謝相關的基因，這些結果可為進一步理解多浪羊與小尾寒羊皮下脂肪沉積與脂代謝差異機制提供理論基礎。