基于毛竹重測(cè)序的基因可變剪接研究

王炯亮 趙韓生 高志民

（國(guó)際竹藤中心北京100102）

竹子是生長(zhǎng)在溫帶和熱帶森林中的多年生草本植物，其纖維素和半纖維素含量與木本植物相媲美。竹子在全基因組層面的研究有限，僅在2013年首次報(bào)道了毛竹基因組草圖 （2.05 Gb），Scaffold N50為328 Kb，并預(yù)測(cè)出31 987個(gè)基因。得益于測(cè)序技術(shù)和分析方法的改進(jìn)，獲得高精確度、連續(xù)性更好的染色體水平毛竹基因組將會(huì)促進(jìn)竹子基因功能和進(jìn)化的研究。

在真核生物中，可變剪接 （Alternative splicing，AS）是由有限的基因增加蛋白質(zhì)復(fù)雜性和多樣性的主要機(jī)制。經(jīng)預(yù)測(cè)，人類中超過95%的多外顯子基因可以表達(dá)生成多個(gè)不同的剪接體 （isoform），而且據(jù)報(bào)道，在擬南芥、大豆、二穗短柄草、棉花、玉米、水稻中分別約有 61%、52%、42%、40%、40%、33%的基因存在可變剪接事件。單個(gè)基因可產(chǎn)生不同的剪接體，是基因具有功能可塑性的主要原因，可能在植物生長(zhǎng)發(fā)育、脅迫反應(yīng)、信號(hào)轉(zhuǎn)導(dǎo)和開花時(shí)間等生物過程中扮演著重要角色。在一定程度上，物種特有的可變剪接將使有限蛋白質(zhì)編碼的基因具有功能多樣性。然而，可變剪接調(diào)控植物進(jìn)化的機(jī)制尚不清楚；可變剪接在保守程度不同的基因中的特性同樣仍不清楚。

基因組質(zhì)量和轉(zhuǎn)錄組測(cè)序深度對(duì)可變剪接的鑒定有著至關(guān)重要的作用。由于毛竹第一版基因組組裝不完整、且片段分散，轉(zhuǎn)錄組的樣本組織數(shù)量少覆蓋率低，所以很難得到完整的可變剪接圖譜，全面鑒定毛竹可變剪接需要高質(zhì)量的基因組和更多的轉(zhuǎn)錄組數(shù)據(jù)。因此，本研究利用從中國(guó)6個(gè)主要竹產(chǎn)區(qū)采集的毛竹 （Phyllostachys edulis）26個(gè)不同組織樣本測(cè)序獲得的轉(zhuǎn)錄組數(shù)據(jù)，進(jìn)行全基因組分析，最終獲得毛竹的可變剪接圖譜。轉(zhuǎn)錄組數(shù)據(jù)由Illumina和Pacific Biosciences（PacBio） 平臺(tái)測(cè)序產(chǎn)生。本研究鑒定得到大量可變剪接基因和可變剪接事件，并對(duì)它們進(jìn)行了分類。另外，通過全基因組范圍的研究，確定了氨基酸保守性與可變剪接之間的關(guān)系，并研究木質(zhì)素生物合成相關(guān)基因可變剪接的進(jìn)化規(guī)律。本研究不僅為進(jìn)一步研究竹子基因功能和調(diào)控網(wǎng)絡(luò)提供了全面的可變剪接圖譜，而且從進(jìn)化角度揭示了可變剪接的作用。

1 材料方法與數(shù)據(jù)分析

1.1 樣本采集

為得到完整的可變剪接圖譜，于2015年春從中國(guó)6個(gè)主要竹產(chǎn)區(qū)采集毛竹樣本進(jìn)行測(cè)序。6個(gè)采集地分別是：江蘇宜興 （31°15′08.41″N， 119°43′42.55″E，海拔212 m）、 浙江天目山 （30°19′13.42″N，119°26′55.21″E，海拔 480 m）、湖北咸寧 （29°81′10.02″N， 114°31′21.12″E， 海拔 150 m）、 湖南桃江 （28°28′39.74″N， 112°11′18.62″E， 海拔 320 m）、 廣西桂林（28°28′39.74″N， 112°11′18.62″E， 海拔 216 m）、 貴州赤水 （28°28′15.27″N，105°59′41.43″E， 海拔 120 m）。收集了處于不同生長(zhǎng)階段毛竹的26個(gè)組織，包括地下莖、根、筍、葉、籜片和芽。有關(guān)生物樣本采集的詳細(xì)信息參見原文附表S19。

1.2 基因組測(cè)序、組裝、注釋

使用WGS和Hi-C策略對(duì)毛竹基因組進(jìn)行組裝，并使用前人的研究方法進(jìn)行注釋。詳細(xì)描述見Protocols.io。Hi-C庫(kù)是按照前文所描述的流程進(jìn)行準(zhǔn)備，具體說明見原文附件。

1.3 RNA的分離、Illumina RNA-seq建庫(kù)與測(cè)序、數(shù)據(jù)分析

