基于简化基因组测序筛选白斑狗鱼耐热性状关联的InDel 标记

海萨·艾也力汗，张 钰，杨博文，咸玉兰，高 攀，沈玉帮

(1.新疆维吾尔自治区水产科学研究所，新疆维吾尔自治区 乌鲁木齐 830000；2.上海海洋大学水产与生命学院，上海 201306)

白斑狗鱼(Esox lucius)隶属狗鱼目(Esociformes)狗鱼科(Esocidae)狗鱼属(Esox)。分布于北纬45 °以北的亚洲、欧洲和北美洲的北极圈周边地区[1]，中国境内仅分布于额尔齐斯河流域。该流域是我国唯一的白斑狗鱼种质资源库宝库。白斑狗鱼生长速率快、肉质鲜美、营养价值高，已成为新疆地区，乃至河北等地主要养殖品种。由于额尔齐斯河位于白斑狗鱼自然分布区域的最南端，因此全球气候变暖引起的极端天气对我国境内白斑狗鱼养殖和野生群体的影响较大。已有研究表明，水温超过18 °C 时白斑狗鱼受精卵孵化率下降，水温超过27 °C 幼鱼的摄食和生长受到影响[2]。因此白斑狗鱼热耐受性研究和耐高温性状的遗传改良已成为白斑狗鱼种质资源保护和产业化可持续发展亟待解决的问题之一。

基于高通量测序技术，利用全基因组关联分析(GWAS 分析)方法筛选与目标性状关联InDel标记和候选基因已广泛应用于牛[3]、山羊[4]、鸡[5]、羊[6]等动物的遗传改良研究中。在水产动物中，伞利择等[7]基于卵形鲳鲹(Trachinotus ovatus)基因组重测序的 InDel 标记挖掘与耐低氧性状关联分析，发现3 个InDel 位点接近显著阈值，注释分析获得了9 个候选基因。陈静[8]基于全基因组重测序获得了鲂属(Megalobrama)鱼类高质量的SNP、InDel 等变异位点信息，并构建了鲂属鱼类群体的基因组变异数据库。该数据库在描述基因型-表型图谱以及亲鱼选择中发挥了重要的作用。黄智康等[9]基于 InDel 标记建立了斑节对虾(Penaeus monodon)早期性别鉴定方法。竞争性等位基因特异性PCR (Kompetitive allele specific PCR，KASP)技术是英国政府化学家实验室(Laboratory of the Government Chemist，LGC)开发的检测SNP 和InDel 的第三代标记技术，具有高通量、成本低、准确性高等特点[10]，已成为利用 SNP 进行育种的首选标记[11]。

本研究为了筛选白斑狗鱼耐热性状关联的分子标记和候选基因，通过热胁迫获得2 个极端群体样品，结合简化基因组测序和全基因组关联分析筛选出与目标性状关联的InDel 位点，利用KASP 技术对显著关联位点进行验证，挖掘与白斑狗鱼耐热性状显著关联的InDel 标记和候选基因，为白斑狗鱼耐热性状的遗传改良和分子机制研究提供理论依据。

1 材料与方法

1.1 实验鱼

简化基因组测序材料：实验鱼亲鱼采自额尔齐斯河北屯段、布尔津县段以及乌伦古湖。2021年4 月12 日至4 月14 日通过人工采集和授精获得受精卵。受精卵孵化、苗种培育在同一循环系统的14 个矩形槽(长×宽×高：400 cm×60 cm×80 cm)中进行培育至幼鱼，随机挑选450 尾(平均全长90.45 mm，平均体重5.85 g)作为下一步热胁迫材料。标记验证材料：2022 年4 月通过人工繁育获得受精卵，与上一步相同的方法培育至幼鱼期，随机挑选400 尾(平均全长83.31 mm，平均体重3.96 g/尾)作为下一步的热胁迫材料。本研究在实验过程中，严格按照实验动物管理条例和使用伦理规范执行。

1.2 热胁迫实验

