基于低深度全基因组测序分析内江猪群体结构和遗传多样性

内江猪是我国著名的优良地方猪种，2000年被列入《国家畜禽品种资源保护名录》，原产于四川盆地中部的内江、资阳一带，属于西南型猪种。内江猪体型较大，全身被黑毛，体躯宽而深，背腰宽且长，腹部较大略微下垂，臀部较宽略微后倾，头短额宽，额上有较深的横向褶皱，有“狮子头”和“二方头”之分，四肢坚实，皮较厚[1]。内江猪肉质优良，性成熟早、抗逆性强、耐粗饲，但瘦肉率低、屠宰率不高且背膘较厚[2]。内江猪拥有极佳的遗传配合力，在各类二元杂交试验中，杂交后代均显示出了明显的杂种优势，因而曾被国内各地区大量引种，一度成为国内引种最多的地方品种[3]。随着市场经济发展和外种猪迅速占领国内市场，内江猪群体迅速减少，内江猪保种成为了重点工作，近年来受非洲猪瘟的影响，内江猪保种工作面临严峻挑战[4]。

SNP(single-nucleotide polymorphism)即单核苷酸多态性，具有数量大、分布广、覆盖密度大、遗传稳定性高等特点。通过高密度SNP芯片和重测序方法可以直接获得全基因组水平的SNP，在猪[5-7]、牛[8-11]、羊[12-15]、鸡[16-18]等各类畜禽研究和全基因组育种中被广泛应用。随着高通量测序技术的发展，测序成本急剧下降，但高深度全基因组测序相比SNP芯片仍然较为昂贵。低深度全基因组重测序结合基因填充技术，能够获得覆盖全基因组的高质量SNP标记[19-20]，同时通过降低测序深度有效降低测序成本，已被广泛用于各类研究[21-23]。通过全基因组模拟数据或SNP芯片全基因组填充数据评估基因组选育效果，研究发现全基因组数据可以进一步提高选择准确性，基于基因组填充或全基因组重测序能够将选择准确性提高30%~40%[24-25]。

本研究利用1×低深度重测序技术获取全基因组SNP，通过不同的SNP质控策略获取了两组数据，全面分析了133头内江猪群体结构和遗传多样性，为内江猪遗传资源保护提供依据。

本试验选择内江市国家内江猪保种备份场的内江猪为研究对象，共采集各家系具有代表性的133头种猪样品，其中公猪12头，母猪121头。耳组织样品采集后置于含无水乙醇的离心管中，迅速放入-20 ℃冰箱保存备用。

1.2.1 DNA提取及测序 　　使用CTAB法提取133份样品的基因组DNA。在1%琼脂糖凝胶上检测DNA降解和污染程度，采用分光光度计检测DNA纯度，采用QubitⓇ2.0荧光计(Life Technologies, CA, USA)中的QubitⓇDNA测定试剂盒测定DNA浓度。

获得高质量的基因组DNA后，使用制造商建议的Truseq Nano DNA HT试剂盒构建测序文库，并添加barcode以区分不同样品。利用Illumina NovaSeq 6000平台对上述文库进行测序，测序模式PE 150 bp，插入片段大小约为350 bp，测序深度为1×。

1.2.2 测序数据处理及比对 　　测序原始数据需要进行严格的质控以保证后续分析的准确性，质控包括：1)去除含有≥10%未知碱基(N)的序列；2)去除含50%以上碱基质量值小于5的序列；3)去除＞10个碱基比对到接头的序列，允许≤10%的不匹配；4)去除文库构建过程中PCR扩增的重复序列。随后通过BWA(Burrows-Wheeler Aligner)将每个样本的clean read与参考基因组进行比对，采用GATK默认参数进行SNP检测。

