简化基因组测序在群体遗传进化分析中的应用

研究背景

简化基因组测序技术通过利用限制性内切酶对约1-10%的基因组进行测序分型，能够大幅降低基因组的复杂度，具有适用于大量样本低成本基因分型，可不依赖于参考基因组分析等优势 (Campbell et al., 2018)。近些年来，以RAD (restriction site-associated DNA sequencing) 为代表的简化基因组测序技术已成为广泛应用于动植物遗传领域基因组分型的重要手段。2b-RAD (type IIB restriction endonucleases site–associated DNA) 技术是唯一可对等长标签 (约30-35 bp) 进行测序分型的RAD类技术 (Wang et al., 2012)，该技术具有可对所有酶切片段分型，可通过选择性碱基调节标签密度，建库流程简便快速，可处理高度降解或者痕量样品等优势。该技术已经被广泛应用于高精度遗传图谱绘制、重要性状遗传解析、群体遗传学分析以及全基因组选择等研究方向。Multi-isoRAD (Multi-isolength restriction site–associated DNA) 是在2b-RAD技术基础上发展起来的串联标签测序技术 (Wang et al., 2016)，在一个测序片段 (read) 实现了对多达5个标签同时分析，并首次实现了全基因组SNP分型和DNA甲基化的同步联合分析。该技术解决了2b-RAD技术无法应用于双末端测序平台的局限，使得简并基因组分析成本大大降低 (单标签测序成本仅为技术改进前的十分之一)。Multi-isoRAD技术能够使研究者以较低的成本获得基因组大量的变异信息，为开展群体多样性和遗传结构分析以及群体进化和适应性机制解析提供了重要的技术手段。

材料与方法

1. 0.2 ml PCR管 (Axygen, catalog number: PCR-02-L-C)
2. 1.5 ml EP管 (Axygen, catalog number: MCT-150-C)
3. 各种型号的枪头 (10 μl, 100 μl, 200 μl, 1000 μl, Axygen)
4. BsaXI (2 U/μl) (NEB, catalog number: R0609，-20 °C储存)
5. CutSmart Buffer (10x) (BsaXI试剂盒配套)
6. T4 DNA Ligase (5 U/μl) (NEB, catalog number: M0202，-20 °C储存)
7. T4 DNA ligase buffer (10x) (T4 DNA ligase试剂盒配套)
8. ATP (NEB, catalog number: P0756S，-20 °C储存)
9. Phusion High-Fidelity DNA Polymerases (2 U/μl) (NEB, catalog number: M0530，-20 °C 储存)
10. dNTP Mix(10 mM each) (Vazyme, catalog number: P031-01，-20 °C储存)
11. LguI (5 U/μl) (Thermo Fisher, catalog number: ER1931，-20 °C储存)
12. QIAquick PCR Purification Kit (QIAGEN, catalog number: 28106，-20 °C储存)
13. 链霉亲和素磁珠 (南通维莱纳米科技有限公司, catalog number: BCS8047，2-8 °C储存)
14. 普通DNA产物纯化试剂盒 (天根，DP204)
15. 丙烯酰胺/甲叉双丙烯酰胺 30%溶液 (29:1) (生工，B546017-0500)
    

仪器设备

1. 96 孔磁力架
2. PCR仪 (SensoQuest)
3. 涡旋振荡器 (VORTEX-6)
4. 台式高速冷冻离心机 (Eppendorf 5427R)
5. Qubit Flex荧光计
6. 0.1-2.5 μl, 0.5-10 μl, 2-20 μl, 10-100 μl, 20-200 μl, 100-1000 μl 移液枪 (Eppendorf系列)
7. 水平电泳槽 (Mupid-2plus)
8. ATTO聚丙烯酰胺电泳槽 (ATTO, 6220)
   

实验方法

一、串联2b-RAD简化基因组文库制备

