利用山兔和褐兔的遗传多态性来调查野野兔的本地适应和保护状况(sylvaticus，尼尔森1831)

在瑞典有两个形态上截然不同的山兔亚种，名义上的亚种(Lepus timidus timidus林奈，1758)，广泛分布于北极/亚北极欧亚大陆，以及地方性的野野兔(sylvaticus尼尔森,1831)¹这是一个更适合温带环境的亚种，分布在斯堪的纳维亚半岛南部和爱沙尼亚西部²．最明显的区别是它们冬季皮毛的颜色，名义上的山兔亚种是白色的，而野野兔是灰色的。

名义亚种的地理分布在瑞典北部已逐渐消退。其中一个可能的原因是来自引进的棕兔的竞争(天兔座europaeus帕拉斯，1778年)^th世纪^3.但它们的撤退也可能是由气候变化引起的伪装不匹配或其他栖息地的变化引起的^4，5，6．虽然野野兔在无雪地区伪装得更好，但正如狩猎收获数据持续下降所表明的那样，由于栖息地改变和来自棕野兔的竞争，它也可能面临局部灭绝⁷．

早期的研究也表明了两个物种之间的种间基因流动(即渐渗)，天兔座timidus而且天兔座europaeus⁸．渗透可以促进非本地棕色野兔适应北方的先决条件(例如季节性转移到白色冬季毛皮^3.，7寒冷的气候和/或病原体抗性⁹)．因此，这个两种物种的综合体为探索人类引起的生态变化如何快速影响物种，包括它们的行为、人口和进化相互作用，提供了一个极好的模型系统。除了上述的保护挑战之外，在瑞典，山兔和野野兔都是文化上重要的狩猎物种，数量的下降引起了人们的关注。

这里报告的基因组数据可用于解决复杂的适应模式，评估物种相互作用，并作为开发受威胁种群监测工具的资源。它们还为亚种差异和进化历史提供了有价值的见解。我们描述了全基因组重测序数据中的11,363,883个多态性，以检测差异标记，并提供了一个可能影响冬季毛皮颜色的色素沉着基因的示例分析筛选¹⁰．除了评估两种野兔的当地适应性、健康野兔的保护状况和开发管理工具外，所报告的基因组数据还可用于增加对整个分布的其他接触区山区野兔和棕兔之间相互作用性质的理解^{3.，7，11，12，13}．结合收获数据、颗粒库存和相机捕获，这些具有潜在添加功能的基因数据可用于开发针对山地兔和褐兔的综合监测工具。

样本采集和DNA提取

肌肉组织样本来自2005年至2017年瑞典常规狩猎季节拍摄的7只兔子。在瑞典南部，山地野兔的常规狩猎季节从9月1日持续到2月15日，在瑞典北部，狩猎季节一直持续到2月的最后一天。在瑞典，褐色野兔的常规狩猎季节是从9月1日到2月的最后一天。样本代表两种野兔和一种野兔亚种(表1)在瑞典(图。1):两个人天兔座europaeus,两个Lepus timidus timidus还有三个sylvaticus．所有三个sylvaticus(Le911, Le918和Le919)样本采集于1月和2月，此时蓝色/浅灰色冬季毛皮确认了亚种状态(表9)1；Angerbjörn & Flux 1995)。此外，样品从一个L. t.蒂米达斯(Le921)与l . europaeus在一月份采集，此时该亚种的白色冬季毛皮便于鉴定。最后三个样本的样本，一个L. t.蒂米达斯(Le922)和两个l . europaeus(LeA01和LeV01)，分别采集于不同物种分布的地区(详情见表1)．模板DNA按照制造商的说明使用QIAsymphony DNA Mini Kit提取。

文库制备和测序

这7个样品在Illumina HiSeqX的单个流池上的4个通道上测序。在瑞典斯德哥尔摩(Stockholm, Sweden)的sciilifelab (National Genomics Infrastructure Stockholm, Sweden)为每个兔样本准备了一个Illumina TruSeq PCR免费文库(Illumina, USA)。库的目标插入大小为350 bp。然后将样品以等摩尔比混合，并在Illumina Hiseq X平台上使用2 × 151 Hiseq X SBS化学试剂在4个通道上测序。使用Illumina软件包CASAVA的bcl2fastq_v2.19.1.403执行Bcl到FastQ的转换。

生物信息学

