王海蓉，齐 明，叶金俊，何贵平

（1.浙江省遂昌县林业技术推广总站，浙江 遂昌 323300；2.中国林业科学研究院亚热带林业研究所，浙江 杭州 311400；3.浙江省林木育种技术研究重点实验室，浙江 杭州 311400）

杂交育种是杉木Cunninghamia lanceolata遗传改良的主要途径[1]。在杉木杂交育种研究中，生长性状的杂种优势现象已有报道[2]，但对杂种的生长、材质等性状分离的分子机理研究的较少。杉木是我国重要的工业用材树种，按照短周期工业用材的要求，除了速生优质外，另外就是要收获期一致。杉木遗传改良进入了第三代，杉木育种群体是多世代与高强度选择的产物，杂合体居多。福建省已观察到杉木高世代种子园的后代生长明显分离，这明显不利于杉木工业用材林的收获。杉木无性系和亲本性状互补的亲本双系杂交种子园有着独有的优势，但是无性系对立地条件要求较高，且存在“C”效应（即年龄效应和位置效应）。杉木双系种子园有独有的优势，但是这必须建立在优良杂交组合的子代分离研究的基础上。通过选择性状互补的优良亲本进行控制杂交，再经苗期评价，可以选育出生长节律大致相同的杂交新品种。

近年来，我国对一年生杂种实生苗分离状况的研究报道多在李[3]Prunus salicina，柚[4]Citrus maxima，枣[5]Ziziphus jujuba，苹果[6-7]Malus pumila和梨[8-9]Pyrusspp.等园艺植物上，研究多在1～ 3 年生以上的生长势或亲本与童期的关系方面，结果发现，一年生植物枝干、叶片、针刺等性状可作为早期鉴定、早期选择的主要指标。在用材林研究上，邱有德等[10]发现巨桉Eucalyptus grandis控制优质大径材的EgLBD基因及其作用机制，但杉木杂交组合的子代分离对收获期的影响尚未引起同行的注意。不少学者在杉木研究中借助RNA-seq 技术来研究杉木纤维性状发育的形成机制和木材形成层活动的机理[11-13]。

目前，尚未发现有人采用转录组技术研究杉木优良杂交组合内生长性状分离的分子机理。本研究在齐明等[14-15]对杉木杂交子代试验林生长量前期试验调查的基础上，利用无参转录组测序技术，在基因差异研究的基础上，以基因差异表达分析为切入点，以杉木一个优良组合的超亲、低亲两组杂交子代及其亲本为研究材料，展开杂种子代和亲本两两间的比较分析，来探讨杉木杂种生长性状分离的分子机理，揭示杉木优良杂交组合高生产力的原因。

1 材料与方法

1.1 研究材料

研究材料取自浙江省遂昌县湖山乡大桥村的杉木杂交子代试验林，1996 年重复制种，1997 年育苗，1998年春营造试验林，试验设计为完全随机区组，19 个处理，8 次重复，纵向单行5 株小区，每个家系40 株参加试验，造林株行距为2 m×2 m。2005 年11 月进行生长量调查[14-15]。2017 年6 月完成叶片取样，叶样来自样树顶部当年生嫩枝和针叶，在研究其遗传变异的基础上，选取表现最优的一个杂交组合（龙15×1339）中的超高亲子代（生长超过最优亲本的子代，HF1）和超低亲子代（生长低于中亲的子代，LF2），和同龄双亲（龙15 和1339）。3 个生物学重复，共12 个样本，参与测序分析。

1.2 研究方法

1.2.1 遗传变异及生长性状的分离研究 本试验分析了研究材料的遗传变异、生长等性状的分离情况、以及优良组合龙15×1339 的生产力，并将其与测序结果联系起来，进行研究。

1.2.2 cDNA 文库准备及RNA-seq 测序 文库构建、无参转录组测序，以及随后的unigene 功能注释、reads 富集与分类和基因表达等内容的分析方法参见文献[16-20]。

