綿羊不同發(fā)育階段背最長肌組織中可變剪接的鑒定與分析

鮑晶晶 浦亞斌 馬月輝 趙倩君

（中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193）

在真核生物轉(zhuǎn)錄過程中，pre-mRNA通過多種不同的方法切除內(nèi)含子，重新拼接外顯子產(chǎn)生成熟的mRNA，這個(gè)過程稱為可變剪接（Alternative splicing，AS）［1］，可變剪接是高等真核生物蛋白質(zhì)多樣性和遺傳多樣性的重要調(diào)控機(jī)制［2］。人基因組中95%以上的多外顯子基因能發(fā)生可變剪接［3］。單個(gè)基因通過可變剪接產(chǎn)生的mRNA異構(gòu)體的數(shù)量從2到數(shù)千個(gè)不等［4］?？勺兗艚拥幕绢愋桶ㄍ怙@子跳躍（Skipped exons，SE），5′端可變剪接位點(diǎn)（Alternative 5′ splice sites，A5SS），3′端 可 變 剪接位點(diǎn)（Alternative 3′splice sites，A3SS），互斥外顯子（Mutually exclusive exons，MXE）和內(nèi)含子保留（Retained introns，RI）［5］。通過可變剪接產(chǎn)生的mRNA和蛋白質(zhì)異構(gòu)體在結(jié)構(gòu)、穩(wěn)定性、亞細(xì)胞定位和功能上是不相同的，且一些基因的異?？勺兗艚訒?dǎo)致發(fā)育無法正常進(jìn)行［6］。Blue等［7］在剛出生和成年小鼠心臟中鑒定了由可變剪接產(chǎn)生的3個(gè)特異性基因包涵體重鏈、包涵體輕鏈-a和轉(zhuǎn)運(yùn)蛋白結(jié)合蛋白-1，并表明包涵體重鏈基因的剪接調(diào)控被抑制時(shí)，小鼠的體重和肌肉重量就會增加。PKM基因的一種剪接體SRSF3能調(diào)節(jié)結(jié)腸癌細(xì)胞的生長和腫瘤特異性能量代謝的維持［8］。調(diào)節(jié)可變剪接事件是心肌細(xì)胞表達(dá)功能細(xì)胞骨架和肌節(jié)蛋白的必要步驟。大約3%的特發(fā)性擴(kuò)張型心肌病患者存在RNA結(jié)合蛋白RBM20的突變，RBM20是可變剪接的組織特異性調(diào)節(jié)劑［9］。鋅指蛋白ZNF148的兩種剪接異構(gòu)體在大腸癌細(xì)胞凋亡和侵襲中發(fā)揮不同作用，ZNF148能促進(jìn)CRC細(xì)胞增殖、侵襲和遷移，ZNF148表現(xiàn)出相反的作用，由此得出ZNF148的兩種剪接異構(gòu)體可能對惡性生物學(xué)活性產(chǎn)生相互拮抗作用［10］。有研究表明約1/3的mRNA轉(zhuǎn)錄本會被無義衰變降解，以調(diào)節(jié)正?；虮磉_(dá)或靶向降解異常表達(dá)的基因［11］。這些研究表明AS形成的功能各異的轉(zhuǎn)錄本在肌肉組織發(fā)育及疾病發(fā)生過程中扮演著重要作用。