使用前人的方法對(duì)RNA進(jìn)行分離、純化、濃度測(cè)定、反轉(zhuǎn)錄和cDNA文庫(kù)的構(gòu)建。所有cDNA文庫(kù)的構(gòu)建和標(biāo)準(zhǔn)化過程信息詳見原文附件。在質(zhì)控之后，對(duì)文庫(kù)進(jìn)行光學(xué)檢驗(yàn)，隨后使用HiSeq-4000平臺(tái) （雙端測(cè)序，插入長(zhǎng)度為150 bp）進(jìn)行測(cè)序。然后，計(jì)算了序列的質(zhì)量，并使用FastQC（軟件版本0.11.3）在默認(rèn)參數(shù)下對(duì)低質(zhì)量序列進(jìn)行過濾。最后，使用RNA-SeQC（軟件版本1.1.8）在默認(rèn)參數(shù)下計(jì)算得到RNA-seq中的關(guān)鍵統(tǒng)計(jì)指標(biāo)。

使用AStalavista（軟件版本 4.0） （AStalavista，RRID：SCR 001815）在默認(rèn)參數(shù)下鑒別可變剪接基因和可變剪接事件。隨后，分析比較了主要的4類可變剪接IR （intron retention）、A3SS （alternative 3′splice site donor）、 A5SS （alternative 5′splice site acceptor） 和 ES （exon skipping）。 使用 Ontologizer（軟件版本 2.0） 和從Gene Ontology數(shù)據(jù)庫(kù) （GO，RRID：SCR_002811）得到的注釋信息對(duì)基因進(jìn)行富集分析。最后，基于表達(dá)水平 （每百萬個(gè)比對(duì)上的序列里面每千個(gè)堿基中的片段總數(shù)）計(jì)算了基因在每個(gè)樣本中的組織特異性值 （Ts）。

1.4 Iso-Seq文庫(kù)的構(gòu)建、測(cè)序與數(shù)據(jù)分析

按照測(cè)序儀制造商 （PacBio）提供的方法，對(duì)Iso-Seq文庫(kù)進(jìn)行構(gòu)建和測(cè)序。首先根據(jù)生物信息學(xué)預(yù)測(cè)得到的轉(zhuǎn)錄本長(zhǎng)度分布 （原文附表S22）構(gòu)建了3類SMRTBell文庫(kù) （轉(zhuǎn)錄本長(zhǎng)度在1 Kb到2 Kb之間的3個(gè)文庫(kù)，轉(zhuǎn)錄本長(zhǎng)度為2 Kb到3 Kb之間的2個(gè)文庫(kù)，轉(zhuǎn)錄本長(zhǎng)度大于3 kb的4個(gè)文庫(kù)），然后使用PacBio測(cè)序平臺(tái)對(duì)這9個(gè)SMRT文庫(kù)進(jìn)行測(cè)序。

通過分析PacBio RSII的測(cè)序數(shù)據(jù)，得到了合理的全長(zhǎng)轉(zhuǎn)錄本異構(gòu)體，隨后合并不同文庫(kù)得到的轉(zhuǎn)錄本異構(gòu)體，并在對(duì)序列進(jìn)行插入、分類、聚類等處理后，去除冗余的數(shù)據(jù)，最終得到了合理的轉(zhuǎn)錄本。組裝得到的轉(zhuǎn)錄本使用PASA（軟件版本2.0.2） （PASA，RRID：SCR_014656） 在默認(rèn)參數(shù)下比對(duì)到參考基因組上。最后，與短片段序列一樣，使用AStalavista在默認(rèn)參數(shù)下對(duì)得到的gtf文件進(jìn)行分析，從而鑒別可變剪接。

1.5 進(jìn)化分析

依據(jù)前人的研究方法，進(jìn)行基因家族的鑒定、系統(tǒng)發(fā)育樹的構(gòu)建和分化時(shí)間的預(yù)測(cè)。在原文附件和protocol.io提供了詳細(xì)的信息。

1.6 木質(zhì)素生物合成相關(guān)基因的全基因組鑒定與分析

從ENSEMBL數(shù)據(jù)庫(kù) （Ensembl，RRID：SCR_002344）下載了擬南芥 （基因組版本 10）、小麥（基因組版本3.1）、水稻 （基因組版本 7.0）、毛果楊 （基因組版本JGI 2.0.31）和高粱 （基因組版本3.1）5個(gè)物種基因組的序列，收集得到了經(jīng)過實(shí)驗(yàn)驗(yàn)證的140個(gè)木質(zhì)素生物合成相關(guān)的基因 （原文附表S28）。通過局部比對(duì)搜索 （BLAST） （美國(guó)國(guó)家生物技術(shù)信息中心 [NCBI]BLAST，RRID：SCR_004870）和結(jié)構(gòu)分析2種方法鑒定毛竹中木質(zhì)素生物合成相關(guān)的基因。BLAST搜索閾值設(shè)置為：E-value＜1e-10，identity＞40%，coverage＞95%。然后使用 hmm-search（軟件版本 3.1b2）和 Pfam-A.hmm數(shù)據(jù)庫(kù) （2017.5.31發(fā)布的版本）對(duì)所有由BLAST搜索得到的序列進(jìn)行分析和人工矯正后，進(jìn)行系統(tǒng)發(fā)育分析。最后，使用yn00（PAML軟件包中預(yù)測(cè)同義替換率和非同義替換率的子軟件包）對(duì)木質(zhì)素生物合成相關(guān)的13個(gè)基因家族進(jìn)行同義替換率的計(jì)算。Ks率通過公式 T＝Ks/2r（r＝6.5×10-9）計(jì)算分化時(shí)間。