1.2.3 SNP数据质控 　　将获取的SNP数据通过Plink(v1.90)软件进行质控。质控筛选包括：1)按照SNP检出率大于90%和95%获得两组SNP数据，分别编号为NJ90和NJ95；2)分别剔除两组数据中次等位基因频率小于5%的SNP；3)根据哈代-温伯格定律，分别剔除两组数据中不符合哈温平衡的偏离位点(P < 1×10-6)；4)去除性染色体SNP。后续分析使用质控后的NJ90和NJ95两组SNP数据。

1.2.4 遗传多样性分析 　　杂合度是反映群体遗传多样性的重要指标，包括观测杂合度和期望杂合度。观测杂合度(observed heterozygosity, Ho)是指群体中杂合子所占比例，期望杂合度(expected heterozygsity, He)是指群体内随机抽取两个等位基因各不相同的概率，其计算公式为He=$\sum_{i=1}^n$Pi2，其中Pi为等位基因频率，n为等位基因总数。有效等位基因数(effective number of alleles, Ae)是指理想群体中，一个位点上产生与实际群体杂合度相同时需要的等位基因数，其计算公式为Ae=$\frac{1}{\sum_{i=1}^n P i^2}$，其中Pi为等位基因频率，n为等位基因总数。多态性标记比(proportion of polymorphic marker, PN)是指具有多态性标记的位点比例，其计算公式为$P N=\frac{N}{M}$，其中N指的是有多态性表现的位点数目，M为总的位点数。多态性标记比与有效等位基因含量也是评估遗传多样性的常用指标。遗传多样性分析通过Plink(v1.90)软件[26]完成。

1.2.5 群体结构分析 　　主成分分析(principal component analysis, PCA)常用于高纬数据降维，使具有相同数据特征的数据聚集在一起，简化复杂变异。本研究使用Plink(v1.90)软件利用SNP标记构建亲缘关系G矩阵计算个体间的亲缘系数，并计算前2个主成分，结果通过Rstudio进行可视化。

随后利用Plink(v1.90)软件进行连锁不平衡分析，以50个SNPs为一个窗口，每次步移5个SNPs，r2阈值为0.6，去除强连锁位点。再基于状态同源(idendity of state, IBS)计算个体间的遗传距离(D)，计算公式为D=1-DST，DST=$\frac{0.5 \times I B S 1+I B S 2}{m}$，其中m为SNP标记位点总数，IBS1为两个个体所有位点中只有一个观测值相同的等位基因数量，IBS2为两个个体所有位点中两个观测值均相同的等位基因数量。随后使用MEGA X软件利用D值通过邻接法(Neighboring-Joining Tree, NJ-tree)构建进化树。其中公猪依据个体间遗传距离，以三代以内无亲缘关系，即亲缘系数小于6.25%为依据划分家系。

1.2.6 群体近交系数 　　基因组近交系数计算常用方法包括基于SNP纯合性计算的Fhom和基于连续性纯合片段(runs of homozygosity，ROH)计算的FROH。Fhom可利用Plink(v1.90)软件计算，其公式为Fhom$\frac{O(\text { hom })-E(\text { hom })}{N-E(\text { hom })}$, 其中O(hom)表示观测纯合子，E(hom)表示期望纯合子，N表示位点总数。连续性纯合片段是指基因组某一段区域内，一定数量一定密度的SNPs表现为纯合时，可以判定该区域存在连续性纯合片段。利用ROH计算的近交系数为FROH，其公式为: FROH=$\frac{\sum L_{R O H}}{L_{{auto }}}$，其中LROH为常染色体上ROH片段的总长度，Lauto为常染色体总长度(本研究参考猪V10.2版本基因组，常染色体长度约为2 450 462.292 kb)。ROH可以通过Plink(v1.90)软件，使用滑动窗口的方法对常染色体进行检测，具体检测参数为：滑动窗口大小为50个SNPs；每个窗口最多含有1个杂合子；最多含有1次缺失；ROH至少包含连续50个SNPs；长度大于1 000 kb；连续SNP间的距离小于1 000 kb。