可變剪接事件在畜禽動物中的研究也陸續(xù)被報(bào)道。在木川烏骨雞中鑒定了TYR基因的5個(gè)可變剪接異構(gòu)體，并證明其在黑色素形成機(jī)制中扮演不同角色［12］。對小尾寒羊和陶塞特羊的卵巢組織中差異AS事件基因功能注釋到與繁殖力、發(fā)育相關(guān)的通路［13］。Agouti是一個(gè)與毛色相關(guān)的基因，在白色山羊中發(fā)現(xiàn)了5種可變剪接異構(gòu)體，而在黑色山羊中只發(fā)現(xiàn)了2種異構(gòu)體［14］。許多涉及性腺發(fā)育、激素代謝、晝夜節(jié)律的基因通過選擇性剪接被差異調(diào)控［13］。但對多浪綿羊肌肉不同發(fā)育過程中可變剪接的時(shí)序變化及其作用機(jī)制未見報(bào)道。本研究基于多浪綿羊不同發(fā)育階段背最長肌的RNA-seq測序數(shù)據(jù)，通過rMATS軟件［15］分析首次鑒定了不同發(fā)育階段背最長肌組織中剪接事件和差異剪接基因（Differential spilicing gene，DSG），通過GO和KEGG分析這些DSG的功能，旨為研究肌肉不同發(fā)育階段AS事件及其調(diào)控機(jī)制提供理論依據(jù)。

1 材料與方法

1.1 材料

本研究所用多浪綿羊來自中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所昌平試驗(yàn)基地。隨機(jī)選擇健康的妊娠90日齡的綿羊胎兒（F90）3只，出生后30日齡羔羊（L30）3只和3歲成年羊（A3Y）3只，屠宰后分別采集9只多浪綿羊的背最長肌組織立即保存于液氮中。

1.2 方法

1.2.1 組織總RNA的提取和質(zhì)量檢測 根據(jù)RNeasy? Plus Universal Mini Kit（Qiagen）試劑盒的操作說明進(jìn)行綿羊背最長肌組織總RNA的提取。用1%的瓊脂糖凝膠電泳檢測總RNA是否降解以及是否有污染，用Nanodrop ND-1000檢測RNA的純度，Agilent 2100 bioanalyzer精確檢測RNA的濃度和完整度（RIN），將RNA濃度≥100 ng/μL，總量≥3 μg，RIN值≥7.0的RNA樣品用于轉(zhuǎn)錄組建庫測序。

1.2.2 cDNA文庫的構(gòu)建和測序 9個(gè)綿羊背最長肌組織的cDNA文庫構(gòu)建和測序均由北京貝瑞和康生物技術(shù)有限公司完成。測序平臺是Illumina Hiseq2500 V4，測序模式為125PE。用NGSQC Toolkit v2.3.3軟件［16］過濾掉含有接頭的reads，過濾掉N含量大于10%的reads，過濾掉低質(zhì)量的reads，得到 clean reads。

1.2.3 RNA-seq數(shù)據(jù)的比對與組裝 從NCBI數(shù)據(jù)庫中下載綿羊的參考基因組（Oar_v4.0）以及基因組注釋GTF文件（http：//www.ncbi.nlm.nih.gov/genome/）。將clean reads 通過HISAT［17］軟件比對到綿羊的參考基因組上，用StringTie［18］軟件將比對上的reads進(jìn)行轉(zhuǎn)錄本的組裝并計(jì)算每個(gè)基因和isoform的表達(dá)水平，然后通過merge［19］將所有樣品的轉(zhuǎn)錄本組裝在一起再次呈遞給StringTie軟件重新估算轉(zhuǎn)錄本豐度，然后將數(shù)據(jù)載入ballgown包［20-21］進(jìn)行差異基因分析。

1.2.4 可變剪接的鑒定及差異分析 使用rMATS［15］軟件對每個(gè)樣品存在的可變剪接類型及相應(yīng)表達(dá)量進(jìn)行檢測。rMATS是一款可利用RNA-Seq數(shù)據(jù)自動檢測和分析差異可變剪接的軟件。rMATS軟件的統(tǒng)計(jì)模型是根據(jù)計(jì)算P值和錯(cuò)誤發(fā)生率（False discovery rate，F(xiàn)DR）來鑒定差異的可變剪接。本實(shí)驗(yàn)采用FDR＜0.05作為標(biāo)準(zhǔn)來篩選差異剪接基因。

1.2.5 差異剪接基因的功能富集分析 對差異剪接基因采用 g：profiler（https：//biit.cs.ut.ee/gprofiler/）在線分析軟件［22-23］進(jìn)行GO富集分析，采用DAVID（https：//david.ncifcrf.gov/）在線分析軟件［24-25］進(jìn)行KEGG富集分析，采用P＜0.05作為顯著富集標(biāo)準(zhǔn)。

