肉雞肌肉與脂肪組織基因組差異剪接基因分析

張 敏，王 杰，孫艷發(fā)*，鐘禮鵬，李 焰，許衛(wèi)華，李虹儀，張 茂，黃其春

(1. 龍巖學院生命科學學院，龍巖 364000； 2. 中國農業(yè)科學院北京畜牧獸醫(yī)研究所，北京 100193； 3. 福建省家畜傳染病防治與生物技術重點實驗室，龍巖 364000)

可變剪接(alternative splicing，AS)是mRNA通過不同的方式產生多種不同剪接RNA轉錄產物的過程[1]，是高等真核生物蛋白質多樣性的重要基因調解機制之一[2]。人基因組中有約95%的多外顯子基因通過100 000多種不同的AS事件產生轉錄產物[3]。新一代RNA測序(RNA-seq)技術已經證實RNA轉錄產物的可變剪接是人基因組中編碼蛋白質多樣性的重要機制[4]。通過對基因組中AS事件的鑒定，已經發(fā)現(xiàn)一些基因的異常AS可導致人類多種疾病[5-7]。

畜禽基因和基因組中一些AS事件逐漸被發(fā)現(xiàn)和證實。研究人員通過設計跨越Lmbr1基因外顯子4的引物，發(fā)現(xiàn)與絲羽烏骨雞趾型呈現(xiàn)顯著相關的Lmbr1基因中Lmbr1-α為其主要轉錄產物，而Lmbr1-β為次要轉錄產物[8]。經重疊延伸RT-PCR結合巢式PCR法擴增，在東北民豬中發(fā)現(xiàn)TLR4基因的3個可變剪接體[9]，鴨TLR4基因中也發(fā)現(xiàn)了兩種不同的剪接體[10]。研究人員通過PCR擴增和5′RACE 方法獲得白色和黑色山羊皮膚組織Agouti基因部分基因組序列及mRNA序列，經過拼接比對，在白色山羊中發(fā)現(xiàn)5種Agouti基因可變剪接體，而黑色山羊中發(fā)現(xiàn)2種可變剪接體[11]。運用RNA-seq技術，畜禽基因組的AS事件得到了全面的鑒定和分析。通過對絨山羊肌內脂肪細胞在成熟前后RNA-seq數(shù)據(jù)分析發(fā)現(xiàn)，絨山羊肌內脂肪細胞中存在A3SS、A5SS、SE和RI 4種主要可變剪接類型，其中A3SS最為普遍[12]。牛基因組的21 755個基因中4 567(21%)個基因發(fā)生了可變剪接[13]。豬睪丸基因組中存在92 738個AS事件，以TSS、TTS、SE和AE 4種類型為主；通過分析發(fā)現(xiàn)，與睪丸素代謝相關的基因中63個基因發(fā)生了可變剪接，以TSS和TTS為主[14]。這些AS事件的鑒定與分析有利于全面研究基因組中AS事件和蛋白質的多樣性。由于AS事件大多數(shù)表現(xiàn)出組織特異性變異，其中10%～30%的AS事件易受細胞、組織和條件特異性調節(jié)的影響[3]。上述的研究并未對組織間差異剪接，特別是差異剪接基因(differential splicing gene，DSG)進行深入的分析。差異剪接分析描述了兩個樣本間AS位點使用情況和剪接異構體的差異，這對基因組中AS及其調控機制的研究至關重要，同時揭示了差異基因表達分析錯過的功能多樣性[15]。DSG指不同組織樣品間因差異剪接而導致基因亞型的相對豐度有差異的基因。因此，DSG分析對于研究組織間AS及其調控機制具有重要的意義。

本研究以AA肉雞RNA-seq測序數(shù)據(jù)為基礎，通過分析鑒定肉雞肌肉和脂肪組織基因組中DSG，分析這些DSG功能和其信號通路，為進一步解析肉雞基因組中組織間AS事件及組織間基因表達的特異性奠定理論基礎。

1 材料與方法

1.1 試驗樣品

AA肉雞母雞苗由福建正大食品股份有限公司提供，于龍巖學院動物飼養(yǎng)房飼養(yǎng)至42日齡。肉雞飼養(yǎng)管理和屠宰參考肉雞生產技術規(guī)范(DB11/T 328-2005)與肉雞屠宰質量管理規(guī)范(NY/T 1174-2006)。3只肉雞禁食12 h后進行屠宰，取胸肌和腹脂組織樣品轉入2 mL RNase-free離心管中，標記后迅速置于液氮中保存?zhèn)溆谩?/p>