测序数据经过整体质控后，按照SNP检出率高于90%和95%分为NJ90和NJ95两组数据用于后续分析。遗传多样性评估结果如表 1所示。NJ90中133个个体包含135 760个SNPs标记，其中等位基因频率为0.87，有效等位基因数为1.27，多态性标记比为0.76，观测杂合度和期望杂合度分别为0.15和0.21，整体来看内江猪该群体拥有较为丰富的遗传多样性。NJ95中133个个体包含32 266个SNPs标记，等位基因频率为0.79，有效等位基因数为1.44，多态性标记比为0.74，观测杂合度和期望杂合度分别为0.30和0.31。相较两组数据结果，NJ95的SNPs位点与NJ90完全重合且NJ90包含更多的SNPs位点，各类指标计算结果有差异。NJ90等位基因频率高于NJ95，有效等位基因数和杂合度低于NJ95，且杂合度差异较大。

对内江猪133个个体进行主成分分析，结果如图 1所示。NJ90主成分分析结果(图 1A)表明, 群体内个体较为分散，主成分1能解释38.65%的变异，主成分2能解释32.31%的变异，且被分为较远的3个群。NJ95主成分分析结果(图 1B)同样说明, 群体内个体较为分散，主成分1能解释36.00%的变异，主成分2能解释32.57%的变异，群体内分层不明显。

基于IBS遗传距离和亲缘关系矩阵分析群体结构，NJ90结果显示内江猪群体内遗传距离平均值为0.20，最大值和最小值分别为0.21和0.15，亲缘系数平均为0.9%(图 2A)。NJ95结果显示群体遗传距离平均值为0.25，最大值和最小值分别为0.27和0.20，亲缘系数平均为0.7%(图 2B)。NJ90和NJ95结果比较，群体平均亲缘系数差异不大，NJ90评估遗传距离较NJ95更小，NJ90和NJ95的群体分析结果表明内江猪群体大部分个体间亲缘关系较远。

由于公猪在保种群中具有重要地位，因而采用邻接法首先对12头公猪进行聚类分析从而划分家系。NJ90结果显示，结合亲缘系数，12头公猪可以分为6个家系(图 3A)。NJ95聚类结果与NJ90差异较大(图 3B)，结合亲缘系数可以分为5个家系，且每个家系中的个体与NJ90不同。

随后基于遗传距离构建群体进化树。NJ90结果显示根部将133个个体分为3个主要分枝，结果与PCA分析一致(图 4A)。NJ95结果有所不同，根部同样将个体分为3枝，但其中有两枝个体较少，个体聚类结果与NJ90相比差异较大(图 4B)。例如公猪129在NJ90中与母猪94、112、106、30、29、16聚为一枝，而在NJ95中却与112、106、119、10、6聚为一枝，只有112、106、119三个个体一致。聚类结果表明，NJ90的公猪聚类结果和群体聚类结果一致性更高，对群体分群和保种工作更具有参考价值。

根据NJ90公猪聚类结果和亲缘关系，可以将内江猪群体划分为6个家系。家系1包含公猪1头，母猪13头；家系2包含公猪3头，母猪15头；家系3包含公猪1头，母猪6头；家系4包含公猪3头，母猪34头；家系5包含公猪3头，母猪34头；家系6包含公猪1头，母猪4头。剩余15头母猪与公猪亲缘关系较远, 因此可以单独分类(表 2)。

运用全基因组方法利用Plink(v1.90)软件计算群体近交系数，NJ90结果显示Fhom平均近交系数为0.27，其中最大Fhom为0.67，最小Fhom为0.05。NJ95结果显示Fhom平均近交系数为0.01，其中最大Fhom为0.55，最小Fhom为-0.23。两组数据近交系数评估结果差异较大，但结果显示该群体存在一定程度的近交。