1. 酶切基因组DNA
   本研究中采用BsaXI限制性内切酶对基因组DNA进行酶切实验，用于酶切的样本DNA用量建议在100-200 ng。可以同时设置不加BsaXI酶的对照组。
   1.1

   配制下述反应体系。
   
   表1. 酶切反应体系
   
   
   
   1.2

   将上述体系至于37 °C PCR仪中反应3 h。酶切结束后将体系放置于-80 °C使酶酶活。
   1.3

   随后利用1%琼脂糖凝胶电泳检测酶切是否完全。使用BsaXI内切酶酶切基因组DNA检测结果如图1所示，对照组主带仍然清晰可见 (泳道1，3，5，7，9)，实验组主带DNA消失，泳道内呈现smear，表明已酶切完全 (泳道2，4，6，8，10)。
   
   
   图1. 酶切组与对照组检测 (对照组的泳道为1、3、5、7、9，实验组的泳道为2、4、6、8、10)
   
   
2. 接头连接
   2.1

   双链接头的制备
   分别将接头的长链和短链稀释至10 μM，等体积 (各加50 μl) 加入到0.2 ml PCR管中，混匀后置于95 °C PCR仪中，变性5 min后放置室温自然冷却, 制备成双链接头。接头长链和短链序列如表2所示，例如Ada1a接头的长短链分别是Ada1a-L和Ada1a-S，Ada1b接头的长短链分别是Ada1b-L和Ada1b-S。
   
   表2. 接头和引物序列表
   
   X代表index序列
   
   
   2.2

   以五份样品为单位进行串联文库的构建，向5份酶切产物中分别加入5组不同的接头，连接体系见表3，根据五组样品所在的串联标签位置加入不同的接头组合，在下一步扩增反应中需使用对应的引物组合 (见表4)。
   
   表3. 连接体系
   
   
   表4. 各串联标签使用的接头和引物对组合
   
   
   
   2.3

   将连接体系放置于PCR仪中，设置4 °C过夜连接
3. PCR扩增，富集标签
   3.1

   参考表5配制扩增体系： 将上述步骤所得到的五份连接产物分别利用不同的引物对组合 (见表4) 进行PCR扩增，富集连接有接头的酶切片段。反应体系如下：
   
   表5. 扩增体系
   
   
   PCR程序如下 [98 °C 5 s，60 °C 20 s，72 °C 10 s] × 16 cycles，72 °C 10 min，4 °C ∞
   
   3.2

   切胶回收产物： PCR产物用8% 非变性聚丙烯酰胺琼凝胶电泳检测，扩增产物大小约为100 bp, 将含有目标片段的胶块切下，随后用研磨棒碾碎，12000 rpm离心3 min后，加入45 μl ddH₂O, 室温静置溶胶2 h。
   3.3

   PCR产物回收
   将吸附柱置于新的1.5 ml离心管中，将上一步得到的溶胶液加入到柱中央，13000 rpm离心5 min，收集切胶产物。
   
4. 再次扩增富集产物
   4.1

   以胶回收产物为模板，再次进行扩增，如步骤3.1。
   4.2

   取3.5 μl PCR产物，PAGE胶电泳检测。观察到清晰均一目的条带后，将得到的5份扩增产物 (标签1，2，3，4，5) 进行混合后 (体积为200 μl)，利用QIAGEN QIAquick PCR purification kit进行纯化，最后用30 μl水洗脱。
5. 酶切连接
   5.1

   对上述纯化后的混合样本PCR产物进行酶切，酶切体系如下表所示：
   
   表6. 酶切体系
   
   37 °C PCR仪 (或金属浴) 中放置3 h。
   
   
   5.2

   平衡磁珠：将磁珠轻轻摇匀，吸出10 μl至0.2 ml PCR离心管中，放在磁力架上静置2 min，弃掉上清；用 20 μl 1x CutSmart buffer 清洗磁珠，反复轻柔吹吸5-10次，尽量避免产生气泡，每次洗涤结束时在磁力架上静置2 min，弃掉上清。
   5.3

   将步骤5.1得到的酶切产物加入管中，室温放置30 min，期间用枪头反复吹吸混匀，随后放置在磁力架上静置3 min，将上清转移至新的微量离心管中，加入1 μl的T4 DNA 连接酶和2 μl的10 mM ATP，16 °C反应45 min。
   5.4

   利用8%聚丙烯酰胺凝胶电泳检测串联标签产物，连接产物大小约为244 bp（如图2所示），切胶回收连接产物。(回收方法同步骤3.2和3.3，溶胶时加50 μl ddH₂O溶胶)
   
   
   图2. 标签串联的电泳检测结果
   
   