1.2 組織總RNA提取

采用TRIzol試劑(Invitrogen，Carlsbad，美國)提取組織樣品總RNA，使用0.8%的瓊脂糖凝膠電泳檢測總RNA是否降解和受污染，并通過Agilent 2100 生物分析儀使用RNA 6000 Nano試劑盒(Agilent Technologies，美國)檢測總RNA的濃度、完整度(RIN)、28S/18S比值和片段大小。用于構建RNA-seq文庫和測序的總RNA均滿足以下條件：總RNA量≥6 μg，RIN≥8.0，28S/18S≥1.0。

1.3 肉雞胸肌和腹部脂肪組織RNA-seq文庫的構建和測序

RNA-seq文庫構建和測序均由深圳華大基因科技有限公司完成。分別取3只肉雞的胸肌和腹脂組織RNA樣品各3份，構建6個文庫，使用Illumina Hiseq 2500測序平臺進行測序，測序獲得長度為150 bp的序列(reads)。去除帶接頭的reads、未知堿基N含量大于5%的reads、去除低質量的reads(質量值低于 15 的堿基占該reads總堿基數(shù)的比例大于 20% 的reads為低質量的reads)，得到clean reads。

1.4 肉雞胸肌和腹部脂肪組織參考基因組比對及AS事件和DSG檢測

得到clean reads之后，使用HISAT軟件[16]將clean reads比對到雞參考基因組序列。使用StringTie軟件[17]對每個樣品進行轉錄本重構，然后用Cuffmerge軟件[18]將所有樣品的重構信息整合在一起，再使用Cuffcompare軟件將整合后的轉錄本與參考注釋信息進行比較，計算轉錄本數(shù)量。

使用rMATS軟件[19]檢測樣品中5種AS事件，分別為外顯子跳躍(skipped exon，SE)、可變供體位點(alternative 5′ splicing site，A5SS)、可變受體位點(alternative 3′ splicing site，A3SS)、互斥外顯子(mutually exclusive exons，MXE)和單內含子保留(retained intron，RI)，并檢測肌肉和脂肪組織樣品間的DSG。rMATS軟件是一款用于RNA-seq數(shù)據(jù)檢測DSG的軟件，它可以計算inclusion亞型和skipping亞型的相對豐度(圖1)，同時計算相應的P-value和FDR值來衡量樣品間剪接差異的顯著性。本研究中使用FDR≤0.05為判斷標準，小于該閾值的DSG定義為顯著差異的DSG。

I. 內含子異構體讀長; S. 跳躍異構體讀長；j. 連接長度；e1，e2. 外顯子長度；r. 讀長；lI. 內含子有效長度；ls. 外顯子異構體有效長度I. Reads of the inclusion isoform; S. Reads of the skipping isoform; j. Junction length; e1，e2. Exon length; r. Read length; lI. Effective length of the inclusion isoform; ls. Effective length of the skipping exon isoform圖1 AS事件中各種剪接類型的相對豐度計算Fig.1 Calculation of relative abundance of different splicing types in AS

1.5 肉雞胸肌和腹部脂肪組織DSG的GO功能富集和KEGG分析

使用KOBAS 3.0軟件[20-21]進行DSG的注釋、GO功能富集分析、信號通路(pathway)分析，并對相關信號通路和基因進行進一步分析。進行GO功能富集分析時，以基因組為背景，統(tǒng)計每個GO分析中基因與背景基因之間的差異，對P值進行校正，定義校正后的P<0.05為這些基因的功能在該GO類別中顯著富集。進行KEGG分析時，以基因組為背景，根據(jù)注釋結果以及功能分類，將DSG進行生物通路分類，校正后的P<0.05時，表明DSG富集在一條通路中。同時對通路中DSG的表達量進行檢測。具體方法為：與參考基因組比對之后，使用RSEM軟件[22]計算基因和轉錄本的表達水平，采用校正的P≤0.05來篩選肌肉和脂肪組織基因的差異表達基因。

1.6 數(shù)據(jù)統(tǒng)計分析

使用IBM SPSS Statistics 19.0軟件(SPSS Inc., Chicago, IL, USA)對肌肉和脂肪組織基因組中的可變剪接數(shù)量進行獨立樣本T檢驗和多重比較(SLD法)，數(shù)據(jù)結果以“平均數(shù)±標準差”表示。

2 結 果

2.1 肉雞肌肉和脂肪組織AS事件鑒定

