插穗年龄对毛白杨不定根形成影响的转录组分析

Transcriptome Analysis of the Effect of Cutting Age on Adventitious Root Formation of Populus

摘要: 为了探究插穗年龄对不定根生根的影响，本研究以1、2、3年生毛白杨杂交三倍体新品种“北林雄株2号”的嫩枝插穗为研究材料，扦插五天后进行转录组测序分析，共检测到5153个差异表达基因。GO分析显示，上调基因主要富集在线粒体呼吸、氧化磷酸化和类黄酮合成等过程，下调基因主要富集在油菜素内酯合成、水杨酸信号转导和赤霉素合成等途径。KEGG结果显示，差异基因在信号转导、能量代谢、翻译、糖代谢、其他次生代谢产物合成、氨基酸代谢和运输等通路富集。最后，对三个比较组中5个共下调基因和4个共上调基因与不定根形成的关系进行了讨论。转录组分析说明随着母株年龄的增长，不定根形成相关基因表达谱发生变化，导致插穗不定根生根受到抑制。

Abstract: In order to explore the effect of mother plant age on adventitious rooting, the five-day-old stem cuttings of a 1-, 2- and 3-year-old hybrid triploid cultivar Populus “Beilinxiongzhu 2” were used as materials for transcriptome sequencing analysis. Significant differentially expressed genes were screened, and GO enrichment analysis and KEGG enrichment analysis were carried out for differentially expressed genes. A total of 5153 differentially expressed genes were identified. Through GO analysis of three comparative groups, the up-regulated genes were mainly concentrated in mitochondrial respiration, oxidative phosphorylation and flavonoid biosynthesis. The down-regulated genes were mainly concentrated in the synthesis of brassinolide, salicylic acid-mediated signaling pathway and gibberellin biosynthesis. KEGG analysis showed that differentially expressed genes were enriched in signal transduction, energy metabolism, translation, sugar metabolism, biosynthesis of other secondary metabolites, amino acid metabolism and transportation. Finally, the relationship between adventitious root formation with 5 down-regulated genes and 4 up-regulated genes was discussed. The transcriptome analysis showed that the expression profiles of adventitious root formation related genes changed with increased age of the mother plant, which led to the inhibition of adventitious root formation in Populus cuttings.

1. 引言

毛白杨是我国北方的主栽树种之一，在全国木材供应、生态环境建设、城乡绿化等方面具有重要意义。北京林业大学林木遗传育种国家工程实验室科研团队选育出的雄性杨树杂种三倍体新品种“北林雄株2号”，父本为银腺杨(Populus alba × P. glandulosa)，母本为毛新杨(P. tomentosa × P. bolleana)，具有雄株不飞絮、生长迅速、抗风能力强等特点，是解决杨树飞絮问题和用于速生丰产林建设的适宜品种 [1] [2] [3]。

通过无性繁殖的技术进行培育，并将无性繁殖与林木育种相结合是目前杨树栽培的主要道路 [4]。然而无性繁殖容易受到成熟效应(age effect)的影响，从而导致品种退化。成熟效应有时也被称为年龄效应，是指无性繁殖材料的发育阶段对无性繁殖效果的滞后影响，例如随着母株年龄的老化，插穗枝条的生根能力逐渐丧失。在杨树等成熟效应比较明显的树种中，苗木幼年期的速生性削弱，很快进入发育的衰老阶段，生长受到抑制，年生长量降低 [5]。本研究通过对1、2、3年生“北林雄株2号”扦插五天的插穗进行转录组测序分析，分析了年龄对不定根生根状态的影响，为杂交三倍体新品种“北林雄株2号”的插穗繁育提供新的理论基础。

2. 材料与方法

2.1. 实验材料

杂交三倍体新品种“北林雄株2号”无性繁殖实验材料来自山东省聊城市冠县国营苗圃。

2.2. 实验处理