通过Plink(v1.90)软件对全基因组连续性纯合片段(ROH)进行检测，NJ90共发现11 062个ROHs，平均ROH长度为1 959.29 kb，长度范围在1 000.05~12 498.5 kb之间，主要集中在1~2.5 Mb之间，1~1.5 Mb之间占比39.44%，1.5~2.5 Mb之间占比40.39%(图 5A)。从染色体上分布来看，1号染色体上ROH数量最多(1 207个)，12号染色体最少(261个)，染色体2、3、4、5、6、7、8、9、13、14、15上较多且长度差异较小(图 5B)。基于该ROH结果评估群体近交系数FROH为6.65%。

NJ95采用同样的参数检测ROH，共发现52个ROHs，平均ROH长度为1 382.76 kb，长度范围在1 005.04~3 104.19 kb之间，主要集中在1~1.5 Mb之间占比76.92%(图 5C)。从染色体分布上来看，2号染色图上ROH最多(17个)，1、6、8、10、11、12、16、17、18号染色体上均未检测到ROH，同时部分个体也未检测到ROH(图 5D)。基于该ROH结果评估群体近交系数FROH为0.02%。

内江猪作为西南特色的地方猪种，以极佳的遗传配合力闻名，是非常宝贵的遗传资源。同其他许多地方品种一样，自20世纪80年代开始，内江猪在政府支持下开始了研究、开发和选育工作，其规模从3万余头一度发展到10万头，生长性能、繁殖性能、肉品质等多种性状研究也取得了一定成果[2]。然而在多种因素影响下，如今的内江猪保种工作仍然面临诸多挑战。一方面，内江猪经济性状相比外种猪不具竞争力，加之饲养环境复杂，使得内江猪种群规模得不到有效保障；另一方面，对于内江猪的开发和利用进展缓慢，虽然具有一定的市场前景但目前仍不能有效发挥其种质优势。本研究利用全基因组SNPs信息，对内江猪保种群体情况做了较为全面的分析，为内江猪保种工作开展提供了依据。

基于全基因组SNPs数据分析群体遗传结构和遗传多样性，是调查评估地方品种群体情况的有效手段。青峪猪[27]、淳安花猪[28]、通城猪[29]、安庆六白猪[30]、里岔黑猪[31]、枣庄黑盖猪[32]等地方品种群体遗传研究均使用商业化SNP芯片，而低深度重测序技术在牛的群体结构和遗传多样性分析中报道较多[9-10]。高深度测序可以通过位点的多层覆盖大大提高基因分型的准确性，并检测到真实的SNP变异，而低深度测序主要依靠基因型填充技术来实现。目前，研究者已提出多种基因型填充的方法，都可以对低覆盖率基因组测序数据进行填充并获得较好的准确性[33-34]。同时研究发现，使用低深度测序数据填充后进行育种值估计相较SNP芯片能够有效提高准确性，其准确度也与遗传力高低相关[11]。因此，低深度重测序相比SNP芯片能更加真实全面的反映基因组SNP变异，可以更准确的进行群体结构和遗传多样性分析。

SNP质控往往首先需要根据检出率剔除部分检出率低的SNP和个体。本研究SNP位点质控中，对SNP检出率按照90%和95%分别过滤，得到NJ90和NJ95两组数据，NJ90总位点数是NJ95的4倍，包含了NJ95全部位点，可以发现低深度重测序获取的SNP检出率在90%~95%之间的位点数量非常多，而这部分数据在NJ95的分析中并未有效利用，造成群体结构和遗传多样性分析结果相比NJ90有较大差异，尤其对群体遗传距离、亲缘关系和近交系数的估计影响较大。此外，商业化SNP芯片检测位点经过筛选后数量有限，常用猪SNP芯片获取SNP数据与NJ95位点总数相当，因此本研究通过数据分组处理可以模拟低深度重测序和SNP芯片两种技术获取的SNP数据，比较其在群体遗传结构和遗传多样性分析中的差别。本研究由于样本量较少，未对个体检出率进行过滤，经统计发现NJ90个体SNP平均缺失率为2.4%，最大缺失率为10.3%；NJ95个体SNP平均缺失率为8.4%，最大缺失率为23%。显然，NJ95个体SNP缺失率更大，这可能也是后续各项分析结果差异较大的原因之一。