苗种培育至幼鱼期(2 月龄)后随机挑选450尾在一个长×宽×高为400 cm×60 cm×80 cm 的矩形槽中暂养7 d 后开展热胁迫实验。热胁迫实验期间利用3 个1 000 W 的智能PID 变频加热棒控温，实验全过程持续充气，提供足够的活体饵料，每天换1/3 相同温度的水。升温策略：暂养期间水温(26.0±0.5) °C，在30.0 °C 之前以1 °C/12 h 的速率升温，在30.0 °C 之后以1 °C/24 h 的速率升温至实验结束。每2 小时测量一次水温和溶解氧。我们对实验鱼进行了24 h 连续监测，及时捞出死鱼，记录体长、全长、体重和累计存活时间，收集肌肉组织，浸泡于无水乙醇，保存于-20 °C。最早死亡的110 尾和最后死亡的110 尾实验鱼分别被认为是热敏感组和耐高温组，进行简化基因组测序。KASP 验证材料的热胁迫实验同上，实验鱼在暂养期死亡13 尾，对387 尾实验鱼开展热胁迫实验。最先死亡的24 尾视为热敏感组，最后死亡的24 尾视为耐高温组，采集鳍条组织，用酚-氯仿法抽提基因组DNA，-20 °C 保存备用。

1.3 简化基因组测序

热胁迫实验筛选出的热敏感组和耐高温组肌肉组织送上海欧易生物医学科技有限公司，通过随机测序式基因型检测(genotyping-by-sequencing，GBS)技术进行简化基因组测序。采用限制性内切酶PstI-HF/MspI 对 DNA 进行酶切，酶切后的片段两端用 T4 连接酶加接头和 barcode；使用改进的磁珠回收系统，通过调整磁珠溶液与连接产物的体积比来回收 300～700 bp 的片段；对回收片段使用高保真酶进行 PCR 扩增；使用 Qubit 测定 PCR产物浓度，浓度需大于 5 ng/μL；将混好的文库上机(Illumina Nova,PE150)测序。

1.4 数据分析

使用软件Stacks 过滤原始测序条目，以去除接头序列和低质量碱基。去除接头序列，去除N(非AGCT)碱基大于或等于5 的读数，去除低质量读数(质量值Q ≤ 20)。利用 BWA 软件将 clean reads 比对到白斑狗鱼基因组Eluc_v4 版本(PRJNA221548)上，根据比对结果统计样品的测序深度和参考基因组覆盖度。使用 GATK4 软件的 Haplotypecaller 模块进行InDel 检测。为了降低InDel 检测的错误率，选用 QD≥2.0 的标准进行过滤，只保留同时满足该条件的突变位点。其中 QD 是突变质量值(quality)除以覆盖深度(depth) 得到的比值，实际上就是单位深度的突变质量值，大部分假阳性突变的 QD 值都小于2。使用 vcftools 进一步过滤：保留 reads 支持深度不小于 4 的位点，删除最小等位基因频率(MAF)小于 0.01 的位点，保留个体检出率大于 80% 的位点。使用 SnpEff 软件对得到的 InDel 进行注释，以确定 InDel 在基因元件的位置、对氨基酸的变化影响等。过滤后的InDel 数据利用TASSEL5.2.80 软件进行亲缘关系分析，利用R 语言可视化。使用GAPIT3 程序包进行关联分析。关联分析时白斑狗鱼耐高温性状(UTT)记录为二元表型数据，分别为0=耐高温组和1=热敏感组。前2 个PCA 为协变量，以MLM 模型(Masked Language Model)进行GWAS 分析，以假阳性检出率(false discovery rate,FDR)校正的P<0.05 为阈值，筛选显著的SNP 位点。IGV 软件可视化分析显著位点的具体位置，统计InDel 位点在每个个体中的分布情况 (表1)。

表1 MLM 模型GWAS 分析中与耐热性状显著关联的InDel 位点Tab.1 InDel sites significantly associated with heat tolerance traits in GWAS analysis of the MLM models