插穗采集：选取生长良好、健壮无病虫害的1、2、3年生母株当年生嫩枝条(直径为0.8~1.2 cm，长度为10 cm)。扦插基质为沙质土：营养土：珍珠岩按2:1:1比例混合。采取扦插五天(D5)的插穗基部以上5 cm范围内部位(1、2、3年生分别取3株重复)，立即放入液氮速冻，于−80℃超低温保存，用于后续转录组测序分析。

2.3. 实验方法

2.3.1. RNA-Seq 技术和Illumina 测序

对1、2、3年生扦插五天的插穗基部进行转录组测序，对样品总RNA进行提取，使用DNase处理杂质DNA，将富集的mRNA打成短片段。然后将这些片段作为模版，用六碱基随机引物合成一链cDNA，之后合成二链cDNA，双链cDNA采用试剂盒进行纯化；纯化后的双链cDNA再进行末端修复、加A尾并连接测序接头和片段大小选择，最后进行PCR扩增；构建好的文库用Agilent 2100 Bioanalyzer质检合格后，使用Illumina HiSeq X Ten测序仪进行测序，产生150 bp的双端数据。

2.3.2. 数据预处理、质量控制与组装

为了得到可供后续分析的高质量reads，需要对高通量测序中产生的原始数据raw reads进行进一步的质量过滤。使用Trimmomatic [6] 软件对数据进行质控并去除接头，过滤掉低质量碱基以及N碱基得到高质量的clean reads。通过Trinity (version:trinityrnaseq_r20131110) [7] 软件paired-end的拼接方法拼接 clean reads得到Transcript序列，根据序列相似性以及长度，挑选出最长的一条作为Unigene。

2.3.3. Unigene 功能注释

采用diamond [8] 软件，取e < 1e−5，筛选序列相似性最高的蛋白，从而得到NR、KOG、GO、Swiss-Prot、eggNOG、KEGG数据库注释功能注释信息。利用HMMER [9] 软件将Unigene与Pfam数据库进行比较，完成对Pfam数据库的注释。

2.3.4. Unigene 定量和差异Unigene 筛选

FPKM (fragments per kilobase of transcript per million mapped reads)是每一百万个片段中映射到每1千个外显子碱基的片段数。FPKM兼顾了测序深度和Unigene长度对片段计数的影响，是目前最常用的估计表达水平的方法。Unigene的FPKM [10] 和count使用bowtie 2 [11] 和eXpress [12] 软件分析得到。通过eXpress软件获取落到各个样本中Unigene的reads数目，使用DESeq [13] R package的estimate Size Factors函数对数据进行标准化，并使用nbinomTest函数计算差异比较的P value和Fold change值。挑选出P value小于等于0.05，Fold change大于等于2的差异表达基因，对差异表达基因进行GO和KEGG富集分析。

2.3.5. 差异Unigene富集分析

GO (Gene Ontology)富集分析可通过注释结果对差异Unigene进行统计学分析和功能描述。统计每个GO term中所包括的差异mRNA个数，并用超几何分布检验方法计算每个GO term中差异Unigene的P value和FDR值，定位差异Unigene最可能相关的GO term。P value小表示差异Unigene在该GO term中出现了富集。

超几何分布检验计算公式为：

$p=1-{\displaystyle \underset{i=0}{\overset{m-1}{\sum }}\frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{(Nn)}}$

Enrichment score计算公式为：

$\text{Enrichmentscore}=\frac{m}{n}/\frac{M}{N}$

KEGG (Kyoto Encyclopedia of Genes and Genomes)是一个整合了基因组、化学和系统功能信息的公共数据库。统计每个pathway中所包括的差异mRNA个数，并用超几何分布检验方法计算每个pathway中差异Unigene的P value和FDR值，定位差异Unigene最可能相关的pathway。P value小表示差异Unigene 在该pathway中出现了富集。相关计算公式参考GO富集分析。通过对差异Unigene进行GO和KEGG富集分析，来研究影响杨树不定根生根过程的生物学功能或者通路。

3. 结果与分析

3.1. 1、2、3年生“北林雄株2号”插穗不定根生根