6. 二轮PCR扩增，引入文库特异性barcode序列
   6.1

   按照表7所示配制反应体系40 μl，每个样品平行扩增2管
   
   表7. 串联标签的扩增体系
   
   
   PCR程序：[98 °C 5 s，60 °C 20 s，72 °C 10 s] × 7 cycles，72 °C 10 min，4 °C ∞
   
   6.2

   取3.5 μl PCR产物PAGE胶电泳检测，观察到清晰的目的条带 (300 bp) 后，将两管PCR产物合并后利用QIAGEN QIAquick PCR purification kit试剂盒纯化，使用30 μl ddH₂O 洗脱。
   6.3

   取1.5 μl纯化产物利用PAGE胶电泳检测，扩增产物大小约为299 bp，图3电泳图上显示目的条带清晰单一，背景干净，无引物二聚体。同时，Qubit定量结果在3 ng/μl以上即可认为获得质量合格的文库。
   
   
   图3. 串联文库的检测结果
   
   

二、串联文库的测序数据分析

1. 基因分型
   参考Wang et al. (2016) 发表的串联RAD标签数据分析方法进行数据的质量控制及数据拆分分型等步骤，其中SNP分型步骤是利用Fu et al. (2013) 发表的RAD-typing软件进行分析。
   1.1

   数据质量过滤
   对Illunima测序得到的原始数据进行质量过滤，去除含有N的序列以及大于5个碱基的质量值小于10的reads。
   $perl extractN.pl –d rawdata
   $perl goodquality.pl –d extractNdata
   
   1.2

   测序Reads拼接及标签序列的拆分
   对获得的高质量双端测序reads利用PEAR软件进行拼接，获得五标签串联的完整reads，然后根据酶切位点所在的位置及五个标签的长度对串联序列进行拆分并提取含有酶切位点的高质量reads，获得所有个体的单标签2b-RAD数据。使用软件PEAR及脚本common-find.pl、f_cut.pl、Extract_cut_site.pl 实现，具体命令如下：
   a.  将pair-end数据进行对齐处理
        $perl common-find.pl -f Lb1_R1.fastq.ATCG.GoodQuality.fq -b Lb1_R2.fastq.ATCG.GoodQuality.fq 
   b.  利用pear软件对双端序列进行拼接
        $PEAR_installed_directory/pear-0.9.6-bin-64 -f Lb1_R1.fastq.ATCG.GoodQuality.fq-filter -r Lb1_R2.fastq.ATCG.
     GoodQuality.fq-filter –m 192 –n 162 -o Lb1 
   c.  根据5个标签所在的位置信息提取标签序列
        $perl f_cut.pl –i 5tagData/Lb1 –s 1 33 66 99 132 –l 50 50 50 50 60
   d.  从标签序列中提取还有酶切位点的序列
        $perl Extract_cut_site.pl -r 5tagData/Lb1_T1.fasta –c [AGTC]{9}AC[AGTC]{5}CTCC[ATGC]{7}-[ATGC]{7}GGAG[AGTC]
     {5}GT[ATGC]{9} –l 27-27 –o 5tagData/Lb1_T1wCS.fasta
   1.3

   基于RAD-typing程序进行SNP分析
   a.  利用脚本从基因组中提取BsaXI的酶切标签，构建高质量Reference。
        perl Extract_cut_site_v1.pl -r genome.fasta -c [AGTCN]{12}AC[AGTCN]{5}CTCC[ATGCN]{10}-[ATGCN]{10}GGAG[AGTCN]
     {5}GT[ATGCN]{12} -l 33-33 -o Ref_BsaXIsites
   b.  序列比对
        调用SOAP软件将群体数据与上步所得高质量Reference进行比对，允许两个碱基错配，使用脚本reads_map.pl实现。
        $ perl reads_map.pl -l 27
   c.  共显性标记分型。根据比对结果，利用最大自然法对每个个体的共显性标记进行分型。使用脚本codom_calling.pl实现。
        $ perl codom_calling.pl -a 0.05 - p 0.8
   d.  群体SNP标记筛选条件
        为筛选群体内可靠的多态性标记位点，根据分型结果，依据下列条件筛选群体多态性SNP标记：(a) 只留取标签内最多有3个SNP的标签位点；(b) 在个体内标签深度在500以上的位点不进行分型；(c) 剔除群体中基因型无多态性的位点；(d)要求群体中至少80%的个体可以分型；
        $ perl codom_filter.pl -max_snp 3 -depth 500 -maf 0.05 - p 0.8
   e.  通过上述分析流程获得群体多态性SNP标记信息可以应用于群体多样性分析、群体结构分析以及系统发生树的构建等。
   

