新疆核桃‘新露’与‘新萃丰’全基因组重测序分析

武鹏雨，李冬，吴翠云，包建平，金强，虎海防 ，张锐*

（1塔里木盆地生物资源保护利用兵团重点实验室，新疆 阿拉尔 843300）

（2南疆特色果树高效优质栽培与深加工技术国家地方联合工程实验室，新疆 阿拉尔 843300）

(3塔里木大学园艺与林学学院，新疆 阿拉尔 843300)

（4新疆佳木果树学国家长期科研基地，新疆 温宿 843100）

核桃（Juglans regia L.）隶属于胡桃科（Juglandaceae）、核桃属（Juglans），是世界4大坚果（核桃、扁桃、板栗、腰果）之一，是不可或缺的坚果和木本油料树种［1］。我国核桃栽培历史悠久，是核桃属（Juglans）植物的起源和分布中心之一［2］，新疆作为我国核桃重要产区之一，核桃栽培面积达3.5×105hm2，其中塔里木盆地周围绿洲产量最高，位居全国前列［3］。新疆栽培核桃资源丰富，拥有野核桃和纸皮核桃等种质资源，也有薄皮核桃如‘扎343’‘新新2号’‘新萃丰’等品种，壳薄如纸，透光可见壳内果仁形状，薄皮核桃硬壳的平均厚度一般小于1.5 mm［4］，在薄皮不露仁栽培种中壳厚最厚的是‘新萃丰’，厚度为1.45 mm［5］；而‘新露’是天然内果皮发育不全的露仁核桃品种，其内果皮的壳面在生长发育过程中发生自动破裂使核桃仁露出，果仁裸露虽方便食用，但不耐运输储藏。对拥有早实、薄壳、抗逆性强的‘新萃丰’和天然果壳裸露的‘新露’进行研究，可以使新疆核桃内果皮发育的遗传背景清晰化，为新疆核桃的保护与科学利用提供理论依据。

随着测序成本的降低和已知基因组序列物种的增多，基因组测序已成为研究植物分子育种、群体进化中最为迅速有效的方法之一［6］。基因组测序技术为不同物种提供了高质量的参考基因组序列，基因组重测序技术根据供参考的高质量序列进行整个基因组、区域或基因的重新测序，每个个体的测序数据与高质量的参考数据进行比对，不同样本之间的遗传变异可以被检测为高度可信的序列差异，并在此基础上对个体或群体进行比较基因组学分析，可以从基因组水平对重要性状的候选功能基因或QTL进行定位与分析［7］，通过与参考基因组比较获得多态性插入与缺失（InDel）标记，使用这些分析标记用于基因鉴定［8］，构建杂交群体使用基因组重测序单核苷酸多态性位点（SNP）数据库筛选与目的形状相关基因［9］。本研究通过Illumina NovaSeq 6000测序平台对新疆栽培核桃中的‘新露’和‘新萃丰’进行全基因组重测序分析，对其SNP和InDel等变异位点进行深度挖掘，从基因组层面分析新疆核桃内果皮发育过程中出现的果仁裸露这一性状的机理，通过检测变异位点所在区域，为后续内果皮发育等性状的相关基因的挖掘、品种选育等研究提供一定的基础。

1 材料与方法

1.1 试验材料

试验所选‘新露’与‘新萃丰’品种材料种植于新疆佳木果树学国家长期科研基地，2019年选取生长情况一致的4年生核桃植株各4株，在5月份新出嫩叶时采摘叶片，放入液氮中保存带回实验室。

1.2 试验方法

‘新露’与‘新萃丰’的DNA提取方法按照DNA提取试剂盒的使用说明操作（TianGen）。DNA样品的质量控制通过0.8%琼脂糖凝胶电泳检测，同时采用紫外分光光度计对DNA进行定量分析。

1.2.1 测序

对提取的DNA进行超声随机打断，随后进行末端损伤修复及连接接头，在与测序接头连接后进行PCR扩增，对测序文库模板进行富集并纯化文库产物，连接产物通过2%琼脂糖凝胶电泳进行最终片段选择与纯化，选择插入片段大约为400 bp回收后根据制造商建库流程构建测序文库。在Illumina Nova-Seq 6000测序平台上对构建的文库进行双端（PE 150 bp）高通量测序。

1.2.2 数据分析

获得原始数据后采用FastQC将数据进行质量控制，进而获得高质量数据。使用Burrows-Wheeler-Alignment Tool软件将高质量数据比对到已测序完成的核桃参考基因组，参考基因组数据为：JrSerr v1.0(GCA_004785585.1)［10］，利用 picard 软件将比对得到的sam结果文件转化为bam文件，使用GATK软件对SNP和InDel进行检测，并用Unified Genotyper对变异位点位置进行获取，使用ANNOVAR软件对获取的SNP和InDel进行功能注释。获取‘新露’与‘新萃丰’相同位置的SNP及InDel后，再根据非同义SNP（nonsynonymous，nsSNP）/InDel获取相关基因。候选基因的富集分析递交于AgriGO软件用于富集Gene ontology terms。