1.2.3 序列比对及差异表达基因与杂种生产力间的关系分析 通过Illumina Hiseq 4000 测序获得的转录组测序数据，需要经过生物信息学处理与分析[16-20]，方能获得有意义的结果。除了样本针叶的采集和研究方案的制定外，所有测序和初步分析（包括GO 富集与分类、KEGG 的富集等分析）委托杭州联川生物技术股份有限公司完成。在此基础上，对注释的基因，GO 和KEGG terms 富集结果，考查terms 间和terms 内基因的平衡状态，并追踪几个显著KEGG terms 内，四个比较组间基因上调或下调规律，以揭示杂种生产力的高产和低产的原因。unigene 的表达量FPKM（Fragments per kb per Million fragments），计算公式：

式中，FPKM是某个基因（A）的表达量，C是唯一比对到基因A 的片段数，N是唯一比对到所有unigene 的总片段数，L为unigene A 的碱基数。

样本比较组的GO-term 或KEGG-term 间基因处于平衡或不平衡的状态指数k=sum（第i个下调term 的基因数/第i个上调term 的基因数）/n，上式中i取1～n，标准差按常规公式计算。当k等于或接近于1 时，样本组间的差异表达基因处于平衡状态，反之则处于不平衡状态。

2 结果与分析

2.1 研究材料的生长变异性和生长分离情况

对9 年生时试验林生长性状的遗传变异进行调查分析[14-15]，结果见表1。

表1 研究群体和优良杂交组合（龙15×1339）9 年生的生长性状大小和分离信息[14-15]Table 1 Growth traits and segregation of C.lanceolata in a 9-year cross hybrid stand and cross combination 15× 1339 in Suichang

由表1 可知，优良杂交组合的平均生长性状水平比群体的平均水平要高；优良杂交组合的生长性状的变异幅度略小于全林试验结果，这与参试子代样本数的多少有关（组合样本数/全林子代样本数=36/602）；除材积的表型变异性，群体水平略高于杂交组合外，在树高、胸径上均为优良组合龙15×1339 超过群体平均变异水平，说明本试验中的优良组合龙15×1339 中子代的分离变异还是较大的。

另外，HF1 的平均胸径是LF2 胸径的134.81%～ 138.32%；HF1 平均树高是LF2 的134.60%～ 150%。采样的子代植株、亲本龙15 和1339 生长发育均正常。

2.2 测序数据质量和基因表达概况

杉木4 个处理12 个样品，测序测得的原始序列（Raw reads）介于（4.3E+07）～（5.2E+07）nt；12 个样本的Clean reads 分别介于（4.3E+07）～（5.2E+07）nt。Clean reads 有效数据占原始Raw reads 的比例在98%以上。二代测序中，每测一个碱基会给出一个相应的质量值，这个质量值是衡量测序准确度的。Q20 代表错误率为1%；Q30 代表错误率为0.1%。测序中的Q20 与Q30 则表示质量值≥20 或30 的碱基所占百分比。12 个样本的Phred数值大于Q20 和Q30 的碱基占总体碱基的百分比，Q20 和Q30 的百分比分别介于98.20%～ 98.84%和95.07%～96.4%。原始测序序列中，碱基G 和碱基C 的数量总和占总碱基数的百分比（GC%）介于44.20%～ 45.08%之间。

综合以上几个测序质量评价指标，说明12 份样品的测序质量较高，能保证后续研究和满足后续数据分析的要求。

2.3 测序数据的拼接和基因注释结果

测序数据的拼接结果是，平均GC%达40.78%；当读长数达50%时，该读长的长度为1 214nt。综合其它项的结果，可以得出测序数据的组装拼接的结果很成功。

对Clean reads 在六个数据库进行BLASTX 分析，将获得的基因分别在Swiss-prot，Nr，Pfam，KEGG，KOG和GO 中注释，结果杉木亲代和子代共注释有80 171 个基因。

杉木为无参测序，测得的数据通常采用blastX 序列比对，各数据库所注释总的基因比率大于100%，这表明有的基因同时在不同的数据库得到了重复注释。