1.5 KASP 验证候选InDel 位点

与白斑狗鱼耐高温性状显著关联的InDel 位点前后100 bp 的序列，利用Primer premier 5.0 软件设计引物 (表2)。每个位点设计两条正向引物和一条反向引物。两条正向引物3′为InDel 位点，5′段分别加FAM(5′-GAAGGTGACCAAGTTCATGCT-3′)和 HEX (5′-GAAGGTCGGAGTCAACGGATT-3′)荧光序列标签。引物由生工生物工程(上海)股份有限公司合成。PCR 反应体系按LGC公司KASP TF V4.0 2×Master Mix 96 Std Rox 试剂盒说明配置，反应程序包括94 °C 预变性15 min；94 °C 变性20 s，63～57 °C 复性/延伸60 s (-0.6/循环)，循环10 次；57 °C 复性/延伸60 s，循环25 次；37 °C 荧光数据读取60 s。根据分型结果再次运行一下程序：57 °C 复性/延伸60 s，循环1～6 次；37 °C 荧光数据读取60 s。反应使用CFX96 Real-Time System 完成。

表2 9N del 位点KASP 引物序列Tab.2 KASP primer sequence of 9N del site

2 结果

2.1 测序结果分析

热敏感组中第1 尾实验鱼在33 °C 存活8 h 0 min 后死亡，最后1 尾实验鱼在34 °C 存活15 h 40 min 后死亡，时间跨度为31 h 40 min。耐高温组中的第1 尾实验鱼35 °C 存活6 h 15 min 后死亡，最后1 尾实验鱼在35 °C 存活13 h 56 min 后死亡，时间跨度尾6 h 15 min。

本研究利用Illumina 平台获得白斑狗鱼热敏感组109 尾和耐高温组103 尾个体的简化基因组数据。全部样本比对率范围为 91.32%～99.79%，平均测序深度为 7.98×，覆盖度范围为 10.52%～9.9%。经过突变位点检测并过滤，共获得 12 821个InDel 位点，其中插入3 683 个，缺失为9 138个。通过对InDel 位点在基因组上的位置分析，发现大部分InDel 分布在内含子，占所有InDel 位点的63.69%，基因间InDel 位点占9.80%，外显子只有1.30% (图1)。获得的InDel 平均密度为36.67 InDel/Mb，NC_047579 号染色体密度最大55.41 InDel/Mb，NC_047593 号染色体密度最小22.59 InDel/Mb(图2)。

图1 InDel 在基因组上的分布Fig.1 The distribution of InDel across the genome

图2 InDel 标记在25 条染色体上的分布密度图Fig.2 The distribution density diagram of InDel markers on 25 chromosomes

2.2 全基因组关联分析

亲缘关系分析结果显示，大多数个体的亲缘关系介于0.1～0.4，群体中亲缘关系值较均匀(图3-a)，表明群体结构对GWAS 分析结果的影响小。QQ 图检验本研究分析的变异位点的有效性，结果显示，在横坐标呈现的均匀分布近似于基因组上的随机漂变，当P-value 小于e-3时，真实的GWAS 结果与均匀分布出现分离(图3-b)，说明遗传变异与耐热性状相关性，同时表明元分析的GWAS 数据的可靠性和代表性。

图3 InDel 标记亲缘关系分析图(a)和QQ 图(b)Fig.3 Phylogenetic analysis diagram (a) and QQ diagram (b) of InDel markers

曼哈顿图呈现了GWAS 分析结果，即12 821个InDel 位点所在的染色体位置和FDR 值大小，虚线分隔为FDR 校正后的阈值FDR <0.05，实线分割线为FDR 校正后的阈值FDR <0.01，5 个InDel 位点通过阈值检验，表现出与耐热性状的显著相关性(图4)。

图4 InDel 标记的曼哈顿图数字1～25 依次代表NC_047569.1 至NC_047593.1 的染色体。Fig.4 Manhattan diagram of InDel markersThe numbers 1-25 represent chromosomes from NC_047569.1 to NC_047593-1.