2 结果与分析

2.1 测序数据预处理

以Illumina Nova Seq测序平台提供的初始测序数据为原始数据（Raw Data），即本次测序中获得的各样本得到的短序列数（Reads）和碱基总数（Bases），对2个核桃样本产出的数据进行汇总（见表1）。2个品种共得到232 715 892个短序列，其中‘新露’核桃短序列数目127 572 756个，Q30比为91.30%；‘新萃丰’核桃中Reads数目105 143 136个，Q30比为91.40%。为剔除Illumina平台错误率对结果的影响，需对原始数据进行质量控制，包括去除低质量序列，去除接头，质量过滤和长度过滤等获得纯化数据（Clean Data），以进行后续工作。质量控制后‘新露’与‘新萃丰’分别获得121 633 052个和100 378 460个短序列，测序总有效短序列与初始下机测序的比例分别为95.34%，95.46%。

表1 测序数据汇总

2.2 短序列匹配统计

将过滤后获得的高质量数据比对到核桃参考基因组，统计结果表明：‘新露’中总匹配的短序列数为121 332 322个，占所有短序列数的98.1%，重复序列占比15.15%，覆盖全基因的深度为26.76 x，当覆盖深度≥20时，覆盖全基因组的百分比为75.32%；‘新萃丰’中总匹配的短序列数为99 895 457个，占所有短序列数的97.96%，重复序列占比14.02%，覆盖全基因的深度为22.45 x，覆盖深度≥20时，覆盖全基因组的百分比为59.92%。

2.3 SNP统计分析

‘新露’与‘新萃丰’核桃样品检测后获得SNP位点统计信息（见表2），得到变异位点所在区域及变异影响。‘新露’样本的总SNP数量为4 700 220个，杂合SNPs数量为2 645 454个，纯合SNPs数量为2 054 766个，编码区非同义突变数量（nsSNP）为81 580个，密度占比1.74%，SNP位于内含子区域内（Intronic）总数为535 307个，在基因间区域内总数（Intergenic）为3 469 438个。

表2 ‘新露’与‘新萃丰’SNPs数目与类型检测结果

‘新萃丰’样本经质量过滤后获得4 675 346个高质量的SNP位点，编码区内非同义突变总数为83 811个，SNP位于内含子区域内总数为541 875个，位于基因间区域的SNP占比73.2%。

2.4 InDel统计分析

在基因组中检测长度小于50 bp的小片段插入或缺失（InDel）并进行变异位置及类型的注释，对‘新露’与‘新萃丰’的InDels变异位点信息进行统计（见表3）。‘新露’的总InDel数量为648 603个，编码区内移码插入总数（Frame shift insertion）和移码缺失（Frameshift deletion）数量分别为3 560个、2090个，InDel在基因5’UTR内总数（5’UTR），3’UTR内总数（3’UTR）和InDel位于不同基因的5’UTR和3’UTR内总数（5’UTR/3’UTR）为3 095个、5 912个、2个，InDel位于内含子区域内总数为106 672个，在基因间区域内总数为420 118个。

表3 ‘新露’和‘新萃丰’InDels数目与类型检测结果

‘新萃丰’鉴定出606 188个InDel位点，其中插入297 690个，缺失311 822个，编码区内移码插入总数和移码缺失数量为3 449个、2 015个。InDels大部分位于基因间区域，有389 525个（64.2%），次之为内含子 202 381个（16.89%），外显子 8 048个（1.33%）。

2.5 DNA水平变异的基因分析

发生在外显子区的变异可能会引起基因功能的异常改变，通过与参考基因组比对，检测到‘新露’和‘新萃丰’基因组间发生非同义突变、同义突变和移码突变等信息，分别获得16 686个、19 829个变异SNPs以及1 687个、2 143个InDels，将这些变异位点比对到Gene ontology数据库中。

2.5.1 ‘新露’核桃变异基因的GO富集分析

‘新露’核桃的SNP和InDel变异位点共注释到4 865个基因，富集分析注释结果分为生物过程、细胞组分和分子功能三类（如图1）。注释结果中，草酸代谢过程、囊泡运输与锰离子结合分别在生物过程、细胞组分和分子功能类别中被显著富集。在这些被显著富集的三类反应中草酸代谢过程在植物体内广泛存在，作用于植物抗逆过程中；高尔基体的囊泡运输参与到植物细胞壁的形成，说明露仁性状可能和次生细胞壁形成有关；而锰离子则是植物体细胞中许多酶的催化剂。

