转录组学对转基因杨树表型变化的代谢调控机制解析

和成成，吴照群，许 楠，张秀省，杜 鑫，杨海灵

(北京林业大学 生物科学与技术学院，北京 100083)

细胞分裂是所有细胞生物的基本特征之一，在动、植物和微生物生长发育过程中发挥着极其重要的作用。对细胞周期的研究一直是科学研究的热点，它对动、植物的生长发育以及人类疾病研究都至关重要。与动物不同，大部分植物属于胚后发育生物，这说明植物器官的形成与细胞分裂、分化以及细胞生长都联系紧密。细胞周期的正常进行需要许多周期调控因子的共同参与。其中，Cyclin作为一种随细胞周期的进行而不断合成和降解的蛋白，在维持细胞周期正常进行的过程中发挥十分重要的作用。一般认为，D型细胞周期蛋白(CYCD)可以与CDKA形成复合物，控制着细胞周期G1/S的过渡[1]。在拟南芥中，共鉴定出10个D型Cyclin[2]；而杨树中鉴定出24个D型Cyclin[3]，其中D2/4亚类在拟南芥和杨树中分别有3个和2个成员[3-4]。这表明木本植物中的D2/4类Cyclin在进化过程中更加保守。

在草本植物的研究中，已经发现CYCD对植物生长发育至关重要。叶片受伤或边缘被去除会促进CYCD;1的表达，通过与CDKA相互作用来增强CDKA激酶活性，进而造成叶片边缘细胞重新进入细胞周期[5]。ICK2/KRP2可与CYCD2;1相互作用，ICK2/KRP2通过保持CYCD2;1失活来抑制根分枝，调节中柱鞘细胞对生长素波动的响应[6]。拟南芥KRP2抑制CDK活性，过表达KRP2导致叶片呈现锯齿状，叶片表皮细胞和叶肉细胞增大，同时抑制叶片核内周期[7]。拟南芥CYCD3参与调控维管发育，cycd3;1-3缺失突变体植物的茎和下胚轴明显变细，形成层活性降低，次生生长受阻[8]。过表达AtCYCD3;1使得叶片无法形成正常的栅栏组织和海绵组织，同时叶表皮细胞由大量不完全分化的小的多边形细胞组成[9]。CYCD3的缺失会损害芽分生组织功能并导致对细胞分裂素的应答降低[10]。AtCYCD2;1的转基因品系提高了CDK激酶活性，缩短了G1期持续时间并减小了在S期和分裂阶段的细胞大小。转基因品系的叶细胞较小且数量更多，表皮细胞和栅栏组织细胞大小分别为对照的59%和69%，而叶片则是正常大小[11]。在木本植物中，过表达PtaCYCD1;2会造成杨树产出波纹状表型，叶片和茎明显卷曲，且维管组织所占的比例增大[12]。过表达PtoCYCD3;3造成杨树叶片变大且褶皱，茎增粗，同时出现分枝[3]。无论是在草本还是木本植物中，CYCD的表达量变化(缺失和过表达)会造成植株出现一系列的表型变化，如叶片形态、株高和茎粗等。

课题组前期研究发现，OE-PtoCYCD2;1造成杨树生长减缓，株高降低，叶片向远轴面卷曲且茎粗变细[13]。本研究主要以OE-PtoCYCD2;1杨树为实验材料，通过转录组学测序来研究PtoCYCD2;1对杨树生长发育的影响，为研究木本植物D型细胞周期蛋白功能提供理论基础。

1 材料和方法

1.1 实验材料

选取1年生的经过qRT-PCR鉴定过的野生型(WT)和过表达(OE-PtoCYCD2;1) ‘741杨’(OE)进行扦插扩繁[13]。扦插后的植株生长于温室中，栽培条件为20～25 ℃，光照强度3 500 Lx，16 h光照/8 h黑暗。选取生长3个月的幼苗幼叶进行转录组测序，WT和OE-PtoCYCD2;1均设置3个生物学重复，共6个样本。将幼叶取下后迅速置于液氮中速冻，-80 ℃保存。

1.2 方 法

1.2.1 总RNA提取和文库构建转录组测序委托北京诺禾致源科技股份有限公司完成。采用标准提取方法从组织中提取RNA，随后对RNA样品进行严格质控，质控方式主要是通过Agilent 2100 bioanalyzer精确检测RNA完整性。

建库起始RNA为总RNA，通过Oligo(dT)磁珠富集带有polyA 尾的mRNA，随后在Fragmentation Buffer中用二价阳离子将得到的mRNA 随机打断。以片段化的mRNA 为模版，随机寡核苷酸为引物，在M-MuLV 逆转录酶体系中合成cDNA 第一条链，随后用RNaseH 降解RNA 链，并在DNA 聚合酶 I 体系下，以dNTPs 为原料合成cDNA 第二条链。纯化后的双链cDNA经过末端修复、加A 尾并连接测序接头，用AMPure XP beads 筛选370～420 bp 左右的cDNA，进行PCR 扩增并再次使用AMPure XP beads 纯化PCR 产物，最终获得文库。将样本库检合格后的样本进行Illumina测序。