测序数据质量、基因表达测序、数据的拼接和基因注释结果更详细的情况参见文献[15]。

2.4 差异表达基因的GO 富集与分类分析和杂种子代的生产力

本研究选择4 个样本组，HF1 VS P1（HF1 与龙15 相比），HF1 VS P2（HF1 与1339 相比），LF2 VS P1（LF2 与龙15 相比），LF2 VS P2（LF2 与1339 相比），对所有的DEGs（differentially expressed genes）进行GO 功能富集分析，从中挑选出基因表达量显著的GO terms 进行GO 分类分析作图。GO 功能分类体系中有参与生物过程（BP）、细胞组分（CC）以及分子功能（MF）3 个大类，43 或45 个小类别。HF1 VS P1 中，GO的功能富集，前10 名的结果是：胁迫反应（k=3）>糖结合（k=9/0）>赤霉素介导的信号通路（k=2）>受体活性（k=20）>防御性反应（k=2）>脱落酸活性（k=5）>蛋白丝氨酸/苏氨酸磷酸酶活性（k=9/0）>转录调控区域DNA结合（k=9/0）>节律过程（k=1）>过氧化物酶活性（k=0.33）。这10 个GO terms 内，表达基因不均衡。

图1 HF1 VS P1 差异表达基因的GO 分类Figure 1 GO classification of differentially expressed genes in HF1 VS P1

比较组LF2 VS P1 的GO 的功能富集结果，前10 名排列顺序为：细胞壁组织（k=0.127 7） >细胞外区域（k=0.145 2） >不可或缺的膜（k=1.213 7）>受体活性（k=2.387 1）>蛋白激酶活性（k=2）>纤维素分解过程（k=0.142 9）>纤维素酶活性（k=0.142 9）>碳水化合物代谢过程（k=0.347 8）>转移酶活性、转移酰基以外的氨基酰基（k=0.428 6）>纤维素生物合成的过（k=0.2）。

图2 是LF2 VS P1 的GO 富集与分类分析结果，与图1 一样，仍以龙15（P1）为参照物，LF2 VS P1 中不同的GO terms 中参与新陈代谢的基因数目不同；针对大多数GO terms，基因下调和上调有着相似的态势。同时在所有的GO terms 上存在显著的上调基因，LF2 VS P1 中GO terms 上基因分布趋于均匀、平衡。

图2 LF2 VS P1 差异表达基因的GO 分类Figure2 GO classification of differentially expressed genes in LF2 VS P1

比较图1 与图2 发现，HF1 VS P1 和LF2 VS P1 的分析结果不一致，其一表现在GO terms 层面，富集分析前10 名的GO terms 中，仅有活体受性一个GO term 相同，其它terms 均不相同；其二表现在GO terms 内上调基因和下调基因数目的不均衡性上，HF1 VS P1 的不平衡结果，比LF2 VS P1 比较组的结果更甚。在HF1 VS P1中，差异表达基因在不同GO terms 上分布处于不均匀、不平衡，而LF2 VS P1 中，整个系统间差异表达的基因数相差不大，平衡系数接近于1.0。由此得出，HF1 生长优于LF2 正是HF1 的基因系统处于非均匀、非平衡状态；而LF2 生长慢则是由于LF2 的基因系统趋向于均匀、平衡状态。这是杉木杂种生长分离的内在遗传基础。

本研究还研究了HF1 VS P2 和LF2 VS P2 比较组GO 的富集与分类研究结果。HF1 VS P2 和LF2 VS P2 的分析结果与HF1 VS P1 和LF2 VS P1 基本一致，尽管参照对象变了，但HF1 和LF2 其系统内差异表达的基因数目所处的状态并没有改变。HF1 和LF2 基因所处的状态、平衡参数的统计分析结果列于表2。

表2 杉木转录组组内差异基因表达平衡与否的状态分析Table 2 Differentially expressed genes in different group

