白桦BpNAC012基因调控白桦木质部发育表达谱分析

郭依萍 刘佳欣 于 颖 王 超 杨传平

（林木遗传育种国家重点实验室，东北林业大学，哈尔滨 150040）

NAC 基因是植物特异性转录因子基因家族之一，其命名是由矮牵牛基因NAM［1］以及拟南芥基因ATAF 和CUC2 的首字母缩略词构成［2］。NAC 转录因子具有多种功能，如参与植物次生长，在细胞分裂和植株衰老中发挥作用，参与激素调控和信号转导［3～4］。植物的次生生长包括维管组织形成、次生细胞壁形成、木质化、细胞程序化死亡以及心材形成等过程［5～6］。近年研究发现多个NAC 基因对细胞次生壁的形成起着正调控作用。例如拟南芥SND1（Secondary wall-associated NAC domain 1），NST1/2/3（NAC secondary wall thicken promoting factor 1/2/3）和VND6/7（Vascular-related NAC domain 6/7）。在拟南芥次生壁正常形成的过程中，SND2/3、MYB103/85/52/54/69/42/43/20 和KNAT7等11 个SND1 转录因子是必需的。抑制SND2/3、MYB103/85/52/54 和KNAT7 的表达，能显著减少纤维细胞次生壁增厚；而SND2/3 和MYB103 的过量表达，则能促进纤维中次生壁的增厚。SND1 在茎的维管束间纤维和木质部纤维中特异表达，过表达SND1促进非厚壁细胞中次生壁的沉积；而抑制SND1 的表达，纤维中缺失次生壁［7～8］。NAC 转录因子对次生壁合成的调控是通过对下游靶基因的调控完成的。例如，VND7 能激活下游转录因子基因及一些参与次生壁形成、细胞壁化学修饰和细胞凋亡等非转录因子基因的表达［9］。

前期研究结果显示，BpNAC012 基因调控白桦（Betula platyphylla Suk.）次生细胞壁的合成［10］。本研究主要以次生生长相关的白桦BpNAC012 基因为研究对象，对其在白桦木质部中调控的下游基因进行研究。通过对过表达、抑制表达转基因株系及对照株系茎部转录组进行分析，鉴定差异表达基因，初步确定BpNAC012基因调控白桦木质部发育的基因表达途径，为白桦的遗传改良奠定基础。

1 材料与方法

1.1 实验材料

分别取3 个月苗龄的转基因白桦和野生型白桦，去除叶片和根，作为转录组测序材料，每个处理3 次生物学重复，每次生物学重复5～7 株（抑制表达株系茎小，7 株可满足建库要求）。WT为野生型株系，OE为基因过表达株系，RS为抑制表达株系。

1.2 建库测序流程

在RNA 样品检测合格后，利用带有Oligo（dT）的磁珠富集总RNA 中的mRNA。随后将mRNA 打断成短片段，以mRNA 为模板，利用random hexamers 合成cDNA 第一链，利用PCR 合成cDNA 第二链，AMPure XP beads 纯化双链cDNA。PCR 扩增获得测序所用文库。为保证文库质量符合测序要求，使用Qubit2.0 对文库进行初步定量，随后使用Agilent 2100 对文库的插入长度进行检测，插入长度符合要求后，使用Q-PCR 方法准确定量文库的有效浓度后（有效浓度＞2 nmol·L-1），进行Illumina HiSeq测序。

1.3 生物信息学分析

利用Corset程序对转录本进行层次聚类，以聚类后序列为参考，进行后续分析。

1.3.1 测序质量评估

通过测序错误率分布检查来检测测序数据的质量，用Qphred=-10log10（e）来衡量，其中测序错误率用e表示。将测序得到的原始测序序列（Sequenced Reads）或者raw reads去除带接头（adapter）的reads；去除N（N表示无法确定碱基信息）比例大于10%的reads；去除低质量reads（质量值Qphred ＜=20的碱基数占整个reads的50%以上的reads）。

1.3.2 转录本拼接

对于无参考基因组的项目，获得clean reads后，采用Trinity 进行转录本拼接。利用比对上转录本的reads数和表达模式对转录本进行Corset层次聚类（https：//code.google.com/p/corset-project/）。对转录本及聚类序列长度分别进行拼接转录本长度分布统计。

1.3.3 基因功能注释及CDS预测