1.2.2 测序结果的统计和功能注释为了保证数据分析的质量及可靠性，需要对原始数据进行过滤。主要包括去除带接头(adapter) 的reads、去除含N(N 表示无法确定碱基信息)的reads、去除低质量reads(Qphred ≤20 的碱基数占整个read长度的50%以上的reads)。同时，对clean data 进行Q20、Q30和GC 含量计算。后续所有分析均是基于clean data进行的高质量分析。

选择NCBI上的毛果杨(Populustrichocarpa)基因组作为参考。使用HISAT2 v2.0.5 构建参考基因组的索引，并使用HISAT2 v2.0.5 将配对末端clean reads 与参照基因组比对。采用StringTie(1.3.3b)进行新基因预测。采用subread软件中的featureCounts工具用于计算比对上的reads映射到每个基因的读数。同时根据基因的长度计算每个基因的FPKM。

1.2.3 数据的统计与分析使用DESeq2 软件(1.20.0)进行差异表达基因分析。作为有生物学重复的样本，分析过程以|log2FC| ≥ 0 且 FDR(Q值) <0.05 作为差异基因筛选条件。通过clusterProfiler 3.4.4软件实现差异表达基因的GO、KEGG、PFAM和转录因子富集分析，其中转录因子富集分析使用PlantTFDB数据库对转录因子进行预测分析。

1.2.4 qRT-PCR验证RNA提取使用天根生化科技(北京)有限公司的RNAprep Pure植物总RNA提取试剂盒；cDNA合成使用天根生化科技(北京)有限公司的FastKing cDNA第一链合成试剂盒(去基因组)；qRT-PCR使用天根生化科技(北京)有限公司的SuperReal PreMix Plus (SYBR Green)试剂盒进行实验。用Primer 5进行引物设计，引物序列见表1。使用2-ΔΔCT计算所选差异表达基因的相对表达水平。

1.2.5 生理指标测定植物叶片叶绿素含量测定采用分光光度法，浸提法提取叶绿素[14]，提取液溶液为丙酮∶无水乙醇∶水(4.5∶4.5∶1)；植物叶片可溶性糖含量检测使用索莱宝生物科技有限公司的植物可溶性糖含量检测试剂盒；植物叶片和茎段木质素含量检测使用索莱宝生物科技有限公司的木质素含量检测试剂盒；以上3种生理指标测定时均分别取3株WT和OE-PtoCYCD2;1植株作为生物学重复，并且测定3次作为技术重复。

表1 qRT-PCR引物序列

2 结果与分析

2.1 测序数据质量分析

所有样品在经过llumina平台测序后生成原始的FASTQ数据，对原始测序数据进行质控统计。由表2可知，所有样品的clean reads均达到6G以上，Q20的百分比均在98%以上，Q30的百分比在94%以上，GC含量均在44%左右，这些数据表明测序数据良好，可以进行下一步分析研究。

2.2 差异表达基因分析

通过计算FPKM，6个样本中共统计得到42 816个基因。将OE-PtoCYCD2;1的所有基因表达量与WT相比，以|log2FC| ≥ 0且FDR(Q值)<0.05为筛选标准进行差异表达基因筛选，共筛选得到1 269个差异表达基因，其中上调基因700个，下调表达基因569个(图1)。

2.2.1 GO分析对统计到的差异表达基因进行GO功能富集分析，主要分为生物过程(BP)、分子功能(MF)和细胞组成(CC) 3大类，气泡图中展示了在三大类中富集最显著的30个小类(图2，A)。其中细胞壁生物发生、植物型细胞壁生物发生、细胞碳水化合物代谢过程、几丁质的响应和植物型细胞壁组织或生物发生是差异基因富集最显著的生物学过程；葡萄糖基转移酶活性、UDP-葡萄糖基转移酶活性、UDP-糖基转移酶活性、天冬氨酸型内肽酶活性和天冬氨酸型肽酶活性是差异基因富集最显著的分子功能；Cul4A-RING E3泛素连接酶复合物、Cul4-RING E3泛素连接酶复合物、光系统Ⅱ析氧复合体、胞外区和类囊体膜是差异基因富集最显著的细胞组分。

2.2.2 KEGG分析对筛选到的差异基因进行KEGG富集分析，气泡图中展示了富集最显著的20个代谢通路(图2，B)，其中碳代谢、半胱氨酸和甲硫氨酸代谢、甘氨酸、丝氨酸和苏氨酸代谢、光合生物中的碳固定、植物昼夜节律、氨基酸的生物合成、苯丙烷生物合成、乙醛酸和二羧酸代谢、磷酸肌醇代谢和淀粉和蔗糖代谢是富集最显著的代谢通路。