測序原始數(shù)據(jù)經過質量控制后，肉雞胸肌3個組織樣分別獲得6.69、6.71和6.75 G，腹脂3個組織樣分別獲得6.68、6.69、6.73 G可用于注釋分析的序列片段(clean reads)。得到clean reads之后，使用HISAT軟件[16]將clean reads比對到雞參考基因組序列。肌肉組織樣品比對率分別為69.18%、69.49%和68.47%，脂肪組織樣品對比率分別為74.60%、73.10%和73.75%。相同組織樣品間均勻的比對率表明樣品之間的數(shù)據(jù)具有可比性。

肉雞肌肉和脂肪組織基因組中AS事件鑒定結果見表1。由表1可知，肉雞肌肉和脂肪組織基因組中分別檢測到(5 966.00±111.66)和(6 757.00±156.51)個AS，脂肪組織中AS數(shù)量比肌肉組織中高13.26%(P<0.05)。共檢測到的5種AS類型中，除MXE外，脂肪組織中4種類型剪接事件(SE、A5SS、A3SS和RI)的數(shù)量均顯著高于肌肉組織(P<0.05)；且以SE為主，分別占肌肉和脂肪組織中總AS事件數(shù)量的54.92%和52.67%。

表1肉雞肌肉和脂肪組織基因組中AS事件鑒定

Table1Theresultofalternativesplicing(AS)eventsidentifiedinmuscleandfattissuesgenomeinbroilers

項目Item肌肉組織Muscle tissue脂肪組織Fat tissueP值P-valueAS事件總數(shù)Total number of AS 5 966.00±111.666 757.00±156.510.002AS類型AS typeSE3 276.33±65.623 558.67±72.110.008A5SS579.67±20.23700.67±15.950.001A3SS989.00±13.861 205.33±39.460.001MXE276.33±11.06265.67±18.230.435RI844.67±22.501 026.67±35.300.002

2.2 肉雞肌肉和脂肪組織DSG鑒定

經rMATS軟件[19]檢測分析發(fā)現(xiàn)，AA肉雞肌肉和脂肪組織中檢測到513個顯著的DSGs，其中SE、A5SS、A3SS、MXE和RI 5種類型的DSG數(shù)量分別為314、34、17、54和94個(圖2)。DSG中發(fā)生SE可變剪接的基因數(shù)量最多；其次為RI和MEX；發(fā)生這3種類型的DSG占總AS事件數(shù)量的90.06%。

2.3 肉雞肌肉和脂肪組織DSG的GO注釋

為了進一步了解肉雞肌肉和脂肪組織DSG的功能，使用KOBAS 3.0在線軟件[20-21]對所有DSG基因進行了基因注釋、基因本體(gene ontology，GO)注釋和功能分析。選取校正P≤ 0.05的DSG作為GO富集分析結果(圖3)。通過GO富集分析發(fā)現(xiàn)，肉雞肌肉和脂肪組織中513個DSGs中的93個DSGs顯著富集在細胞組分(cellular component)和分子功能(molecular function)類別中的11個GO項內，同一基因多富集在多個GO項中。影響肉雞肌肉和脂肪組織DSG涉及的細胞組分包括細胞(GO：0005623)、細胞部分(GO：0044464)、細胞內(GO：0005622)、細胞外(GO：0044424)、細胞器(GO：0043226)、胞內細胞器(GO：0043229)、細胞質(GO：0005737)、胞漿部分(GO：0044444)8個GO項；還有少量DSG分類到分子功能中，包括蛋白質結合(GO：0005515)、磷蛋白結合(GO：0051219)和蛋白質磷酸化氨基酸結合(GO：0045309)3個GO項。

圖2 肉雞肌肉和脂肪組織差異可變剪接基因數(shù)量Fig.2 The number of DSG between muscle and fat tissues in broilers

圖3 肉雞肌肉和脂肪組織基因組DSG的GO富集結果Fig.3 The result of GO analysis for DSG in genome between muscle and fat tissues in broilers

2.4 肉雞肌肉和脂肪組織DSG的通路和相關基因分析