對(duì)木質(zhì)素生物合成相關(guān)基因的編碼序列正選擇作用分析：首先使用PROBCONS（軟件版本1.12）（ProbCons，RRID：SCR_011813）對(duì)每個(gè)家族的蛋白質(zhì)序列進(jìn)行比對(duì) （迭代次數(shù)設(shè)為1 000，其他參數(shù)默認(rèn)）。隨后將比對(duì)后的蛋白質(zhì)翻譯為相應(yīng)的編碼序列。接著使用Gblocks（軟件版本0.91b）軟件從比對(duì)的結(jié)果提取保守序列塊，并使用jModelTest（軟件版本2.1.6）依據(jù)貝葉斯信息準(zhǔn)則分析保守域得到最佳的核酸替換矩陣。然后，使用PhyML（軟件版本 3.0） （PhyML， RRID： SCR_014629） （1 000次檢驗(yàn)）選擇最佳的核酸替換矩陣重新構(gòu)建了系統(tǒng)發(fā)育樹。最后，使用 PAML（軟件版本 4.8） 的branch-site模型對(duì)系統(tǒng)發(fā)育樹中的一些分支進(jìn)行正選擇作用分析 （詳細(xì)的信息見protocols.io）。

2 結(jié)果與分析

2.1 數(shù)據(jù)描述

經(jīng)不同測(cè)序方法產(chǎn)生用于毛竹基因組組裝的基因組數(shù)據(jù)總共約有603.3 Gb。首先，全基因組序列（Whole-genome sequence，WGS）由新得到的約154 Gb數(shù)據(jù)和已有的約220 Gb數(shù)據(jù)進(jìn)行組裝。之后使用由157 Gb原始數(shù)據(jù)經(jīng)質(zhì)量控制得到的約17.58 Gb Hi-C有效序列進(jìn)行Hi-C輔助組裝 （原文附表S1）。此外，由Illumina和PacBio平臺(tái)分別產(chǎn)生的約379 Gb和5 Gb數(shù)據(jù)用于轉(zhuǎn)錄組分析 （原文附表S2-S7）。本研究通過使用染色體水平的參考基因組和大量轉(zhuǎn)錄組數(shù)據(jù)，在毛竹中鑒定出了25 225個(gè)可變剪接基因和266 711個(gè)可變剪接事件。

2.2 數(shù)據(jù)分析

2.2.1 毛竹染色體水平基因組的組裝和基因注釋

為提高毛竹基因組的質(zhì)量，本研究對(duì)61個(gè)文庫(kù)進(jìn)行測(cè)序 （原文附表S1），總共得到讀長(zhǎng)在76 bp到250 bp之間的約603.3 Gb的基因組數(shù)據(jù)。隨后為獲得高質(zhì)量的基因組，使用了不同的組裝策略。首先，WGS組裝得到的基因組大小達(dá)到1.91 Gb，其Contig N50和Scaffold N50的長(zhǎng)度分別達(dá)到55 Kb和894 Kb（原文附表S8）。由WGS組裝得到的新基因組的質(zhì)量和各統(tǒng)計(jì)量與之前的版本相比都有明顯的提高（原文附表 S9，S10），比如 Contig N50和 Scaffold N50分別增加172%和358%，且未知序列減少43%。其次，使用Hi-C數(shù)據(jù)進(jìn)行輔助組裝，得到了總長(zhǎng)度為1.91 Gb的基因組，其Contig N50和Scaffold N50長(zhǎng)度分別達(dá)到53.29 Kb和79.90 Mb（原文圖1A，1B）。由WGS組裝得到的約93.17%Scaffold可以錨定到24條染色體上 （原文附表S10），而且Scaffold N50的長(zhǎng)度增加了約89倍 （原文表1）。根據(jù)關(guān)聯(lián)圖和組裝結(jié)果，可以明顯區(qū)分開24條染色體。將毛竹的染色體比對(duì)到水稻基因組上面，發(fā)現(xiàn)平均覆蓋率約為59.77% （原文附圖S2，附表S11）。此外，使用人工細(xì)菌染色體 （BAC）序列、全長(zhǎng)cDNA和一些已知的毛竹基因序列對(duì)毛竹染色體水平的基因組進(jìn)行評(píng)估 （原文附圖S3，附表S12-14），發(fā)現(xiàn)染色體水平基因組較第一版基因組的覆蓋范圍更廣，且準(zhǔn)確性更高。