2.2.3 Pfam分析对筛选到的差异基因进行Pfam富集分析，气泡图中展示了富集最显著的10种结构域(图3，A)。其中AP2结构域、N 端木聚糖酶抑制剂、C 端木聚糖酶抑制剂、真核天冬氨酰蛋白酶和糖基转移酶家族富集最显著。

通过使用PlantTFDB对转录因子进行预测分析，共预测出2 422个转录因子，其中包括145个差异表达基因，这些差异表达基因分别属于31个转录因子基因家族。图3，B中展示了数目最多的10种转录因子，其中数目最多的为ERF转录因子且均上调表达，这与结构域富集结果相同。植物AP2/ERF类转录因子参与植物生长发育过程[15]，这可能暗示OE-PtoCYCD2;1主要通过调控ERF转录因子来调控植物的生长发育。

表2 样本测序数据质量汇总

2.2.4 转录组数据验证为验证转录组数据可靠性，选取9个差异表达基因进行qRT-PCR实验验证(图4)。结果显示，9个差异表达基因的qRT-PCR结果与RNA-seq测定的差异基因表达水平变化趋势基本一致，表明了本研究所用RNA-seq结果具有极高的可靠性。

2.3 PtoCYCD2;1对杨树的影响分析

2.3.1 木质素合成OE-PtoCYCD2;1与WT相比，明显显示出植株株高降低，叶片卷曲[13]。同时GO富集分析显示，在植物细胞壁合成相关的生物学过程中的相关基因表达明显发生改变。进一步分析发现，有8个差异表达基因富集到木质素合成通路中，且均下调表达(图5)。其中有5个漆酶(Laccase)基因均下调表达，这可能表明PtoCYCD2;1主要通过影响LAC表达来影响植物木质素合成。这表明OE-PtoCYCD2;1会影响杨树木质素合成途径中相关基因表达，这可能是造成OE-PtoCYCD2;1植株株高降低，茎变细的原因。

2.3.2 卡尔文循环KEGG富集分析显示，29个差异表达基因富集到碳代谢通路中，上调表达基因2个，下调表达基因27个。卡尔文循环作为植物中光合碳同化的基本途径，在植物生长发育中至关重要。对碳代谢通路中富集的差异表达基因进一步分析发现，其中8个与卡尔文循环有关的被富集到，且均下调表达(图6)。该通路中有3个果糖-1,6-二磷酸醛缩酶(fructose-1,6-bisphosphate aldolase,FBA)基因均下调表达，可能表明PtoCYCD2;1主要通过影响FBA表达来影响植物木质素合成。总之，OE-PtoCYCD2;1会影响杨树卡尔文循环途径中相关基因表达，对植物生长发育产生影响。这可能是造成OE-PtoCYCD2;1植株整体生物量降低的原因。

2.4 PtoCYCD2;1对杨树生理指标的影响

通过观察WT和OE-PtoCYCD2;1的生长过程，发现与WT相比，OE-PtoCYCD2;1的叶片颜色更绿，由此猜测可能由于叶片的叶绿素含量差异造成的。通过检测WT和OE-PtoCYCD2;1幼叶和成熟叶的叶绿素含量，结果发现，幼叶和成熟叶的叶绿素a(Chla)分别增加52.98%和57.22%；幼叶和成熟叶的叶绿素b(Chlb)分别增加71.44%和138.16%，较Chla增加的比例更大，尤其是成熟叶；幼叶和成熟叶的总叶绿素含量分别增加57.36%和78.22%(表3)。结果说明：OE-PtoCYCD2;1会提高单位质量下叶片(无论是幼叶还是成熟叶)的叶绿素含量，尤其是Chlb的增加更加明显。

植物体内可溶性糖含量通常作为植物碳代谢变化的指标之一。通过检测WT和OE-PtoCYCD2;1的幼叶和成熟叶的可溶性糖含量，结果发现，与WT相比，OE-PtoCYCD2;1的幼叶可溶性糖含量变化不大，但成熟叶的可溶性糖含量下降了12.72%(表4)。结果表明，OE-PtoCYCD2;1会降低单位质量下叶片的可溶性糖含量，且对成熟叶的影响更大。

转录组结果表明，OE-PtoCYCD2;1植株中与木质素合成相关的基因下调表达，因此我们检测了WT和OE-PtoCYCD2;1的叶片和茎段中木质素含量。结果发现，与WT相比，OE-PtoCYCD2;1幼叶中木质素含量下降了4.48%，而成熟叶中木质素含量下降了8.03%；OE-PtoCYCD2;1幼茎中木质素含量下降了20.03%，成熟茎中木质素含量下降最多，达到了31.63%(表4)。该结果表明，OE-PtoCYCD2;1会降低单位质量下叶片和茎部的木质素含量，在成熟茎中下降更加明显。