本研究所用的插穗来源于山东冠县国营苗圃的杨树杂交三倍体新品种 “北林雄株2号”。图1显示来自不同年龄(1、2、3年生)母株的茎插穗在经过5天的栽培后不定根的生根情况。1年生插穗下部观察到明显的较长的不定根，而2年生插穗上观察到较短的不定根，3年生插穗只有根原基出现。此外，1年生插穗呈现绿色，而2、3年生插穗与1年生的相比显示出黄褐色，并且较为干枯。形态学观察表明1年生插穗生根能力最强。

Figure 1. Adventitious roots growth in cuttings from 1-, 2- and 3-year-old hybrid triploid cultivar Populus “Beilinxiongzhu 2” mother plants after five days of cultivation

图1. 1、2、3 年生不同年龄母株杂交三倍体杨树“北林雄株2号”插穗扦插五天后的生根情况

3.2. 数据产出统计

样本双端测序数据来自Illumina平台。鉴于数据错误率对结果的影响，我们采用Trimmomatic [6] 软件对原始数据进行质量预处理，并对整个质控过程中的reads数进行统计汇总。质量预处理步骤：去接头 (Adaptor)，去除低质量Reads，从3’端及5’端以不同方式去除低质量碱基，统计原始测序量、有效测序量、Q30、GC含量，并进行综合评价。本研究共有3个样品，分别为1、2、3年生的扦插五天的插穗，简称为D5_1、D5_2、D5_3，各样本的原始数据Q30分布在 96.06%~96.60%，有效数据量分布在7.03~7.07 G，平均GC含量为44.16 % (表1)。

Table 1. The results of sequencing data quality preprocessing

表1. 测序数据质量预处理结果一览表

3.3. 常见功能数据库注释

在本研究中，3个样品共拼接出Unigene 38,321条，总长度为43,432,678 bp，平均长度为1133 bp。Unigene的数据库注释结果如下：28,984个(75.63%)基因注释到NR库，27,039个(70.56%)基因注释到 eggNOG库，21,268个(55.50%)基因注释到Swissprot库，19,163个(50.01%)基因注释到GO库，15,673个(40.90%)基因注释到KOG库，10,592个(27.64%)基因注释到KEGG库。此外，还有55个(0.14%)基因注释到Pfam库(表2)。其中，注释到NR和eggnog数据库的基因最多，注释到NR、eggNOG、Swissprot和GO数据库的基因均超过了总基因数的50%。

Table 2. Summary for the BLASTx results of Populus “Beilinxiongzhu 2” transcriptome against six databases

表2. “北林雄株2号”转录组数据注释结果统计

图2对各数据库共同注释的基因进行了统计。结果显示，被3个及以上数据库同时注释的基因达到22,760个(59.39%)，其中14,860个(38.78%)基因同时注释到5个及以上数据库，表明基因注释的可信度较高。

Figure 2. Each database annotates the Venn diagram. The number on the bar chart above represents the result of the intersection of databases with corresponding black dots in the matrix below, and the column on the left represents the total number of genes annotated in the databases

图2. 各数据库注释统计图。上方柱状图上的数字代表下方对应矩阵中加黑点数据库的交集，左边的柱子代表各数据库注释到基因总数

3.4. 主成分分析

利用Unigene的表达量进行主成分分析(Principal Component Analysis, PCA)，考察样品分布情况，对样本间关系进行探究或者对实验设计进行验证。PCA可以从不同维度展现样品间的关系。样本聚类距离或者PCA距离越近，说明样本越相似。不同组样本分布在二维的不同区域，同组的样品在空间分布比较集中。由图3可知，3个数据点处于不同象限，可清晰地被PC1 (81%)分离，并进一步被PC2 (19%)分离，说明实验设计合理。此外，1、2年生的插穗样品和3年生的插穗样品之间差距较大，而1年生和2年生的插穗样品间较为相似。

Figure 3. PCA map. PCA figure shows the largest components of variance. PC1 accounts for 81% of the variance, and PC2 accounts for 19% of the variance

