烟草生物碱性状的QTL定位

刘颖超 方敦煌 徐海明 童治军,* 肖炳光,*

1浙江大学农业与生物技术学院作物科学研究所, 浙江杭州 310058; 2云南省烟草农业科学研究院国家烟草基因工程研究中心, 云南昆明 650021

烟草是世界上最重要的经济作物之一, 烟草行业为国家经济发展、就业稳定、脱贫攻坚和乡村振兴做出了巨大贡献。生物碱是烟草在长期进化过程中形成的内源防御性物质, 当烟草受到昆虫、病原菌侵害时, 烟草生物碱是抵御侵袭的武器。总植物碱(total plant alkali, TPA)是烟草产品的物质基础,对消费者有强烈的生理刺激[1]。烟碱(nicotine, NIC),音译为尼古丁, 占烟草生物碱的90%左右, 是最重要的一种生物碱, 烟碱含量关系到烟草生长发育及工业可用性[2]。在工业生产中, 烟碱含量过高或过低都会导致烟叶化学成分含量不协调, 降低卷烟工业可用性。除了烟碱, 烟草生物碱还包括降烟碱、假木贼碱和新烟草碱等其他种类生物碱。降烟碱(nornicotine, NOR)较不稳定, 容易发生氧化、酰化和亚硝化反应, 产生的物质对烟叶的香气和安全性有不利影响[3]。相较于不稳定的降烟碱, 假木贼碱(anabasine, ANAB)在烟叶中的含量较低且更稳定,对烟叶香气浓度有显著影响, 新烟草碱(anatabine,ANAT)对香气质、香气量、刺激性以及烟草特有亚硝胺含量等均有影响[4]。

了解烟草生物碱的遗传机制, 有助于改善烟草的品质。受限于烟草栽培品种间狭窄的遗传基础和极低的遗传多样性, 烟草的QTL定位研究落后于许多其他作物[5]。目前烟草QTL定位主要集中于烟草抗病性的研究, 如烟草青枯病[6]、赤星病[7]和白粉病[8]等病害抗性的QTL分析, 也有部分关于株高和叶数等农艺性状的研究[9], 关于烟草化学成分的探讨相对较少。早期有针对烟叶化学成分(糖类、氮、酚类化合物、氨基酸、生物碱和石油醚提取物)进行QTL定位的研究。2006年, Julio等[10]基于138个标记的遗传图谱, 在12个连锁群上鉴定出75个与农艺性状、化学成分相关的QTL, 对表型贡献率为8.0%～41.5%, 为烟草化学成分QTL定位研究奠定基础。最近, Tong等[11]基于45,081个SNP标记的遗传图谱, 鉴定出17个与烟草干叶化学成分相关的区域,LOD值为3.2～4.0。有研究针对烟叶的化学成分(总糖、还原糖、总氮、尼古丁和总钾含量)进行全基因组关联分析[12], 以确定与烟叶化学成分显著关联的基因座, 筛选含有烟草品质改良所需有利等位基因的优良材料。迄今, 关于烟草化学性状的研究多关注糖类、尼古丁、钾、氮等物质, 针对多种生物碱性状相关QTL定位分析相对较少。此外, 有研究表明在烤烟烟碱的遗传控制中除加性效应外, 上位性效应也具有重要价值[13], 上位性效应对杂种优势具有潜在作用[14], 然而大多数烟草QTL研究主要集中在QTL单位点效应, 上位性及其与环境互作几乎没有被关注。烟草生物碱是受生态环境、栽培技术以及基因型等多种因素影响的复杂性状[2], 在生物碱性状的QTL研究中, 分析QTL的环境特异性效应对培育特定环境的优异品种至关重要。

本研究利用烤烟品种K326和Y3为亲本构建RIL群体, 基于多年份与地点的表型数据, 对5种生物碱性状进行QTL定位。定位到的主效应QTL、上位性QTL及其与环境的互作效应, 为烟草生物碱性状遗传结构的解析提供了更多资源, 为烟草分子育种提供了新的理论依据。

1 材料与方法

1.1 研究材料与试验设计