图1 ‘新露’变异基因的GO注释分类图

2.5.2 ‘新萃丰’核桃变异基因的GO富集分析

‘新萃丰’核桃共注释到6 131个基因，GO terms富集分析如图2所示，在生物过程和细胞组分中微管过程和囊泡运输被显著富集，用于维持细胞骨架和信号传导。在分子功能中富集到了天冬氨酸和苯丙氨酸的2种转移酶，天冬氨酸转氨酶，苯丙氨酸转氨酶在植物体内分别作用于苯丙烷途径前后，最终目的为合成木质素的单体。为进一步对变异基因的生物学功能进行了解，对富集到的天冬氨酸转氨酶和苯丙氨酸氨基转移后的裂解反应的基因进行挖掘，筛选出7个相关基因，结果见表4。

表4 木质素合成相关变异基因

图2 ‘新萃丰’变异基因的GO注释分类图

3 讨论

随着测序技术的快速发展，苹果［11］、柑橘［12］、梨［13］等多种果树的全基因组序列陆续公布，不同果树物种的基因组信息逐渐被完整破译。通过全基因组重测序可以找到大量的SNP、拷贝数变异（CNV）、InDel、结构变异（SV）等变异信息，基于检测到的变异位点能进一步研究该物种的特性［14］、群体进化［15］和定位目标性状基因位点［16］。

SNP在基因组中广泛存在，并含有丰富的遗传信息，随着测序技术的发展以及生物信息学分析能力的提升，SNP作为第三代分析标记应用于相关性状的基因定位研究中［17］。SNP的最早研究用于人类基因组，最终绘制成142万个SNP的人类基因组图谱，这种高密度的SNP图谱为整个基因组单倍型变异提供了依据［18］。将2个橄榄品种杂交鉴定获得10 941个SNP，并使用其中多态性的SNP位点构建覆盖基因组3 049 cM，25个连锁群体的高密度遗传连锁图谱，遗传连锁图谱的构建可更有效地定位数量性状和其他重要经济性状［19］。使用基因测序方法进行基因分型，对180种枣进行全基因组关联分析（GWAS），鉴定出4 651个高质量SNP和45个与SNP性状关联基因，为标记辅助育种提供了数据基础［20］。将SNP作为新型分子标记也应用于核桃分子育种中，将14 761个SNP进行全基因组关联，确定了叶芽萌发和花芽分化等性状的标记位点，为后续QTL鉴定奠定了基础［21］。在自然适应和选择性育种过程中，最常见且研究最多的SNP外还有InDel，InDel变异是不同个体间基因组同一位点的序列发生不等大小片段的插入与缺失，从而使同源序列比对产生空位（gap）的现象［22］。有研究者将SNP、InDel、SV以及简单序列重复(SSR)组装成核桃基因组变异图谱，利用遗传信息丰富的标记位点对中国五个核桃品种进行系统发育研究，将5个核桃种以100%的bootstrap(BS)值划分为两个已知区段核桃/胡桃组(Juglans/Dioscaryon 和核桃楸组 Cardiocaryon)［23］。在核桃上使用454焦磷酸测序技术开发鉴定获得48 165个SNPs和1 037个InDels，转换/颠换（ts/tv）比例为2.79:1，略高于本次新疆核桃测序的比例［24］。本研究选用具有代表性的新疆露仁品种与薄皮中的厚壳品种的‘新露’与‘新萃丰’，利用基因组重测序技术检测核桃内果皮果壳发育不全，果仁裸露这一性状的各种遗传变异。‘新露’与‘新萃丰’的总SNP数目为4 700 220个，4 675 346个，总InDel数目为648 603个，606 188个，变异位点在基因编码区，基因间区域和内含子区域都有分布，在全基因组水平上有着显著的遗传变异，2者的SNP数量显著多于InDel的变异数量，说明外界环境对基因组的变化主要以单核苷酸多态性变异为主，这与梨［25］、马铃薯［26］等测序结果一致。本研究将检测到的‘新露’与‘新萃丰’核桃基因组间所发生的非同义突变、同义突变等基因比对到GO数据库中，发现在薄皮不露仁系列中的厚壳性状品种‘新萃丰’富集到与天冬氨酸转氨酶和苯丙氨酸转氨酶的相关过程，通过与基因组序列比对分析，筛选到7个与木质素生物合成有关的基因。本试验通过全基因组重测序技术对新疆2个核桃栽培种进行变异检测分析，通过比对变异位点，挖掘变异的关键基因，为基因定位和克隆功能验证以及为核桃分子育种的利用奠定了基础。