多鳞(Sillago sihama Forsskål)又称沙锥鱼, 隶属鲈形目、科、属, 为热带印度-西太平洋浅海鱼类, 广泛分布于我国沿海地区(Lin et al, 2021)。多鳞肉质细嫩, 味道鲜美, 营养价值高, 具有重要的经济价值。该物种曾在我国近海渔业中占有重要地位, 近年来由于过度捕捞, 其自然个体低龄化、群体趋小化和产量逐年减少等问题突出(黄洋等, 2013), 市场售价已从2003年的20~30元/kg持续上升到目前的100~110元/kg。本课题组前期突破了多鳞全人工繁育技术, 实现了其全人工规模化育苗, 为下一步多鳞良种选育奠定了基础(黄洋等, 2013)。然而, 多鳞养殖苗种缺乏定向选育, 存在因近亲繁殖导致的生长速度减慢、抗逆性能下降等性状退化风险, 严重限制了其养殖业的健康发展, 亟需开展分子遗传育种研究, 解析重要经济性状相关遗传机制。近年来, 研究人员在多鳞低氧应激(Tian et al, 2020; Pan et al, 2021)、营养饲料(Huang et al, 2020)、群体遗传(Liu et al, 2012)、人工繁育(余家旺等, 2022)以及基因组资源(Lin et al, 2021)等方面进行了探索研究, 但尚未见分子标记开发及生长性状相关QTL定位研究报道。

基于分子标记技术构建遗传连锁图谱是实现重要经济性状QTL定位和分子标记辅助育种的有效途径(桂建芳等, 2016)。高通量测序技术简化了遗传标记分型方式, 降低了实验成本, 使得在众多水产养殖物种中通过构建高密度遗传图谱进行性状QTL精细定位成为可能(叶华等, 2011)。鱼类生长性状是重要的经济性状, 对生长性状进行遗传改良是最有价值的选育目标之一。借助简化基因组测序技术, 截至2020年已有40多种鱼类构建了高密度遗传连锁图谱, 并开展了经济性状QTL定位研究(You et al, 2020), 如鳙(Hypophthalmichehys nobilis) (Liu et al, 2016)、鲤(Cyprinus carpio) (Peng et al, 2016)、褐石斑鱼(Epinephelus bruneus) (Kessuwan et al, 2016)、大口黑鲈(Micropterus salmoides) (Dong et al, 2019)、锦鲤(Cyprinus carpio haematopterus) (Feng et al, 2018)以及斑点叉尾 (Ictalurus punctatus) (Zhang et al, 2019), 鉴定出了一批生长性状相关的遗传标记及其候选基因(陈军平等, 2020)。

目前, 多鳞生长相关性状的定位研究尚属空白, 本研究的目的在于利用基因分型测序(genotyping-by-sequencing, GBS)技术构建多鳞高密度连锁图谱, 并进行6个生长性状的QTL定位分析, 以期为多鳞基因组组装、比较基因组学分析以及分子标记辅助选择育种提供有力支撑, 并为该物种生长调控机制研究奠定基础。

1 材料与方法 1.1 作图家系及DNA提取

试验所用多鳞亲本及其F₁代全同胞家系养殖于广东省湛江市广东海洋大学东海岛繁育基地。2018年6月, 挑选体型良好的性成熟雌雄个体为亲本, 一对一繁殖构建全同胞家系。孵化后的幼苗在水泥池(5.8 m×4.8 m×1.8 m)中进行养殖, 在整个养殖过程中, 水温(27±0.5) ℃、pH 7.2±0.6及溶解氧(7.1±0.5) mg/L。2019年6月采集161尾F₁代个体及父母本鳍条, 用于基因组DNA的提取; 同时, 采集161尾个体的体重、体高、体厚、眼径、体长和背鳍前长6个生长性状[方法见余家旺等(2022)]。随后, 将样本送至广州基迪奥生物科技有限公司使用CTAB (十六烷基三甲基溴化铵, 广州化学试剂厂)法提取基因组DNA (李进波等, 2014); 并通过Nanodrop 2000 (Thermo Scientific, 美国)和1.5%琼脂糖凝胶电泳对提取的DNA进行质量检测。将符合试验要求的各样品DNA在广州基迪奥生物科技有限公司建库。