2 結(jié)果

2.1 可變剪接事件的鑒定

RNA-seq原始數(shù)據(jù)經(jīng)過質(zhì)控后，每個(gè)背最長肌組織樣品均產(chǎn)生不低于10 GB的clean data，GC含量平均為48.6%且Q30＞87%，測序數(shù)據(jù)滿足后續(xù)試驗(yàn)要求。通過HISAT軟件將clean data比對到綿羊參考基因組上，9個(gè)樣品83.8%-88.5%的clean reads能比對到綿羊的參考基因組上。將經(jīng)過Stringtie軟件組裝的轉(zhuǎn)錄本與綿羊參考基因組進(jìn)行比較，綿羊妊娠90日齡胎兒時(shí)期（F90）、出生后羔羊時(shí)期（L30）和成年3歲時(shí)期分別鑒定出13 923、11 959和12 164個(gè)可變剪接事件（圖1）。其中SE占各個(gè)發(fā)育階段鑒定出的可變剪接類型的比例最高，達(dá)到69.19%-72.19%，A5SS、A3SS、MXE、RI分別占各個(gè)階段可變剪接事件總量的6.80%-7.75%，1.04%-1.19%，5.80%-6.12%，4.49%-5.18%。從圖中可以看出不同發(fā)育階段的綿羊背最長肌組織中的可變剪接事件的數(shù)量存在差異，且在胚胎時(shí)期的數(shù)量最多，出生后30 d和成年3歲時(shí)期的變化不大，這也與綿羊出生前后肌肉的發(fā)育相一致。

圖1 綿羊背最長肌組織不同發(fā)育階段可變剪接事件的數(shù)量

2.2 背最長肌組織的差異剪接基因鑒定

對rMATS軟件分析結(jié)果以FDR＜0.05為標(biāo)準(zhǔn)進(jìn)行差異篩選，共篩選到6 935個(gè)顯著的差異剪接基因（圖2），F(xiàn)90與L30比較組在SE、A5SS、A3SS、MXE、RI事件中鑒定到的差異剪接基因數(shù)分別為1 856，155，310，199，25；F90與A3Y比較組在SE、A5SS、A3SS、MXE、RI事件中鑒定到的差異剪接基因數(shù)分別為1 958，165，302，236，28；L30與 A3Y比 較 組 在 SE、A5SS、A3SS、MXE、RI事件中鑒定到的差異剪接基因數(shù)分別為1 216，120，184，152，29。差異剪接基因中SE事件的占比最高（約72.41%），RI的占比最少。

圖2 綿羊背最長肌組織不同發(fā)育階段差異可變剪接基因數(shù)量

2.3 差異剪接基因的功能富集分析

對獲得的差異剪接基因進(jìn)行GO富集分析，多浪綿羊胎兒組與羔羊組背最長肌組織差異剪接基因富集到的前20個(gè)GO terms（圖3），多浪綿羊胎兒組與成年羊組背最長肌組織差異剪接基因富集到的前20個(gè)GO terms（圖4），多浪綿羊羔羊組與成年羊組背最長肌組織差異剪接基因富集到的前20個(gè)GO terms（圖5）。3個(gè)比較組都富集到了與肌肉顯著相關(guān)的GO terms，如橫紋肌發(fā)育（GO：0014706）、肌肉結(jié)構(gòu)發(fā)育（GO：0061061）、肌細(xì)胞分化（GO：0042692）、肌膜（GO：0042383）、發(fā)育進(jìn)程（GO：0032502）。

圖3 綿羊胎兒組與羔羊組肌肉組織差異剪接基因的GO富集結(jié)果

