转录组变异图提供了对遗传基础的见解Pinus Massoniana.羊肉。进化与油树脂屈服的关系

背景

马龙松(马尾松羊)是中国南方主要的本土针叶树种，是木材和油树脂的重要经济来源。然而，对马尾松种质遗传变异的认识仍然有限。本研究利用来自10个主要分布区的204份野生马尾松种质资源，利用转录组测序数据获得的94,194个核心单核苷酸多态性(SNPs)，对马尾松种质资源的遗传多样性和群体结构进行了分析。

结果

平均期望杂合度为0.2724，表明马尾松种质具有丰富的遗传多样性。分子方差分析(AMOVA)表明，3.29%的变异来源于遗传分化。结构分析确定了两个地理上截然不同的群体。主成分判别分析(DAPC)结果表明，其中一个类群进一步分为两个类群。川渝种源是地理起源，沿两条不同的方向向外扩散。松脂的产率反映在两组的演化过程中，并沿两条扩散线呈现出两种不同的趋势。基于SNPs和表达的分析表明，几丁质酶、CYP720B、细胞色素P450、ABC转运体和AP2/乙烯应答转录因子(ERF)基因可能与油脂产量有关。

结论

来自转录组测序的SNP标记能够在不同物种中评估遗传多样性，以及对象特征的遗传控制。将在未来的研究中核实这些基因的功能，并且那些与油树脂产量强烈相关的那些基因将通过早期基因型选择和基因工程来改善产量。