本试验选用烤烟品种Y3和K326为亲本构建RIL群体。2012年, 在云南省烟草农业科学研究院试验基地种植亲本, 杂交得到F1籽粒。同年, 种植F1, 从F2世代开始采用单粒传法, 于2018年构建含有271份F7的重组自交系群体。RIL群体分别于2018年、2019年和2020年种植于云南省玉溪市研和试验基地(E1: 2018_YH)和云南省昆明市石林试验基地(E2: 2018_SL、E3: 2019_SL和E4: 2020_SL)。每个株系种植1行, 每行15株, 株距0.6 m, 行距1.2 m。每个株系随机取5株, 按烟草行业标准(YC/T 160-2002[15]和YC/T 383-2010[16])测定烟叶的总植物碱(TPA)、尼古丁(NIC)、降烟碱(NOR)、假木贼碱(ANAB)和新烟草碱(ANAT), 取5棵单株均值作为该株系的性状表型值。

1.2 数据统计与分析

应用Microsoft Excel和R软件, 对271份烟草RIL材料在4种环境组合(E1、E2、E3和E4)下的5种生物碱表型数据进行整理与分析。方差分量和遗传率估计基于以下混合线性模型:

式中, y是性状表型值向量; β是固定效应向量; ug和ue分别是基因型效应和环境效应的随机效应向量,ug～MVN(0, Iσg2), ug～MVN(0, Iσe2); X、Zg和Ze是对应的系数矩阵; ε是残差效应向量, ε～MVN(0,Iσε2)。利用R包sommer[17]的mmer函数估计随机效应的方差并获得基因型效应的最佳线性无偏预测(BLUP)值uˆg。遗传率通过以下公式估计:

式中,是估计的基因型效应方差,是估计的残差效应方差。利用R包ggcor (https://github.com/hannet91/ggcor)进行相关性分析及可视化。对于表型相关, 按环境组合分别计算成对性状间的皮尔逊相关系数; 而遗传相关, 则用预测的基因型效应值计算性状间的皮尔逊相关系数。

1.3 DNA提取与测序

利用Qiagen Dneasy Plant Mini kit (Qiagen,Hilden, 德国)从7周龄烟草幼苗的叶片中提取274个样本(271个F7代个体, 父母本和F1代个体)的基因组DNA, 利用超声波高性能样品处理系统(Covaris, Woburn, 美国)随机打断技术, 得到长度约350 bp的DNA片段。通过LM-PCR和滚环复制生成DNA纳米球, 质控合格后上机测序。测序文库根据BIGSEQ-500平台的要求构建, 测序策略为PE100[18]。

1.4 READS过滤与比对

使用SOAPnuke (v1.5.6)[19]对原始数据进行过滤,除“-l 40 -q 0.2”外, 其余参数均为默认。以烟草K326为参考基因组[20], 过滤掉长度小于500 bp的Scaffold,并将剩余的942,190 (～4.6 Gb)条Scaffold连成40条假染色体。利用BWA软件[21]“mem”比对模块将前一步过滤得到数据比对到参考基因组上(−t 8 -k 19 -M-R), 生成的比对结果sam文件通过Picard软件(http://broadinstitute.github.io/picard/)的SortSam.jar工具转变为sortbam文件。

1.5 SNP/INDEL检测与过滤

采用GATK (v 3.3.0)流程[22]对群体进行SNP和InDel的检测, 共鉴定了15,731,714个SNP位点和2,327,758个InDel位点, 过滤后剩余10,057,282个SNP位点和569,946个InDel位点。

SNP标记过滤标准为: (1) 两亲本和F1代个体的过滤条件为: DP (reads支持数)≥10且GQ (质量值)≥25, 否则标记该位点为缺失; (2) 271个F7代个体的过滤条件为: DP (reads支持数)≥3, 否则标记为缺失; (3) 过滤两亲本以及F1中有缺失的SNP位点; (4)过滤271个F7代个体中缺失率大于20%的位点; (5)过滤非二型的SNP位点。

InDel标记过滤标准为: (1) 两亲本和F1代个体的过滤条件为: DP (reads支持数)≥10且GQ (质量值)≥40, 否则标记该位点为缺失; (2) 271个F7代个体的过滤条件为: DP (reads支持数)≥3且GQ (质量值)≥10, 否则标记为缺失; (3) 过滤两亲本以及F1中有缺失的InDel位点; (4) 过滤271个F7代个体中缺失率大于20%的位点; (5) 过滤非二型的位点;(6) 过滤超过5 bp的InDel。

多态性分析共检测到两亲本(aa和bb)中表现出2种不同的纯合基因型, 在F1中表现出杂合基因型(ab)的SNP共1,626,811个, InDel共71,130个。

1.6 遗传连锁图谱的构建

由于构建遗传图谱的软件对标记数目限制, 获得的可用标记数目远大于软件分析上限, 因此对一个Scaffold上存在多个标记的情况进行过滤, 具体过滤方法为: 在同一个scaffold上, 100 kb的物理距离内仅保留1个标记(标记间的物理距离越近, 理论上发生重组的概率就越低, 标记冗余的可能性就越高)。过滤后剩余46,967个标记(SNP为44,804个,InDel为2163个), 用LepMAP3软件[23]的ParentCall2模块和SeparateChromes2模块, 最终保留了46,129个标记, 其中SNP和InDel数目分别为44,012个和2117个, 包含7292个bin标记(定义为在相同遗传距离下有多个标记时基因型缺失率最小的标记)。使用R包LinkageMapView[24]可视化遗传连锁密度图, 考虑到RIL群体的大小和标记的分布密度, 剔除了相邻距离小于1 cM的标记, 保留2187个bin标记用于后续的QTL定位分析。

1.7 QTL定位分析

利用新构建的遗传连锁图谱, 结合5种生物碱性状的表型数据, 采用基于混合线性模型的复合区间作图法(Mixed-linear-model-based composite interval mapping, MCIM)[25]及配套软件QTLNetwork 2.0[26](http://ibi.zju.edu.cn/software/qtlnetwork/)进行QTL定位。基于Henderson III的F统计量, 对基因组进行一维和二维扫描, 检测单位点效应显著的候选位点, 互作效应显著的成对候选位点, 后对候选位点进行逐步回归, 得到最终的QTL全模型; 基于全模型, 利用马尔可夫链蒙特卡洛(Markov Chain Monte Carlo, MCMC)方法估算QTL的各项效应并进行显著性检验。QTL扫描时采用置换检验方法检验其显著性, 试验水平的显著性概率值设为0.05。QTL命名方式参照McCouch (https://archive.gramene.org/newsletters/qtl_nome.html)。

1.8 候选基因预测

根据QTL定位的结果, 基于K326参考基因组[20],从K326参考基因组序列注释(https://solgenomics.net/ftp/genomes/Nicotiana_tabacum/edwards_et_al_20 17/annotation/Nitab-v4.5_gene_models_Scf_Edwards 2017.gff) (https://solgenomics.net/ftp/genomes/Nicotiana_tabacum/edwards_et_al_2017/annotation/Nitab-v4.5_proteins_Edwards2017.fasta)中选择QTL区域的候选基因及其对应的蛋白序列。利用NCBI的blastp模块, 从UniProtKB/Swiss-Prot (ftp://ftp.ebi.ac.uk/pub/databases/uniprot/knowledgebase/uniprot_sprot.fasta.gz)数据库中检索候选基因的功能, 查阅可能与生物碱性状相关的基因, 推断潜在的候选基因。

2 结果与分析

2.1 性状表型的基本统计分析及其遗传率

对烟草RIL群体在4种环境组合下生物碱数据的平均值、标准差、最小值、最大值和变异系数进行统计分析, 数据的偏度、峰度指标表明5种生物碱性状在RIL群体中基本呈正态分布(表1)。各生物碱性状在不同个体间存在差异, 变异系数范围为20.59%～33.14%。性状地点间差异较年份间差异大,表明这些性状对不同生长地区更敏感, 不同年份间相对稳定。5种生物碱性状遗传率范围为24.04%～33.09% (表1), 新烟草碱遗传率最高, 降烟碱遗传率最低。4种不同环境组合下, 不论是表型相关性(图1-a～d)还是遗传相关性(图1-e), 各生物碱的正负相关趋势基本一致, 且性状间显著的正相关都较强(r> 0.75)。使用预测的基因型值计算遗传相关系数,整体结果与4种环境组合下表型相关性的模式基本一致, 即烟草中总植物碱和烟碱含量越高, 降烟碱、新烟草碱和假木贼碱含量相应会越多。

图1 4种环境组合下生物碱性状的表型相关与遗传相关Fig. 1 Phenotypic correlation and genetic correlation of alkaloid traits in four types of environment

表1 4种环境组合下生物碱性状的统计Table 1 Statistics of alkaloids in four types of environments

2.2 遗传连锁图谱

本研究构建的图谱总计46,129个标记, SNP和InDel数目分别为44,012个和2117个, 包含7292个bin标记(定义为在相同遗传距离下有多个标记时基因型缺失率最小的标记), 分布在24个连锁群上, 全长3499.88 cM, 平均遗传距离0.480 cM; 最长连锁群LG6长度201.11 cM, 平均遗传距离0.457 cM; 最短连锁群LG24, 长度55.75 cM, 平均遗传距离0.526 cM (图2和表2)。

图2 连锁图谱的标记密度分布Fig. 2 Marker density distribution of linkage map

表2 24个连锁群中2类标记的概况Table 2 Summary of two types of markers in 24 linkage groups

2.3 生物碱性状的QTL定位

2.3.1 烟草生物碱相关QTL的加性效应及其与环境互作 对5种烟草生物碱性状进行全基因组QTL定位, 共检测到15个具有显著加性效应的QTL,除了降烟碱的部分QTL, 表型变异解释率相对较小,其他大多数QTL解释了2%以上的表型变异, 平均贡献率为5.12% (表3)。其中4个QTL对生物碱性状遗传变异有较大贡献, 分别为总植物碱的qTPA14(ha2=10.97%)、烟碱的qNIC14(ha2=10.04%)、假木贼碱的qANAB14(ha2=11.57%)和新烟草碱的qANAT14(ha2=10.77%), 说明这些QTL对调控烟草对应表型起主要作用, 可推测上述生物碱性状由一个主效QTL控制。降烟碱的4个QTL贡献率较小,表明降烟碱性状由微效多基因控制(表3)。

此外, 多个不同性状的QTL均位于连锁图谱上大致相同的遗传位置(表3),qTPA6、qNIC6、qNOR6、qANAB6和qANAT6均位于连锁群LG6的36.1～38.4 cM处;qTPA14、qNIC14、qNOR14、qANAB14和qANAT14均位于连锁群LG14的111.0～113.6 cM处。控制总植物碱的qTPA6、qTPA9-1和qTPA14, 控制假木贼碱的qANAB6、qANAB9-1和qANAB14, 以及控制新烟草碱的qANAT6、qANAT9和qANAT14基本处于对应连锁群的同一位置(LG6: 37.7～38.3 cM; LG9: 0～1.0 cM; LG14: 111.5～112.5 cM), 说明控制这几种性状的基因可能处于同一位置, 存在一因多效的可能性, 进一步解释了总植物碱、新烟草碱和假木贼碱间遗传相关性较强的现象(图1)。

虽然近一半的QTL具有显著的加性与环境互作效应(ae), 但它们对表型变异的贡献较小(ha2e<2%),贡献率最高的QTL (qANAT6)与环境互作效应也仅占新烟草碱表型方差的1.81%。部分QTL的ae效应表现出与主效应相同的方向, 如qTPA6、qNIC6、qNOR6、qANAB6和qANAT6在环境E1中与其主效应均正向促进对应生物碱含量的增加。然而与不同环境的互作也会产生与主效应方向相反的结果, 如控制烟碱的qNIC6在环境E2中与其主效应相反, 抑制烟碱的表达, 同样情况也包括控制降烟碱的qNOR6与控制假木贼碱的qANAB6。

2.3.2 烟草生物碱相关QTL的上位性效应及其与环境互作 二维基因组扫描检测到6对具有加性-加性上位性(aa)效应和/或环境互作上位性(aae)效应的QTL, 其中4对仅具有aa效应, 1对仅具有aae效应, 1对同时具有aa和aae效应(表4)。虽然这6对QTL存在的aa效应和/或aae效应显著, 但它们对表型变异的贡献较小, 大多低于2%。此外, 不同的环境中,aae的表现可能相反。例如, 总植物碱的qTPA1和qTPA9-2易受环境影响, 在环境E1中, 上位性与环境互作效应为-0.1093, 会抑制其表达; 在环境E4中, 效应为0.1711, 会促进其表达, 从而增加烟叶中总植物碱的占比。

2.4 QTL候选基因的预测

根据QTL定位分析的结果, 主效位点qTPA14、qNIC14、qANAB14和qANAT14分布在8号染色体6,624,235 ～7,447,394 bp处(表5)。qTPA14、qANAB14和qANAT14对应区间内共有6个编码基因(表5), 注释功能为ADP核糖基化因子GTP酶激活蛋白AGD8、泛素羧基末端水解酶和含有锚蛋白重复序列的蛋白At5g02620等。其中Nt08g00269编码与拟南芥ADP核糖基化因子GTP酶激活蛋白AGD8功能相似的蛋白, 与蛋白质运输、植物生长发育有关[27];Nt08g00273编码含有锚蛋白重复序列的蛋白, 参与植物激素调节过程[28]。qNIC14对应区间内共有20个编码基因(表5), 注释功能为蛋白磷酸酶2C 29、类叶绿体蛋白ABCI7、铝激活苹果酸转运蛋白和MND1-相互作用蛋白1等。其中,Nt08g00266编码转录因子MYC2, 为根瘤菌介导的诱导系统抗性(ISR)过程中增强茉莉酸(JA)响应基因表达的潜在调节剂[29]。当烟草遭受伤害时, 茉莉酸可激活烟草植株系统性应急反应, 促进烟草根系生物碱的合成,并加速向地上部运输与聚集[30], 从而有效对抗不同类型的病原体和昆虫等食草动物侵害。因此, 可推测Nt08g00266为控制烟草烟碱的关键候选基因。

?

3 讨论

本研究分别在云南省玉溪市研和、云南省昆明市石林2种地点3个年份下探索烟草RIL群体的生物碱性状遗传规律。构建包含46,129个标记的遗传连锁图谱, 基于此图谱对生物碱性状进行QTL定位分析, 共检测到4个与生物碱性状相关的主效QTL,分别是控制总植物碱的qTPA14、烟碱的qNIC14、假木贼碱的qANAB14和新烟草碱的qANAT14, 均解释10%以上的表型变异, 可推测这几种生物碱性状受主效应QTL控制, 在培育相应品种时应予以重视。烟草化学成分遗传机制的探讨相对较少, 尤其是降烟碱、新烟草碱和假木贼碱。本研究定位到的部分QTL (qNOR14、qTPA14、qANAB14、qANAT14和qNIC14)与前人定位到的位点位置相同[11], 在LG14的86.40～135.30 cM处, 部分QTL (qNIC6、qTPA6、qANAB6、qANAT6、qNOR6、qNOR7、qTPA9-1、qANAB9-1、qANAT9和qNOR11)在前人研究中未见报道, 是新发现的QTL位点, 对烟草化学成分遗传机制的剖析、功能基因的挖掘具有重要价值。

此外, 考虑到烟草生物碱受遗传及环境等多种因素影响[2], 本研究利用软件QTLNetwork 2.0, 将加性、上位性及基因与环境互作效应整合到一个QTL定位系统中, 这是大多数烟草QTL研究未考量的。基因与环境间的互作强弱会影响品种的稳定性,对其效应的准确估计是评价品种稳定性的基础[31]。关注上位性及其与环境互作效应, 探究基因间的互作关系及其对性状表达的影响, 有助于分子标记辅助改良育种[32]。本研究的5种性状中共检测到6个具有环境互作效应的QTL, 6对具有上位性和/或环境互作上位性效应的QTL, 为今后选育环境特异性或适应性强的品种提供良好理论基础。

烟草生物碱是烟草化学成分的重要组分。卷烟工业生产方面, 烟草生物碱是卷烟商品价值、口感风味的物质基础; 植物生长发育阶段, 烟草生物碱在自身保护中起重要作用。本研究定位到的烟碱主效QTL (qNIC14)对应区间内基因编码类似转录因子MYC2的蛋白, 通过参与茉莉酸的调控表达, 促进烟草根系生物碱的合成, 进而影响植物的抗病虫性[29]。后续可利用分子生物学技术对相关候选基因进行功能验证, 这对理解生物碱性状的遗传机制以及烟草品质的分子育种具有重要意义。

4 结论

本研究对烟草品种Y3、K326构建的RIL群体进行测序, 构建包含46,129个标记的遗传连锁图谱,检测5种生物碱表型, 利用基于混合线性模型QTL定位方法及软件QTLNetwork 2.0, 进行QTL定位分析。共检测到15个具有显著加性效应的QTL, 加性效应对表型贡献率为0.58%～11.57%。其中4个主效QTL, 分别为控制总植物碱的qTPA14、烟碱的qNIC14、假木贼碱的qANAB14和新烟草碱的qANAT14, 均可以解释相应性状10%以上的表型变异, 均位于14号连锁群上。6个具有显著加性与环境互作效应的QTL, 对表型贡献率为0.80%～1.81%。6对QTL具有显著加性-加性上位性效应和/或环境互作上位性效应, 其中4对仅具有aa效应, 1对仅具有aae效应, 1对同时具有aa和aae效应, 上位性效应对表型的贡献率为0.15%～2.31%, 上位性与环境互作效应对表型的贡献率为0.81%～1.16%。本研究挖掘控制生物碱的QTL, 为今后进一步研究相关的候选基因和烟草分子改良奠定了基础, 为培育优质新品种提供了新的见解。