在注釋重復(fù)序列后，染色體水平基因組將更有利于后續(xù)的基因注釋 （原文附表S15）?；诖罅康霓D(zhuǎn)錄組 （原文附表S16）、全長(zhǎng)cDNA和同源蛋白數(shù)據(jù)，預(yù)測(cè)出了51 074個(gè)具有完整結(jié)構(gòu)的蛋白編碼基因 （原文附表S17），這些基因的內(nèi)含子和外顯子平均長(zhǎng)度分別為668 bp和284 bp（原文附表S18）。通過單分子實(shí)時(shí)測(cè)序數(shù)據(jù)和人工校驗(yàn)對(duì)不合理的注釋進(jìn)行校正，通過添加非編碼區(qū)域 （UTR）注釋和內(nèi)部結(jié)構(gòu)調(diào)整，改進(jìn)了約17%的基因模型 （原文附表S19）。基因組注釋完整性評(píng)估 （單拷貝同源物分析）顯示：毛竹 （95.2%）的注釋比玉米（92.2%）更完整，與水稻 （95.6%）接近 （原文圖1D，附表S20）。與前一版本的注釋相比，本版本中97.23%的基因模型在公共數(shù)據(jù)庫(kù)被鑒定，這有助于精確檢測(cè)可變剪接事件 （原文附表S21）。原文附表S22-S24和原文附圖S3-S9提供了基因模型預(yù)測(cè)和基因組進(jìn)化研究相關(guān)的詳細(xì)信息。此外，最新版本的基因組和基因注釋已經(jīng)通過GigaScience的GiGaDB資源庫(kù)發(fā)布。包括新發(fā)布的基因組序列、基因集、重復(fù)元件、tRNA、miRNA和基因簇等最新毛竹基因組的數(shù)據(jù)為基因組學(xué)研究、遺傳學(xué)研究、分子生物學(xué)實(shí)驗(yàn)等提供了可靠的數(shù)據(jù)資源。

2.2.2 轉(zhuǎn)錄組數(shù)據(jù)分析

為促進(jìn)毛竹全基因組水平的可變剪接圖譜研究，分析轉(zhuǎn)錄后水平影響可變剪接的因素，本研究利用Illumina HiSeq-4000平臺(tái)進(jìn)行高通量轉(zhuǎn)錄組測(cè)序（RNA sequencing，RNA-seq）?？偣矊?duì) 26個(gè)獨(dú)立RNA樣本進(jìn)行雙端測(cè)序 （讀長(zhǎng)為150 bp，原文附表S2，附圖S10-11），測(cè)得每個(gè)樣本高質(zhì)序列平均約為9 000萬條 （約為 13.6 Gb），占原始序列的92.78%。高質(zhì)量序列中約80.57%可以比對(duì)到參考基因組唯一的位置上，并將其標(biāo)記為唯一序列（Unique reads，原文附表S3-4）。大部分序列比對(duì)到外顯子區(qū)域，外顯子的平均比對(duì)率為81.94%。余下的序列中，有8.46%比對(duì)到內(nèi)含子區(qū)域，有9.6%比對(duì)到基因間區(qū)域 （原文附表S5，附圖S12-13）。每個(gè)樣本的平均外顯子覆蓋率約為2 521倍 （原文附圖S14）。因此，大規(guī)模、高深度、高質(zhì)量的轉(zhuǎn)錄組測(cè)序和染色體水平的參考基因組，將有助于精準(zhǔn)識(shí)別基因組中的可變剪接。

為準(zhǔn)確鑒定全長(zhǎng)的剪接體，使用PacBio平臺(tái)對(duì)由毛竹26個(gè)樣本混合得到的RNA進(jìn)行全長(zhǎng)可變剪接體 （FL-cDNA sequencing of alternatively spliced isoforms，Iso-Seq）的測(cè)序。根據(jù)全部樣本的轉(zhuǎn)錄本長(zhǎng)度分布，由混合樣本構(gòu)建了三類單分子實(shí)時(shí)Bell文庫(kù) （3個(gè)轉(zhuǎn)錄本長(zhǎng)度在1 Kb到2 Kb之間的文庫(kù)，2個(gè)轉(zhuǎn)錄本長(zhǎng)度在2 Kb到3 Kb之間的文庫(kù)，4個(gè)轉(zhuǎn)錄本長(zhǎng)度大于3 kb的文庫(kù)），并對(duì)9個(gè)文庫(kù)進(jìn)行測(cè)序，共產(chǎn)生約5 Gb的原始數(shù)據(jù)和214 372條序列（read-of-insert， ROI）。 在所有的ROI中，有133 599個(gè)是全長(zhǎng)ROI（包含了5′、3′和 poly（A） ），剩下的ROI則是非全長(zhǎng)序列 （原文附表S7，附圖S15）。本研究通過將ROI比對(duì)到新基因組進(jìn)行精確度評(píng)估，評(píng)估顯示每一個(gè)核苷酸的錯(cuò)誤率大約是2.05%，其中0.32%的錯(cuò)配，0.98%的插入，0.75%的刪除。

2.2.3 毛竹中有大量的基因存在可變剪接

借助高質(zhì)量基因組和高通量轉(zhuǎn)錄組數(shù)據(jù)，使用前人的分析流程在全基因組范圍內(nèi)鑒定毛竹的可變剪接。結(jié)果顯示，總共在25 225個(gè)可變剪接基因中鑒定了266 711個(gè)特異的可變剪接事件。在所有鑒定的可變剪接基因中，有12 653個(gè)基因在基因注釋環(huán)節(jié)中被注釋為可變剪接基因，剩下的12 572個(gè)基因被認(rèn)定為新鑒定的可變剪接基因 （原文附圖S16）。