图3. 主成分分析图。主成分分析图显示方差最主要的组成部分。PC1表示主成分1 (81% 的方差贡献率)，PC2表示主成分2 (19%的方差贡献率)

3.5. 差异Unigene筛选

将三组数据进行两两比较，并将差异基因以柱形图、维恩图和火山图进行统计和分析(图4)。在D5_2/D5_1、D5_3/D5_1、D5_3/D5_2三组中，差异上调表达基因数目分别为972，1670，1418个，差异下调表达基因数目分别为1292，918，877个(图4(a))。维恩图显示，三个比较组的共同差异表达基因有52个(图4(b))。火山图可以显示差异基因的整体分布情况。图中，红色为显著上调的差异基因，绿色为显著下调的差异基因，灰色为无差异基因。X轴的绝对值越大，差异倍数越大；Y轴数值越大表明差异越显著。如图4(c)~(e)所示，D5_2/D5_1组下调表达基因最多，而D5_3/D5_1组上调表达基因最多。三组中只有D5_2/D5_1组下调表达基因多于上调表达基因。与D5_2/D5_1相比，D5_3/D5_1和D5_3/D5_2 中差异倍数较大(X > 5)、显著性较大(Y > 5)的基因数目明显增多，说明3年生样品与1年生、2年生插穗之间差异更大，与PCA结果一致。

(a) (b) (c) (d) (e)

Figure 4. (a) The number of genes was significantly up-regulated and down-regulated in the comparison groups; (b) Venn Map of Differentially Expressed Genes; (c)~(e) Volcano maps of differentially expressed genes. The differences are reflected in volcano maps, where the gray is non-differential Unigene, the red is up-regulated significant difference Unigene, the green is down-regulated significant difference Unigene, the X axis is log2 Fold change, and the Y axis is log10 P value

图4. (a) 各比较组中显著上调和下调基因数统计；(b) 差异表达基因维恩图；(c)~(e) D5_2/D5_1、D5_3/D5_1、D5_3/D5_2差异基因表达火山图。灰色为无差异基因，红色为显著上调的差异基因，绿色为显著下调的差异基因；X轴为log₂ Fold change，Y轴方向为log₁₀ P value的展示

3.6. 差异Unigene GO富集分析

使用Gene Ontology数据库对三组差异基因进行功能显著性富集分析。分析发现，在D5_2/D5_1、D5_3/D5_1、D5_3/D5_2三个比较组中，差异上调表达基因主要有105、83和73个显著富集条目，差异下调表达基因主要有106、126、101个显著富集条目。

与1年生杨树插穗相比，2年生样品中差异上调表达基因在氧化磷酸化、线粒体电子传递、呼吸链、细胞色素c氧化酶活性等功能富集，差异下调表达基因在油菜素内酯的合成、乙烯激活信号通路、水杨酸介导的信号通路等功能富集(图5(a))。3年生与1年生插穗相比，差异上调表达基因在木葡聚糖代谢、类黄酮生物合成、植物型细胞壁等功能富集，差异下调表达基因在乙烯激活信号通路、赤霉素生物合成、丝氨酸型内肽酶抑制剂活性、蛋白质水解调节等功能富集(图5(b))。3年生插穗与2年生插穗相比，差异上调表达基因在木聚糖代谢过程、乙烯激活信号通路、过氧化氢分解过程、植物型细胞壁等功能富集，差异下调表达基因在有氧呼吸、呼吸电子传递链、蛋白质水解的调节、细胞色素c氧化酶活性等功能富集(图5(c))。综上所述，不同年龄的杨树插穗材料在不定根生根时，呼吸作用、植物激素信号转导、细胞壁合成和蛋白质代谢可能是影响不定根生根的主要因素。

Figure 5. Gene ontology enrichment top 30 term map. (a) D5_2/D5_1; (b) D5_3/D5_1; (c) D5_3/D5_2. In the figure, the X coordinate is GO entry name and the Y coordinate is-log₁₀ P value

图5. GO富集top 30条形图。(a)：D5_2/D5_1；(b)：D5_3/D5_1；(c) D5_3/D5_2。图中X坐标为GO条目名称，Y坐标为-log₁₀P value

