基于简化基因组测序分析云南小桐子遗传多样性

郑周宜　刘雨婷　陈春　王莎莎　唐唯

关键词：小桐子；遗传多样性；分子标记；全基因组关联分析；单核苷酸多态性

中图分类号：S565.9 文献标识码：A

小桐子（Jatropha curcas L.）属大戟科麻疯树属的多年生木本油料植物，原产美洲，现广泛分布于全世界热带及亚热带地区，在我国尤以云南分布面积最广、资源数量最大[1-2]。小桐子具有生物产量大、适应范围广、耐贫瘠等特性，是维持生态[3]、医药开发、生物防治及生产有机肥料的理想资源[4-5]。由于其种子含油率高达60%以上，也被认为是具有巨大开发潜力的多用途能源植物树种[6]。尽管近年来国内外有关小桐子的研究报道数量大增，但大多集中在药理毒理活性及新药开发[4-5]、生理適应机制[7-8]、生物柴油产业化[6]、经济用途开发[3]等方面，关于其种质资源的遗传多样性研究起步相对较晚。另一方面，种质资源不清、品种良莠不齐、基因资源库保存不完善等问题已成为限制国内小桐子产业开发的重要瓶颈[9]。

迄今为止，来自中国、巴西、印度、墨西哥等多国的研究人员已先后采用RAPD、AFLP 等多种分子标记对国内外小桐子种质进行了遗传多样性的鉴定和评价[10-11]，发现我国范围内收集的小桐子种质资源遗传背景虽普遍低于它国（南美、东南亚地区等）[12-13]，但云南地区的小桐子居群具有较为丰富的遗传多样性，且其居群内的遗传分化高于居群间[14-16]。上述研究在一定程度上印证了向振勇等[17]提出的云南小桐子种质来源较为复杂的推测。

基于上述溯源，为了进一步了解近年云南地区小桐子种群结构特点，本研究首先从种子和幼苗性状特点入手，结合简化基因组测序和分子标记技术解析云南各地区小桐子的遗传多样性特点；利用GWAS 分析方法，挖掘地下部分鲜重这一重要性状差异基因及地区独特性基因，预期研究结果可为小桐子种质资源的保护利用、遗传改良以及杂交育种亲本选配提供重要遗传基础。

1 材料与方法

1.1 材料

供试材料来自云南省不同地区的9 个居群，种源地信息见表1。每个地区根据植株分布情况随机选取50 个具有代表性的样株，每个样株至少间隔5 m 以上，每个样株分别取3 颗种子，编号整理。

1.2 方法

1.2.1 幼苗培养 在各种源地的种子材料中挑选颗粒饱满、大小均一的种子50 粒，清洗净后以1% CuSO4 溶液浸泡消毒20 min，经蒸馏水充分漂洗、浸泡、吸涨后于26℃萌发。将萌发14 d 后的种子移栽至装有等量营养土的花盆中，置于人工气候箱中继续培养14 d，用于幼苗生物学性状的测定。

1.2.2 性状测定 （1）百粒重及种子形状测定。

在各种源地随机选取种子，取3 个重复，每个重复30 粒，以电子天平分别称重（精确至0.01 g），计算种子百粒重。以游标卡尺测定长径、宽径及侧径（精确至0.01 mm）。

（2）生长指标及生物量测定。每个种源随机抽样10 株苗，于移栽后14 d 进行株高测定（精确至0.01 cm）。生物量测定采用全称重法。样株在流动自来水下充分洗净、吸干表面水分，分别对地上部分、地下部分进行鲜重称量。

（3）数据处理。以SPSS 软件（v26.0）进行数据统计、方差分析、多重比较（Duncans 法或非参数检验，以不同小写英文字母表示显著差异，P<0.05）及相关性分析。

利用变异系数（CV）分析种子性状变异性，计算公式为：CV=（S/X）×100%（S 为标准差，X 为平均值）。计算根冠比=地下部分（鲜重）/地上部分（鲜重）×100%。

1.2.3 RAD-seq 基因组测序及分子标记开发 （1）基因组DNA 的提取。综合考虑幼苗性状（鲜重、根冠比）及平均变异系数，取差异显著的2 个居群材料（各10 株）为BSA 分离群体，其中版纳勐海BNMH 为BSA1，楚雄元谋CXYM 为BSA2，采用CTAB 法提取基因组DNA（以第2 片真叶为材料），经琼脂糖凝胶电泳（0.8%～1.0%）和超微量分光光度计检测质量、纯度和完整性。