利用生物信息学数据库进行基因功能注释，包括：Nr，Nt，Pfam，KOG/COG，Swiss-prot，KEGG，GO。从比对结果中提取出转录本的ORF 编码框信息，并按照标准密码子表将编码区序列翻译成氨基酸序列（按照5′-＞3'的顺序）；对于没有比对结果或者比对未预测出结果的序列，则采用estscan（3.0.3）软件预测其ORF。

1.3.4 基因表达水平分析

采用RSEM 软件将每个样品的clean reads 往参考序列（ref）上做mapping。进行参考序列比对。比对结果进行统计，得到样品比对的readcount 数目，并对其进行FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced）转换，来进行基因的表达水平分析。

1.3.5 差异表达分析

利用基因表达水平分析中得到的readcount数据进行基因差异表达分析。采用DESeq2 软件，筛选阈值为padj＜0.05 且|log2FoldChange|＞1，对原有假设检验得到的p-value进行校正。

1.3.6 差异基因富集分析

利用GOseq 进行GO 富集分析。以KEGG（Kyoto Encyclopedia of Genes and Genomes）Pathway 为单位，应用超几何检验，找出差异基因相对于所有有注释的基因显著富集的pathway，进行Pathway 显著性富集分析。绘制KEGG 富集散点图，KEGG 富集程度通过Rich factor、qvalue 和富集到此通路上的基因个数来衡量。其中Rich factor指差异表达的基因中位于该pathway 条目的基因数目与所有有注释基因中位于该pathway 条目的基因总数的比值。qvalue 是做过多重假设检验校正之后的Pvalue，qvalue 的取值范围为［0，1］，越接近于零，表示富集越显著。

2 结果与分析

2.1 RNA质量分析

RNA 质量决定着文库的质量，RNA 提取质量分析结果显示如图1，样品质量满足建库测序要求，且总量满足2次以上建库需要。

2.2 表达谱测序结果评估

2.2.1 测序质量评估

测序过程存在错误率，测序错误率分布检查可以反映测序数据的质量。一般情况下，单个碱基位置的测序错误率应该低于1%。从表1可以看出，本文库测序错误率均低于0.1%，说明测序质量符合要求。

从表2 可以看出，6 个文库的G/C 含量符合正常规律。接头序列、未知序列、低质量序列等不合格序列的百分比低于2%，cleanread 应占原始数据的90%以上，这样的数据为合格数据。从表1 和附图1 可以看出，该文库cleanread 占原始数据的比例均符合后续分析的要求。

2.2.2 转录本拼接

对转录本及聚类序列长度分别进行统计，结果见表2、附表1 和附图2。从图表中可以看出，本次分析转录组的基因长度分布符合进一步分析的要求。

表1 数据产出质量情况一览表Table 1 List of data output quality

表2 拼接长度频数分布情况Table 2 List of frequency distribution of splice length

2.3 基因功能注释

为获得基因功能全面信息，对本研究文库获得的数据进行7 大数据库基因功能注释，包括：Nr，Nt，Pfam，KOG/COG，Swiss-prot，KEGG，GO 数据库。本项目数据在7 大数据库中的注释成功率情况见表3。从注释成功率可以看出，获得功能注释的基因比率符合后续基因功能分析的要求。

表3 基因注释成功率统计Table 3 Gene annotation success rate statistics

2.4 基因表达水平分析

2.4.1 参考序列比对

参考序列比对统计结果见附表2。在RNAseq 技术中，FPKM 分析是目前最为常用的基因表达水平估算方法。所以我们将readcount数进行了FPKM转换用于下一步差异表达分析。

2.4.2 基因差异表达分析及差异基因筛选

本研究3个处理间共获得差异表达基因1 476条。从图2 中可以看出，与对照相比，过表达株系OE 中上调的基因有627 条，下调的基因有229 条，抑制表达株系中上调的基因有299条，下调表达的基因有207 条。说明在转基因株系中，BpNAC012基因的表达影响一系列的基因表达变化，进而调控木质部发育性状的变化。

2.5 差异基因GO富集分析

