转录组分析二马牙和长脖类型林下参表型差异

李舒欣, 张浩, 郑厚胜,2, 郑培和, 逄世峰, 许世泉*

(1.中国农业科学院特产研究所， 长春 130112; 2吉林农业大学中药材学院， 长春 130118)

人参(PanaxginsengC.A. Mey.)被誉为“百草之王”，已有几千年的药用历史，几乎应用于中医各科临床，现代药理学也证明人参具有抗疲劳[1-3]、抗肿瘤[4-5]、降血糖[6-8]和提高免疫力[9-10]等诸多功效。由于人参生长所需的生态环境遭到破坏，加之人们过度采挖，野生人参资源几近枯竭，当前来源主要为人工栽培(园参)和仿生栽培(林下参)2种途径。与园参相比，林下参经济价值更高，不仅对生长年限有要求，更注重五形六体[11]，其品质也可与野生人参媲美。但林下参种植存在投资多、回报周期长等问题，因此林下种植人参的选种问题极为重要。

由于人参驯化时间短，群体间遗传多样性高，野山参的形态特征也呈现多样化。判断林下参外观品质包括芦(根茎)、艼(芦上不定根)、体(主根)、纹(主根肩上的环纹)、须(支根上生长较细的根)[12]，其中芦头包括雁脖芦、二马牙、堆花、缩脖芦等，艼有枣核艼、毛毛艼、艼变等，体有灵体、疙瘩体、顺体、过梁体、笨体、横体等不同形态，纹讲究细而深，须讲究柔韧性，有弹性及珍珠点。我国栽培的人参在培养驯化过程中形成了马牙型和长脖型两大系统七个品系[13]，马牙系统包括大马牙、二马牙，长脖系统包括长脖、圆膀圆芦、竹节芦、线芦、草庐。二马牙和长脖类型差异较大，芦、艼、体、纹、须均存在差异[14]：二马牙型林下参芦较短粗壮，芦碗大而明显，主根较粗长，多顺体，产量较长脖高；长脖型林下参芦细长，芦碗较小，体长而顺直，多顺体，少笨体。马小军[15]利用RAPD和AFLP技术分析5种农家品种，与形态学分类结果有一致性，二马牙与长脖亲缘关系较远。虽然已有研究表明这2个类型人参皂苷含量存在一定差异，但是其化学成分和形态差异的分子机制尚不明确[16]。

高通量测序可以高效、快捷地获得基因序列，覆盖全部转录本，在挖掘功能基因中得到广泛应用[17]。目前采用转录组测序技术已经在人参皂苷生物合成[18-22]、抗病基因筛选[23-24]、连作障碍[25-26]等研究中取得丰硕成果，但关于不同类型人参种质资源的转录组分析尚未见报道。本文通过比较二马牙与长脖的表型特征以及转录水平的差异，筛选调控林下参形态特征的差异基因，为林下参的品种选育以及揭示林下参表型差异的分子机制提供理论参考。

1 材料与方法

1.1 试验材料

供试材料选自同一林下环境13年生二马牙(E)和长脖(C)农家品种，移栽至吉林省吉林市左家镇，由中国农业科学院特产研究所许世泉研究员鉴定。试验地为阔叶林，北坡，坡度5°左右，土壤类型为黑色腐殖土。

样品采集：第二年人参生长至14年生开花后期，采集叶片(LD)、芦头(T)和芦头以下5 cm主根(G)用于RNA-seq，设置三个生物学重复。再选择三株长势一致的人参，组间随机选取三株根部样品，测量9种人参皂苷Rg1、Re、Rf、Ro、Rb1、Rc、Rb2、Rb3、Rd的含量，45 ℃烘干，磨成粉，参照逄世峰等[31]方法测定人参皂苷含量。

1.2 转录组测序