原始测序读数的质量¹⁴采用FastQC 0.11.9¹⁵(无花果。1 c)．然后使用cutadapt版本2.10从原始测序读取中删除适配器序列¹⁶接下来是进一步的阅读质量检查。然后将修剪后的读取映射到伪参考Lepus timidus基因组¹⁷(NCBI接入:GCA_009760805.1)与Burrows Wheeler Aligner, BWA版本0.7.17¹⁸的默认设置bwa mem算法。每个个体在不同的测序车道上进行测序，因此每个个体的所有文件都被用作Picard版本2.21.4 markduplduplicate的输入¹⁹步骤，每个个体得到一个单独的bam文件，去除PCR和光学重复。使用freebayes版本1.3.1在100,000 bp区域的基因组中调用变体^20.．VCFtools版本为0.1.16²¹用于对原始变量文件进行四轮变量过滤。首先，我们允许最大缺失基因型率为50%，最小等位基因数为3，最小变异phred质量评分为30。其次，我们应用了一个较小的等位基因频率截断值0.05和最小平均深度20。第三，我们筛选等位基因平衡(AB)的变体:AB > 0.25 & AB < 0.75 | AB < 0.01。最后，应用了一个允许最大平均深度为30的硬过滤器。这种过滤产生了具有高质量调用的变体文件。命令和工具的详细列表可以在表中找到2．最终VCF文件中的基因型分期是使用beagle 4.1版本实现的²²．

最终VCF文件的详细信息²³根据样本种或亚种对VCF文件进行划分，计算共享和唯一snp的数量。然后使用vcffixup 1.0.0版本重新计算变量统计数据，并过滤以只保留变量位点(AC > 2和AF < 1)。²³．最后，将私有VCF文件转换为床文件，并利用干涉算法计算床文件的交点²⁴．R版本4.0.4估计种群分化(2021-02-15)²⁵使用glpca函数包含在包adegenet 2.1.5版中^26，27．

为了生成与其他野兔物种的线粒体系统发育树，我们组装了每个个体的线粒体基因组。使用NovoPlasty 4.3版本组装线粒体基因组²⁸．NovoPlasty需要一个种子序列来构建组装，因此我们使用了整个线粒体基因组天兔座timidus(Genbank登录NC024040)作为种子序列。然后将组装好的线粒体基因组与小鼠的基因组进行比对Oryctolagus cuniculus(注册编号AJ001588.1)，l . capensis(注册编号GU937113.1)，l .也(注册编号NC_024043.1和KJ397613.1)，l . townsendii(注册编号NC_024041.1和KJ397609.1)，l . coreanus(注册编号KF040450.1)，l . granatensis(注册编号NC_024042.1, KJ397610.1)，l . tolai(注册编号KM609214.1)，l . yarkandensis(注册编号MN450151.1和MG279351.1)，l . sinensis(注册号KM362831.1)，l . hainanus(注册编号JQ219662.1)，l . capensis(注册编号GU937113.1)，l分析(注册编号LC073697.1)，l . europaeus(登录号AJ421471.1和KY211025.1)l . timidus(登录号NC_024040.1, KJ397605.1, KR030072.1, KR030070.1, KR030069.1, KR013248.1和KR019013.1)使用MAFFT v 7.490²⁹．然后使用智商树web服务构建系统发育树^30.，由ModelFinder确定最佳模型³¹并使用2000次超快速自举迭代。第二种系统发育树是利用已确定的全基因组变异生成的。系统发育树估计为RaXML³²版本8.2.12，带有GTR + Γ + ASC模型和200个快速自举复制。然后在Ape 5.4版本中使用chronos功能对整个基因组ML树进行超测量³³的估计散度时间天兔座europaeus而且天兔座timidus．

验证

使用FastQC对测序质量进行技术验证，然后使用MultiQC编译数据³⁴．为了验证数据集的来源，我们使用了两种方法，首先确认冬季表型，其次确认系统发育关系及其祖先系数。系统发育关系是根据上述包括其他野兔物种在内的线粒体基因组进行评估的，然后使用来自封闭分析的全基因组变异。使用Poppr版本2.9.3确定祖先系数³⁵基于DAPC方法。

使用的例子

为了进一步说明数据集的潜在用途，我们对59个已知色素沉着基因中的snp进行了示例分析¹⁰，这可以用来了解皮毛颜色的遗传学(图。3.)．首先，下载表中每个基因的兔或小鼠参考序列2Hoekstra，(2006)从NCBI数据库。然后这些基因被映射到伪参考Lepus timidus基因组使用minimap2³⁶拼接感知转录算法。然后，我们计算了基于BEDtools映射所识别的皮毛颜色基因区域重叠的snp的数量³⁷并给出该数据(图;3.)．