GO 功能显著性富集分析给出与基因组背景相比，在差异表达基因中显著富集的GO 功能条目，从而给出差异表达基因与哪些生物学功能显著相关。结果如表4 所示。差异表达基因主要富集在转录因子相关的GO 分类中。说明转录因子在本实验样本的基因表达差异及性状差异中起到重要作用。统计被显著富集的各个GOterm 中的基因数，以柱状图的形式展示，如图3 所示。结合差异表达基因数量分析，从图中可以看出，过表达BpNAC012 能够调控更多的基因表达变化。而抑制表达BpNAC012 更多的是影响蛋白修饰和转运类基因的表达。

表4 样品中差异基因富集的Gene Ontology分类Table 4 Gene Ontology classification of differential gene enrichment in samples

2.6 差异基因KEGG富集分析

差异基因KEGG 富集散点图（见图4）图形化展示了KEGG 富集分析结果。从图4中可以看出，涉及受体信号通路、营养代谢、氨基酸合成及苯丙烷生物合成相关代谢通路基因比较富集。其中苯丙烷生物合成，是次生细胞壁合成的关键代谢通路，因此说明BpNAC012 的表达，调控了次生细胞壁合成的代谢通路。

2.7 表达差异最大的前十位基因分析

对过表达株系相对于野生型表达量上调最大的10个基因和抑制表达株系相对于野生型表达量下调最大的10 个基因进行分析（见表5～6）。可证明这些基因强烈受BpNAC012 的调控，并在过表达或抑制表达株系发育的生理程序中起着重要作用。

2.8 BpNAC012 调控木质部发育及次生壁合成相关基因的表达

在BpNAC012过表达株系中，鉴定到一些和纤维素合成相关的基因（见表7），和一些细胞壁相关酶类。同时，鉴定了一些与木质素合成及沉积相关的基因，和木质部发育相关的基因。这些基因可能受BpNAC012 的诱导，参与过表达转基因白桦木质部发育及次生细胞壁的合成。

在抑制表达株系中，鉴定到一些木质部发育及细胞壁合成相关基因下调表达（见表8），如一些与纤维素合成相关的基因下调表达。

2.9 BpNAC012 调控其它转录因子基因的差异表达

在过表达和抑制表达株系中，分别鉴定到了与野生型对照差异表达的转录因子（见表9～11）。其中，与野生型对照相比，在过表达株系中鉴定了97 个上调表达的转录因子，15 个下调表达的转录因子，在抑制表达株系中，鉴定了42个上调表达的转录因子，15个下调表达的转录因子。

表6 抑制表达株系中差异表达倍数最高的前10位基因Table 6 The top 10 genes with the highest differential expression folds in the suppressed expression lines

3 讨论

BpNAC012 对白桦次生细胞壁合成的调控是通过下游基因表达变化实现的，与对照相比，过表达株系OE 中上调的基因有627 条，下调的基因有229 条，抑制表达株系中上调的基因有299 条，下调表达的基因有207 条。两个转录组中共有140条基因差异表达。说明在转基因株系中，Bp⁃NAC012 基因的表达影响了一系列的基因表达变化，进而调控次生壁的合成，这些差异表达基因可能参与白桦次生细胞壁及木质部形成的调控。

在转基因株系中，表达变化最明显的基因，可能对于BpNAC012 转基因造成的变化起着重要作用，对过表达株系相对于野生型表达量上调最大的10个基因和抑制表达株系相对于野生型表达量下调最大的10个基因进行分析。UTP—glucose-1-phosphate uridylyltransferase，EC number 2.7.7.9，同义名称UDP glucose pyrophosphorylase，UDPG pyrophosphorylase，UGPase［11］，是 一 种参 与 糖 代 谢 的酶。它以葡萄糖-1-磷酸和UTP 为底物合成UDP

葡萄糖；UTP-葡萄糖-1-磷酸尿苷酰基转移酶是在所有三个领域（细菌、真核生物和古细菌）中发现的酶，因为它是糖原生成和细胞壁合成的关键参与者［12］。UDP-葡萄糖焦磷酸化酶基因过表达可提高黄麻（Corchorus capsularis L.）的纤维素含量［13］。说明UGPase 在营养生长、细胞壁形成及纤维素合成中起重要作用。在BpNAC012抑制表达株系中，该基因显著下调表达，可能与抑制表达株系纤维素含量下降有关。值得注意的是No lysine kinase 1 isoform 1，同时出现在过表达株系上调倍数前十中，和抑制表达株系下调倍数前十中，说明该基因与BpNAC012 的表达密切协同，受BpNAC012 的强调控。No lysine kinase 1（WNK）具有丝氨酸/苏氨酸激酶的催化赖氨酸的独特特征，因而得名。WNK1 在植物体内离子转运的调节、激素对渗透和离子胁迫的调节、昼夜节律和开花时间、蛋白质转运中、离子通道调控等过程中起重要作用［14］。因此说明BpNAC012 表达的改变，在某种程度上可能影响了转基因株系的信号传导途径。

