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

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

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

1.2.4 可變剪接的鑒定及差異分析 使用rMATS［15］軟件對每個樣品存在的可變剪接類型及相應表達量進行檢測。rMATS是一款可利用RNA-Seq數據自動檢測和分析差異可變剪接的軟件。rMATS軟件的統計模型是根據計算P值和錯誤發生率（False discovery rate，FDR）來鑒定差異的可變剪接。本實驗采用FDR＜0.05作為標準來篩選差異剪接基因。

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

2 結果

2.1 可變剪接事件的鑒定

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

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

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

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

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

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

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

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

對3個比較組的差異剪接基因的KEGG通路分析也顯著富集到與肌肉顯著相關的通路，見表1-表3，如胰島素信號通路、Wnt信號通路、MAPK信號通路、肌動蛋白細胞骨架調控。其中胞吞作用、胰島素信號通路、焦點粘連、泛素介導的蛋白水解是這3個比較組中都顯著富集到的共同信號通路。細胞的胞吞作用參與分子物質運輸轉運；胰島素信號通路參與肌肉的發育；焦點粘連參與細胞膜受體和肌動蛋白骨架之間的結構連接；泛素介導的蛋白水解調控細胞周期；肌動蛋白細胞骨架調控維持細胞形態和調控細胞粘連。

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

3 討論

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

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

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

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

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

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

4 結論

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