以雞基因組為參照，使用KOBAS 3.0在線軟件[20-21]對肉雞肌肉和脂肪組織DSG進行了通路分析。分析發(fā)現(xiàn)，31個DSGs富集在肌動蛋白細胞骨架調節(jié)、焦點粘連、丙酮酸代謝和細胞內吞作用4條信號通路中。肌動蛋白細胞骨架調節(jié)與焦點粘連兩條信號通路關系密切，其中10個基因均在兩條信號通路之中(表2)。肌動蛋白細胞骨架在維持細胞形狀、調控細胞粘連等方面具有重要功能；焦點粘連則參與細胞膜受體和肌動蛋白骨架之間的結構連接；丙酮酸代謝作為三羧酸循環(huán)的重要組成部分，在調解物質代謝中具有重要作用；細胞內吞作用在分子運輸轉運方面具有重要作用。

為了驗證DSG在肌肉和脂肪組織中的表達情況，以脂肪組織為參照，經RSEM軟件[22]分析，信號通路中的31個DSGs有13個在肌肉組織中上調表達，4個下調表達；其中14個DSGs(占信號通路基因中48.39%)表達量在肌肉和脂肪組織中未發(fā)生顯著的表達量變化(表3)。

表2肉雞肌肉和脂肪組織基因組中DSGs的通路分析結果

Table2TheresultofpathwayanalysisforDSGsingenomebetweenmuscleandfattissuesinbroilers

代謝通路PathwayDSG數(shù)量Number of DSGsKEGG ID校正的P值Corrected P-valueDSG名稱DSG name肌動蛋白細胞骨架調節(jié)Regulation of actin cytoskeleton13048100.00CRK,ENAH,FGFR1,ITGA3,ITGA7,ITGB1,MRAS,MYL9,P38C-CRK,PPP1R12A,PPP1R12B,PXN,ROCK2焦點粘連Focal adhesion11045100.02CCND3,CRK,ITGA3,ITGA7,ITGB1,MYL9,P38C-CRK,PPP1R12A,PPP1R12B,PXN,ROCK2丙酮酸代謝Pyruvate metabolism5006200.03ACSS2,HAGH,LDHD,PDHA1,PKM細胞內吞作用Endocytosis12041440.04AP2B1,CAPZB,CCDC53,CLTCL1,CYTH1,EPN3,FAM21A,GIT2,H3KBP1,SPG21,USP8,VPS29

表3肉雞肌肉和脂肪組織基因組中DSGs差異表達結果

Table3TheresultofdifferentialexpressionlevelsforDSGsingenomebetweenmuscleandfattissuesinbroilers

序號N基因名稱Gene name基因IDGene ID差異表達倍數(shù)Fold change校正的P值Corrected P-value上/下調Up/Down-regulation1ACSS2419158-0.411.59×10-38Down2AP2B1417525-2.420.01*3CAPZB3964180.535.78×10-28Up4CCDC53418094-5.210.42*5CCND34199280.862.26×10-13Up6CLTCL14167650.364.86×10-53Up7CRK107054794-2.370.08*8CYTH1422085-1.132.78×10-9*9ENAH374180-246.060.99*10EPN31008592900.304.38×10-46Up11FAM21A423772-0.739.97×10-14Down12FGFR1396516-61.850.96*13GIT23740350.992.78×10-5Up14HAGH4165370.456.14×10-37Up15ITGA3373946-0.501.69×10-21Down16ITGA71017490980.414.32×10-21Up17ITGB1374058-0.994.90×10-8Down18LDHD4156890.473.38×10-10Up19MRAS3951490.720.00Up20MYL9396215-1.103.46×10-6*21P38C-CRK417553-1.580.60*

(轉下頁 Carried forward)

使用RSEM軟件計算基因和轉錄本的表達水平，采用校正的P≤0.05來篩選肌肉和脂肪組織的差異表達基因；Up表示基因在肌肉組織中上調，脂肪組織中下調，Down表示基因在肌肉組織中下調，在脂肪組織中上調；*表示基因表達量差異不顯著

The expression level of genes and transcripts were calculated using RSEM software. Genes with correctedP-value ≤ 0.05 were considered DEGs. Up indicted gene was up-regulated expression in muscle tissue and down-regulated expression in fat tissue; Down indicted gene was down-regulated expression in muscle tissue and up-regulated expression in fat tissue; * indicted the genes expression level were not different between muscle ant fat tissues

3 討 論