1.2.1文库构建及质控 文库构建委托百迈客生物科技有限公司完成。用Trizol试剂盒(Invitrogen,USA)分别提取18个样品的总RNA，用Nanodrop 2000 与 Agilent 2100 检测总 RNA 质量、纯度和完整性。样品检测合格后构建文库：①用带有OligodT的磁珠富集mRNA；②加入Fragmentation Buffer将mRNA随机打断；③以mRNA为模板，用六碱基随机引物(random hexamers)反转录第一条cDNA链，合成第二条cDNA链，利用AMPure XP beads纯化cDNA；④纯化的双链cDNA再进行末端修复、加A尾并连接测序接头，然后用AMPure XP beads进行片段大小选择；⑤最后通过PCR富集得到cDNA文库。

文库构建完成后，使用Qubit 2.0荧光计和 Agilent 2100检测插入片段在基因上的分布以及插入片段长度，评估mRNA片段化的随机性、mRNA的降解情况、插入片段长度的离散程度，使用QPCR方法准确定量文库的有效浓度(>2 nmol·L-1)，以保证文库质量。

1.2.2测序及数据分析 库检合格后用Illumina Hiseq 2000平台进行高通量测序。

测序数据首先用q20质控对原始Reads进行过滤，去除接头N比例大于10%以及质量值Q≤10的碱基数占整条Read50%以上的Reads，得到Clean Reads。然后将各样品的Clean Reads用HISAT2软件[27]与指定的参考基因组[28]进行序列比对，利用String Tie[29]对比对上的reads进行组装得到Mapped Data。

1.2.3测序数据分析 使用BLAST软件将Mapped Data在 NR、GO、KOG、KEGG 和 Swiss-Prot 等数据库进行比对注释，用FPKM方法计算基因的表达量。利用百迈克云平台对数据进行分析，设置错误率FDR≤0.1、绝对值 log2FC≥2，将两类型对应部位进行比较，筛选差异表达基因(differentially expressed gene，DEGs)，对其进行GO及KEGG等富集分析，深入挖掘关键通路数据。

1.3 实时荧光定量PCR

筛选7个芦头部位代表性差异表达基因，用实时定量聚合酶链(qRT-PCR)技术检测其表达量检测。使用scaffold128323为内参基因，根据差异表达基因核苷酸序列设计荧光定量PCR引物(表1)。采用反转录试剂盒(TUREscript 1stStand cDNA SYNTHESIS Kit，Aidlab公司)得到cDNA，-80 ℃保存。使用荧光定量PCR仪(analytikjena-qTOWER2.2型，德国)进行荧光定量PCR(10 μL反应体系：5 μL 2×SYBR©Green Supermix、正反引物各0.5 μL、1 μL cDNA模板、3 μL ddH2O；反应条件：95 ℃ 3 min；95 ℃ 10 s，58 ℃ 30 s，39个循环；熔解曲线60～95 ℃，升温程序1 ℃，持续4 s)。按2-△△Ct相对定量法[30]计算相对表达量，设计 3 个生物学重复。

表1 定量PCR引物信息Table 1 Primers used for qRT-PCR analysis

2 结果与分析

2.1 两种类型林下参的形态特征分析

2.1.1根部表型差异 两种类型林下参的根部性状如图1所示：两类型主根长无显著差异，二马牙的芦头直径及根质量均极显著高于长脖，而芦头的长度极显著低于长脖。因此，14年生同环境的二马牙根部较长脖更加粗壮，产量更高，较同年生长脖短，长脖的芦头细长。观察所有样参发现，二马牙的须根较长脖更加发达。

注：**和*表示两个品种间差异在P<0.01和P<0.05水平具有显著性。Note：** and * indicate significant difference at P<0.01 and P<0.05 levels, respectively.图1 两种类型林下参根部性状差异比较Fig.1 Comparison of root traits of two types of forest cultivated Panax ginseng

2.1.2地上部分差异 分析两种类型林下参地上部分的农艺性状差异(图2)发现，二马牙的果柄直径极显著大于长脖，茎高显著大于长脖，其余4项指标差异不显著。说明二马牙较长脖地上部分也更粗壮，利于光合产物的积累。