作為Illumina RNA-Seq分析的平行實(shí)驗(yàn)，本研究還用Iso-Seq數(shù)據(jù)在相同的流程下檢測(cè)可變剪接。結(jié)果顯示，共鑒定出2 218個(gè)可變剪接基因和4 246個(gè)可變剪接事件。對(duì)兩種數(shù)據(jù)得到的可變剪接結(jié)果進(jìn)行的重疊分析 （評(píng)估可變剪接基因預(yù)測(cè)結(jié)果的可靠程度）顯示，由Iso-Seq數(shù)據(jù)鑒定得到的81.21%可變剪接事件和97.34%可變剪接基因分別與由RNA-seq數(shù)據(jù)鑒定得到的完全重疊。隨后，對(duì)可變剪接進(jìn)行了分類，其中主要的可變剪接類型分別是IR、A3SS、A5SS和ES。平均80.37%的可變剪接事件和95.59%的可變剪接基因?qū)儆谶@主要的4類可變剪接類型 （原文附圖S17）。由PacBio與Illumina平臺(tái)數(shù)據(jù)分別鑒定得到的兩組可變剪接基因之間的覆蓋率很高，證明由計(jì)算機(jī)預(yù)測(cè)得到的可變剪接具有較高的可靠性。

可變剪接事件數(shù)量與可變剪接基因數(shù)量、主要的4類可變剪接的基因數(shù)量之間存在強(qiáng)烈的正相關(guān)（Mann-Whitney U 檢驗(yàn)， R2＞0.91， P＜0.05） （原文圖2C）。毛竹中主要的4類可變剪接由經(jīng)典的剪接模型 （GT-AG、GC-AC和AT-AC剪接位點(diǎn)）從可變剪接事件中鑒定得到。如原文圖2B所示，在主要的4類可變剪接事件中，IR（38.22%）數(shù)量最多，接下來分別是 A3SS（20.20%）和 A5SS（10.48%），最少的是ES（2.92%）。

基因功能富集分析顯示，所有樣本共有的885個(gè)可變剪接基因顯著富集于RNA代謝、mRNA加工、RNA修飾和RNA剪接等功能 （原文附表S25）?？勺兗艚佑兄@著的組織和發(fā)育特異性，鑒定出了181 105個(gè)組織特異的可變剪接事件，是所有可變剪接事件的2/3（67.57%，標(biāo)識(shí)為 “組織間可變剪接”）。剩下的1/3的可變剪接事件通過比較單個(gè)組織內(nèi)的異構(gòu)轉(zhuǎn)錄本得到 （標(biāo)識(shí)為 “組織內(nèi)可變剪接” ）（原文附圖S18）。

轉(zhuǎn)座子 （TE）分析顯示，在26 366個(gè)基因 （占所有基因的51.62%）中存在著轉(zhuǎn)座子的插入，這些插入轉(zhuǎn)座子的總長(zhǎng)度約為46 Mb。通過分析轉(zhuǎn)座子插入內(nèi)含子的位置，發(fā)現(xiàn)插入轉(zhuǎn)座子的內(nèi)含子主要集中于基因的頭尾部分 （原文附圖S19）。此外，使用率最高和分布最多的剪接位點(diǎn)是GT-AG（占所有可變剪接事件的 97.31%），接下來的是 GC-AG（2.33%）和 GT-AT（0.32%）位點(diǎn) （原文附圖S20）。除經(jīng)典的剪接位點(diǎn) （GT-AG、GC-AG和ATAC）之外，將剩下的2 406個(gè)剪接位點(diǎn)標(biāo)記為非經(jīng)典剪接位點(diǎn)，這些非經(jīng)典剪接位點(diǎn)包括2 373個(gè)GTAT類型的剪接位點(diǎn)和33個(gè)其他類型的剪接位點(diǎn)。

2.2.4 毛竹中可變剪接的進(jìn)化分析

通過在8個(gè)物種 （無油樟、擬南芥、油棕、二穗短柄草、水稻、浮萍、高粱和毛竹）中鑒定全基因組同源基因和系統(tǒng)發(fā)育樹的構(gòu)建 （原文圖3A，3B），根據(jù)同源基因起源的時(shí)間不同，定義了8個(gè)特異的同源基因數(shù)據(jù)集 （D8-D1）。比如，同源基因數(shù)據(jù)集7（D7）只包含在1.649億年 （Mya）到2.136億年前之間起源的同源基因 （原文圖3A）。此外，分別從上述的8個(gè)數(shù)據(jù)集中提取單拷貝基因集，標(biāo)記為D8s-D1s。毛竹特有的基因 （D1，包含4 023個(gè)同源基因）較不保守，而8個(gè)物種中都存在的基因 （D8，包含18 997個(gè)同源基因）則高度保守，基因保守程度由D8到D1單調(diào)下降。在所有數(shù)據(jù)集中都檢測(cè)到了可變剪接，但可變剪接基因的比例由D8到D1逐漸降低 （Mann-Whitney U檢驗(yàn)，P＜0.05）。同樣的趨勢(shì)也出現(xiàn)在單拷貝基因集 （D8s-D1s）中。因此，毛竹保守基因集中包含更多的可變剪接基因。