从表2 中可以发现，HF1 的差异表达基因的分布处于不均匀、不平衡状态，LF2 的差异表达基因的分布趋于均匀、平衡状态。就KEGG terms 的富集结果，我们使用了更多的KEGG terms（上调/下调基因）资料，从系统的角度，探讨了全部资料的KEGG terms 的富集结果与杂种生长性状分离间的关系，计算其平衡系数，结果也列于表2。由表2 可知，255 个KEGG terms 中，相对优良亲本，HF1 的平衡系数k>1.2；而LF2 的平衡系数k十分接近于1.0。

2.5 差异表达基因的Pathway 显著性富集分析与子代的生长优势

通过KEGG 数据库中的Pathway 富集分析，确定DEGs 参与的主要生化代谢途径和信号转导途径等，结果显示，子代/亲代4 个比较组中，差异表达的基因分布在74～ 130 条Pathway 中。表3 为处理比较组HF1 VS P1，HF1 VS P2，LF2 VS P1 和 LF2 VS P2 富集结果中最显著的前10 条Pathway。由表3 可见，HF1 VS P1 和HF1 VS P2富集结果接近；LF2 VS P1 和LF2 VS P2 的富集结果接近。但在相同的参照物下，HF1 与LF2 的KEGG 富集结果差异很大：LF2 VS P1 与HF1 VS P1 仅有两个KEGG terms（苯丙氨酸代谢和亚油酸的新陈代谢）相同，其它八个KEGG terms 不相同；同时相同的KEGG term 内，上调基因数与下调基因数不同。LF2 VS P2 与HF1 VS P2间的比较结果：也有两个KEGG terms（苯丙素的生物合成和苯丙氨酸代谢）相同，其它八个KEGG terms 不相同，同时在相同的KEGG terms 内，上调基因数与下调基因数相同或不同。

表3 4 个比较组的KEGG 富集结果中前10 个代谢途径（表中p<0.01）Table 3 KEGG enrichment and top 10 metabolic pathways of four group samples of parents and cross seedlings

综合以上Pathway 富集分析结果可以发现，HF1 和LF2 比较组的遗传差异，表现在两个层面上：其一是KEGG-term 间的不平衡；其二是KEGG-term 内，上调/下调的基因不均衡。这与GO 的富集分析结果基本一致，这正是HF1 生长优于LF2 的原因。但要揭示杉木杂交组合生长性状分离的内在机制，最好在HF1 和LF2 比较组间，选择相同的富集显著的terms，然后再比较terms 内基因是上调表达，还是下调表达。

按照这一思路，我们比较了LF2 VS P1 与HF1 VS P1，HF1 VS P2 与LF2 VS P2 几个相同的KEGG terms 基因表达情况，列于表4。

表4 相同terms 下基因的差异表达情况Table 4 Differential expression of genes within same terms

由表4 可知，在LF2 中，基因一般仍表现为下调表达或沉默，而在HF1 中上调表达，只有亚油酸代谢中，在LF2 中下调表达，在HF1 中上调表达，这是解释杂种生长分离的原因之一，但它属于terms 间的差异。

3 结论与讨论

杉木杂种高生产力的原因是其内在基因系统中基因分布不均匀、不平衡的原因。杉木优良组合中子代生长会出现分离：HF1 生长比LF2 生长快，是由于HF1 GO terms 和KEGG terms 层面不平衡，以及GO terms 和KEGG terms 内的基因表达数量，以及其内上调/下调的基因分布处于不均匀、不平衡状态；LF2 GO terms 和KEGG terms的差异基因的数量，以及其内上调/下调的基因分布趋于均匀、平衡状态（图1 和图2）。杉木杂种生产力的高低符合耗散结构理论：不平衡体系，能量大，信息多，生产力也最高。

杉木杂种生长性状存在分离现象，对收获期不利。GO 和KEGG 富集分析，揭示有大量的基因参与了生长过程。标记辅助选择MAS 研究工作量太大，成本高，技术复杂。利用分子生物学技术，对育种群体展开遗传多样性研究，选择表达基因互补的亲本，进行杂交组配，形成杂合体子代，这样降低杂种生长性状的分离，以达到收获期的大体一致，是一条可行的技术路线[21-22]。