马尾松是本地的优势树种(Pinus Massoniana.是一种具有重要商业价值的针叶树，为我国提供木材和油树脂。自然分布于北纬21°41′~ 33°56′，东经102°10′~ 123°14′，种植面积200万公顷[1]．广东、广西、湖南、四川、重庆、贵州、浙江、福建和江西等省是我国马尾松的主要自然分布区[2]．因为树具有快速生长和贫瘠土壤的特点，因此通常被认为是污稀覆盖山区内造林的先驱物种。遗传多样性对于物种的长期存活至关重要，推动物种以适应各种非生物和生物应力，以避免灭绝[3.]．对大多数树种而言，根据种源或家系分析，在自然居群之间或居群内部可以观察到生长、萜类产量、抗病能力等方面的巨大遗传变异[4.那5.那6.]．为了揭示马尾松主要经济性状之间的遗传变异，自1978年以来在中国开展了大规模的种源试验。中国已经建立了两条完整的原生源区和许多局部原生源区[2]的研究为揭示导致表型多样性的各种进化力量之间的相互作用和重要性提供了良好的材料，并制定了捕捉片段内自然遗传多样性的基因保护策略。马尾松胸径(胸径)随纬度的地理变化呈经典模式[7.]．

作为由Masson Pine获得的二级物质，大黄素是一种重要的天然产物，用作化学工业中许多不同化合物的来源[8.那9.]，以防止昆虫和疾病[10那11，用作高级液体生物燃料[12]．马尾松不同科间与油松脂产量相关的遗传变异显著，为14.12 ~ 50.55 g / d [13那14]．Zeng等人和Westbrook等人也报道了火炬松松脂产量的变异是可遗传的(P. Taeda.)，通过选择，可在一代内增加1.5 ~ 2.4倍。

分子标记对于鉴定种质，评估生物多样性以及描述遗传变异的地理模式非常有用。单核苷酸多态性（SNP）通常用于遗传研究。利用下一代测序（NGS）技术，以较低的成本迅速发展数百万SNP [9.]．这些高吞吐量SNP已成功用于评估遗传多样性并推导人口结构[15那16和亲属关系[17]．作为中国南方的重要森林树，高密度SNP地图对于遗传创新至关重要，并在未来的繁殖计划中提高马龙松的特征。然而，迄今为止，没有关于开发SNP标记的完整基因组序列的报道，以研究Masson Pine的遗传多样性和结构。使用随机扩增多晶型DNA（RAPD）进行分析仅部分Masson Pine Germplasms [18，间简单序列重复(ISSR) [19]，简单的序列重复（SSR）[20.[和反横向跨越多态性（IRAP）[21]．

缔合转录组织对鉴定与表型变异有关的序列多态性和转录性丰富，特别是在非模型物种中，显着贡献22那23]．此外，高质量的全长转录本对于功能分析和理解遗传多样性至关重要[24]．在这项工作中，我们首先通过全长转录组和基于NGS的unigenes组合构建高质量的转录参考序列。然后采用204个代表性载体的RNA测序（RNA-SEQ）进行DE Novo SNP发现，以产生基因组范围的变异图。我们研究的目的是：（1）评估Masson Pine的遗传多样性，人口结构和地理来源;（2）揭示与油树脂产量相关的基因。我们的研究结果将有助于管理该物种，并阐述形成高产大肠素的机制。

测序和变异发现

马尾松的基因组和近缘物种的基因组一样大而复杂，缺乏高质量的基因组序列，阻碍了马尾松的研究和育种计划P. Taeda.[25]．为了克服这种障碍，我们通过PACBIO单分子，实时（SMRT）测序平台，从次级Xylem转录组中构建了一种高质量的全长转录物组。从18 GB的PACBIO亚克斯获得了总共81,837个高质量和非冗余的全长转录物。为了探讨遗传多样性的起源和模式，我们还为204个地理上多样的Masson Pine型，我们设计了对Illumina Hiseq™2000序列测序平台的人口转录组实验，从中国的主要栖息地收集了他们的主要栖息地;总共组装了341,714个非冗余的unigenes。随着全长转录物和未成年人组合，将423,288个非冗余转录物视为进一步分析的参考序列（附加文件1：表S1）。平均而言，每个样本的85.02％的读数成功映射到参考序列，表明参考转录物的完整性高达（附加文件2：表S2）。

使用GATK包裹从204种基因型的转录组中检测到总共1,326,230个SNP和153,459个插入/缺失（诱导）[26]，平均SNP密度为每份转录版3.13（附加文件3.:表S3)。其中次要等位基因频率(MAF)≥0.05和缺失基因型呼叫< 5%的核心SNPs 94,194个，占总数的7.1%。这些核心SNPs包括23,864(25.33%)非同义(nsSNPs)(附加文件4.：表S4）。该转录组变异图将有益于核心种质鉴定，遗传变异研究和人工育种。

马龙松的遗传多样性

遗传多样性P. Massoniana基于94,194个SNP，研究了来自各种分布的主要区域的种质。观察到的杂合性（H_0.）值低于预期的杂合性（H_E.）每个人口的价值，从0.2211（广西）到0.2358（四川和重庆）。的H_E.不同群体中的值相似，范围为0.3011（江西）至0.3124（四川和重庆）（表1）.近交系数(F索引从0.2242（四川和重庆）到0.2714（广西），平均值为0.2731，表明四川和重庆人口中的SNP具有最高的多态性。通过基于94,194个SNP的Amova测试九种种群之间的推定差异（表2）.结果表明，群体中的分化解释了总方差的3.29％。在不同的亚步骤中发现了0.01％的变异，暗示了他们内部的封闭性血缘关系。总之，我们的变体数据集提供了各种尺度的基因组多样性的全面概述，并代表了学术和农业研究社区的剥削遗传信息的丰富来源。

建设的P. MassonianaSNP高多样性核心种质

等位基因多样性P. MassonianaSNP标记可以最大限度地增加选材。冗余曲线显示，马尾松的等位基因多样性可被40多个核心种质所代表。这40个具有代表性的基因型仅占样本总数的20%，占等位基因多样性的90.7%以上。1和表3.）.因此，这40份具有代表性的材料中，浙江9份、贵州9份、四川和重庆7份都可以构建出最小的核心种质5.：表S5）。据我们所知，这是基于大量人群的高密度SNP地图的核心种质的首先综合鉴定，这是有价值的P. Massoniana育种实践。

人口结构P. Massoniana种质

为了进一步了解Masson Pine的进化历史，我们使用了混合物软件[27]估计每次加入的祖先比例。遗传分配分析显示了最佳价值K. = 2, which clearly separated the accessions of Chongqing and Sichuan from those of the other wild genotypes (Additional file6.：图。S1和附加文件7.:图S2A)。第一组主要是来自重庆和四川的无性系，种群间混合的信号水平较高(图1)。2一种）。为了k = 3.，来自中国中部和中国东南部的两个新群体，从重庆和四川的加入中出现（图。2b）。值得注意的是，包括重庆和四川省克隆的群体也显示出高水平的混合物。II集团II，其中包括来自中国中部的主要克隆，包括贵州，广西，广东和湖南省份，表现出了群体中的速降信号。这可能是由于在分离后由于动物或风的分散而发生的自然杂交。．中国东南部，包括福建，江西，浙江和安徽省的克隆被分配到第三组。有趣的是，III群体保持均匀的遗传背景，可能是由于它们的地理孤立，阻止了杂交种类的杂交（附加文件7.:图开通)。结果表明，马尾松基因型的居群结构与地理分布基本一致。

主成分（DAPC）的判别分析显示了三种遗传簇，驱动了我们面板内的多样性分配（图。2c)。群集我只包括四川和重庆省（94.4％）（额外文件）8.:表S6);聚类II主要包括江西(100%)、福建(100%)、浙江(100%)和安徽(100%)的克隆，占聚类II的88.3%;III类主要种质来自贵州(100%)、广西(96.6%)、广东(61.1%)和湖南(88.2%)。这三个基因群在地理上是孤立的。集群I由主要生活在中国西部的品种组成;聚类II主要分布在中国东南部;第三类则包括来自中国中南地区的加入。DAPC分析结果与种群结构分析结果一致k =3.

我们进一步估计了不同群体的遗传多样性。聚类I的全基因组核苷酸多样性(π) (2.91 × 10⁻²）高于群体II（2.77×10⁻²）和群集III（2.83×10⁻²）并展出了最高的多样性。该结果也得到了支持H_E.值，显示了基于杂合位点的序列多样性(表1)4.）.的Nei's遗传距离显示0.135（簇II与簇III）至0.303（群集I与群集II）的值，而成对固定指数F_圣范围从0.024（集合II与群集III）到0.110（群集I与集群II），以及Nei's和F_圣群集I与群集II的遗传距离高于群集I与群集III（附加文件）9.：表S7）。这些观察结果表明，SICHANG盆地的Masson Pine保持了高遗传多样性，并且比中国东南部更大的分化。

地理上的起源与扩散P. Massoniana种质

为了进一步阐明马尾松的进化图谱和传播途径，我们构建了204个马尾松基因型的系统发育树(图1)。2d)。此外,P. Taeda.被分配为最大可能性树的超级，以确定最早的分歧人口，被视为现代的祖先P. Massoniana（附加文件10：图。S3）。系统发育树表明，四川盆地（四川，重庆）的基因型最接近P. Taeda.然后是其他林，暗示四川盆地是马龙松的地理起源。四川盆地是在最后一代冰川冰川期间多种物种的几种冰川难民之一[28那29]，可能已经从消灭事件中救出了马龙松。马龙松在冰川时期结束后逐渐迁移到贵州高原，并逐渐适应高原栖息地。湖南基因型从贵州思潮分支组成的亚洲，随后是中国南部（广东和广西）和中国东南部（江西，福建，浙江和安徽）的其他马龙松基因型。值得注意的是，这些基因型清楚地根据其地理分布分为两个亚基。这种观察允许我们提出了Masson Pine Evolution Map中的两种不同取向扩展线的假设（图。3.一种）。一个迁徙线从四川和重庆到贵州，然后到湖南，然后散布在广东和广西。另一条线从四川和重庆到贵州，然后到湖南，然后传播到江西，福建，安徽和浙江。人口结构证据强烈支持这一假设（图。3.b).四川/重庆与贵州之间的种群分化显著更大(F_圣 = 0.13) than between those of other populations, implying a strong variation in the genome for the new natural adaptation when first transferring from the basin to the plateau habitat. The nucleotide diversity is slightly higher for the progenitors of the Sichuang/Chongqing population (Fig.3.一种）。TREEMIX计划检测到两种传播线的群体之间钝化信号。检测到广东/广西人群的四川盆地人口的杂交信号（图。3.c)。

联想转录组和油树脂产量

木瓜中的油树脂产量P. Massoniana在204个马尾松无性系中差异较大，水平在0.00 ~ 6.07 g·cm之间⁻¹·D.⁻¹（附加文件11:表S8)。各供试材料中树脂产量分布呈正偏态(图1)。4.一种）。四川盆地和贵州省的屈服的产量明显低于湖南进程（图。4.(2)在中南扩散路径上，马尾松向广东和广西扩散时，松脂产量略有下降，但仍高于四川盆地和贵州的产量。在中国东南部扩散路径上，马尾松向中国东南部扩散时，油松脂产量显著增加，尤其是向安徽、浙江和江西等地扩散。

结合转录组学分析从109个转录本中鉴定出121个SNPs，这些SNPs与油树脂产量显著相关P. < 10⁻⁶意义水平（图。4.c,额外的文件12：图。S4和附加文件13:表S9)。最显著的SNP (c51955_f1p3_1546，R.² = 0.51,P.= 3.74E−19)定位于几丁质酶I类注释的转录本(表5.）.突变SNP位于编码区上游，但c51955_f1p3_1546的表达量与油脂产量无显著相关性。

属于细胞色素P450单氧基酶的CYP720B系列（P450），是参与二萜树脂酸生物合成的重要类别，作为大黄树脂的主要含量[30.]．一个CYP720B (c19795_f1p0_1763)和一个细胞色素P450 (c9591_f1p0_1663)被发现与油脂产量(9.60E−07,1.85E−07)序列相关。CYP720B突变SNP导致非同义突变，密码子CTC向TTC过渡。细胞色素P450突变的SNP (c9591_f1p0_1663)属于同义突变。这两个转录本(c19795_f1p0_1763, c9591_f1p0_1663)的表达与树脂产量显著相关(P.= 3。6.1E− 08,P.= 2.13 e−08年)。对10份高产和低产松脂材料进行实时定量聚合酶链反应(qRT-PCR)，结果显示这两种转录本在高产松脂马尾松中均有较高的表达水平。4.d)。

