簡州大耳羊早期胎兒到新出生階段背最長肌中可變剪接的鑒定與分析

郭家中，陶海溪，李鵬飛，李　利，張紅平

(四川農業大學 動物科技學院，成都　611130)

可變剪接(alternative splicing，AS)是指真核生物體內多外顯子基因通過選擇不同的剪接位點，從而形成多個不同RNA轉錄本的生物學過程。可變剪接普遍存在于各種真核生物中。例如，基于人類大腦、肝、心臟等多個組織的轉錄組比較分析表明，90%以上的蛋白編碼基因發生過可變剪接[1-2]。在低等動物中，Gibilisco等[3]在4種果蠅的多個組織中發現，大約有20%～37%的多外顯子基因發生過可變剪接，而且不同性別的同一組織發生的可變剪接模式也有顯著不同。在植物中，Shen等[4]發現在大豆的生長發育過程中大約有63%的基因發生了可變剪接，而且早期發育階段發生的可變剪接多于晚期階段。此外，已有研究表明，在動物體內比例最高的可變剪接類型是外顯子跳躍[1,3,5]，而在植物體內比例最高的可變剪接類型則是內含子保留[4,6]，表明動植物體內可變剪接發生的分子機制具有差異。

已有研究認為，可變剪接形成的部分轉錄本并不具有功能，是生物體不精密調控產生的噪音[7-8]；但越來越多的證據表明，可變剪接形成的大量功能不同的成熟轉錄本[2,9]，在高等動物的組織發育和疾病的發生發展過程中扮演著重要作用。在肌肉發育方面，Singh等[10]通過RNA-Seq數據發現敲弱Rbfox2的表達后， Mef2d 和 Rock2基因的可變剪接體喪失活性，導致小鼠成肌細胞融合過程受阻，證明單個可變剪接因子所調控的可變剪接事件是小鼠成肌細胞融合過程的必要條件之一。最近，Du等[11]在強直性肌營養不良的小鼠骨骼肌中發現了由CUG重復引起的可變剪接紊亂。另外，可變剪接在動物心肌的發育過程中也扮演著重要作用。例如，ASF/SF2基因的敲除影響CaMKIIδ基因的轉錄本類型，從而導致小鼠生長過程中的細胞死亡[12]。

截至目前，盡管有關畜禽骨骼肌生長發育過程中的轉錄組變化已經被陸續報道[13-15]，但是對骨骼肌發育過程中可變剪接的時序變化特征及其作用機制鮮有研究，而近幾年迅猛發展的高通量測序數據為在全基因組水平上研究可變剪接提供了豐富的數據資源。本試驗采用轉錄組測序數據分析簡州大耳羊妊娠期第45天、60天和105天胎兒以及出生后第3天羔羊背最長肌中的可變剪接事件，從而為山羊骨骼肌生長發育的分子機制研究提供基礎數據和理論參考。

1　材料與方法

1.1　組織樣品的收集和高通量測序

本研究從四川簡州大耳羊原種場選取一批成年簡州大耳羊空懷母羊，進行同期發情處理并配種。根據配種記錄，分別采集母羊妊娠后第45天(E45-1、E45-2和E45-3)、第60天(E60-1和E60-2)、第105天(E105-1、E105-2和E105-3)胎兒和出生后第3天(B3-1、B3-2和B3-3)羔羊的背最長肌樣本(n=11)。分別提取總RNA，構建高通量測序文庫，采用HiSeq 2500測序儀進行讀段長度為125 bp的雙末端測序(SRR3567144)。

1.2　序列比對和轉錄本組裝