三、群体遗传进化分析

群体遗传分析主要使用POPGENE、PLINK1.9、Structure和Mega等软件工具包实现。

1. 基本参数计算
   利用POPGENE软件统计分析多态性位点比例 (P)、Nei's基因多样性指数 (H)、Shannon氏指数 (I)、群体间基因分化系数 (G_ST)、基因流 (Nm)、Nei's 群体间遗传距离 (D) 等。采用Arlequin 2.0 软件群体内核群体间的遗传变异进行分子变异分析 (AMOVA)，根据计算出的参数，可以分析群体的遗传变异规律。使用PLINK1.9软件可以计算群体之间的遗传分化系数 (Fst)、群体内分化系数 (Fis) 以及哈温平衡检验等。
   
2. 群体结构分析
   2.1

   PCA分析
   基于群体内获得的SNP位点数据使用plink进行主成分分析，计算出群体中PC1~PCn的数值，使用主成分1 (PC1) 和主成分2 (PC2) 使用R软件绘制散点图对样本群体进行区分。
   
   
   图4. 两个群体结构主成分分析显示
   
   
   2.2

   群体结构分析
   用structure软件 (verison 2.3.4) 进行群体结构聚类分析。首先预设群体亚群数K = 1~5，基于贝叶斯模型的计算方法，根据模拟算法的到的K值分布，选择最优K值。下图所示为k值为4时，两种不同壳色群体的海湾扇贝 (100只) 的structure群体结构。
   
   
   图5. 来自两个群体的100只海湾扇贝的Structure群体结构
   
   
3. 系统发生分析
   选取待分析的SNP位点所在的标签序列，将每个个体的序列首尾相连，如果缺失相应的位点，则用"-"代替。获得的序列生成fasta文件导入MEGA软件中，分别构建NJ (Neighbor-joining) 树和ML (Maximum likelihood) 树，树的可靠性通过bootstrap 法进行检验 (重复1000次)。下图所示为100只海湾扇贝2b-RAD数据所构建的NJ树，整棵树分成主要的两大支，其中海湾扇贝红壳群体50个个体 (r1-r50) 主要聚为一支，黑壳群体大多数个体 (b1-b50) 聚为另外一支。
   
   
   图6. 基于SNP位点序列构建的两个群体的NJ系统树
   
   

方法运用

本实验方案介绍了一种串联2b-RAD标签测序技术，该技术兼容多种测序平台、串联的标签组合可灵活设置。其中，标签类型可采用不同样品、不同内切酶以及不同应用 (SNP分型或DNA甲基化水平检测) 的组合。串联标签测序技术是一种高效、灵活的全基因组遗传变异和表观遗传变异筛查和检测的手段，在大规模遗传变异筛查、群体遗传进化分析以及分子育种等研究领域具有广泛的应用前景。

致谢

本研究受山东省泰山学者计划项目和三亚崖州湾科技城管理局项目 (SKJC-KJ-2019KY01) 资助。

参考文献

1. Campbell, E. O., Brunet, B. M. T., Dupuis, J. R. and Sperling, F. A. H. (2018). Would an RRS by any other name sound as RAD? Methods Ecol Evol 9(9): 1920-1927.
2. Wang, S., Meyer, E., McKay, J. K. and Matz, M. V. (2012). 2b-RAD: a simple and flexible method for genome-wide genotyping. Nat Methods 9: 808-810.
3. Wang, S., Liu, P., Lv, J., Li, Y., Cheng, T., Zhang, L., Xia, Y., Sun, H., Hu, X. and Bao, Z. (2016). Serial sequencing of isolength RAD tags for cost-efficient genome-wide profiling of genetic and epigenetic variations. Nat Protoc 11(11): 2189-2200.
4. Fu, X., Dou, J., Mao, J., Su, H., Jiao, W., Zhang, L., Hu, X., Huang, X., Wang, S. and Bao, Z. (2013). RADtyping: an integrated package for accurate de novo codominant and dominant RAD genotyping in mapping populations. PLoS One 8(11): e79960.