自1978年Gilbert[23]發(fā)現(xiàn)AS現(xiàn)象以來，Ensembl數(shù)據(jù)庫[24]已經收錄了人、小鼠、斑馬魚、雞、鴨、豬、牛、羊等80多種生物的可變剪接事件?；诟咄縍NA-seq技術挖掘基因組AS的應用[25]，使得人、小鼠、豬、牛、果蠅等多個物種基因組的可變剪接事件得到了全面的鑒定和分析[14]。冉茂良等[14]研究發(fā)現(xiàn)，AS事件具有發(fā)育時期特異性。徐鐵山等[26]在北京鴨肌肉和皮脂中分別發(fā)現(xiàn)了5 958.00和7 528.67個AS。本研究中，在肉雞肌肉和脂肪組織基因組中分別檢測到5 966.00和6 757.00個AS，其中脂肪組織中AS數(shù)量顯著高于肌肉組織。由此可見，可變剪接事情具有組織特異性，脂肪組織發(fā)生AS事件的頻率遠遠高于肌肉組織。本研究發(fā)現(xiàn)，肌肉和脂肪組織中DSG發(fā)生SE可變剪接的基因數(shù)量最多。SE是選擇性剪接最簡單和最普遍存在的一種方式，這與Li等[27]對豬肌肉組織中差異剪接類型的研究相一致，說明SE可變剪接類型在肉雞肌肉和脂肪組織中發(fā)生頻率較高。

差異剪接分析對于理解細胞和器官的生理過程十分重要[28]。為了探索肌肉和脂肪組織生理過程的差異，本研究進行了肉雞肌肉和脂肪組織間差異剪接分析。肉雞肌肉和腹脂組織基本組成差異較大，前者由肌纖維構成，后者則屬于白色脂肪，由脂肪細胞構成[29]。本研究肌肉和脂肪組織間的DSG的功能大部分與細胞基本結構相關。由此可見，差異剪接在調控肌肉和脂肪組織細胞的基本組成上具有重要作用。信號通路分析發(fā)現(xiàn)，DSG富集在肌動蛋白細胞骨架調節(jié)、焦點粘連、丙酮酸代謝和細胞內吞作用4條信號通路中，這些通路與肌肉組成、脂肪代謝等密切相關。由此可見，肌肉和脂肪組織細胞組成和行使功能的差異，部分源于基因組中差異可變剪接?？勺兗艚诱{控的基因通常不同于差異表達所調控的基因。因此，在基因組水平描述基因調控時，差異剪接與差異基因表達通?；ハ嘌a充[15]。本研究發(fā)現(xiàn)，約50%的DSG在肌肉和脂肪組織中的表達水平未發(fā)生顯著變化，說明肌肉和脂肪組織中DSG和差異表達基因發(fā)揮各自的功能，共同調控肌肉和脂肪組織的基本組成、生理功能等。

對于DSG可以通過剪接因子對整個基因組的可變剪接進行研究，也可以針對某一基因進行具體的結構和功能的研究。剪接因子在mRNA前體剪接過程中發(fā)揮重要作用。研究發(fā)現(xiàn)，可以通過剪接因子SRSF10對雞基因組的剪接作用來研究整個基因組中的可變剪接[30]。在與本試驗相關的研究中發(fā)現(xiàn)，對脂肪代謝具有重要調控功能的Apelin信號通路對肉雞腹脂的沉積具有重要的調控作用(未發(fā)表)。其中Apelin前體由 77 個氨基酸組成，編碼基因通過可變剪接可指導合成多種成熟的Apelin活性多肽片段：Apelin-36、Apelin-17、Pyr-apelin-13、Apelin-13和Apelin-12等[31]。因此，對Apelin基因的可變剪接、時空表達規(guī)律以及功能驗證等研究，能夠加深對肉雞中該基因DSG、基因功能的深入理解。

4 結 論

4.1本研究在肉雞肌肉和脂肪組織中分別檢測到(5 966.00±111.66)和(6 757.00±156.51)個可變剪接事件，脂肪組織發(fā)生可變剪接的頻率顯著高于肌肉組織(P=0.002)。

4.2本研究鑒定的5種可變剪接類型中以外顯子跳躍為主，分別占肌肉和脂肪組織中可變剪接事件的54.92%和52.67%。

4.3本研究在肌肉和脂肪組織基因組中檢測到513個顯著的DSGs，其中93個DSGs顯著富集于細胞組分和分子功能中；31個DSGs富集在肌動蛋白細胞骨架調節(jié)、焦點粘連、丙酮酸代謝和細胞內吞作用信號通路中，14個DSGs在肌肉和脂肪組織基因組中的表達水平未發(fā)生顯著變化。

4.4本研究對肉雞肌肉和脂肪組織基因組中DSG的鑒定與分析，為研究可變剪接及其調控的機制、組織間差異剪接的多樣性提供理論依據(jù)。

注：DSGs的詳細信息需要者可向作者索取。