對3個(gè)比較組的差異剪接基因的KEGG通路分析也顯著富集到與肌肉顯著相關(guān)的通路，見表1-表3，如胰島素信號通路、Wnt信號通路、MAPK信號通路、肌動蛋白細(xì)胞骨架調(diào)控。其中胞吞作用、胰島素信號通路、焦點(diǎn)粘連、泛素介導(dǎo)的蛋白水解是這3個(gè)比較組中都顯著富集到的共同信號通路。細(xì)胞的胞吞作用參與分子物質(zhì)運(yùn)輸轉(zhuǎn)運(yùn)；胰島素信號通路參與肌肉的發(fā)育；焦點(diǎn)粘連參與細(xì)胞膜受體和肌動蛋白骨架之間的結(jié)構(gòu)連接；泛素介導(dǎo)的蛋白水解調(diào)控細(xì)胞周期；肌動蛋白細(xì)胞骨架調(diào)控維持細(xì)胞形態(tài)和調(diào)控細(xì)胞粘連。

圖4 綿羊胎兒組與成年羊組肌肉組織差異剪接基因的GO富集結(jié)果

3 討論

1978年，Gilbert第一次發(fā)現(xiàn)可變剪接以來［26］，已證明可變剪接參與高等真核生物的基本生命調(diào)節(jié)過程并已成為真核基因調(diào)控領(lǐng)域的研究熱點(diǎn)?？勺兗艚邮寝D(zhuǎn)錄后基因表達(dá)調(diào)控的關(guān)鍵步驟，它通過使單個(gè)基因產(chǎn)生不同的RNA異構(gòu)體而有助于蛋白質(zhì)和功能的多樣性。可變剪接通過細(xì)胞和組織特異性調(diào)節(jié)細(xì)胞分化和機(jī)體發(fā)育，可變剪接的失調(diào)會導(dǎo)致人類各種遺傳疾病，包括早衰癥、強(qiáng)直性營養(yǎng)不良、額顳葉癡呆、囊性纖維化和癌癥［27］，并且可變剪接最近已成為預(yù)后和治療癌癥的標(biāo)志［28］。可變剪接通過控制特定細(xì)胞類型和特定時(shí)間來表達(dá)特定的RNA異構(gòu)體［29］。LPIN1基因有2種RNA剪接體，LPIN1α主要作用于前脂肪細(xì)胞的分化，而LPIN1β則會促進(jìn)前脂肪細(xì)胞增殖［30］。在陜北絨山羊睪丸組織中首次鑒定了CTNNB1的4個(gè)轉(zhuǎn)錄本異構(gòu)體，命名為Ctnnb1-A/-B/-C/-D，Ctnnb1-C在睪丸中表達(dá)最高，表明其在雄性生育中起重要作用［31］。NFIX是NFI家族的成員，在肌肉和大腦發(fā)育中發(fā)揮重要作用，有研究在奶山羊中鑒定并驗(yàn)證了NFIX基因的2個(gè)新轉(zhuǎn)錄本（NFIXa和NFIXb），正常NFIX轉(zhuǎn)錄本在肺中高表達(dá)，NFIXa異構(gòu)體在胰腺中高表達(dá)，而NFIXb異構(gòu)體在肺和胰腺中高表達(dá)。此外，NFIXa異構(gòu)體在肝臟、脾臟、脂肪、腸道和睪丸中的表達(dá)水平顯著高于正常NFIX和NFIXb異構(gòu)體，NFIXa異構(gòu)體在大腦中的表達(dá)顯著高于NFIXb異構(gòu)體［32］。說明可變剪接事件具有時(shí)空特異性和組織特異性。高通量轉(zhuǎn)錄組測序技術(shù)已成為研究RNA轉(zhuǎn)錄產(chǎn)物的可變剪接的可靠工具。冉茂良等［33］對豬60胚齡、90胚齡、30日齡和180日齡4個(gè)不同發(fā)育時(shí)期的睪丸進(jìn)行轉(zhuǎn)錄組測序鑒定出92 738個(gè)可變剪接事件，篩選到63個(gè)發(fā)生可變剪接的基因與睪丸素代謝密切相關(guān)。張敏等［34］利用RNA-seq技術(shù)在肉雞肌肉與脂肪組織中分別鑒定到5 966和6 757個(gè)AS事件，并檢測到513個(gè)參與脂肪和肌肉發(fā)育的顯著差異的剪接基因。郭家中等［35］在簡州大耳羊早期胎兒（妊娠45 d、60 d、105 d）到新出生（出生后3 d）階段背最長肌中鑒定了137 308個(gè)可變剪接事件，4個(gè)階段均發(fā)生可變剪接事件的基因經(jīng)GO富集分析表明顯著富集在細(xì)胞表達(dá)調(diào)控和物質(zhì)合成等過程，并且發(fā)現(xiàn)肌肉早期發(fā)育階段發(fā)生的剪接事件多與晚期時(shí)期。黃艷群等［36］研究發(fā)現(xiàn)與絲羽烏骨雞趾型顯著相關(guān)的是Lmbr1基因轉(zhuǎn)錄本Lmbr1-α。