通過對(duì)每個(gè)數(shù)據(jù)集中主要的4類可變剪接分布的研究，發(fā)現(xiàn)不同類型可變剪接的分布趨勢(shì)相同（原文圖3C），但不同類型可變剪接在不同數(shù)據(jù)集中所占的比例不同 （IR ＞A3SS＞A5SS＞ES，Chisquare檢驗(yàn)，P＞0.86）。D8中IR所占的比例為57.76%，約是D1中IR （16.95%）的3.4倍。其他類型可變剪接的比例隨著基因保守度的降低而增高。在兩類數(shù)據(jù)集中，可變剪接事件的數(shù)量從D8到D1、D8s到D1s逐漸降低 （原文圖3C）。此外，對(duì)組織表達(dá)特異性不同的基因 （maxTs，maxTs＝1和maxTs＝0分別表示基本表達(dá)和組織特異性表達(dá)，詳細(xì)見方法部分）的可變剪接事件比較發(fā)現(xiàn)，maxTs與D8-D1這8個(gè)數(shù)據(jù)集中基因的起源時(shí)間存在著負(fù)相關(guān)關(guān)系 （R2＞0.86，P＜0.01），即組織特異性隨著基因保守程度的降低而增大 （原文圖3D）。綜上所述，保守基因集傾向于包含更多的可變剪接基因和事件，且具備更低的組織特異性。

為從整體了解可變剪接，認(rèn)識(shí)可變剪接與基因特征之間的關(guān)系和分析影響可變剪接的因素，研究了不同數(shù)據(jù)集中可變剪接的分布與基因特征之間的關(guān)系 （原文附圖S21）。分析表明不同數(shù)據(jù)集 （D8-D1）的基因與基因長(zhǎng)度、對(duì)應(yīng)編碼序列大小、內(nèi)含子大小、外顯子數(shù)量呈正相關(guān) （R2＞0.9，P＜0.05），與外顯子長(zhǎng)度、內(nèi)含子長(zhǎng)度呈負(fù)相關(guān) （R2＞0.81，P＜0.05）。除此之外，同樣分析了轉(zhuǎn)座子基因在8個(gè)數(shù)據(jù)集 （D8-D1）中的分布，大體上呈負(fù)相關(guān) （R2＞0.77，P＜0.05），表明越保守的基因存在越多的轉(zhuǎn)座子插入。

2.2.5 木質(zhì)素生物合成相關(guān)基因家族的擴(kuò)張和其對(duì)功能多樣性的影響

在擬南芥、二穗短柄草、水稻、楊樹、毛竹和高粱的基因組序列中系統(tǒng)地鑒定出了13個(gè)與木質(zhì)素生物合成相關(guān)的基因家族，這13個(gè)基因家族中的大多數(shù)都存在擴(kuò)張 （原文附表S26）。在毛竹基因組中每個(gè)基因都有著多個(gè)拷貝，其中木質(zhì)素生物合成相關(guān)基因家族的規(guī)模是最大的 （每個(gè)家族平均有19個(gè)成員）。過氧化物酶 （包含77個(gè)成員）和香豆酸-3-羥化酶 （包含3個(gè)成員）基因家族分別是毛竹中成員數(shù)量最多和最少的基因家族。此外，參與木質(zhì)素生物合成的基因在500萬年到1 600百萬年前之間分化，這與毛竹在700萬年到1 200萬年前之間發(fā)生的全基因組復(fù)制 （WGD）事件相對(duì)應(yīng)。

之后，對(duì)木質(zhì)素生物合成相關(guān)基因的可變剪接進(jìn)行了分析 （原文圖4）。除了阿魏酸-5-羥化酶（F5H）基因家族中的可變剪接比例過低，查耳酮合成 （CHS））基因家族和咖啡酸鄰O-甲基轉(zhuǎn)移酶（COMT）基因家族沒有檢測(cè)到可變剪接基因之外，總計(jì)13個(gè)基因家族中的10個(gè) （超過總數(shù)的一半）存在可變剪接基因。在肉桂酸4-羥基肉桂酰輔酶A連接酶 （4CL）、羥基肉桂酰轉(zhuǎn)移酶 （HCT）、肉桂醇脫氫酶 （CAD）基因家族中觀察到高比例 （＞75%）的可變剪接事件。最后，使用branch-site模型對(duì)參與木質(zhì)素生物合成的相關(guān)基因家族進(jìn)行正選擇作用分析表明，HCT和CAD這兩個(gè)基因家族中的一些基因檢測(cè)到正選擇作用。使用的最適模型得到的系統(tǒng)發(fā)育關(guān)系和對(duì)數(shù)似然比 （log likelihood ratio）等信息詳見原文附表S27。

3 討論

3.1 毛竹基因組的完善與質(zhì)量提升

目前，隨著新技術(shù)的發(fā)展和更多高通量數(shù)據(jù)的出現(xiàn)，高通量基因組測(cè)序和改進(jìn)的組裝技術(shù)普遍應(yīng)用于植物基因組研究中。在2013年，通過對(duì)毛竹基因組的初步分析，得到了毛竹的基因組、基因結(jié)構(gòu)、關(guān)鍵功能基因等數(shù)據(jù)，還了解了毛竹中發(fā)生的WGD事件。通過本次研究，提高了毛竹基因組的精確度和完整性，更新了基因組注釋，且通過不同物種之間的比較研究精準(zhǔn)定位了毛竹在進(jìn)化上的位置。此外，利用從最新版本毛竹基因組得到的信息對(duì)毛竹不同的生物特性進(jìn)行了詳細(xì)的研究。染色體水平的參考基因組和準(zhǔn)確的注釋將有利于今后毛竹和其他近緣物種的基因組學(xué)研究。