AP2结构域转录因子和ABC转运体的序列与油树脂产量有关P. Taeda.[14]．在该研究中，来自AP2 /乙烯响应转录因子（ERF）（C24091_F1P1_1286，C8825_F1P0_1733）和来自ABC转运蛋白的一个SNP（C189021.Graph_C0）的两个SNP也与Masson Pine中的产油素产量显着相关。来自AP2 / ERF（C8825_F1P0_1733）的SNP导致非同义编码和从胱氨酸转换为精氨酸的编码氨基酸。

此外，来自微管蛋白α链（C20772_F1P4_1467）的转录物的一个SNP与序列中的油树脂产量显着相关（P. = 8.73E− 16) and expression level (P. = 4.83E− 08) simultaneously. The SNP resulted in non-synonymous mutation, with the transition of codon CTC to TTC. However, the function of the tubulin alpha chain during the biosynthesis of oleoresin is unclear.

SNP标记已被用于评估许多物种的多样性，包括杨树trichocarpa[31]，vitis Vinifera[32),而银杏毕洛巴巴[33]．源自转录组测序的SNP是一种更有效的策略，用于表征非模型或大规模基因组物种中的多样性，因为在编码区域而不是整个基因组上检测到序列。在本研究中，通过转录组测序获得的94,194个SNP用于研究来自中国十个省和市的马龙松的多样性。的H_0.值（约0.22）略低于H_E.价值观（群体中的0.30），暗示人群中发生了频繁的近亲繁殖事件（表1）.任何一个H_E.或H_0.可以用来评估基因变异，但是H_0.价值往往受到人口内近亲繁殖水平的影响。所以，H_E.更常用于比较不同物种之间或同一物种内种群的遗传多样性[3.]．