3.7. 差异Unigene KEGG富集分析

KEGG是有关Pathway的主要公共数据库，利用KEGG数据库对差异Unigene进行Pathway分析，通过Pathway分析可以找到显著富集差异Unigene的Pathway条目，寻找不同样品的差异Unigene可能和哪些细胞通路的改变有关。

在D5_2/D5_1、D5_3/D5_1、D5_3/D5_2三个比较组中，差异表达基因在KEGG Level 2水平主要注释到以下几个代谢通路：信号转导、能量代谢、翻译、糖代谢、其他次生代谢产物的生物合成、氨基酸代谢、运输和分解代谢。这些结果表明，随着母株年龄的增长，插穗的内源激素、总糖、次生代谢产物、蛋白质等的含量发生了变化，对不定根生根产生了影响(图6(a)~(c))。

图7中，X轴为富集分值，气泡越大的条目包含的差异Unigene数目越多，气泡颜色由紫-蓝-绿-红变化，表明其富集P value值逐渐变小，显著程度逐渐增大。由KEGG pathway富集分析可知，植物激素信号转导(ko04075)、淀粉和蔗糖代谢(ko00500)、苯丙烷类生物合成(ko00940)、类黄酮生物合成(ko00941)、氨基糖和核苷酸糖代谢(ko00520)、核糖体(ko03010)、碳代谢(ko01200)是共同显著富集的代谢途径(图7(a)~(c)、表3)。

在D5_2/D5_1组中，差异基因富集的KEGG通路下调基因明显多于上调基因，而D5_3/D5_1、D5_3/D5_2组中差异基因富集的KEGG通路上调基因明显多于下调基因(表3)。这些结果表明，2年生的插穗和1年生插穗相比，植物激素信号转导、淀粉和蔗糖以及氨基酸代谢、次生代谢产物合成等代谢途径起到了抑制作用，而3年生插穗和1年生插穗相比、3年生插穗和2年生插穗相比，这些代谢途径起到了促进作用。

3.8. 差异共表达Unigene分析

在D5_2/D5_1、D5_3/D5_1、D5_3/D5_2三个比较组中，差异共表达的基因共有52个(图4(b))。其中共下调表达的基因有5个(表4)，共上调表达的基因有4个(表5)，其他基因呈现出先上调后下调或先下调后上调的趋势。

Table 3. KEGG differential gene enrichment pathway

表3. KEGG差异基因富集pathway

Figure 6. Up-regulated and down-regulated differential expressed Unigene in KEGG Level 2 distribution maps. (a): D5_2/D5_1; (b) D5_3/D5_1; (c) D5_3/D5_2

图6. 上下调差异表达Unigene在KEGG Level 2水平分布图。(a) D5_2/D5_1；(b) D5_3/D5_1；(c) D5_3/D5_2

Figure 7. KEGG enrichment top 20 bubble diagrams. (a) D5_2/D5_1; (b) D5_3/D5_1; (c) D5_3/D5_2

图7. KEGG富集top 20气泡图。(a) D5_2/D5_1；(b) D5_3/D5_1；(c) D5_3/D5_2

Table 4. Differentially co-down-regulated genes

表4. 差异共下调表达基因

Table 5. Differentially co-up-regulated genes

表5. 差异共上调表达基因

共下调表达的基因中，TRINITY_DN21087_c0_g1_i1是响应乙烯的转录因子。APETALA 2/ETHYLENE RESPONSE FACTORS (AP2/ERFs)家族是植物最大的转录因子调控家族之一 [14]，家族成员可以分为AP2、ERF、DREB、RAV和Soloist 5个亚家族。其中ERF可通过参与乙烯介导的信号转导途径参与基因表达调节，进而参与植物发育的多个进程。在杨树中，PtaERF003基因表达的上调促进了不定根生长 [15]。共上调表达的基因中，TRINITY_DN20139_c0_g2_i1基因功能为过氧化物酶，应答氧化应激胁迫并参与生长素的代谢。过氧化氢可通过二胺氧化酶增加不定根的生成 [16]，过氧化物酶的增加可导致过氧化氢含量下降，从而抑制不定根的形成。

