基于毛竹重测序的基因可变剪接研究

王炯亮 赵韩生 高志民

（国际竹藤中心北京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个代表性植物物种的比较分析表明，保守基因有着高表达和低组织特异性的趋势。对毛竹中木质素生物合成相关的基因进行了可变剪接分析和进化分析，观察到木质素生物合成相关的基因家族发生了扩增，这些基因家族包含了大量的可变剪接和正选择作用。综上所述，本研究为毛竹独特的材性研究和从进化角度探索可变剪接提供了重要的资源。