实验对25 个染色体中的InDel 标记，以前两个PCA 为协变量，利用MLM 模型与白斑狗鱼的耐热性状进行了关联分析。结果发现，4 个位点与目标性状极显著关联(FDRP-value<0.01)，1 个位点为显著关联(FDRP-value<0.05)。这些位点位于未知基因LOC117593903、CLSTN2、safb基因的内含子中。其中，9N del、4N del-1 和4N del-2三个位点位于NC_047571.1 号染色体的CLSTN2基因第二个内含子中，9N del 和4N del-1 间有40个碱基，4N del-1 和4N del-2 间有246 个碱基。9N del 和4N del-1 位点在热敏感组中的DI 基因型个体为37 尾、DD 为72 尾，耐高温组中DI 基因型个体为83 尾、DD 为20 尾。4N del-2 位点DI基因型个体为34 尾、DD 为75 尾，耐高温组中DI 基因型个体为80 尾、DD 为23 尾。2N del 位点在热敏感组中的DD 基因型个体为21 尾、DI为35 尾、II 为53 尾，耐高温组中DD 基因型个体为20 尾、DI 为87 尾、II 为1 尾。3N ins 位点在热敏感组中DD 基因型个体为89 尾、DI 为20 尾，耐高温组中DD 基因型个体为94 尾、DI为7 尾、II 为2 尾。

2.3 KASP 验证

IGV 软件可视化分因型为20 尾，4N del-2 位点DI 基因型为80 尾、DD 基因型为23 尾。9N del 和4N del-1 位点在212 尾个体析显示，热敏感组中9N del 和4N del-1 位点DI 基因型个体为37 尾、DD 基因型为72 尾，4N del-2 位点DI 基因型为34 尾、DD 基因型为75 尾。耐高温组中9N del 和4N del-1 位点DI 基因型个体为83 尾、DD基中的基因型分布无差异，4N del-2 位点与前两者的差异仅为2.83%，212 尾个体中3 个位点基因型分布高度连锁。因此我们针对9N del、2N del和3N ins 位点设计引物，以验证群体DNA 为模板，利用KASP 技术进行了验证，分型结果见图5。结果显示，2N del 和3N ins 位点未通过验证，9N del 位点在热敏感组DD 基因型占37.5%(9 尾)，DI 基因型为25%(6 尾)，未知占37.5%(9 尾)，在耐高温组中DD 基因型占8.3%(2 尾)，DI 基因型为58.3%(14 尾)，未知占33.3%(8 尾)。9N del 位点DD 基因型个体在热敏感组中占优势，DI 基因型个体在耐高温组中占优势，未发现插入纯合(II)，验证结果与简化基因组测序结果基本一致。

图5 KASP 基因分型图橘色圆圈为缺失纯合FAM 基因型(DD)，绿色三角为杂合基因型(DI)，黑色棱形为对照或未知基因型。Fig.5 Genotyping graph based on KASPThe orange circle represents the missing homozygous FAM genotype(DD),the green triangle represents the heterozygous genotype (DI),and the black prisms represent control or unknown genotypes.

3 讨论