测序数据

测序是成功的，结果总共有32.4亿次读取。测序阅读质量较高，平均phred值大于Q30(图3)。1 c)．对象的平均映射速率假参考Lepus timidus93.59%，平均79.9%的reads映射为合适的配对(表2)3.)．原始测序数据¹⁴可以从NCBI短读档案下的项目登录和个人登录号码可以在表1．

变量调用和过滤

变体调用是在CSC Puhti服务器上使用freebayes软件的并行版本执行的2(CSC -科技中心，芬兰)。Freebayes是一个单倍型敏感的变体调用者，在整个基因组中一共识别了78,968,182个变体。应用各种过滤器后(表2)，我们回收了11,363,883个高质量变异(表4)分布在整个假参考基因组中。

线粒体基因组组装

我们为每个个体组装了接近完整的线粒体基因组序列。线粒体基因组组装长度从16,648 bp到17,612 bps，有4个个体产生环状contigs (LeA01, LE911, LE919和LE921)。所有线粒体基因都存在于每个个体中。在本研究中产生的线粒体基因组的对齐包括作为phylip格式的Dryad数据集的一部分²³．

主成分分析(图;1 b)和complot(图;2摄氏度)，进一步证实了不同的物种分组。类似的分组模式已在天兔座在芬兰各地取样的物种³⁸的遗传多样性l . europaeus较小，可能是由于方正效应导致这些样本的紧密聚类，而更广泛的聚类l . timidus样品(图。1 b)．种间渗透也可能影响观察到的模式³⁸．在11,363,883个高质量的变种中，1,23,035个是独一无二的l . europaeus596,291到l . timidus(无花果。1 c)．

为了确认野兔样本之间的关系，我们基于全基因组snp生成了最大似然(ML)系统进化树(图2)。2)．这些样本的系统发育聚类证实了样本的起源，因为第一个分支分离了天兔座europaeus而且天兔座timidus，其次是分离的l . timidus亚种L. timidus sylvaticus．该亚种在来自大陆的样本和瑞典西南部岛屿上收集的样本之间产生了进一步的分裂。同样，基于线粒体基因组的系统发育(图。2 b)，首先显示l . europaeus而且l . timidus样本。然后L. timidus而且L. timidus sylvaticus样品被进一步分割。有趣的是，还有其他的l . timidus来自远东的基因组L. timidus sylvaticus与这些样本聚在一起。观察到的全基因组和线粒体基因组数据之间的差异可能受到渐渗模式以及不同的遗传方式的影响。

所描述的基因组数据可用于评估野兔在局部、区域或大范围分布中的适应性，例如色素沉着、耐寒性和病原体抗性。一个特别的价值是开发一个综合监测工具，以评估瑞典南部和中部的野野兔的保护状况。这可以通过整合基因数据、收获数据、颗粒库存和摄像机捕获来实现。除了一般的保护之外，从资源的角度来看，长期保护野野兔，以及理解山兔和棕兔之间相互作用的本质是很重要的，因为野兔在瑞典和其他地方是一种狩猎物种，因此是肉类和娱乐的来源。

尽管样本集很小，数据表明在两个山兔亚种中有显著数量的私有等位基因(图2)。1 d)．由此导致的两个演化支之间的分歧(图。2 a, B)在系统发育和祖先合集分析中清晰可见(图。2摄氏度)，支持亚种代表生物实体，而不仅仅是单一基因的冬季毛皮颜色差异³⁹．此外，这些私有等位基因可用于通过SNP面板检测野兔及其亚种的多样化。亚种在无花果中分裂。2可能从最后一个冰期末期开始就一直保持着同域关系⁴⁰尽管时间还需要进一步验证⁴⁰．维持遗传分裂的机制值得进一步研究，并提出了一个关于遗传分裂的有趣问题善意的野野兔的亚种状况。

为了说明如何将这些基因组数据与未来的分析相结合，我们计算了在已知色素沉着基因中发现的变异的数量。虽然，这项分析的样本量本身太小，无法提供任何结论，但我们能够识别出基因内部和周围的一些变异(图2)。3.)，对色素沉着有已知影响^10，39．如果与更大的数据集集成，我们相信这里提供的基因组数据将为种群遗传学、保护生物学和进化研究提供有价值的资源。