（2）基因组测序及分子标记挖掘。构建插入片段350 bp 的文库，测序在ILLUMINA Hiseq2000 平台进行（Novogenes Inc.），测序读长为双端150 bp，测序序列用GenomeScope 2.0[18]进行K-mer 分析确定染色体倍性。SNP 标记挖掘由CLC Genomics Workbench 12.0.1（QIAGEN， Inc.）完成，工具为Basic Variation Detection（参数默认）。

（3）SSR 引物设计。取上述2 个居群的小桐子材料进行RAD-seq 基因组测序（5×覆盖度），以测序基因组为基础，利用基于Perl 的MISA 程序挖掘SSR 引物[19][（X）n；（XX）n-2；（XXX）n-3；（XXXX）n-4；（XXXXX）n-6；（XXXXXX）n-7，n=10] ， 挖掘到的SSR 标记以软件primer3（https：//primer3.org/）设计引物。引物位于SSR标记2 个侧翼，长度在20～22 bp 之间，正向引物的5端加上荧光标记FAM，PCR 目标产物预期大小在50～150 bp 范围内。引物用QC filter 进行筛选，去掉纯合位点，引物由昆明擎科生物科技有限公司合成。

（4）分子标记的电泳检测（琼脂糖凝胶、毛细管电泳）及遗传多样性分析。对检测和纯化好的不同种群小桐子基因组DNA 进行筛选，配制反应体系后进行PCR 扩增，Tm 以引物合成报告单为准。对琼脂糖凝胶电泳初步确定扩增效果理想的引物进行毛细管电泳分析（ABI 3500），用GeneMapper 5.0（https：//tools. thermofisher.com/content/sfs/manuals/4476603A.pdf）读取条带，用MEGA 5.0[20] 构建UPGMA 系统发育树， 用PopGene 1.31（https：//sites. ualberta.ca/～fyeh/popgene.pdf）进行遗传多样性计算。

（5）基于重要性状差异的GWAS 分析。以地下部分鲜重统计数据为性状表型，以全基因组SNPs 为基因型，使用R 包GAPIT 混合线性模型（MLM）进行GWAS 分析（全基因组关联分析，Genome Wide Association Study），以CMplot 包进行GWAS 结果输出，分别生成SNPs 密度分布图和柱状曼哈顿图。

（6）基因组注释及比较基因组分析。对上述重测序2 份材料的clean data 进行全基因组比对（CLC Genomics Workbench 12.0.1）得到共线性结果。对2 个基因组进行de novo 组装（SOAPdenovo2[21]），得到的组装序列用Augustus（http：//bioinf.uni-greifswald.de/augustus/）进行CDS 预测。对预测到的基因以NCBI 中protein databank 和swiss-prot protein databank 进行同源注释和比较分析，找出同时具有以下2 个特征的特异蛋白：①仅在上述2 份材料中存在的；②与地下部分鲜重性状紧密连锁。

2 结果与分析

2.1 云南地区不同种源小桐子种子性状及其变异分析

百粒重测量结果（表1）显示其变幅在46.06～80.40 g，不同种源的差异可达显著差异水平（P<0.05）。长径、宽径、侧径变幅依次为13.00～19.50、7.50～12.00、5.20～10.10 mm，最大平均值均为双江大文乡种源，最小平均值为丽江华坪种源。对种子性状的变异性分析发现，供试种子变异系数差别较为明显， 其中百粒重变异系数最大（13.93%），长径、宽径与侧径的变异系数依次为3.04%、5.34%和6.20%。

2.2 云南地区不同种源小桐子幼苗性状及其变异分析

从表1 可以看出，株高性状表现最好的是临沧耿马种源，平均苗高为11.46 cm，版纳勐海种源具有相对最高的生物量积累（ 全株鲜重6.60 g），丽江华坪、楚雄禄丰和楚雄元谋种源在苗高和生物量积累上均表现最差，尽管楚雄元谋种源的根冠比表现最佳（0.25）。方差分析表明不同种源的苗高、地上部分鲜重、全株鲜重及根冠比均存在显著差异。不同种源地小桐子表型性状变异程度依次为：地上部分鲜重＞全株鲜重＞地下部分鲜重＞根冠比＞株高＞百粒重＞种子侧径＞种子宽径＞种子长径（表1），尤以生物量积累的性状变异系数相对较高，暗示这些表型性状存在遗传不稳定性，易受环境压力的选择。