1.2 GBS文库构建及高通量测序

多鳞亲本及161个F₁代基因组DNA测序采用GBS方法(Elshire et al, 2011)进行, 由广州基迪奥生物科技有限公司完成。质检合格的基因组DNA被两种限制性内切酶EcorI和NIaIII (New England Biolabs, Ipswich, MA)消化, 然后进行末端修复, 加A-尾, 并使用NEBNext^® ΜLtra^TM DNA Library Prep Kit库(NEB, USA)添加Illumina测序接头。选取300~500 bp的DNA片段进行PCR扩增富集, 最后使用AMPure XP系统(Beckman Coulter, Brea, CA, USA)对PCR产物进行纯化, 使用Agilent 2100生物分析仪(Agilent, Santa Clara, CA)对测序文库进行检测, 并使用实时PCR进行文库定量。测序在Novaseq 6000测序仪上进行, 采用PE150测序策略。

1.3 SNP标记筛选与过滤

使用FASTP(0.18.0) (Chen et al, 2018)对测序后的raw data进行SNP过滤。筛选标准为: (1) 去除含有未知核苷酸(N)≥10%的reads; (2) 去除phred质量评分≤20及碱基≥50%的reads; (3) 删除含接头的reads。过滤后的clean reads用于组装分析。使用Burrows-Wheeler Aligner (BWA) (Li et al, 2009) (0.7.12; 比对参数为-k 32 -M)采用mem算法将过滤后的reads比对到参考基因组(SRA PRJNA642704); 比对后结果使用软件picard (http://sourceforge.net/projects/picard/.) (1.129)进行标记。使用变异检测软件GATK (Van Der Auwera et al, 2013) (3.4~46) (设置参数: -Window 4, -G_filter “QD < 2.0 || FS > 60.0 || MQ < 40.0”)进行群体SNP检测, SNP标记过滤标准为: 去除分型比例低于30%的位点, 去除杂合比例大于75%的位点, 按理论比对标记位点的基因型比例进行卡方检验, p值小于0.001的位点视为严重偏分离位点并去除, 保留分离类型为母本杂合型lmxll、父本杂合型nnxnp和双亲杂合型hkxhk的标记。使用ANNOVAR软件(Wang et al, 2010)进行功能注释SNP。

1.4 遗传图谱构建

利用筛选出来的SNP标记, 使用Joinmap4.1划分连锁群, 构建遗传图谱。对划分好的连锁群采用最大似然法进行排序, 对排序后的结果进行校正, 然后再使用回归算法进行排序, 使用Perl SVG模型获得整合图谱(Van Ooijen, 2011)。整合图谱构建中, 以多鳞基因组为参考, 根据标记的物理位置信息将其归并成对应的连锁群, 去除非染色体上的标记, 在每个遗传位置的标记中取完整度最高(缺失率最低)的标记作为作图标记。

1.5 生长相关性状QTL定位及候选基因预测

使用R/qtl (Broman et al, 2003)的scanone()函数对6个生长性状进行复合区间定位(CIM), 扫描步长为1 cM, 使用LOD=3.0作为阈值筛选QTL。利用多鳞参考基因组注释信息进行候选基因预测及功能注释。通过BLASTN将确定的显著性SNP位点比对到多鳞参考基因组, 并根据基因组注释结果筛选区间的基因, 并进行基因功能注释。

2 结果 2.1 生长相关性状的表型参数

本研究使用161尾多鳞F₁代全同胞个体进行连锁图谱构建及QTL定位分析。课题组前期对该家系的生长性状相关性和通径分析显示, 多鳞体高、体厚、眼径、体长和背鳍前长对多鳞的体质量有显著影响, 初步确定为多鳞的生长性状选育指标(余家旺等, 2022)。形态学分析显示, 6个主要形态性状的变异系数(CV)值在11.83%~33.51%, 表明性状变异丰富, 适合于进行QTL分析。6个表型性状的频率分布符合正态分布(图 1), 可用于后续的QTL初步定位分析。

2.2 多鳞作图家系GBS测序数据及SNP筛选

对多鳞父母本和161个F₁代进行GBS测序。两个亲本共获得3.40 G bp Raw data (测序深度: ~6.5×; 多鳞基因组为521.63 Mb), 过滤后亲本的平均Clean data为1.60 G bp, 161个子代样本共获得880.55 G bp Raw data (平均每个个体550 M bp), 过滤后平均每个子代个体的Clean data为534 M bp, 测序数据Q20≥98%, Q30%≥94%, 样本GC分布正常(平均43%)。过滤后共获得657.88 M Clean Reads, 平均每个个体reads长度为142 bp。将测序数据与参考基因组进行比对, 比对结果显示, 所有reads与参考多鳞基因组(SRA PRJNA642704)的平均比对率为98%以上, 用于后续的变异检测及相关分析。

群体变异检测共获得205 471个SNP位点, 这些位点中转换标记有120 302个, 占58.54%, 颠换标记有85 169个, 占41.46%。根据亲本的基因型确定标记的分离类型(亲本测序深度不低于4), 保留分离类型为母本杂合型lmxll、父本杂合型nnxnp和双亲杂合型hkxhk的标记, 共筛选获得143 886个多态性SNP位点(表 1); 随后, 对标记进一步过滤, 去除分型比例低于30%、或杂合率大于75%的位点、或严重偏分离的位点, 共保留107 406个高质量的SNP标记用于后续作图分析。

2.3 多鳞遗传连锁图谱构建

以筛选获得的107 406个SNP位点为基础, 根据标记的物理位置信息划分连锁群, 去除非染色体上的标记, 并在每个遗传位置的标记中保留完整度最高的标记, 最终获得4 735个SNP标记, 被分配到24个连锁群, 与多鳞单倍体染色体数目一致(表 2, 图 2)。该遗传连锁图谱全长2 154.803 cM, 平均遗传距离0.455 cM。

如表 2所示, 24个连锁群长度介于65.059 cM (LG23)与114.276 cM (LG08)之间, 平均长度89.783 cM, 各连锁群标记平均遗传距离介于0.383 cM (LG04)至0.554 cM (LG22)之间。其中, LG01连锁群分布有最多的SNP标记(232个), 其连锁群长度为112.766 cM, 标记平均距离为0.486 cM; 而LG23上分布有最少的SNP标记(131个), 其长度为65.059 cM, 平均遗传距离为0.497 cM。

2.4 生长相关性状QTL的定位

基于上述遗传图谱, 利用R/qtl软件对多鳞体重、体高、体厚、体长、眼径、和背鳍前长6个生长性状进行QTL复合区间定位分析, 共定位到20个QTLs (LOD > 3.0), 分布在LG04、LG05、LG07、LG08、LG16、LG17、LG18以及LG24 8个连锁群(图 3)。单个QTL的LOD峰值范围为3.02 (qGRT-15) ~ 4.23 (qGRT-13), 可解释的表型变异(PVE)范围为0.14% (qGRT-7) ~ 8.42% (qGRT-12), 平均为3.42% (表 3)。连锁群LG08上聚集了8个生长性状显著相关的QTL, 其中79.908~114.276 cM区间的6个显著QTL与体重、体高、体长以及背鳍前长4个生长性状显著相关。

2.5 候选基因的预测

本研究通过与多鳞参考基因组序列比对, 对QTL置信区间进行了基因查找, 共鉴定得到19个生长相关候选基因(表 3, 表 4)。相关基因包括组蛋白去甲基化酶(kdm5c)、跨膜蛋白9 (tmem9)、Tel2相互作用蛋白1 (tti1)、RNA结合基序蛋白10 (rbm10)、胰岛素样生长因子1 (igf1)、胰岛素样生长因子2 (igf2)、转化生长因子B受体2 (tgfbr2)、生长抑制特异性蛋白1 (gas1)、类胰岛素生长因子酸不稳定亚基(igfals)、生长抑素受体5 (sstr5)、生长抑素1a (sst1a)、成纤维细胞因子6 (gfg6)、成纤维细胞因子20 (gfg20)、哺乳动物雷帕霉素靶蛋白(mtor)、葡萄糖转运蛋白4 (slc2a4)、骨形态发生蛋白7 (bmp7)和肿瘤抑制蛋白p53 (p53)等基因。上述候选基因为研究多鳞生长调控提供了重要的基因资源, 其功能有待深入研究。

3 讨论

作图群体的选择是构建遗传连锁图谱的先决条件。不同的作图家系类型各有利弊, 主要包括重组近交系(recombined inbred lines, RIL)、单倍体(haploid, HAP)和双单倍体(doubled haploid, DH)、子一代(F₁)、子二代(F₂)和回交(back crossing, BC)等家系群体(Peng et al, 2016)。通常在鱼类中构建RIL、HAP和DH家系较为困难, 而F₂和BC家系的构建耗费时间长。基于养殖鱼类繁殖周期较长、基因组复杂、遗传多态性高、繁殖率高等特点, 目前大部分水产物种采用拟测交策略的F₁作图家族构建遗传连锁图谱。除了作图家系的选择对遗传连锁图谱质量产生影响外, 作图标记类型的选择也是直接影响遗传图谱质量的原因之一。与其他遗传标记相比, SNP标记是基因组中最丰富的标记类型, 具有高多态性。然而, 在具有相对较大作图群体的非模式水产物种中, 很难开发足够数量的SNP (Tsigenopoulos et al, 2014; Zhu et al, 2014)。随着高通量测序技术的发展, 各种基于基因组序列的SNP基因分型技术(包括RAD、2b-RAD、GBS和SLAF等)被开发出来, 为构建高密度连锁作图提供了快速便捷的手段。GBS-seq是在第二代测序基础上发展而来的一种基于全基因组酶切位点的基因分型测序技术(Elshire et al, 2011)。GBS-seq不仅实验操作简单, 性价比高, 更为重要的是它一次测序即可获得数以万计的多态性遗传标记, 已经广泛用于群体遗传学、基因组学和遗传连锁图谱构建等研究领域(Guo et al, 2021; Weng et al, 2021)。本研究基于拟测交策略, 对多鳞F₁代作图群体进行基因分型测序, 获得了大量分离比为1:1的SNP位点, 成功构建了该物种高密度连锁图谱。多鳞高密度连锁图谱全长2 154.803 cM, 分为24个连锁群, 包括4 735个SNP标记, 标记间平均遗传距离为0.455 cM。近几年公布的鱼类遗传图谱中, 鳜(Siniperca chuatsi)遗传图谱全长1 694.3 cM, 上图标记数2 328个, 标记间平均遗传距离为0.7 cM (Guo et al, 2021)。黄姑鱼(Nibea albiflora)图谱全长4 300.20 cM, 上图标记数6 219个, 标记间平均遗传距离为0.69 cM (Xu et al, 2021)。丝尾鳠(Hemibagrus wyckioides)图谱全长2 067.35 cM, 上图标记数2 369个, 标记间平均遗传距离为0.87 cM (Zhou et al, 2021)。与其他鱼类遗传图谱相比, 多鳞遗传图谱上图标记数量合适、标记密度高。此外, 当遗传连锁图谱构建时的上图标记数量和分布合适时, 遗传连锁图谱连锁群的数量应该等于该物种单倍体的染色体数目。多鳞染色体水平基因组分析结果表明, 多鳞单倍型染色体数目N=24 (Lin et al, 2021)。在本研究中获得的遗传连锁图谱连锁群数量也为24条, 与多鳞单倍体染色体数目一致, 这表明本研究构建的多鳞遗传连锁图谱质量高。综上, 本研究构建了高密度和高分辨率和的遗传图谱, 为后续QTL精细定位提供了重要工具。

QTL分析有助于发现性状相关的连锁标记和预测候选基因。现有研究表明, 鱼类生长性状多受微效多基因的调控。本研究中, 在多鳞的全同胞家系中定位到多个QTL位点且分布在不同的连锁群上, 未能定位到生长性状相关主效位点, 这表明多鳞生长相关性状可能受到微效多基因的调控, 跟已知的多数鱼类生长性状QTL定位结果相似。Zhang等(2018)在布氏鲳鲹(Trachinotus blochii)中成功检测到了23个与体长、体高及体重等生长相关性状的QTL。Zhou等(2021)通过QTL分析显示, 在10个连锁群上鉴定了丝尾鳠4个形态性状(体长、头长、体高和体重)的17个QTL, 其可解释的表型变异为7.4%~13.3%。Zhang等(2020)对瓦氏黄颡鱼(Pelteobagrus vachelli)体重相关的8个性状进行QTL分析, 鉴定到4个生长性状显著相关的QTL, 其可解释的表型变异分别为11.2%~12%。本研究获得的生长相关QTL的可解释的表型变异比值相比其他水产物种较低, 平均为3.42%, 可能与选取的作图家系类型(F₁全同胞家系)及其个体数目不足有关, 在后续研究中, 选用不同的作图群体类型(如F₂代群体)或增加作图群体数量或选用不同发育时期的个体可以直接和间接地改善可解释的表型变异, 从而提高检测QTL的准确性, 降低因假阳性而导致的错判。

位于连锁群L08上79.908~114.276 cM区间的6个显著QTL与体重、体高、体长以及背鳍前长4个生长性状显著相关, 定位到的生长相关QTL具有成簇分布的特点, 这表明连锁群LG08的上述区段可能与多鳞的生长发育过程显著相关。鱼类生长相关QTL往往呈现成簇分布的特点, 如Wang等(2018)在黄河鲤(Cyprinus carpio haematoperus)中发现, 定位到生长相关QTL呈现出成簇分布的特点; Yang等(2020)在金钱鱼(Scatophagus argus)中研究发现其体长和全长两个生长相关性状的大多数QTL均分布在连锁群LG2的重叠置信区间内。现有研究表明, 鱼类生长性状多受微效多基因的调控, 生长相关QTL在同条染色体上的成簇分布, 区域内的功能基因对生长相关性状的影响可能存在一因多效或不同的功能基因间互作影响。对多鳞生长性状相关QTL的深入研究, 可明确各区域内基因位点的相互作用和生长性状间的调控机制, 有助于指导多鳞的选育工作。

现有研究结果表明, 鱼类生长与发育受多基因调控, 已有物种的研究成果尚不足以解释鱼类生长调控机理, 有必要在更多的物种中开展鱼类生长调控相关基因的发掘研究。本研究中, 在显著QTL区域内鉴定出22个生长相关候选基因, 其中包括胰岛素样生长因子1 (igf1)、胰岛素样生长因子2 (igf2)、类胰岛素生长因子酸不稳定亚基(igfals)、生长抑制特异性蛋白1 (gas1)、生长抑素1a (sst1a)、生长抑素受体5 (sstr5)、成纤维细胞因子6 (gfg6)、转化生长因子B受体2 (tgfbr2)、成纤维细胞因子20 (gfg20)以及骨形态发生蛋白7 (bmp7)等已知参与调控鱼类生长发育的功能基因。此外, 在连锁群L08 79.908~114.276 cM区间注释到组蛋白去甲基化酶(kdm5c)、Tel2相互作用蛋白1 (tti1)及RNA结合基序蛋白10 (rbm10)等候选基因。KDM5家族蛋白是激活大量促增殖细胞周期基因所必需的, 在神经发生、DNA损伤反应和癌症发展等方面发挥重要作用(Shen et al, 2021)。tti1基因参与促进哺乳动物组装、稳定和维持mTORC1和mTORC2复合物的活性, 这些复合物调节细胞生长和存活, 以响应营养和激素信号(Kaizuka et al, 2010)。rbm10基因参与调节前体mRNA的选择性剪接和mRNA稳定, 该基因的过度表达已被证明可抑制肺腺癌的恶性增殖(Li et al, 2020)。上述候选基因为研究多鳞生长调控提供了重要的基因资源, 为该物种的为分子遗传育种研究提供有价值的参考信息。

4 结论

本研究基于多鳞高密度遗传连锁图谱进行生长性状的QTL初步分析及候选基因的预测, 获得20个生长相关性状的显著QTL和多个生长调控相关的候选基因。多鳞的QTL初步定位及候选基因的初步鉴定将为该物种遗传育种研究提供有用信息和奠定基础, 而在之后的F₂群体中可利用不同养殖环境的多个家系进行全面而精确的QTL分析, 以获取更稳定、精确的生长相关性状QTL。