本研究对白斑狗鱼109 尾热敏感组和103 尾耐高温组个体进行简化基因组测序，经过滤后获得12 821 个InDel 位点。基因组上的位置分析发现，大部分InDel 分布在内含子(63.69%)，外显子分布的InDel 位点较少(1.30%)。这一结果符合非编码区的InDel 变异多于编码区的生物学特性[12]。GWAS 分析，未发现与白斑狗鱼的耐高温性状显著关联的编码区InDel 位点，与耐热性状显著关联的5 个InDel 位点均分布于内含子。由于内含子不是编码区序列，一度认为基因组中的“垃圾序列”。但近几年研究证实，内含子与基因表达、细胞骨架构建和动态变化密切相关[13-14]。在饥饿条件下酵母中大多数内含子可促进细胞的生存[15]，内含子保留的mRNA 不仅在细胞的正常发育中调控基因表达，还在胁迫的应答中起着重要的作用[16-17]。随着内含子生物学功能的不断证实，内含子区域的变异也被重视。虽然鱼类InDel 标记研究较少，但在牛、绵羊、山羊、鸡等经济动物研究中发现了与经济性状关联的基因内含子区InDel 位点，并影响了所在基因的表达。Shi 等[3]在发现中国黄牛品种的SMAD3 基因内含子中检测到1 个17 bp 的插入片段与SMAD3 转录水平显著相关；大连雪龙牛SREBP1 基因第5 内含子InDel变异与背膘厚具有显著相关[18]；Li 等[6]在绵羊PRNP基因中发现了4 个内含子多态性，通过关联分析结果显示，4 个内含子InDel 均与绵羊13个不同生长性状显著相关；陕西白绒山羊KDM6A基因内含子的一个16 bp 插入/缺失显著影响KDM6A基因的表达，并与陕西白绒山羊生长相关性状之间的关联[4]；Wei 等[19]在鸡LDB2 基因第2内含子区域中发现了一个31 bp InDel 突变位点，通过关联分析发现该InDel 与多个生长性状和胴体性状显著相关，并影响LDB2 基因在肌肉组织中的表达。由此推断，本研究中发现的5 个InDel 突变可能会对白斑狗鱼的耐热性状产生显著的影响，可作为白斑狗鱼耐热性状改良的候选分子标记。

对所发现的5 个与白斑狗鱼耐热性状显著关联InDel 位点注释和IGV 可视化结果显示，2N del 位点位于未知基因LOC117593903 内含子中，3N ins 位于safb基因内含子中。与目标性状极显著关联(FDR <0.01)的9N del、4N del-1 和4N del-2 三个InDel 位点均位于CLSTN2 基因第3 内含子，并9N del 和4N del-1 间仅有40 个碱基，4N del-1和4N del-2 间有246 个碱基。上述3 个位点在212 尾个体中的基因型分布高度连锁。在验证群体中KASP 基因分型结果显示，9N del 位点DD基因型个体在热敏感组中占优势，DI 基因型个体在耐高温组中占优势。有意思的是我们[20]之前发现的白斑狗鱼耐热性状显著关联的SNP 位点HT1位于4N del-1 和4N del-2 间。钙离子作为细胞内重要的第二信使，参与包括盐胁迫、冷和热刺激等非生物胁迫应答，并胁迫信号能改变细胞内钙离子水平[21]。过高的体温具有直接的细胞毒性，影响细胞膜的稳定性和跨膜转运蛋白的功能。因此，跨膜离子转运被中断，导致胞内钙离子的浓度增加。蛋白质和DNA 的合成也在该途径中不同阶段被中断[22]。CLSTN2(Calsyntenin-2,CLSTN2)是钙结合蛋白家族(Calsyntenins，CLSTNs)成员之一，是I 型跨膜蛋白，分布于脑、心脏、骨骼肌、肝脏、胰腺、胎盘和肺组织中[23]。CLSTNs对Ca2+有高度亲和性的蛋白质，其能够通过与Ca2+结合参与Ca2+的运输、吸收与分泌[24]。此外，白斑狗鱼CLSTN2 基因第2 内含子为超长内含子(34 365 bp)。王玲平[25]认为，超长的内含子通常采用外显子限定的剪接模式，因此更容易发生可变剪接；此外，由于内含子越长，伴随的基因转录时间也就越长，转录时受到干扰而终止的可能性也就越高，因此越长的内含子，其处理的不可预测性就越大。由此推断，CLSTN2 可能是白斑狗鱼耐高温性状相关的重要候选基因。与白斑狗鱼耐热性状显著关联的3 个InDel 和1 个SNP 位点均位于CLSTN2第二内含子，而CLSTN2 第二内含子为超长内含子，容易发生可变剪切和转录终止从而影响CLSTN2 基因的转录。

（作者声明本文无实际或潜在的利益冲突）