本研究分别利用NJ90和NJ95两组数据分析了内江猪遗传多样性，其中等位基因频率、有效等位基因数和多态性标记比结果差别不大，但杂合度估计中NJ90的观测杂合度和期望杂合度都远低于NJ95，这可能是因为NJ90包含更多的SNPs位点，其中杂合位点的比例相对来说低于NJ95。NJ95的结果与其他通过SNP芯片分析地方猪种遗传多样性的结果更加相似，但两组结果都显示出了内江猪较为丰富的遗传多样性。在亲缘关系分析中，NJ90和NJ95两组数据结果差异不大，都说明了群体内大部分个体亲缘关系较远。NJ90和NJ95的主成分分析和遗传距离聚类分析结果分别相互印证，但两相比较结果差异较大，尤其是聚类分析。在公猪遗传距离聚类分析中，个体129、125、133、124虽然都被聚到一起，但NJ90能够更准确的评估其中的差异，126、122、132聚类结果也不相同，相较而言NJ90可能更能反映基因组的真实情况，更具有参考价值。

近交系数用于衡量个体的近交水平，表示个体随机选择两个等位基因是血缘同源的概率。在有限的群体规模中，近交水平增加会导致有害纯合子增加，从而产生近交衰退，后代生产性能降低。传统近交系数计算基于完整准确的系谱资料，但实际生产中很难保证系谱资料完整、准确，地方猪保种工作中情况可能更加严峻。基于基因组数据计算群体近交系数能有效规避系谱资料缺失错误的问题，同时也能较为准确的反映基因组的真实情况。常用方法主要包括基于SNP纯合位点的方法和基于连续性纯合片段的方法。基于SNP纯合位点的方法计算结果可能出现负值[35]，其结果区间为-∞~1，如NJ95结果显示最小近交系数为-0.23，负值的出现可能是因为当前研究样本群体具有密切的相关性或较高水平的近交，则远交亲本后代的近交系数将远低于参考群体，出现负值[36]。连续性纯合片段(ROH)是亲代将来自共同祖先的同一段单倍型传递给后代形成，基于ROH估算近交系数也被广泛应用[37-38]，其结果一般在0~1区间内，且结果跟ROH的长度具有一定的相关性[39-40]。本研究针对两组数据同时使用了上述两种方法，结果差异较大。NJ90和NJ95计算Fhom平均值分别为0.27和0.01，结果说明群体具有一定程度的近交。而FROH分别为0.066 5和0.000 2，该结果又显示群体近交程度较低。两组数据不同结果可能是由于质控后SNP数量差异较大，覆盖基因组程度不同，导致相同方法检测的ROH数量和长度不同，从而导结果差异。同时也可能受基因型填充的影响而不能反映基因组真实情况[38]。有研究人员提出，在有效规模较小、群体关系更紧密、选择相对不太强烈的群体中，基于SNP纯合估计近交系数出现明显的高估，基于ROH估计近交系数更加准确[41]。因而本研究中基于NJ90，通过ROH评估近交系数更具参考价值。

本研究基于低深度全基因组测序，分析了133头内江猪的遗传多样性和群体结构，结果表明内江猪群体具有较为丰富的遗传多样性。同时根据公猪聚类分析和亲缘系数可以将该群体划分为6个家系。该群体整体亲缘关系较远，近交程度较低，但仍需要严格规划选配，控制群体近交。同时低深度重测序能够更大范围覆盖基因组信息，对群体遗传多样性分析和遗传结构分析更具参考价值。