注：**和*表示两个品种间差异在P<0.01和P<0.05水平具有显著性。Note：** and * indicate significant difference at P<0.01 and P<0.05 levels, respectively.图2 两种类型林下参地上部分农艺性状比较Fig.2 Comparison of aerial portion traits between two types of forest cultivated Panax ginseng

2.2 两种类型林下参人参皂苷含量差异

人参皂苷含量如图3所示，二马牙Rg1、Re、Ro、Rb1、Rc、Rb2、Rb3人参皂苷含量均高于长脖，而Rf、Rd低于长脖，总人参皂苷含量二马牙高于长脖，但两种类型间无显著差异(P>0.05)。

图3 两种类型林下参人参皂苷含量Fig.3 Ginsenoside content between two types of forest cultivated Panax ginseng

2.3 两种类型林下参基因表达水平

对18个样品进行高通量测序，共获得123.70 Gb Clean Data，各样品Clean Data均达到6.09 Gb，Q30碱基百分比在92.07%及以上，与指定的参考基因组比对，效率94.29%到95.37%不等。基于比对结果，发掘新基因24 835个，其中22 087个得到功能注释。差异基因表达量主成分分析(图4)显示，两类型叶片聚合在一起，二马牙的主根芦头聚合在一起，长脖的主根芦头聚合在一起，而这两种类型的主根芦头部分又有所重叠，说明在根部与叶片基因表达不同，而两个品种在基因表达中也存在差异。筛选两类型差异表达基因(图5)：叶片中有124个差异表达基因，其中上调基因48个，下调基因76个；芦头中有604个差异表达基因，其中上调基因283个，下调基因321个；主根中有89个差异表达基因，其中上调基因31个，下调基因58个，两种类型林下参差异表达基因主要存在于芦头中，这些信息对于分析两品种根部形态差异提供必要线索。

高中物理中有如图1所示的斜面模型，物体放在斜面上，斜面分为两种情况，一是光滑斜面，二是粗糙斜面。物体在光滑斜面沿斜面加速下滑，斜面无摩擦力；而对于粗糙斜面，下滑的物体做受力分析，物体受重力mg、支持力FN、摩擦力f，由于支持力FN=mgcosθ，摩擦力f=μFN=μmgcosθ，而重力平行斜面向下的分力为mgsinθ，所以当mgsinθ=μmgcosθ时，物体沿斜面匀速下滑，由此得sinθ=μcosθ，亦即μ=tanθ。其中μ为物体与斜面的摩擦系数，θ为斜面倾角。

注：CG—长脖根部；CLD—长脖叶片；CT—长脖根茎；EG—二马牙根部；ELD—二马牙叶片；ET— 二马牙根茎。Note: CG—Changbo root；CLD—Changbo leaf；CT—Changbo rhizome;EG—Ermaya root；ELD—Ermaya leaf；ET—Ermaya rhizome. 图4 两类型林下参叶片、芦头及主根主成分分析Fig.4 PCA of two types of forest cultivated Panax ginseng

图5 两种类型林下参的差异表达基因Fig.5 Genes significantly differentially expressed in two types of forest cultivated Panax ginseng

2.4 两种类型林下参差异表达基因的GO分析及KEGG富集分析

2.4.1叶片差异表达基因GO分析及KEGG富集分析 对叶片筛选的124个差异表达基因进行GO功能的分类统计(图6)，包括生物过程(biological process)、细胞组成(cellular component)及细胞功能(molecular function)三类：按基因参与的生物过程可分为22类，其中，新陈代谢过程(metabolic process)、细胞过程(cellular process)基因总数最多，单生物过程(single-organism process)、生物调节(biological regulation)次之，发展过程(developmental process)、信号传递(signaling)、多生物过程(multi-organism process)的基因最少；按基因在细胞内的组成可以分为15类，其中，细胞(cell)和细胞部分(cell part)这两类基因数量所占比例最大，膜(membrane)及细胞器(organelle)次之，膜封闭腔(membrane-enclosed lumen)最少；按基因在细胞中的功能可分为15类，其中催化剂活性(catalytic activity)及结合(binding)最多，转运体活性(transporter activity)次之，结构分子活性(structural molecule activity)、信号传导活性(signal transducer activity)、分子传导活性(molecular transducer activity)最少。说明两种类型林下参叶片在许多生理生化过程中表达量均不相同。