2.3 RAD-seq 基因组测序及结果分析

测序结果显示BNMH（BSA1）及CXYM（BSA2）的染色体为2 倍体（图1A），得到具有EcoRⅠ酶切位点的Reads 数分别为64 734 412 条和60 187 766 条。平均Q20 为98.6%，平均Q30为93.9%，对reads 进行有参组装（参考基因组RJC1）后分别得到2783 和3033 条高质量的contigs（表2）。

2.4 SSR 引物开发及遗传多样性分析

利用MISA 随机开发出8 对SSR 引物，PCR扩增后利用琼脂糖凝胶电泳确认有效扩增引物为4 对（表3）；以4 对引物、9 个居群共95 个DNA样本进行毛细管电泳检测，结果发现引物多态性扩增条带最多为6 条，最少为2 条，各居群间利用4 个毛细管电泳引物扩增遗传多样性水平均较高（表4）。以每个地区所有个体多态性条带为居群进行UPGMA 聚类分析，结果表明保山地区和西双版纳地区可能作为云南省小桐子的亲本或起源，各居群间亲缘关系和地理位置有一定相关性（图1B）。

2.5 基于地下部分鲜重差异的GWAS分析

基于在种子及幼苗性状调查分析中发现，地下部分鲜重（underground fresh weight， UFW）为9 个居群间差异最大，因此本研究针对UFW，首先在全基因组水平挖掘到相关InDel 标记8182个，SNPs 标记82 155 个。因小桐子为杂合二倍体，因此经纯合和等位基因频率（≥50%）过滤后，共挖掘到有效SNPs 22 756 个。GWAS 分析进一步发现，SNPs 分布于11 对染色体（图2A），UFW 高频SNPs 集中于9 号染色体（图2B）。

2.6 注释及比较基因组分析

本研究通过对基因组测序材料BNMH（版纳勐海）和CXYM（楚雄元谋）进行注释和功能预测，分别得到27 133 和26 282 个候选功能蛋白。Blast2GO 分析对预测蛋白进行GO 功能注释，发现在生物学途径（biological process）中，获得注释条目最多的是蛋白修饰、跨膜运输和转录调节；在分子功能（molecular function）和细胞组件（cellular component）中分别为催化活性、结合功能、膜整合组分和胞质（图3）。

为了初步明确云南本地小桐子的独特性，进一步比较分析发现，仅在BNMH 和CXYM 两份测序材料中存在的候选蛋白共40个，其中未注释的功能蛋白有16 个，和UFW 性状连锁的有3 个（表5）。

3 讨论

小桐子（J. curcas L.）是一种在全世界（亚）热带地区广泛种植的多用途油料作物，其作为一种新兴能源植物的开发潜力而引起了全世界的广泛关注，然而，现有的育种种质在遗传多样性上的匮乏成为小桐子完全商业化的主要阻碍，而种质资源的收集和评估对于培育高产、优质和高适应性的小桐子新品种至关重要。迄今为止在研究多样性的各种分子标记中，基于PCR 的标记如RAPD（随机扩增多态性DNA）、ISSR（简单序列间重复）、SSR（简单序列重复）、EST-SSR 和AFLP（扩增片段长度多态性）均已被用于探究小桐子的遗传多样性和系统发生。PECINAQUINTERO等[22]借助AFLP 技术对来自墨西哥中部和东南部9 个州的175 份种质（包括有毒性、无毒性、蛋白质含量和种油含量对比大的基因型）展开全面的遗传多样性分析，发现墨西哥小桐子具有很高的遗传多样性（明显高于GRATIVOL 等[23]对巴西来源小桐子的相关报道），并伴随着不同州的种群结构。SUDHEER 等[24]使用RAPD 和AFLP技术对42 个印度小桐子种质的分子多样性分析，指出其多态性的平均百分比为26.47，显著低于墨西哥（33.18）。这些结果都支持了目前的主流观点，墨西哥和中美洲可能是麻疯树起源和多样化中心[23-24]，而非洲和亚洲的小桐子种质存在低的遗传变异性[11， 25-26]。WEN 等[13]早期利用36 个EST-SSR 和20 个G-SSR 从来自南美洲、中国和印度尼西亚的45 个种质中检测到183 个多态性等位基因，平均遗传多样性指数为0.5572，发现中国云南的种质依然显示出比其他地理区域（如格林纳达和中国海南）更高水平的遗传变异，可作为丰富小桐子遗传背景的重要种质资源。在本研究中，小桐子种子采集自云南省的不同地域，并首先用于多项形态学或农艺性状分析，包括种子径围、百粒重、幼苗株高和生物量积累等。这些性状对作物育种改良研究很重要，但应用于小桐子遗传多样性的系统研究还很少见到报道。本研究还发现小桐子种群内存在显著性状差异，尤以生物量相关性状的变異系数最大（ 21.90%～25.87% ）， 这与PESTANA- CALDAS 等[27] 和LAVIOLA 等[28]的研究相似，均认为单株植物种子产量（包括种子百粒重）、植株高度对小桐子遗传变异性贡献最小，而生长有关的特征则表现出最大的表型变异，这可能和地区适应性和驯化有密切关系。

