基于毛竹重測序的基因可變剪接研究

王炯亮 趙韓生 高志民

（國際竹藤中心北京100102）

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

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

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

1 材料方法與數據分析

1.1 樣本采集

為得到完整的可變剪接圖譜，于2015年春從中國6個主要竹產區采集毛竹樣本進行測序。6個采集地分別是：江蘇宜興 （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）。收集了處于不同生長階段毛竹的26個組織，包括地下莖、根、筍、葉、籜片和芽。有關生物樣本采集的詳細信息參見原文附表S19。

1.2 基因組測序、組裝、注釋

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

1.3 RNA的分離、Illumina RNA-seq建庫與測序、數據分析

使用前人的方法對RNA進行分離、純化、濃度測定、反轉錄和cDNA文庫的構建。所有cDNA文庫的構建和標準化過程信息詳見原文附件。在質控之后，對文庫進行光學檢驗，隨后使用HiSeq-4000平臺 （雙端測序，插入長度為150 bp）進行測序。然后，計算了序列的質量，并使用FastQC（軟件版本0.11.3）在默認參數下對低質量序列進行過濾。最后，使用RNA-SeQC（軟件版本1.1.8）在默認參數下計算得到RNA-seq中的關鍵統計指標。

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

1.4 Iso-Seq文庫的構建、測序與數據分析

按照測序儀制造商 （PacBio）提供的方法，對Iso-Seq文庫進行構建和測序。首先根據生物信息學預測得到的轉錄本長度分布 （原文附表S22）構建了3類SMRTBell文庫 （轉錄本長度在1 Kb到2 Kb之間的3個文庫，轉錄本長度為2 Kb到3 Kb之間的2個文庫，轉錄本長度大于3 kb的4個文庫），然后使用PacBio測序平臺對這9個SMRT文庫進行測序。

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

1.5 進化分析

依據前人的研究方法，進行基因家族的鑒定、系統發育樹的構建和分化時間的預測。在原文附件和protocol.io提供了詳細的信息。

1.6 木質素生物合成相關基因的全基因組鑒定與分析

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

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

2 結果與分析

2.1 數據描述

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

2.2 數據分析

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

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

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

2.2.2 轉錄組數據分析

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

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

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

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

作為Illumina RNA-Seq分析的平行實驗，本研究還用Iso-Seq數據在相同的流程下檢測可變剪接。結果顯示，共鑒定出2 218個可變剪接基因和4 246個可變剪接事件。對兩種數據得到的可變剪接結果進行的重疊分析 （評估可變剪接基因預測結果的可靠程度）顯示，由Iso-Seq數據鑒定得到的81.21%可變剪接事件和97.34%可變剪接基因分別與由RNA-seq數據鑒定得到的完全重疊。隨后，對可變剪接進行了分類，其中主要的可變剪接類型分別是IR、A3SS、A5SS和ES。平均80.37%的可變剪接事件和95.59%的可變剪接基因屬于這主要的4類可變剪接類型 （原文附圖S17）。由PacBio與Illumina平臺數據分別鑒定得到的兩組可變剪接基因之間的覆蓋率很高，證明由計算機預測得到的可變剪接具有較高的可靠性。

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

基因功能富集分析顯示，所有樣本共有的885個可變剪接基因顯著富集于RNA代謝、mRNA加工、RNA修飾和RNA剪接等功能 （原文附表S25）。可變剪接有著顯著的組織和發育特異性，鑒定出了181 105個組織特異的可變剪接事件，是所有可變剪接事件的2/3（67.57%，標識為 “組織間可變剪接”）。剩下的1/3的可變剪接事件通過比較單個組織內的異構轉錄本得到 （標識為 “組織內可變剪接” ）（原文附圖S18）。

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

2.2.4 毛竹中可變剪接的進化分析

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

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

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

2.2.5 木質素生物合成相關基因家族的擴張和其對功能多樣性的影響

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

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

3 討論

3.1 毛竹基因組的完善與質量提升

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

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

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

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

3.3 毛竹一些生物學特性可能與不同類型可變剪接的比例有關

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

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

目前，基因保守性和可變剪接之間的關系仍不清楚。為了探究這個問題，在基因組層面分析兩類同源基因集中的可變剪接，這2類基因集按保守程度的不同又各自細分出8個子集，分別是D8-D1和D8s-D1s。分析表明可變剪接基因傾向于集中在高保守基因集中，且高保守基因集的可變剪接基因包含著更多的可變剪接事件。本研究在2類同源基因集中發現了相同的趨勢，說明這個發現是可靠的。有研究表明復制是功能分化和新基因產生的主要原因之一，且保守基因在基因互作網絡中有著更高連接度的傾向。新基因因為新加入到基因互作網絡，所以其基因連接度較低，但隨著時間推移，新基因的連接度和重要程度逐漸增加。本研究表明，保守度高的基因比保守度低的基因傾向于擁有更多的可變剪接事件，這與保守度高的基因在基因互作網絡中連接度更高的趨勢是一致的。因此，推測在進化過程中基因連接度的增加可能與可變剪接有關。此外，與保守度低的基因集相比，保守度高的可變剪接基因集的組織特異性低，表明在基因互作網絡中有著更高的連接度的基因可能在基礎功能中發揮著重要作用。因此，推測一些基因之所以重要可能是因為它們包含著更多的可變剪接事件。可變剪接作為一個重要的生物學過程，在生物功能多樣化的過程中扮演著重要的角色，這可能解釋了保守度高的可變剪接基因為什么包含著更多的可變剪接事件。推測這個現象可能不僅存在于毛竹中，而且還存在于其他植物甚至是動物中。

3.5 毛竹中木質素生物合成相關基因家族擴增和可變剪接可能與WGD事件有關

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

4 結論

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