3.2 毛竹中可變剪接具有普遍性，且在不同組織中表現(xiàn)出多樣性

RNA-seq和Iso-Seq高通量數(shù)據(jù)有助于精準(zhǔn)檢測(cè)低表達(dá)水平的轉(zhuǎn)錄本、鑒定完整的基因結(jié)構(gòu) （特別是在可變剪接的研究中），通過分析獲得了毛竹全基因組范圍的可變剪接圖譜。通過對(duì)毛竹中可變剪接的研究，加深了對(duì)轉(zhuǎn)錄后調(diào)控層面上的可變剪接的認(rèn)知，包括可變剪接基因與事件的鑒別、不同種類可變剪接的分布、不同種類剪接位點(diǎn)的使用率、可變剪接體外顯子的長(zhǎng)度分布等?？勺兗艚颖徽J(rèn)為是基因數(shù)量有限的生物產(chǎn)生多樣性的主要機(jī)制。例如，果蠅中細(xì)胞黏附分子基因通過組合分別含有12、48、33和2個(gè)外顯子的4個(gè)可變剪接基因，最多可產(chǎn)生38 016（12×48×33×2） 種蛋白質(zhì)異構(gòu)體。在毛竹中，共鑒定出266 711個(gè)特異的可變剪接事件和25 225個(gè)可變剪接基因，平均每個(gè)樣本檢測(cè)到15 071個(gè)可變剪接事件和9 080個(gè)可變剪接基因。因此，可變剪接可能存在組織特異性，而且毛竹中實(shí)際的可變剪接比例可能被低估了。隨著測(cè)序深度的增加，可以從表達(dá)水平低的轉(zhuǎn)錄本中檢測(cè)出更多的可變剪接事件。

根據(jù)觀察，毛竹中地下莖組織比根組織包含更多的可變剪接事件，這可能是因?yàn)檫@兩個(gè)組織在毛竹的生長(zhǎng)過程中扮演的角色不同。處于高速生長(zhǎng)期的竹筍沒有葉子不能進(jìn)行光合作用，所以竹筍生長(zhǎng)所需的大量營(yíng)養(yǎng)物質(zhì)和能量需要通過地下莖從成熟竹子運(yùn)輸?shù)缴L(zhǎng)部位。因此，作為地下莖植物，毛竹的地下莖在營(yíng)養(yǎng)物質(zhì)和能量運(yùn)輸中扮演著核心角色，這可能解釋了為什么會(huì)在地下莖組織樣本中檢測(cè)到更多的可變剪接事件。

3.3 毛竹一些生物學(xué)特性可能與不同類型可變剪接的比例有關(guān)

在全基因組范圍上鑒定可變剪接時(shí)，經(jīng)常需要分析不同種類可變剪接之間的差異，因?yàn)椴煌N類可變剪接的頻率或比例的差異可能反映了pre-mRNA剪接的差異。不同種類可變剪接的分布結(jié)果表明：IR類型的可變剪接在毛竹中數(shù)量最多，且IR普遍存在于其他處在不同進(jìn)化位置上的植物中，這可能與IR的重要性有關(guān)。然而，與其他植物相比，毛竹中IR類型 （38.22%）和其他類型 （28.18%）的可變剪接比例更高。這可能與毛竹的特性或轉(zhuǎn)錄組測(cè)序深度有關(guān)，也可能是兩者共同造成的，未來可以通過比較分析進(jìn)行辨別。此外，毛竹中主要的4類可變剪接的分布與擬南芥、大豆、玉米中的一致。然而，在動(dòng)物和酵母中不同種類可變剪接的分布與植物中的不同。在動(dòng)物和酵母中，可變剪接事件數(shù)量最多的一類是ES，接下來分別是A3SS和A5SS，最少的是IR?？梢娭参锖蛣?dòng)物在不同種類可變剪接的分布上存在差異，這暗示了植物和動(dòng)物的基因組結(jié)構(gòu)和剪接位點(diǎn)識(shí)別機(jī)制存在差異。此外，對(duì)單一基因剪接位點(diǎn)的鑒定提供了理解可變剪接和剪接體結(jié)構(gòu)的重要數(shù)據(jù)。與各類可變剪接的分布情況不同的是，毛竹中主要的3種可變剪接位點(diǎn) （GT-AG，GC-AG和AT-AC）分布和之前在動(dòng)物和其他植物中觀察到的規(guī)律是一致的。

3.4 包含更多可變剪接事件的保守基因可能在進(jìn)化和功能中扮演著重要的作用