圖5 綿羊羔羊組與成年羊組肌肉組織差異剪接基因的GO富集結(jié)果

表1 綿羊胎兒組與羔羊組肌肉組織差異剪接基因的KEGG通路分析

表2 綿羊胎兒組與成年羊組肌肉組織差異剪接基因的KEGG通路分析

表3 綿羊羔羊組與成年羊組肌肉組織差異剪接基因的KEGG通路分析

基因組中差異可變剪接會導(dǎo)致細(xì)胞組成和功能的差異。本研究采用rMATS軟件鑒定可變剪接事件及差異剪接基因，在妊娠90日齡與出生后30日齡、妊娠90日齡與成年3歲綿羊以及出生后30日齡與成年3歲綿羊背最長肌組織比較組中分別鑒定到了2 545，2 689和1 701個(gè)差異剪接基因。從研究結(jié)果可以看到出生前后背最長肌組織的可變剪接事件存在差異，由此可見可變剪接具有發(fā)育階段特異性，且可變剪接在出生前的行動更活躍，這是與哺乳動物出生前主要是肌纖維數(shù)量的增加，出生后則是肌纖維的增長和變粗有密切關(guān)系［37］。本研究可變剪接事件中SE事件的數(shù)量是最多的，這與前人的研究一致［38］，說明SE是生物體中最普遍發(fā)生的可變剪接事件。對差異剪接基因進(jìn)行功能富集分析，GO顯著富集到了與肌肉發(fā)育相關(guān)通路上且共同富集到細(xì)胞內(nèi)部分（GO：0044424），參與細(xì)胞內(nèi)物質(zhì)合成，KEGG共同富集到了胞吞作用、胰島素信號通路、焦點(diǎn)粘連、泛素介導(dǎo)的蛋白水解4條信號通路上，這些信號通路與肌肉發(fā)育密切相關(guān)，表明可變剪接在綿羊背最長肌的發(fā)育過程中發(fā)揮重要作用。差異可變剪接基因與差異表達(dá)基因不一樣，本研究中約26%的差異可變剪接基因的表達(dá)水平是差異的，說明差異可變剪接基因與差異表達(dá)基因各自發(fā)揮功能共同維持細(xì)胞生命活動［39］。

4 結(jié)論

多浪綿羊背最長肌組織在妊娠90日齡胎兒時(shí)期（F90）、出生后羔羊時(shí)期（L30）和成年3歲（A3Y）時(shí)期分別鑒定出13 923、11 959和12 164個(gè)可變剪接事件，其中SE占各個(gè)發(fā)育階段鑒定出的可變剪接類型的比例最高，約占70.69%，A5SS、A3SS、MXE和RI分別約占7.28%、11.18%、5.96%和4.84%。綿羊背最長肌組織在F90_vs_L30、F90_vs_A3Y和L30_vs_A3Y比較組中鑒定到差異剪接基因數(shù)分別為2 545、2 689和1 701個(gè)。本研究中3個(gè)比較組的差異剪接基因的GO和KEGG富集分析都顯著注釋到與肌肉發(fā)育顯著相關(guān)的通路上。