通过对两种类型林下参叶片中差异表达的基因进行KEGG富集(图6)，将124个差异基因富集到了16个代谢通路中，包括植物激素信号转导通路(plant hormone signal transduction)以及一些糖类、脂类、蛋白和次生产物的代谢，其中有两个代谢通路富集达到了显著水平(P<0.05)，氨基糖和核苷酸糖(amino sugar and nucleotide sugar metabolism)、蛋白酶体(proteasome)、激素信号及糖类等代谢的差异表达为后续分析两品种芦头膨大差异提供了重要信息。

图6 两种类型林下参叶片中差异表达基因的GO分类及KEGG富集的前16个通路Fig.6 GO categories and top 16 enriched KEGG pathways in leaf of two two types of forest cultivated Panax ginseng

2.4.2芦头差异表达基因GO分析及KEGG富集分析 对芦头筛选的604个差异表达基因进行GO功能的分类统计(图7)发现，芦头参与生物过程的基因分为20类: 新陈代谢过程、单生物过程最多，细胞过程次之，多细胞生物过程(multicellular organismal process)、多生物过程、增长(growth)最少；按基因在细胞内的组成可以分为15类，膜最多，细胞(cell)、细胞部分、膜部分(membrane part)次之，超分子络合物(supramolecular complex)、病毒(virion)、病毒部分(virion part)、胞外区(extracellular region part)最少；按基因在细胞中的功能可分为15类，催化剂活性最多，结合次之，信号传导活性、核酸结合转录因子活性(nucleic acid binding transcription factor activity)、电子载体活性(electron carrier activity)、分子传导活性、转录因子活性(transcription factor activity), 蛋白结合(protein binding)、养分贮液囊活性(nutrient reservoir activity)最少。说明两种类型林下参芦头在许多生理生化过程中表达量均不相同。

图7 两种类型林下参芦头中差异表达基因的GO分类图及KEGG富集的前20个通路Fig.7 GO categories and top 20 enriched KEGG pathways in rhizoma of two two types of forest cultivated Panax ginseng

在604个差异基因中，用160个差异基因共富集出70个代谢通路(图7)，其中有26个代谢通路富集达到显著水平(P<0.05)，17个代谢通路富集达到极显著水平(P<0.01)。首先富集到的是苯乙醇苷生物合成通路(phenylpropanoid biosynthesis)，有38个差异基因，占所有差异表达基因的23.75%，说明两种类型林下参芦头的苯乙醇苷含量有所不同；淀粉和蔗糖代谢(starch and sucrose metabolism)、戊糖和葡萄糖醛酸的相互转化(pentose and glucuronate interconversions)、糖酵解/糖异生(glycolysis/gluconeogenesis)、半乳糖代谢(galactose metabolism)中差异基因均显著，说明两种林下参的芦头在干物质的积累和能量代谢中有所不同；在组氨酸代谢(histidine metabolism)、脂肪酸降解(fatty acid degradation)、赖氨酸降解(lysine degradation)等通路中差异基因均达到显著水平，说明两种类型林下参的芦头在氨基酸的代谢方面有所不同；光合作用-天线蛋白质(photosynthesis-antenna proteins)通路中基因差异显著，说明两类型林下参光合速率也有所不同；类黄酮生物合成(flavonoid biosynthesis)、倍半萜类化合物和三萜类化合物的生物合成(sesquiterpenoid and triterpenoid biosynthesis)、异喹啉生物碱生物合成(isoquinoline alkaloid biosynthesis)等在两类型林下参中也可能存在差异，这些信息对于分析人参生长发育、次生代谢产物以及逆境胁迫等提供线索，其中苯乙醇苷生物合成的差异表达为芦头膨大差异研究机制提供了重要信息。