目前，基因保守性和可變剪接之間的關(guān)系仍不清楚。為了探究這個(gè)問題，在基因組層面分析兩類同源基因集中的可變剪接，這2類基因集按保守程度的不同又各自細(xì)分出8個(gè)子集，分別是D8-D1和D8s-D1s。分析表明可變剪接基因傾向于集中在高保守基因集中，且高保守基因集的可變剪接基因包含著更多的可變剪接事件。本研究在2類同源基因集中發(fā)現(xiàn)了相同的趨勢(shì)，說明這個(gè)發(fā)現(xiàn)是可靠的。有研究表明復(fù)制是功能分化和新基因產(chǎn)生的主要原因之一，且保守基因在基因互作網(wǎng)絡(luò)中有著更高連接度的傾向。新基因因?yàn)樾录尤氲交蚧プ骶W(wǎng)絡(luò)，所以其基因連接度較低，但隨著時(shí)間推移，新基因的連接度和重要程度逐漸增加。本研究表明，保守度高的基因比保守度低的基因傾向于擁有更多的可變剪接事件，這與保守度高的基因在基因互作網(wǎng)絡(luò)中連接度更高的趨勢(shì)是一致的。因此，推測(cè)在進(jìn)化過程中基因連接度的增加可能與可變剪接有關(guān)。此外，與保守度低的基因集相比，保守度高的可變剪接基因集的組織特異性低，表明在基因互作網(wǎng)絡(luò)中有著更高的連接度的基因可能在基礎(chǔ)功能中發(fā)揮著重要作用。因此，推測(cè)一些基因之所以重要可能是因?yàn)樗鼈儼嗟目勺兗艚邮录??？勺兗艚幼鳛橐粋€(gè)重要的生物學(xué)過程，在生物功能多樣化的過程中扮演著重要的角色，這可能解釋了保守度高的可變剪接基因?yàn)槭裁窗嗟目勺兗艚邮录?。推測(cè)這個(gè)現(xiàn)象可能不僅存在于毛竹中，而且還存在于其他植物甚至是動(dòng)物中。

3.5 毛竹中木質(zhì)素生物合成相關(guān)基因家族擴(kuò)增和可變剪接可能與WGD事件有關(guān)

木質(zhì)素是一種由木質(zhì)素單體組成的復(fù)雜芳香族聚合物，它與纖維素和半纖維素相互作用共同構(gòu)成次生細(xì)胞壁。木質(zhì)素約占毛竹干重的25%。通過可變剪接分析和進(jìn)化分析，檢測(cè)到了木質(zhì)素生物合成相關(guān)基因家族的擴(kuò)增。結(jié)合木質(zhì)素生物合成相關(guān)基因的分化時(shí)間研究結(jié)果和我們前期的研究，估測(cè)毛竹基因組在前700萬年到1 200萬年前之間發(fā)生了一次WGD，這暗示在其進(jìn)化過程中可能存在一次四倍體事件。之后，四倍體祖先進(jìn)化成為現(xiàn)在的二倍體毛竹。WGD可以提供更多的基因拷貝，并通過產(chǎn)生新功能而加速基因的進(jìn)化。因此，毛竹中發(fā)生的WGD事件可能導(dǎo)致了木質(zhì)素生物合成相關(guān)基因的擴(kuò)增。此外，HCT和CAD這兩個(gè)基因家族被檢測(cè)到正選擇作用，且包含著更多的可變剪接事件。HCT介導(dǎo)對(duì)-香豆酰輔酶A（也被CHS介導(dǎo)產(chǎn)生類黃酮）產(chǎn)生木質(zhì)素。HCT與CHS都可以結(jié)合對(duì)-香豆酰輔酶A，它們之間存在相互競(jìng)爭(zhēng)作用。在毛竹中，HCT基因家族與CHS基因家族相比，有著更多的基因成員和可變剪接事件，說明HCT基因家族可能在與CHS基因家族對(duì)-香豆酰輔酶A競(jìng)爭(zhēng)結(jié)合中占據(jù)著主導(dǎo)地位。CAD可以催化不同的底物合成不同的木質(zhì)素。在毛竹中木質(zhì)素芳香聚合物由以下三個(gè)單體組成：對(duì)羥基苯基丙烷 （H）、香草醛 （G）、丁香醛（S）。前人的研究表明，在毛竹中G、S木質(zhì)素的含量高，而H木質(zhì)素的含量低。毛竹中CAD基因家族的擴(kuò)增和正選擇作用可能解釋了毛竹中不同木質(zhì)素單體比例差異帶來的基質(zhì)偏好性。毛竹木質(zhì)素生物合成相關(guān)基因中存在的大量可變剪接事件、基因擴(kuò)增和正選擇作用，這與毛竹強(qiáng)大的木質(zhì)素合成能力一致。

4 結(jié)論

為從進(jìn)化角度深入探究毛竹可變剪接情況，組裝出了染色體水平的基因組，并重新進(jìn)行了基因組注釋。通過分析染色體水平的基因組和大量轉(zhuǎn)錄組數(shù)據(jù)，得到了毛竹中完整的可變剪接圖譜，鑒定出了25 225個(gè)可變剪接基因中的266 711個(gè)特異的可變剪接事件。此外，毛竹中可變剪接結(jié)果的綜合分析和8個(gè)代表性植物物種的比較分析表明，保守基因有著高表達(dá)和低組織特異性的趨勢(shì)。對(duì)毛竹中木質(zhì)素生物合成相關(guān)的基因進(jìn)行了可變剪接分析和進(jìn)化分析，觀察到木質(zhì)素生物合成相關(guān)的基因家族發(fā)生了擴(kuò)增，這些基因家族包含了大量的可變剪接和正選擇作用。綜上所述，本研究為毛竹獨(dú)特的材性研究和從進(jìn)化角度探索可變剪接提供了重要的資源。