表7 过表达株系木质部发育相关上调表达基因Table 7 Genes related to up-regulated expression of xylem development in overexpressing lines

表8 抑制表达株系木质部发育相关下调表达基因Table 8 Genes related to down-regulated expression of xylem development-inhibited expression lines

表9 过表达株系中差异表达转录因子Table 9 Differentially expressed transcription factors in overexpressing lines

表10 抑制表达株系中差异表达转录因子Table 10 Inhibition of differentially expressed transcription factors in expressing lines

过量表达BpNAC012可促进木质部发育，次生壁增厚。因此本研究进行差异表达基因分析时，关注了和这些性状相关的基因。CesA是细胞壁生物合成过程中催化纤维素合成的关键酶。该蛋白参与纤维素合成代谢途径，参与次生细胞壁的形成，是木质部细胞壁增厚所需［15］。本研究中，Bp-NAC012 过表达株系中鉴定的2 条CesA 相关基因的表达量比对照高，而在抑制表达株系中鉴定的1条CesA 基因的表达量比野生型对照低。表明Bp-NAC012 通过增强过表达株系中CesA 基因的表达促进纤维素生物合成。细胞壁结构的调节在植物的生长和分化中起着关键的作用。Xyloglucan endotransglucosylase/hydrolase family protein（XTHs）催化木葡聚糖链的分解和分子连接，在细胞壁的松弛和重新排列中起作用。本研究在过表达株系中鉴定了2 条上调表达的Xyloglucan endotransglucosylase/hydrolase family protein，综合分析结果说明，该基因受BpNAC012 的诱导，在白桦木质部发育中起重要的作用。endo-1，4-beta-glucanase 家族基因在植物生长和发育的多方面具有功能，能促进植物生长和侧根发育，增强细胞扩展。杨树的endo-1，4-beta-glucanase 是木质部组织特异的，在纤维素生物合成中起作用，能够促进植物的生长［16］。本研究中，在抑制表达株系中，鉴定了1 个下调表达的endo-1，4-beta-glucanase 基因，综合分析，说明endo-1，4-beta-glucanase 随着BpNAC012表达的下调而下调，受其调控，和抑制表达转基因株系的木质部发育及纤维素合成受阻相关。

本研究在过表达和抑制表达株系中，分别鉴定到与野生型对照差异表达的转录因子AP2-EREBP、Tify、WRKY、Orphans、bHLH、NAC、GRAS、C2H2、SWI/SNF、SET、MYB、HSF、AUX/IAA、ARF、Sigma70-like、mTERF、G2-like、C3H、HB、FAR1、CPP、ABI3VP1、TRAF、GNAT、C2C2-Dof和LOB 等。在过表达株系中，差异表达数量最多的转录因子是AP2-EREBP 类转录因子，而且全部是上调表达；其次是Tify 类转录因子，WRKY 类转录因子；此外bHLH、NAC、AUX/IAA 和MYB 转录子表现了不同的表达模式。在抑制表达株系中，上调表达的转录因子数量明显少于过表达株系，其中差异表达数量最多的，依然是AP2-EREBP 类转录因子，并且有1 个下调；上调表达的WRKY、Tify和Orphars类转录因子数量也明显下降；此外，在过表达株系中都上调的转录因子，在抑制表达株系中出现了下调表达，如C3H 和C2H2。在过表达株系中，GRAS、SWI.SNF 和HSF 类转录因子都有大于3条的上调表达，而在抑制表达株系中没有检测到；在抑制表达株系中，TRAF、C2C2-dof 和GNAT类转录因子都下调表达，而在过表达表达株系中没有检测到。以上结果说明，在BpNAC012调控的转基因株系中，存在多个层级的转录调控，BpNAC012能够调控其他转录因子的表达；在过表达和抑制表达株系中，存在不同的、复杂的调控模式，一些特异转录因子，在各自的性状变化中起着相应的调控作用。