2.4.3主根差异表达基因GO分析及KEGG富集分析 对主根筛选的89个差异表达基因进行GO功能的分类统计(图8)，主根参与生物过程的基因中分为20类：新陈代谢过程、细胞过程、单生物过程最多，定位类(localization)次之，刺激响应(response to stimulus)、多细胞生物过程、增长最少；按基因在细胞内的组成可以分为15类，细胞、细胞部分、膜最多，细胞器、膜部分次之，膜封闭腔(membrane-enclosed lumen)最少；按基因在细胞中的功能可分为15类，催化剂活性最高，结合次之，核酸结合转录因子活性最少。说明两种类型林下参主根在许多生理生化过程中表达量均不相同。

图8 两种类型林下参主根中差异表达基因的GO分类及KEGG富集的前17个通路Fig.8 GO categories and top 17 enriched KEGG pathways in root of two two types of forest cultivated Panax ginseng

89个差异基因共富集到17个代谢通路，包括糖、脂肪、氨基酸以及一些次生代谢产物等通路(图8)，说明两种类型林下参主根生命活动也有许多差异，其中有1个代谢通路富集达到显著水平(P<0.05)，为柠檬酸循环(citrate cycle/TCA cycle)，其是三大营养素(糖类、脂类、氨基酸)的最终代谢通路，又是糖类、脂类、氨基酸代谢联系的枢纽，还有许多中间产物参与其他生命活动，说明两品种林下参主根存在差异。

2.5 两种类型林下参关键通路基因转录表达分析

在KEGG通路中，筛选4条关键通路(图9)。在同源重组代谢通路中，二马牙3个部位Rad54基因均表现为下调，该基因在DNA损伤的修复及染色质重塑交换反应中有重要作用[32]；在淀粉和蔗糖代谢代谢通路中，二马牙大部分基因(β-淀粉酶、果糖激酶、海藻糖6-磷酸合酶/磷酸酶、磷酸化酶、葡萄糖磷酸变位酶、葡萄糖-1-磷酸腺苷酸转移酶、1,4-α-葡聚糖分支酶、海藻糖6-磷酸合酶/磷酸酶)表现为下调，芦头中2个基因(β-果糖呋喃糖苷酶、α-海藻糖酶)表现为上调，1个基因(β-葡糖苷酶)表现为既有上调又有下调，最终可能表现为叶片与主根中的淀粉等物质积累，芦头中葡萄糖合成增加。

在植物激素信号转导通路中，二马牙叶片乙烯(ethylene)信号通路中ETR以及EBF1/2上调，可能会促进果实成熟及植株的衰老，芦头植物生长素(auxin)信号通路中响应因子AUX1和AUX/IAA下调，脱落酸(abscisic acid)代谢通路中蛋白磷酸酶2C(PP2C)下调，茉莉酸(jasmonic acid)代谢通路中MYC2下调，激素的差异响应可能会造成两类型林下参生长发育及逆境胁迫响应不同。

2.6 差异表达基因qRT-PCR验证分析

在芦头中筛选出7个差异表达基因：AUX/IAA(生长素响应蛋白，auxin-responsive protein)、AUX1(生长素输入载体，auxin influx carrier)、PP2C(蛋白磷酸酶，2C protein phosphatase 2C)、MYC2-1(MYC2转录因子，MYC2 transcription factor )、MYC2-2(双螺旋DNA家族结合蛋白转录因子，MYC2 basic helix-loop-helix DNA-binding family protein transcription factor )、CCR(肉桂烯醛基-CoA还原酶，cinnamoyl-CoA reductase)、LUP4(β-淀粉合酶，beta-amyrin synthase)，利用qRT-PCR技术对它们的表达量进行验证。差异基因的表达结果均与RNA-seq结果一致(图9)，两品种间无显著性差异(P>0.05)，说明本次转录组测序结果可信度较高。qRT-PCR的结果也进一步验证了两类型林下参根部形态学差异的分子机制，在植物激素代谢通路中，二马牙芦头的AUX1、AUX/IAA、PP2C以及MYC2表现为下调，激素信号基因的表达差异为两类型林下参表型差异提供了重要信息，木质素的合成关键酶CCR基因表现为上调，可能导致木质素合成增加。