有研究表明在小桐子中，分子标记和性状数据连锁不紧密[28-29]。经典的遗传定位依靠基因型和表型的紧密连锁[30]，结合基因组测序的BSA 分析法扩大了群体内一致性数据和群体间差异性数据[31]，本研究选择地下鲜重这一极端差异性状的2 个分离群体，BNMH（版纳勐海）和CXYM（楚雄元谋）以BSA 分离法进行全基因组测序，分别组装获得2783 和3033 条contigs 并通过GWAS对重要SNP 标记的基因组区域进行高精度的关联和标记[32]，在9 号染色体上确定了5 个UFW 高频SNPs 的峰值区域，有3 个候选蛋白与UFW 性状连锁。从基因组测序和比较分析来看，云南省小桐子具有较高的单核苷酸多态性，具有特有蛋白40 个（含未注释功能蛋白16 个），显示了较好的多样性、独特性和地区适应性。此外，地下鲜重主要涉及根系生长和发育，这一性状的差异也明显和地区驯化有关。

本研究利用MISA 开发并随机选出8 对SSR引物（有效扩增引物4 对），对9 个居群共95 个DNA 样本进行毛细管电泳检测，结果发现引物多态性条带最多为6 个，最少为2 个，引物位点的多态信息含量指数（PIC）变幅为0.041～0.711（平均为0.370），高于V?SQUEZ-MAYORGA 等[33]（PIC=0.274±0.165）和KUMAR 等[34（] PIC=0.25±0.16）的报道，可作为小桐子种质亲缘关系鉴定、遗传变异检测等分子生物学研究的重要可用资源。此外，群体结构研究虽表现为中度多态性[35]，但本研究发现云南小桐子的杂合度（ 平均He=0.437）与目前世界范围内的小桐子报道相比十分可观，明显高于WEN 等[13]对南美洲和云南种质资源的报道（He 分别为0.33、0.3473）。多数研究认为小桐子起源于中美洲，但传入我国时间和路线不详[36]。本研究基于UPGMA 聚类分析发现各居群间亲缘关系和地理位置有一定相关性，各居群间遗传多样性较高，这与WEN 等[13]、RAFII等[37]的结论相一致，暗示小桐子可能由多点多分化种传入云南，传入后在各生长地之间基因交流较少，变异主要为独立群体演化事件，其主要结果为充分了解小桐子在中国的起源中心的遗传多样性和种群结构提供参考，也为其作为可持续能源作物的复兴奠定新的理论基础。在后续研究中应进一步结合胞质基因组多态性特点及相关历史考据文献对小桐子在云南的传入、传播及分布特点进行更加深入的探讨。此外，小桐子作为云南省的一个主要的外来物种，在将来的育种工作中有必要加大对云南省（尤其是边境地区）的小桐子种质资源的勘探与发掘，结合InDel 等大片段标记[34， 38]，开发适用于种质资源检测的各类有效标记，进一步通过分子设计育种加速小桐子的种植资源改良。

4 结论

本研究采集了云南省10 个小桐子聚居地种子，调查了9 个重要农艺性状，利用地下鲜重差异显著的2 个BSA 分离群体，通过RAD-seq 测序结合GWAS 分析挖掘到SNPs 22 756 个，分析发现高频SNPs 集中在9 号染色体，进一步定位到来自于版纳和楚雄群体的特有候选蛋白40 个，筛选了3 个UFW 候选蛋白。此外，研究基于4對SSR 标记对9 个地区小桐子种质资源进行了群体结构分析，得到了较好的多态性和推测了可能的起源中心。