原始測序數據通過FastQC軟件[16]進行質量評估，使用rMATS(v3.2.2 beta)軟件[17]剪切尾部低質量序列后，獲得長度為120 bp的高質量讀段。采用STAR(v2.5)軟件[18]將高質量讀段比對到山羊參考基因組(https://www.ncbi.nlm.nih.gov/genome/10731?genome_assemb-ly_id=281266)[19]。利用Stringtie(v1.3.3)軟件[20]組裝轉錄本并評估基因表達量(默認參數)。

1.3　可變剪接的鑒定及基因富集分析

使用Astalavista(v4.0)軟件[21]對每個樣本中的可變剪接事件進行鑒定、分類和統計。根據說明，Astalavista 軟件能夠鑒定可變剪接外顯子跳躍(skipping exon,SE)、5′端可變剪接位點(alternative 5 splice site,A5SS)、3′端可變剪接位點(alternative 3 splice site,A3SS)、內含子保留(intron retention,RI)及其他類型。利用自編的Python程序提取發生可變剪接的基因名稱和數量，然后應用Gene Ontology Consortium(http://geneontology.org/)進行GO功能注釋和富集分析(數據庫版本日期2017-09-26)。

1.4　可變剪接因子基因表達量的分析

真核生物體內的可變剪接過程受到RNA結合蛋白(RBMs)、核內不均一核糖核蛋白 RNA聚合酶Ⅱ轉錄物(HNRNPs)和絲氨酸/精氨酸富集剪接因子(SRSFs)等多種類型轉錄因子的共同調控[22-24]。據此，本研究比較分析RBMs、HNRNPs和SRSFs3個家族的成員在簡州大耳羊背最長肌中的表達變化模式。

2　結果與分析

2.1　讀段的比對和轉錄本組裝

對4個發育階段的簡州大耳羊背最長肌測序原始讀段進行質量控制后，總共獲得526 497 414條高質量讀段。由表1可見，在各樣本中，88.18%～96.62%的高質量讀段成功比對到山羊參考基因組。其中，比對到基因組唯一位置的讀段所占比例為77.36%～83.33%。上述高深度的測序結果為山羊背最長肌發育中可變剪接的鑒定提供了可靠的數據基礎。

表1　11個文庫讀段比對結果Table 1　Mapping results of clean reads in eleven libraries

2.2　全基因組水平的可變剪接事件鑒定

采用Stringtie軟件進行轉錄本組裝，共獲得2 050 999條轉錄本。將上述轉錄本與山羊參考基因組進行比較，共鑒定出137 308個可變剪接事件，577種不同類型的可變剪接組合模式。其中，以跳躍外顯子、內含子保留、可變的3′端剪接位點、可變的5′端剪接位點4種基本類型的可變剪接數量最多。在各個階段鑒定出的可變剪接類型中，SE的比例最高，占各樣本總量的27.80%～34.44%；其次，A3SS、RI和A5SS分別為18.24%～20.78%，12.14%～22.10%和10.49%～11.21%(圖1)。值得注意的是，隨著生長發育的進行，RI類型的比例在持續降低，尤其是在出生后第3天的可變剪接數目較妊娠期顯著地降低(P<0.05)，檢測出的可變剪接總數也顯著減少(P<0.05)。另外，在妊娠期第45天和妊娠期第60天鑒定出可變剪接數量最多的基因是伴肌動蛋白(Nebulin，NEB)，其數量別為71和79；妊娠期第105天和出生后第3天鑒定出可變剪接數量最多的基因是遮掩蛋白(Obscurin，OBSCN)，其數量分別為131和218。

圖1　4個發育階段可變剪接事件的數量及比例Fig.1　Numbers and proportions of alternative splicing events at four developmental stages

基因組注釋結果顯示，鑒定出的可變剪接事件來自于4 562個已注釋基因，占參考基因組中已注釋基因總數(27 199)的16.77%(圖2)。妊娠期第45天、妊娠期第60天、妊娠期第105天和出生后第3天的可變剪接分別來自于2 779、3 388、3 006和2 127個已注釋基因。其中1 319個基因在4個發育階段均發生了可變剪接，表明在出生前的背最長肌發育過程中大量基因發生了可變剪接。各個階段發生可變剪接的基因比較結果顯示，妊娠期第60天的特異性基因最多(615個)，表明此階段轉錄表達較活躍；妊娠期第45天、第105天和出生后第3天依次為361、283和137個。

2.3　可變剪接基因的功能富集分析

功能注釋結果(圖3)表明，4個階段均發生可變剪接的1 319個基因分別顯著富集(P<0.05)在 mRNA加工(GO:0006397)、腺苷親核轉酯反應的RNA剪接(GO:0000377)、RNA調控(GO:0043484)等291個生物學過程。在分子功能方面，上述基因顯著富集在RNA結合(GO:0003723)、轉錄因子結合(GO:0008134)、鈣黏著蛋白結合(GO:0045296)等64個GO類別中。在細胞組成方面，上述基因顯著富集在核斑(GO:0016607)、收縮纖維(GO:0043292)、肌原纖維(GO:0030016)等95個GO類別。上述結果表明，各個發育階段的可變剪接事件主要參與到大分子代謝過程和基因表達調控等過程。

圖2　4個發育階段發生可變剪接的基因數量Fig.2　Numbers of AS genes at four developmental stages

對4個階段發生可變剪接的特異性基因進行GO功能注釋和富集分析，僅有妊娠期第45天和妊娠期第60天的基因富集結果顯著。由圖4可見，妊娠期第45天的特異性基因主要富集在5個GO類別中，包括生物學過程中的細胞組成結構調節(GO:0051128)；細胞組成中的細胞質基質(GO:0005829)、胞內細胞器(GO:0044446)、膜內膜結合細胞器(GO:0043231)；分子功能中的蛋白質結合(GO:0005515)。上述結果表明，妊娠期第45天的背最長肌中發生可變剪接的特異性基因主要功能為調節細胞內物質合成。

如圖5所示，妊娠期第60天特異性基因顯著富集在RNA聚合酶Ⅰ轉錄本延長(GO:0006362)、染色體組織(GO:0051276)、有絲分裂(GO:0000278)、細胞定位(GO:0006996)等24個GO生物學過程類別中。在細胞組成方面，上述基因顯著富集在中心體(GO:0005814)、含磷基團轉移酶復合物(GO:0061695)、細胞核(GO:0005730)、轉移酶復合物(GO:1990234)、微管組織中心(GO:0005815)、微管骨架(GO:0015630)等35個GO類別。在分子功能方面，上述基因顯著富集在嘌呤依賴性解旋酶活性(GO:0070035)、ATP-依賴性解旋酶活性(GO:0008026)、作用于RNA的催化活性(GO:0140098)、轉移酶活性(GO:0016740)、雜環化合物結合(GO:1901363)、有機環狀化合物結合(GO:0097159)等9個GO類別。

2.4　可變剪接因子的表達差異分析

轉錄組測序結果顯示，RBMs、HNRNPs和SRSFs3個家族中58個基因的平均表達量大于5 FPKM。對上述基因的表達量進行單因素方差分析，結果表明，包括 HNRNPA0、 RBM15B、 SRSF9等19個基因在4個階段中呈現顯著性差異表達(P<0.05)。 如圖6所示，15個可變剪接調控因子的表達量在發育過程中呈現降低趨勢。而 RBM7和 RBM48的表達量表現為在出生前逐漸升高，出生后急劇下降。

3　討 論

發生在轉錄后水平的可變剪接是形成高等動植物體內蛋白質多樣性和復雜性的重要機制之一。Wang等[1]報道，在人類大腦、肝臟、肌肉等多個組織中，大約有95%基因發生了可變剪接。盡管如此，在不同物種中發生可變剪接基因的比例具有差異。例如，相關研究表明，豬、牛、鴨基因組分別有80.6%[25]、21%[26]、46.12%[27]基因發生了可變剪接。本研究利用轉錄組測序數據，在妊娠期第45天、妊娠期第60天、妊娠期第105天和出生后第3天共4個發育階段的簡州大耳羊背最長肌共鑒定出137 308個可變剪接事件，發生可變剪接的基因占山羊注釋基因組總數的16.77%。該結果表明，可變剪接廣泛發生在山羊骨骼肌發育過程中；值得注意的是，上述比例低于人、豬和鴨中的結果。一方面，可能是因為不同組織中發生可變剪接的基因比例不同。例如，在人類各組織中，大腦發生可變剪接的基因比例最高，而肌肉、脂肪組織相對較低[1]。另一方面，山羊參考基因組不夠完善也是導致上述結果的一個重要原因。此外，本研究中外顯子跳躍類型的可變剪接比例最高，這與在其他動物中發現的結果相一致。

圖3　 4個發育階段共有可變剪接基因的GO功能富集分析Fig.3　GO Functional enrichment analysis of the overlapped As genes across four stages

圖4　 妊娠期第45天特異性基因GO功能富集分析Fig.4　GO Functional enrichment analysis of unique AS genes at day 45 of gestation

已有研究表明，哺乳動物骨骼肌肌纖維數量在其出生前已經穩定，出生后主要是肌纖維的增長和變粗，故出生前是哺乳動物骨骼肌發育關鍵階段[28]。例如，Greenwood等[29]和Nissen等[30]對豬早期胚胎和出生后胎兒肌肉的冷凍切片觀察發現，初級纖維和次級纖維的增殖在出生前就已經停止，出生后纖維只有長度和體積的變化。在本研究中，4個發育階段均發生可變剪接的基因進行GO功能富集分析后，細胞組成結果主要與肌肉發育相關，如肌原纖維、收縮纖維、肌節、A帶和I帶等，表明肌纖維在上述階段持續發育，可變剪接基因參與到肌肉發育過程中。Wilson等[31]的研究表明，綿羊脛骨前肌在胚胎期第32天出現原發性肌管，胚胎期第38天出現少量的次級肌管，發育至第62天次級肌管形成。本研究中，山羊胎兒背最長肌在妊娠期60天的可變剪接數量和可變剪接調控因子的表達量均高于其他階段，表明背最長肌發育進入活躍期，大量細胞在該階段增殖和分化。同時GO功能注釋和富集分析結果也顯示，可變剪接基因主要富集在細胞增殖和大分子合成相關過程，表明可變剪接在背最長肌肌纖維發育的分化調控過程中具有重要作用。

圖5　妊娠期第60天特異性基因GO功能富集分析Fig.5　GO Functional enrichment analysis of unique AS genes at day 60 of gestation

圖6　19個可變剪接因子基因在4個發育階段中的表達量熱圖Fig.6　Heatmap of expression profiles of 19 alternative splicing factors across four developmental stages

在高等動植物體內，可變剪接是一個受到RBMs、HNRNPs和SRSFs等多種因子協同調節的生物學過程[20-22]。例如， HNRNP A/B蛋白可以結合HIV-1病毒的3′端剪切位點和第3外顯子活化區域，從而抑制其某一種轉錄本的表達[32-33]。Motta-Mena等[34]發現HNRNP 蛋白可以直接抑制 CD45基因的第4外顯子，并且通過抑制鄰近的剪接增強子影響第5外顯子的表達。最近，Liu等[35]在 HNRNPA1基因敲除小鼠中證明， HNRNPA1對肌肉相關基因的表達和選擇性剪接起著至關重要且不可替代的作用。在本研究中， HNRNPA1、HNRNPL和 HNRNPR等基因的表達量都呈現出隨生長發育而逐漸降低的趨勢。該結果表明，在生長發育早期，可變剪接因子的高表達保證了背最長肌增殖分化所需要的豐富轉錄本[36]。綜上所述，可變剪接事件的時序變化特征與山羊背最長肌的發育過程密切相關[37]。

參考文獻Reference：

[1]WANG E T，SANDBERG R,LUO S,etal.Alternative isoform regulation in human tissue transcriptomes[J].Nature,2008, 456(7221):470.

[2]PAN Q,SHAI O,LEE L J,etal.Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing[J].NatureGenetics,2008,40(12):1413-1415.

[3]GIBILISCO L,ZHOU Q,MAHAJAN S,etal.The evolution of alternative splicing inDrosophila[J].Biorxiv,2016,2(21):054700.

[4]SHEN Y,ZHOU Z,WANG Z,etal.Global dissection of alternative splicing in paleopolyploid soybean[J].ThePlantCell,2014,26(3):996-1008.

[5]LIN L,PARK JW,RAMACHANDRAN S,etal.Transcriptome sequencing reveals aberrant alternative splicing in Huntington’s disease[J].HumanMolecularGenetics,2016,25(16):3454-3466.

[6]LI W,LIN W D,RAY Petal.Genome-wide detection of condition-sensitive alternative splicing inArabidopsisroots[J].PlantPhysiology,2013,162(3):1750-1763.

[7]PICKRELL J K,PAI A A,GILAD Y,etal.Noisy splicing drives mRNA isoform diversity in human cells[J].PLoSGenetics,2010,6(12):e1001236.

[8]TRESS M L,ABASCAL F,VALENCIA A.Alternative splicing may not be the key to proteome complexity[J].TrendsinBiochemicalSciences,2017,42(2):98-110.

[9]NILSEN T W,GRAVELEY B R.Expansion of the eukaryotic proteome by alternative splicing[J].Nature,2010,463(7280):457-463.

[10]SINGH R K,XIA Z,BLAND C S,etal.Rbfox2-coordinated alternative splicing of Mef2d and Rock2 controls myoblast fusion during myogenesis[J].MolecularCell,2014,55(4):592-603.

[11]DU H,CLINE M S,OSBORNE R J,etal.Aberrant alternative splicing and extracellular matrix gene expression in mouse models of myotonic dystrophy[J].NatureStructural&MolecularBiology,2010,17(2):187-193.

[12]XU X,YANG D,DING J H,etal.ASF/SF2-regulated CaMKIIδ alternative splicing temporally reprograms excitation-contraction coupling in cardiac muscle[J].Cell,2005,120(1):59-72.

[13]ZHAO X,MO D,LI A,etal.Comparative analyses by sequencing of transcriptomes during skeletal muscle development between pig breeds differing in muscle growth rate and fatness[J].PloSOne,2011,6(5):e19774.

[14]ZHAN S,ZHAO W,SONG T,etal.Dynamic transcriptomic analysis in hircine longissimus dorsi muscle from fetal to neonatal development stages[J].Functional&IntegrativeGenomics,2017(S1):1-12.

[15]WANG Y,ZHANG C,FANG X,etal.Identification and profiling of microRNAs and their target genes from developing Caprine skeletal muscle[J].PloSOne,2014,9(5):e96857.

[16]VESELOVSKA L,SMALLWOOD S A,SAADEH H,etal.Deep sequencing and de novo assembly of the mouse oocyte transcriptome define the contribution of transcription to the DNA methylation landscape[J].GenomeBiology,2015,16(1):209.

[17]SHEN S,PARK J W,LU Z X,etal.rMATS:robust and flexible detection of differential alternative splicing from replicate RNA-Seq data[J].ProceedingsoftheNationalAcademyofSciences,2014,111(51):E5593-E5601.

[18]DOBIN A,DAVIS C A,SCHLESINGER F,etal.STAR:ultrafast universal RNA-seq aligner[J].Bioinformatics,2013,29(1):15-21.

[19]BICKHART D M,ROSEN B D,KOREN S,etal.Single-molecule sequencing and chromatin conformation capture enable de novo reference assembly of the domestic goat genome[J].NatureGenetics,2017,49(4):643-650.

[20]PERTEA M,PERTEA G M,ANTONESCU C M,etal.StringTie enables improved reconstruction of a transcriptome from RNA-seq reads[J].NatureBiotechnology,2015,33(3):290-295.

[21]FOISSAC S,SAMMETH M.ASTALAVISTA:dynamic and flexible analysis of alternative splicing events in custom gene datasets[J].NucleicAcidsResearch,2007,35(2):W297-W299.

[22]MANLEY J L,KRAINER A R.A rational nomenclature for serine/arginine-rich protein splicing factors(SR proteins)[J].Genes&Development,2010,24(11):1073-1074.

[23]ELLIOTT D J,OGHENE K,MAKAROV G,etal.Dynamic changes in the subnuclear organisation of pre-mRNA splicing proteins and RBM during human germ cell development[J].JournalofCellScience,1998,111(9):1255-1265.

[24]HWANG B,LIM J H,HAHM B,etal.hnRNP L is required for the translation mediated by HCV IRES[J].BiochemicalandBiophysicalResearchCommunications,2009,378(3):584-588.

[25]冉茂良,陳斌,李智,等.基于RNA-seq測序數據鑒定和分析豬基因組可變剪接事件[J].中國科學(生命科學),2016(3):274-284.

RAN M L,CHEN B,LI ZH,etal.Identification and analysis of alternative splicing events in sus scrofa using RNA-Seq data[J].ScientiaSinica(Vitae),2016(3):274-284.

[26]CHACKO E R S.Genome-wide analysis of alternative splicing in cow:implications in bovine as a model for human diseases[J].BMCGenomics,2009,3(3):S11.

[27]徐鐵山,顧麗紅,侯水生,等.應用RNA-seq數據開展鴨基因組可變剪接的鑒定與分析[J].中國家禽,2016(17):10-16.

XU T SH,GU L H,HOU SH SHetal.Identification and analysis of alternative splicing in duck genome using RNA-seq data[J].ChinaPoultry,2016(17):10-16.

[28]PICARD B,LEFAUCHEUR L,BERRI C,etal.Muscle fibre ontogenesis in farm animal species[J].ReproductionNutritionDevelopment,2002,42(5):415-431.

[29]GREENWOOD P,HUNT A,HERMANSON J,etal.Effects of birth weight and postnatal nutrition on neonatal sheep:Ⅱ.Skeletal muscle growth and development[J].JournalofAnimalScience,2000,78(1):50-61.

[30]NISSEN P,DANIELSEN V,JORGENSEN P,etal.Increased maternal nutrition of sows has no beneficial effects on muscle fiber number or postnatal growth and has no impact on the meat quality of the offspring[J].JournalofAnimalScience,2003,81(12):3018-3027.

[31]WILSON S,MCEWAN J,SHEARD P,etal.Early stages of myogenesis in a large mammal:formation of successive generations of myotubes in sheep tibialis cranialis muscle[J].JournalofMuscleResearchandCellMotility,1992,13(5):534-550.

[32]BILODEAU P S,DOMSIC J K,MAYEDA A,etal.RNA splicing at human immunodeficiency virus type 1 3′ splice site A2 is regulated by binding of hnRNP A/B proteins to an exonic splicing silencer element[J].JournalofVirology,2001,75(18):8487-8497.

[33]ZHU J,MAYEDA A,KRAINER A R.Exon identity established through differential antagonism between exonic splicing silencer-bound hnRNP A1 and enhancer-bound SR proteins[J].MolecularCell,2001,8(6):1351-1361.

[34]MOTTA MENA L B,HEYD F,LYNCH K W.Context-dependent regulatory mechanism of the splicing factor hnRNP L[J].MolecularCell,2010,37(2):223-234.

[35]LIU T Y,CHEN Y C,JONG Y J,etal.Muscle developmental defects in heterogeneous nuclear Ribonucleoprotein A1 knockout mice[J].OpenBiology,2017,7(1):160303.

[36]GABUT M,S AMAVARCHI-TEHRANI P,WANG X,etal.An alternative splicing switch regulates embryonic stem cell pluripotency and reprogramming[J].Cell,2011,147(1):132-146.

[37]KALSOTRA A,COOPER T A.Functional consequences of developmentally regulated alternative splicing[J].NatureReviewsGenetics,2011,12(10):715-729.