Masson Pine在中国各地区拥有持续的本土分布。Hamrick等人。[34发现，使用丙比中分析，具有广泛分布的树种群中的平均预期杂合性为0.228。黄和张[35报道了H_E.贵州省马尾松六种自然群体的价值为0.27，通过同工酶分析确定。评估了福建省马尾松五个群体的遗传多样性，平均值H_E.被发现为0.22 [36]．在本研究中，马尾松的遗传变异性较高，可能是由于取样地区面积较大，几乎覆盖了马尾松的整个原产于地区。同样的结果也在苏格兰松(P. Sylvestris.)[37那38]．此外，在这些研究中评估的遗传多样性的差异也可能是不同标记类型，采样位置和大小的结果[39]．

结构分析和DAPC都将四川和重庆的样品分开（图。2(a, c).这个区别也符合F_圣和Nei's遗传距离(附加文件9.:表S7)，显示四川和重庆两省的种质值最高F_圣和Nei's遗传距离，分别由DAPC。虽然结构分析显示了最小的交叉验证误差K. = 2, the cross-validation error atK.= 3只比at的略高K.= 2。为了K. = 2, most of the germplasm from the other provinces not including Sichuan and Chongqing were grouped into one cluster (Fig.2a).然而，这个集群被分为两组K.= 3，与DAPC的聚类具有较强的对应关系，尽管每个聚类的成员之间存在较小的差异。集群II和集群III之间的差异相对较小F_圣和Nei's遗传距离值分别为0.024和0.135，暗示来自中国中南部和东南部的马龙松种质比来自中国西南部更密切相关。此外，与另外两个集群相比，我由四川和重庆种质组成的集群，位于中国西南部的遗传多样性最高，有一个H_E.价值0.318。

气候是森林树木适应性进化的主要驱动因素之一[40那41.]．在北半球，始新世的暖亚热带和暖温带气候，有丰富的裸子植物，从渐新世到新生代，在中纬度和高纬度陆块，特别是第四纪有大规模的冰覆盖和冰川，转变为寒冷的强季节性气候[42.那43.那44.]．在此期间严重感冒，许多树种被灭绝。然而，一些树种的装备更好地适应持续的较冷条件。在中国南方，复杂的地形有助于众多温带森林在各种“难民”中的最后冰川最大值存活下来[45.),如银杏毕洛巴巴那metasequoia glyptostroboides.那Glyptostrobus pensilis， 和Liroyendron Chinense.，这些都在中国仍然存在。四川盆地包括四川省的中东部和重庆市，被青藏高原、大巴山、巫山和云贵高原(海拔1000-3000米)所环绕。然而，四川盆地底板的标高仅为250 ~ 750米。虽然在第四纪四川盆地也发生了冰川作用，但较冷的气候并没有导致植物和动物物种的灭绝，如幸存下来的大熊猫(Ailuropoda melanoleuca）.因此，四川盆地是众多物种的主要避难所之一[46.]．这些物种在冰川期间扩展到较低的海拔，并在中间尖阶段期间撤退到更高海拔的避难所[45.那47.那48.]．

为探讨马尾松的进化历史，我们采用火炬松(P. Taeda.）作为参考。我们的研究表明，四川盆地是马龙松的地理起源。这符合Qin报告的先前结果[46.[他对两种物种针的特征的研究。然而，结构分析显示了两个地理上不同的群体，并且DAPC鉴定了本研究中的三种簇，这表明基因一直在改变以适应栖息地。

马龙松最初从四川盆地蔓延到贵州省。虽然邻近的贵州省位于广西和湖南省的北部和西部边界，分别是马龙松，随后遍布贵州湖南。云南 - 贵州高原可能是妨碍来自贵州到广西的马龙松传播的障碍。虽然广西经络与使用结构，DAPC和钢皮图分析的几个其他杂种不高度区别，但广西和广东杂志与增长特征的其他杂种之间的差异是显着的。来自广西和广东杂烩的马龙松的传播更快，与那些起源的热资源有关[7.]．

通过选择与目标性状相关的基因，可加速马尾松的选育。Np据报道，ABC1是第一个参与大豆萜类分泌的转运体[49.]．在针叶树中，油树脂从活细胞运输到树脂管道中，当茎患有非生物刺激时从伤口流动[50.]．Westbrook等。[14]发现位于ABC转运蛋白中的SNP与油树脂产率相关，并推断出ABC转运蛋白参与作业素运输。在该研究中，SNP的结果也表明ABC转运蛋白与大素产率显着相关，这表明ABC转运蛋白通过改变序列来调节油树脂产量的重要作用。

几丁质酶在修改细胞壁的结构方面发挥着关键作用。钟等。[51.发现几丁质酶突变基因(elp1)会导致木质素异位沉积在茎中拟南芥木质素化细胞壁未增厚。几丁质酶的作用可能影响油树脂从活细胞到树脂管道的运输速率。

据报道，大多数细胞色素P450s涉及次生新陈代谢的进展[52.]．Cyp720b的细胞色素P450s的基因家族对于针叶树特异，并且可以催化各种二萜树脂酸的生物合成途径中的连续氧化步骤，作为油树脂的主要成分[30.]．我们发现CYP720B中的一个SNP和细胞色素P450中的一个SNP与油树脂产量显着相关，CYP720B中的SNP导致非同义突变。因此，推断CYP720B中的SNP对通过改变序列和表达水平来确定大素产量的重要影响。

乙烯可以诱导许多针叶树的创伤性树脂导管的生物合成和形成[53.]．AP2/ERF转录因子参与非生物胁迫乙烯信号通路中乙烯应答基因的表达调控。表达的欧塞尔巴普1导致水稻脂质代谢相关基因表达增加[54.]．在P. Taeda.，AP2结构域转录因子中的一个SNP也与油树脂产量相关[14]，通过我们的Masson Pine的结果验证。

了解Masson Pine的遗传建筑以改善遗传育种过程中的大素产量非常重要。虽然Masson Pine的基因组尚未被测序，但我们基于使用全长转录组的94,194个SNP获得遗传多样性，人口结构和特性 - 基因关联的令人满意的结果。Masson Pine显然分为两组，四川和重庆出处被证明是其地理来源，从中马尾松沿着两条不同的线路向外扩散。油树脂产率沿两条扩散呈两种不同的趋势，与几丁质酶，CYP720B，细胞色素P450，ABC转运蛋白和AP2 / ERF的基因相关，其中一些也被证实存在于其他针叶树中。这些基因的功能将在未来的研究中进行核实。

样品采集

浙江西部崂山林场(119°02′e, 29°33′n)马尾松克隆试验海拔为152 m)。本试验包括来自10个省市的400个克隆(基因型)。20世纪80年代，中国成立了马尾松国家技术合作小组，由当地政府授权的各省级技术合作小组鉴定并提供马尾松无性系的接穗。1985年4月，在马尾松上冠处采集了野生乔木的健壮枝条作为接穗。随后，采用髓形成层配对嫁接的方法将接穗嫁接到杭州中国林业科学研究院亚热带林业研究所2年生马尾松本地幼苗上。第二年，以这些嫁接苗为试验材料，采用完全随机设计，10个重复，单株间距为2.0 m × 3.0 m。对这些植物的实验研究，包括收集植物材料，完全遵守机构、国家和国际准则。实地调查是根据当地法例进行的。作者遵守了自1975年起生效的《濒危野生动植物种国际贸易公约》(CITES，又称《华盛顿公约》)。

本研究从中国10个省市随机选取204个马尾松健康无性系14:表S10)。在1997年重庆成为直辖市之前，它是四川省的一个城市。以位于四川盆地的四川和重庆为样本，将重庆和重庆作为单一群体进行分析。

根据Liu的方法测量油树脂产率[55.] 2017年5月至2018年5月之间，每棵树的大黄树脂产量计算为每CM条纹长度的单个树的产量。同时，在去除树皮和韧皮柱后，从样品树上收获5mm的深新鲜次级木耳组织。将这些样品立即置于野外液氮中，然后储存在-80℃以进行RNA提取。这些实验是在杭州中国林科院校亚热带林业研究所进行的。

RNA提取和pacbio基测序

每个样本的总RNA分别提取，并按照Liu的方法进行评估[56.]．简而言之，使用植物RNA Kit RN38 Easyspin Plus（Aidlab Biotech，中国）提取来自每个样品的总RNA。使用Ultraspec TM 2100 Pro UV /可见分光光度计和Agilent 2100生物分析仪检测总RNA的浓度和完整性。高质量的RNA样品用于构建cDNA文库。将来自每个样品的一种微克数量的RNA合并在一起，并且使用Smarter TM PCR cDNA合成试剂盒合成全长cDNA。使用Bluepippin（SiveScience，Beverly，MA，USA）选择全长CDNA的大小，并建立了三种不同大小的cDNA文库（1-2 kB，2-3 kB和> 3 kB）。然后使用qubitfluorometer（IN）量化cDNA的大小分布vitrogen使用Agilent 2100 BioAnalyzer评估了Thermo Fisher Scientific，Waltham，Ma，USA）和三个图书馆的质量。随后，使用太平洋生物学RS II进行SMRT测序（门罗Park, CA, USA)平台，生物标记技术公司，北京，中国。

下一代测序

使用磁珠富集程序从每个样本的高质量总RNA中获得mRNA。用裂解缓冲液随机对mRNA进行裂解。合成了第一和第二链cDNA。所有cdna均用AMPure XP珠纯化。经末端修复、加A、接头连接后，用AMPure XP珠选择纯化cDNA的片段大小。通过PCR扩增扩增cDNA片段，利用qubit荧光仪和Agilent 2100生物分析仪对每个样本的cDNA文库质量进行评估。最后，在Illumina HiSeq™2000测序平台上对每个样本的cDNA文库进行配对测序。

RNA-Seq数据的质量控制

对低质量的读数据按以下4条规则进行过滤:(1)如果对端读数据的一端有> 5%的“N”个碱基，则删除对端读数据;(2)对于每对读取，如果其中一方的碱基平均质量小于20，则将其全部剔除;(3)对于每一篇阅读，如果质量分数小于13，我们会削减3个基数。在基部停止修剪，质量评分≥13。修剪后，如果剩余碱基数小于40，则去除对端读基;(4)去除对端重复读取。干净的数据随后被用来调用SNPs和indel。

SNP和Indel呼叫

然后使用挖掘机轮椅对准器的BWA-MEM算法映射到参考序列（全长转录组）的过滤读取。在Masson Pine的204个样本上使用Gatk中的单倍型来电称为SNP。最后，删除了低质量的SNPS（质量<30，MQ <40.0，FS> 60.0和QD <2.0）。使用与SNP呼叫相同的管道调用indels。过滤原始衬衫以减少使用GTAK变体过滤使用参数：FS> 200，QD <2，读取POS等级< - 20.0的误报。

遗传多样性分析

观察杂合度的遗传参数(H_0.)、期望杂合度(H_E.)、次要等位基因频率(MAF)、近交系数(F）使用PLINK软件（1.9版）估计;http://zzz.bwh.harvard.edu/plink.）.通过使用Arlequin软件的分析（AMOVA）分析（3.5.2版）计算群体内群体和克隆内克隆的种群的变化（AMOVA）计算（3.5.2版;http://cmpg.unibe.ch/software/arlequin35）.

系统发育分析

在系统发育树中，火炬松(P. Taeda.）首先从NCBI数据库（SRX4454630）下载，然后与全长转录组序列对齐。随后，从遗失松树的基因组调用SNP。然后使用ITOL进行系统发育树可视化和编辑分配（http://itol.embl.de.）.马尾松与火炬松之间的发散时间采用在线时间树软件(http://timetree.org.）.最后，通过最大可能性（PML）包（PML）封装，使用MCMCTREE程序计算和可视化每个种质的发散时间（http://abacus.gene.ucl.ac.uk/software/paml.html)[57.]．

人口结构分析

掺合料软件(https://peciationgenomics.github.io/admixture.)[27[用于可视化人口的遗传结构。预设的祖先人口数范围从K. = 1 toK. = 10. The most likely number of ancestral genetic groups was determined by the minimumK.交叉验证曲线上的值。

主成分的判别分析（DAPC）

对于DAPC，首先利用主成分分析(PCA)将遗传数据转化为不相关成分。然后使用k-means聚类定义遗传簇的数量，这是一种寻找值的算法K.这使得不同群体之间的差异最大化。计算了贝叶斯信息准则(BIC)K. = 1–40, and theK.值最小时，选择最优聚类数。然后利用dapp功能对前120个主成分进行判别分析，有效地描述遗传集群。

油脂产量相关基因的鉴定

采用混合线性模型(MLM)、关联、进化和连锁性状分析(TASSEL) (https://www.maizegenetics.net/tassel.)[58.]．的P.计算每个关联对应的-value，当P.-value≤1.06E−5，用1/N命名为bonferroni校正（N是SNP标记的数量）。

定量分析中存在

采用10份高和10份低产油树脂的RNA样本进行qRT-PCR。引物对(附加文件15表S11)对几丁质酶、微管蛋白α链、AP2/ERF、ABC转运体、CYP720、细胞色素P450七个基因进行设计，按Liu的方法扩增cDNA [56.]使用2计算基因的表达水平^−ΔΔCt方法 [59.]．使用伸长系数1-α（EF 1-α）来规范转录物谱。

本研究期间生成或分析的所有数据都包含在本发布的文章中（及其附加信息文件）。本研究中生成的所有测序数据可从SRA存档中获得（https://www.ncbi.nlm.nih.gov/sra.)，生物项目编号:PRJNA636925。在PacBio SMRT测序平台上，一个全长转录组的原始序列在NCBI GeneBank中通过SRA运行获得SRR11912706。在Illumina HiSeqTM 2000测序平台上，204个样本的RNA-seq原始序列在SRR11912579- SRR11912705和SRR11912707- SRR11912783下存入NCBI GeneBank。

Amova：

分子方差分析

AP2：

Apetala2结构域转录因子

DAPC:

主成分的判别分析

ERFS：

乙烯响应转录因子

门店:

下一代测序

QRT-PCR：

定量逆转录 - 聚合酶链反应

SNP：

单核苷酸多态性

我们非常感谢殷恒福博士在手稿写作方面的建议。

本研究得到了中国国家重点研发计划在第13届五年计划期间（授予2017YFD0600300）和项目第十三次五个计划期间，浙江省浙江省新森林树木的主要特殊项目支持的项目（2016C02056）-4）。资金代理商在研究设计，数据收集和分析或数据解释中没有作用。

隶属关系

1. 2 .中国林业科学研究院亚热带林业研究所，浙江杭州311400

   刘庆华，谢怡妮，刘斌，周志春

2. 2 .浙江省林木育种重点实验室，浙江杭州311400

   刘庆华，谢怡妮，刘斌，周志春

3. 四川农业大学，四川成都，611130

   华桓扬

4. 春南春山林农场，春南，311700，浙江，中华人民共和国

   中平冯

5. 生物标记技术公司，北京，101300

   雅通陈

贡献

LQH参与了该作品的设计;数据的获取、分析和解释以及手稿的撰写。XYN有助于数据的获取和分析。LB有助于数据的获取和分析。YHH为数据的获取和分析做出了贡献。FZP有助于植物材料的收集。CYD对数据的获取和分析作出了贡献。ZZC为植物材料的收集、构思和设计工作做出了贡献。作者阅读并批准了最终稿件。

通讯作者

对应到Zhichun周．

伦理批准和同意参与

不适用。

同意出版物

不适用。

利益争夺

提交人声明他们没有竞争利益。

出版商的注意事项

Springer Nature在发表地图和机构附属机构中的司法管辖权索赔方面仍然是中立的。

开放获取本文是基于知识共享署名4.0国际许可,允许使用、共享、适应、分布和繁殖在任何媒介或格式,只要你给予适当的信贷原始作者(年代)和来源,提供一个链接到创作共用许可证,并指出如果变化。本文中的图像或其他第三方材料包括在文章的创作共用许可中，除非在材料的信用线中另有说明。如果材料没有包含在文章的创作共用许可证中，而您的预期使用不被法律法规允许或超过允许的使用，您将需要直接获得版权持有人的许可。如欲浏览本许可证的副本，请浏览http://creativecommons.org/licenses/by/4.0/．“创作共用公共领域”豁免书(http://creativecommons.org/publicdomain/zero/1.0/）适用于本文中提供的数据，除非另有用入数据的信用额度。

重印和权限