4. 讨论

杨树年龄效应引起的插穗生根困难是无性林业系发展的一个难题。本文通过对杂交三倍体新品种“北林雄株2号”不同年龄五天插穗的转录组分析，研究了母株年龄对插穗生根的分子调控机制。

由GO功能分析可知，在D5_2/D5_1组中，涉及线粒体呼吸作用的基因出现显著上调，而差异下调基因主要涉及植物激素。D5_3/D5_1组差异上调的基因主要涉及植物细胞壁和类黄酮，下调表达的基因在植物激素和蛋白代谢相关的功能富集。D5_3/D5_2组差异上调的基因多富集在植物细胞壁和激素和信号分子代谢，差异下调的基因主要富集在线粒体呼吸作用和蛋白质代谢。呼吸作用是产生活性氧(reactive oxygen species，ROS)的主要途径，呼吸作用的增强可导致ROS的增多，ROS在植物的生长、发育和在胁迫下的反应以及细胞程序性死亡(PCD)调控过程中是一个重要的信号分子 [16]。油菜素内酯(brassinolide，BR)是属于甾醇类化合物的植物激素，很低浓度的BR就可以诱导植物细胞的伸长和分裂，BR通过促进细胞膜系统质子泵对氢离子的泵出，导致自由空间酸化，使细胞壁松弛从而促进生长。BR还可以抑制生长素氧化酶的活性，提高植物内源生长素的含量 [17]。乙烯(Ethylene)在不同植物中起的作用不同，包括抑制、促进不定根的形成，也可能不起作用 [18]，但已知参与调节生长素和细胞分裂素的合成 [19]。水杨酸(salicylicacid，SA)可以促进绿豆下胚轴不定根的发生 [20]，并与过氧化氢(H₂O₂)和茉莉酸甲酯存在相互作用 [21]。木葡聚糖是存在于所有高等植物细胞壁的具有生物活性的半纤维素多糖，能够控制植物的生长和代谢。类黄酮(flavonoids)可以清除ROS，研究表明类黄酮的含量与不定根形成呈负相关。赤霉素(gibberellin)在许多研究中被证实在不定根生根的起始阶段抑制不定根生成，而在伸长阶段则促进不定根形成 [21]。丝氨酸型内肽酶主要作用是促进蛋白质降解，其活性增强被认为是植物衰老过程中蛋白质含量下降的主要原因 [22]。通过对D5_2/D5_1、D5_3/D5_1、D5_3/D5_2三个比较组的分析可知，与幼龄插穗相比，老龄插穗中上述不定根形成的正调控基因下调，负调控基因上调。

由KEGG pathway富集分析可知，在1、2、3年生的插穗中涉及糖、蛋白质、植物激素和次生代谢产物的基因发生富集。在插穗衰老过程中，部分富集到代谢通路的基因表达量呈先下调后上调的趋势。

以上结果说明，不定根的形成是一个复杂的同时受到内外源多因素调控的过程 [23]。随着母株年龄的增长，杨树不定根生根受到抑制。本研究获得的不同年龄插穗中不定根形成相关基因表达谱变化，可以较好地解释杨树插穗的年龄效应。同时，转录组分析发现的共差异表达基因，可为“北林雄株2号”扦插不定根生根能力下降的问题提供新的解决思路。

5. 结论

雄性杨树杂种三倍体新品种“北林雄株2号”的插穗易受母株年龄的影响。随着杨树母株年龄的增长，不定根形成的正调控基因下调，负调控基因上调，导致插穗不定根生根受到抑制。

致谢

感谢洪汉辉同学为本文提供的帮助，感谢康老师、汪老师、薛老师为本文提供的宝贵修改意见。

基金项目

“十三五”国家重点研发计划林业资源培育及高效利用技术创新重点专项(2016YFD0600403)。

参考文献

NOTES

^*通讯作者。