表3 2种材料叶片叶绿素含量

表4 2种材料叶片可溶性糖和木质素含量

3 讨 论

植物生长发育离不开细胞周期的正常进行。与动物相比，植物细胞具有全能性，单个植物细胞可以生长成为一个新的个体，这离不开细胞分裂和分化。细胞周期蛋白(Cyclin)、细胞周期依赖性激酶(CDK)和细胞周期依赖性激酶抑制剂(ICK)形成三元复合物[1]，控制着细胞周期的正常进行，其中D型细胞周期蛋白(CYCD)参与G1/S的转换。本研究在转录组中鉴定出所有的CDK和ICK，但发现只有CDKG;4上调表达。在草本植物拟南芥中，CDKG亚家族中有两个成员：CDKG1和CDKG2[16]。CDKG1蛋白在高温下维持花粉育性和产量有重要作用。缺乏 CDKG1的植物表现出正常的营养生长，在高温环境下生长不育，但在较低温度下可育。相比之下，缺乏 CDKG2 的植物在高温和低温下都可以完全繁殖[17]。CDKG1 和 CDKG2 均已显示介导下游基因的 AS。在营养组织中CDKG1 调节剪接因子 U2AF65A 的 AS，而在花粉中它调节 CalS5 mRNA 的 AS[18-19]。在杨树中，CDKG基因家族扩张至5个成员[3]，它们之间的功能尚未有报道。OE-PtoCYCD2;1植株中PtoCDKG;4表达上调，它们之间的联系和造成这种变化的原因仍需进一步研究。

CLR5为AP2/ERF转录因子家族蛋白的一个成员，它受生长素调控表达并可以通过正向调节 A 型 RR、OsRR1 来抑制细胞分裂素信号传导，从而促进冠根启动[20]。PtaERF1在韧皮部组织中表达，可能暗示PtaERF1在维管组织发育和功能中具有潜在作用[21]。BOLITA基因(属于 AP2/ERF 家族)的转基因植物的叶子小于非转基因植物，叶片变小是由于细胞大小和细胞数量均减少所造成[22]。本研究发现有26个AP2/ERF转录因子均上调表达，表明PtoCYCD2;1可能主要通过影响AP2/ERF转录因子的表达来调控植物的生长发育。

木质素作为植物生长发育中的一种次生代谢物质，在植物的生长发育和抗性方面具有重要的生物学功能[23]。木质素是植物细胞壁的主要成分之一，是一种分子量大、组成和结构复杂的天然酚类聚合物。木质素代谢参与植物生长发育，干扰木质素生物合成会导致植物生长发育受抑制和畸形发育[24]。2个拟南芥CCR1敲除突变体的花序茎表现出木质素水平降低，且两者都具有植株矮小和延迟衰老的表型[25]。本研究通过转录组分析发现8个木质素合成通路基因下调表达，且叶片和茎中单位质量下的木质素含量均下降，又因为OE-PtoCYCD2;1植株的叶片更小，茎段更细，表明OE植株的总体质量会有所下降[13]，导致植株总体木质素总量更少，这可能是造成OE植株株高降低、茎粗变细的原因。

拟南芥CYCD2的表达受到蔗糖的诱导[6,26]，这可能表明CYCD2的表达在植物碳代谢上存在一定的相关性。本研究发现，共有29个差异基因富集到碳代谢通路上，2个差异基因上调表达，27个差异基因下调表达。其中，一些重要的碳代谢相关基因，如3-磷酸甘油醛脱氢酶(GADPH)、磷酸甘油异构酶(TIM)、丙酮酸激酶(PK)、苹果酸脱氢酶(MDH)等均下调表达。其中卡尔文循环通路有8个差异表达基因也下调表达。以植物体内的可溶性糖作为检测植株碳代谢变化的指标，发现叶片的可溶性糖含量降低，又由于OE植株的叶片面积变小，整体植株的质量有所下降[13]，推测植株总体可溶性糖含量降低，与转录组分析结果相呼应。表明PtoCYCD2;1可能调控植物的碳代谢水平，降低植物对CO2的利用效率，造成过表达植株的总生物量下降。

综上，根据OE-PtoCYCD2;1的转录组数据，结合已发表的转基因表型和文献分析，认为PtoCYCD2;1的过表达能上调多个AP2/ERF转录因子的表达，并且下调碳代谢通路和木质素通路中一些关键基因的表达，减少植株中蔗糖等碳源和木质素的含量，从而造成过表达植株表型改变，植株总体生物量下降。本研究对木本植物中细胞周期蛋白表达量改变引起植株表型变化方面的功能探讨提供了很好的研究思路和方向。