图9 关键 DEG 序列的定量 PCR 检测Fig.9 Real-Time PCR verification for some important DEGs

3 讨论

根系形态直接影响植物对水分及养分的吸收及运输[33-34]，从而影响产量及有效物质的积累。比较两种类型林下参一些农艺性状指标及人参皂苷含量可以发现，同年生二马牙的芦头较长脖短，但粗壮，根部及地上部分也更粗壮，须根更发达，与前人研究结果一致[14]，而两类型人参皂苷含量无显著差异，这与前人研究两类型人参皂苷含量有差异结果不同，可能与测量人参皂苷样品量较少、林下环境复杂等因素有关。

两种类型林下参在糖类积累与代谢、激素作用以及木质素的合成中存在差异基因的表达，这可能是两类型林下参参形存在差异的主要原因。β-淀粉酶基因参与淀粉降解，与土豆块根发育有关[35]，并且为人参皂苷中合成齐墩果烷型人参皂苷的关键酶[36]，二马牙与长脖相比，与糖类合成与代谢相关的酶在植株的表达中体现在叶片与主根中储存淀粉，芦头合成葡萄糖，可能作用途径为：主根与叶片中的淀粉运输到芦头中转换为葡萄糖为芦头膨大提供能量，因此，二马牙芦头更加粗壮。

激素在许多植物根部生长发育过程中起重要作用，如马铃薯[37]、地黄[38]块根的形成。生长素在主根伸长及侧根发育中有重要调节作用[39-40]，二马牙的芦头在植物生长素信号通路中响应因子AUX1和AUX/IAA下调，而AUX/LAXs蛋白家族是目前已知的最主要的生长素输入载体[41]。Yang 等[42]研究发现，AUX1主要在根尖表皮细胞、中柱细胞、侧根根冠及原生韧皮部细胞中表达，参与生长素由根中向顶式和向基式的运输。何静等[43]研究发现，AUX1定位丰度的减少可以导致生长素在根尖顶部的积累，根部生长素浓度梯度的改变，会引起根生长的抑制及侧根的提前发生。二马牙的植物生长素表达少于长脖，可能造成芦头伸长生长受到抑制，而侧根生长提前，导致芦头短、须根多。

二马牙在木质素的生物合成通路中诸多基因表现为上调。人参根为肉质直根，肉质根膨大是由于薄壁细胞分裂、形成层细胞增生及次生组织发育共同作用的结果。随着次生细胞壁的薄壁细胞不断分裂，需要韧皮部和木质部不断扩大，导致木质部的木质化程度逐渐增加[44-46]，芦头的膨大可能与其木质化程度有关，二马牙的木质化程度强是其芦头膨大较大的结果也是其继续膨大的条件。

转录组结果还显示，两品种在抗逆方面可能存在差异，乙烯响应因子EFR1/2受多种防御因素诱导，可能与植物的木质化程度有关[47]，脱落酸代谢通路在植物抗逆响应中有很大的作用[48]；茉莉酸代谢通路中MYC2在调控植物体积累萜类物质等次生代谢产物和诱导植物防御反应中具有重要的作用[49]，二马牙林下参在这些基因的表达量与长脖均不相同，两品种的抗逆性差异还需要进一步的实验验证。

根据两种类型林下参的叶片、芦头及主根转录组的分析结果可以初步推测出两类型林下参芦头形态差异的分子机制：糖类的合成与代谢不同，为二马牙芦头与根部的生长提供能量，木质素的合成使得芦头木质化与膨大，生长素的作用使得芦头短粗，须根发达。通过对人参中形态差异较大的两种类型林下参差异表达基因的分析发现，其差异基因表达主要注释到淀粉和蔗糖代谢、植物激素信号转导通路以及苯乙醇苷生物合成通路中。为人参形态建成研究、分子育种、种质资源筛选以及优良品种选育提供参考。