混合分子标记扩增子高通量测序方案
A Mixed Molecular Marker Amplicon Sequencing Scheme Based on High Throughput Sequencing Platform   

材料与试剂

1. 200 μl和1.5 ml微量离心管 (BBI Life Sciences，常温储存)
2. PCR八联排管和96孔板 (BBI Life Sciences，常温储存)
3. 各种型号枪头 (10 μl，200 μl，1 ml，BBI Life Sciences，常温储存)
4. 基因组DNA提取试剂盒 (北京天根生化科技有限公司，货号: DP304-03，常温储存)
5. 通用型DNA 纯化回收试剂盒 (北京天根，货号: DP214-03，常温储存)
6. EasyTaq DNA Polymerase (北京全式金生物技术有限公司，货号: AP111-01，-20°C储存)
7. TransTaq DNA Polymerase High Fidelity (HiFi) (北京全式金，货号: AP131-03，-20°C储存)
8. NEBNext dsDNA Fragmentase (New England Biolabs，货号: M0348L，-20°C储存)
9. NEBNext Ultra DNA Library Prep Kit (New England Biolabs，货号: E7645L，-20°C储存)
10. 10x Tango Buffer (Fermentas lnc.，货号: BY-5，-20°C储存)
11. T4 DNA ligase (5 U/μl) (Thermo Fisher Scientific，货号: EL0014，-20°C储存)
12. T4 DNA Polymerase (5 U/μl) (Thermo Fisher Scientific，货号: EP0062，-20°C储存)
13. T4 Polynucleotide Kinase (10 U/μl) (Thermo Fisher Scientific，货号: EK0032，-20°C储存)
14. 10 mM dNTPs (上海生工生物工程股份有限公司，货号: A610056，-20°C储存)
15. 100 mM ATP (上海生工，货号: A600311，-20°C储存)
16. 100 mM dATP (上海生工，货号: B500044，-20°C储存)
17. 0.5 M EDTA， pH8.0 (上海生工，货号: B300599，4°C储存)
18. 1x TE Buffer (上海生工，货号: B548106，4°C储存)
19. 琼脂糖 (Invitrogen，货号: 75510-019，常温储存)
20. 双蒸水 (广州誉维生物科技仪器有限公司Unique超纯水机制)
    

仪器设备

1. 0.5-5 μl，2-20 μl，20-200 μl，100-1000μl移液器 (Thermo Fisher Scientific F2系列)
2. Gene Pro PCR扩增仪 (杭州博日科技股份技术有限公司，catalog number: TC-E-96G)
3. Vortex Genie2漩涡混合器 (Scientific Industries，catalog number: G560E)
4. Cubee迷你离心机 (广州美津生物技术有限公司，catalog number: aqbd-i)
5. 超微量分光光度仪 (Thermo Fisher Scientific，catalog number: NanoDrop2000)
6. 通用电泳仪 (北京百晶生物技术有限公司，catalog number: BC-power6001)
7. 高压蒸汽灭菌锅 (Kagoshima Seisakusyo Inc.，catalog number: SX-500)
8. 服务器 (linux系统，32G以上的内存，500G的存储)
   

软件

1. Trimmomatic v 0.39 http://www.usadellab.org/cms/?page=trimmomatic
2. Seqtk_demultiplex https://github.com/jameslz/fastx-utils/blob/master/seqtk_demultiplex/
3. Trinity_r20140717 https://sourceforge.net/projects/trinityrnaseq/files/PREV_CONTENTS/previous_releases/
4. Java 1.7.0_80 https://www.java.com/zh-CN/download/help/linux_x64_install.html
5. BWA 0.7.10 https://sourceforge.net/projects/bio-bwa/files/
6. Samtools 1.4 https://sourceforge.net/projects/samtools/files/samtools/
7. BLAST 2.6.0 ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/
8. Python 2.7.9 https://www.python.org/ftp/python/2.7.9/
9. Biopython 1.68 https://pypi.org/project/biopython/1.68/
10. Phylopypruner https://gitlab.com/fethalen/phylopypruner/-/wikis/Installation
11. FastTree http://www.microbesonline.org/fasttree/#Install
12. PRANK http://wasabiapp.org/download/prank/
13. RAxML v8.2.0 https://github.com/stamatak/standard-RAxML
14. PartitionFinder2 https://github.com/brettc/partitionfinder/releases/tag/v2.1.1
15. ASTRAL 5.7.4 https://github.com/smirarab/ASTRAL
    

实验步骤

读者首先需要对所选的分子标记进行PCR扩增，获得大量的扩增子 (目标扩增子的选择可参考研究背景中列出的适用于不同生物类群的分子标记工具箱)，随后再按照: 扩增子混合，DNA片段化，样本特异性标签连接，多样本混合和片段筛选，文库构建和高通量测序共五个步骤，实现混合分子标记扩增子的高通量测序，其工作流程图如图1所示。

图1. 混合PCR产物高通量测序方案的实验流程图

一、对所选分子标记进行PCR扩增
在大规模扩增分子标记之前，建议开展一个预实验，随机挑选2-3个样本和5-10种分子标记分别进行PCR扩增，目的是确保PCR反应体系，引物和DNA模板可正常使用。

1. 为保证PCR扩增成功率，推荐使用巢式PCR分别对每个分子标记进行两轮扩增 (标准PCR扩增与巢式PCR扩增的电泳结果比对图如图2)。如果分子标记数量较多, 推荐使用96孔板进行反应。以下PCR程序仅为示例, 需根据自己实验室习惯进行修改。
   第一轮PCR: 按下表建立20 μl的反应体系:
   
   
   PCR反应程序:
   
   
   第二轮PCR: 按下表建立20 μl的反应体系:
   
   
   PCR反应程序:
   
   
   
2. 使用1.2 %的琼脂糖凝胶电泳检测二轮PCR扩增产物 (如图2b)。
3. 第一轮和第二轮PCR产物可保存于-20 °C (建议保存时间不超过1年)。为了提高实验效率，建议读者在获得所有样本的PCR产物之后再进行后续操作。
   
   
   图2. 标准 PCR 扩增 (a) 与巢式 PCR 扩增 (b) 的电泳结果比对图。S1-14对应14个样本。

二、扩增子混合 (参考图1步骤1)
为了简化实验操作流程，降低实验劳动强度，本方法并没有对每个PCR产物进行单独纯化，而是直接等体积混合同一物种不同的PCR产物，对于无可见条带或有非特异性扩增条带的PCR产物也同等对待。

1. 将同一物种的不同PCR产物等体积 (5 μl) 混合。
2. 使用DNA 纯化回收试剂盒纯化混合PCR产物，具体操作参阅试剂盒说明书。https://www.tiangen.com/asset/imsupload/up0609376001604559269.pdf
3. 使用NanoDrop 2000进行浓度测定，纯化产物可长期保存于-20 °C。

三、DNA片段化 (参考图1步骤2)
该步骤的目的是为了获得符合Illumina测序平台的要求的，具有一定长度范围的DNA片段。以HiSeq X-ten测序平台为例，该平台期望的插入片段大小在200-400 bp之间，所以片段化后的DNA片段最好集中在这个区域。片段化的方法有超声处理，酶处理等，对于新手而言，推荐使用后者，因为该过程不依赖特殊仪器 (超声波DNA打断仪)，且入手简单。本方法以NEBNext dsDNA Fragmentase为例，介绍片段化方法，其他注意事项可参阅产品说明书。https://international.neb.com/protocols/0001/01/01/digestion-with-nebnext-dsdna-fragmentase-m0348.

1. 取1 μg纯化的混合PCR产物于0.2 ml离心管中。
2. 按下表建立片段化反应体系:
   
   
   注: 建议片段化体系的准备工作在冰上完成。
   
   
3. 将配制好的酶切反应体系充分振荡混匀，离心，放入PCR仪中，37°C孵育10-25 min。
   注: 强烈建议读者做时间梯度预实验 (如: 10 min，12 min，14 min等) 来确定最优的片段化时间。判断标准是，与片段化之前的混合PCR产物相比，片段化后的DNA片段集中在200-400 bp之间，如图3。
   
   
   图3. 片段化电泳结果示意图。1号泳道代表片段化之前的混合PCR产物， 2号泳道代表片段化之后的混合PCR产物。
   
   
4. 迅速将离心管转移至冰上，加入2 μl 0.5 M EDTA终止片段化反应。充分振荡混匀，离心，65°C孵育15 min，使片段化酶彻底失活。
5. 取2 μl片段化产物，使用2%的琼脂糖凝胶电泳进行检测，如果依然有主带或片段大小不在期望范围内，需要重新进行片段化实验。
6. 使用DNA 纯化回收试剂盒纯化片段化产物，最后使用37.4 μl 0.1x TE缓冲液进行洗脱，建议立即进行下一步操作。
   
   

1. 平末端与磷酸化
   1)

   取新的0.2 ml离心管， 按下表建立末端磷酸化体系:
   
   
   
   2)

   使用移液器吹打混匀 (动作要轻柔，不要产生气泡)，之后简短离心，将离心管置于PCR仪中，12 °C孵育30 min，37°C 孵育30 min，65°C孵育20 min，4°C短暂保存。
   3)

   使用DNA纯化回收试剂盒纯化末端磷酸化产物，最后使用43 μl 0.1x TE缓冲液进行洗脱。
   
   
2. 平末端加A尾
   1)

   将上一步获得的平末端与磷酸化的纯化产物转移至新的0.2 ml离心管中。
   2)

   按下表建立平末端加A尾的反应体系:
   
   
   
   
   
   3)

   充分振荡混匀，之后简短离心，将离心管置于PCR仪中，72°C孵育30 min。
   4)

   待反应体系降至室温后，使用DNA纯化回收试剂盒进行纯化，最后使用34 μl 0.1x TE缓冲液进行洗脱。
3. 特异性标签连接
   为了区分不同样本，本方法设计了16个特异性标签 (见附录)。标签设计有如下要求: (1) 标签序列长度为9个碱基; (2) 标签两两之间的差异大于两个碱基; (3) 标签内部避免连续两个以上的碱基重复; (4) 退火后双链标签的一条链的3’末端突出一个碱基C， 目的是为了防止退火后的双链特异标签发生自连接，另一条链的5’末端突出一个碱基T，目的是与补平加A磷酸化后的DNA片段间进行T/A连接，提高连接效率 (如图4)。若读者需要更多的特异性标签，可参考上述思路自行设计。
   1)

   将上一步获得的平末端加A尾的纯化产物转移至新的0.2 ml离心管中。
   2)

   按下表建立连接体系:
   
   
   *双链特异性标签的制备方法见溶液配方。
   
   
   
   3)

   充分振荡混匀，之后简短离心，将离心管置于PCR仪中，16°C连接过夜 (至少10 h)。
   
   
   图4.样本特异性标签连接反应流程图
   

五、多样本混合和片段筛选 (参考图1步骤4)
连接反应结束后，将具有不同特异标签的连接产物进行混合，并在2%的琼脂糖凝胶中割胶回收200-400 bp的DNA片段。

六、文库构建和高通量测序 (参考图1步骤5)
使用高通量文库构建试剂盒 (如: NEBNext Ultra DNA Library Prep Kit) 对割胶产物进行文库构建 (具体方法参阅产品说明书https://international.neb.com/protocols/0001/01/01/digestion-with-nebnext-dsdna-fragmentase-m0348)，并进行高通量测序。

数据分析

数据分析的主要目的是将高通量原始数据处理成可直接进行系统发育分析的数据，其流程包括: 数据分选与组装，目标基因序列的鉴定与提取，序列比对，数据质量控制和系统发育分析 (如图5)。本方法提供了与上述分析流程相关的Python脚本，供读者参考使用，下载地址为https://datadryad.org/stash/dataset/doi:10.5061/dryad.n21cr.
一、数据分选与组装

1. 使用Trimmomatic (v 0.39) 对测序得到的原始数据 (raw data) 进行质量控制和筛选，去除reads中的测序接头和引物序列，去除平均碱基质量值低于20的reads，去除reads首端和末端碱基质量小于3或N的碱基，得到高质量的测序数据 (clean reads)。
   以处理双端测序数据为例:
   $ java -jar trimmomatic-0.39.jar PE RawReads_R1.fq.gz
   RawReads_R2.fq.gz
   CleanReads_R1_paired.fq.gz
   CleanReads_R1_unpaired.fq.gz CleanReads_R2_paired.fq.gz
   CleanReads_R2_unpaired.fq.gz ILLUMINACLIP:TruSeq3-
   PE.fa:2:30:10:2:keepBothReads AVGQUAL:20 LEADING:3 TRAILING:3 MINLEN:36
2. 依据各个样本的特异性标签 (barcode)， 使用fastq-multx将clean reads分选到对应的实验样本中。
   以分选双端测序数据为例:
   $ fastq-multx -B barcode.txt -b CleanReads_R1_paired.fq.gz
   CleanReads_R2_paired.fq.gz -o %_R1_clean.fq %_R2_clean.fq
   注: %指的是输出文件名与barcode文件中的样本名称相对应。详细使用方法可参阅https://github.com/brwnj/fastq-multx
   
   
   图5.高通量数据分析流程图
   
   
3. 使用转录组拼接软件Trinity对每个样本的clean reads进行组装。
   以组装样本A的clean reads为例:
   $ Trinity --seqType fq --JM 24G --CPU 4 --PasaFly --left
   taxaA_R1_clean.fq --right taxaA_R2_clean.fq --out taxaA
4. 使用BWA和Samtools计算每个contig的测序深度。
   以处理样本A的contigs为例，按顺序运行以下命令:
   $ bwa index taxaA_trinity.fasta
   $ bwa -t 4 taxaA_trinity.fasta taxaA_ R1_clean.fq taxa_R2_clean.fq > BWA_taxaA.sam
   $ samtools view -@ 4 -bS BWA_taxaA.sam > BWA_taxaA.bam
   $ samtools sort -@ 4 BWA_taxaA.bam -o BWA_taxaA_sorted.bam
   $ samtools depth -@ 4 BWA_taxaA_sorted.bam > taxaA_coverage.txt
5. 删除平均测序深度较低的contigs，避免由于测序深度不足而引入的拼接错误:
   $ python 2_compute_contigs_coverage.py (计算contigs平均测序深度)
   $ python 3_delete_contigs_with_low_coverage.py (删除平均测序深度小于10的contigs)

1. 使用BLAST方法从每个物种过滤后的contigs中鉴定候选的目标基因序列。以鉴定102 个脊椎动物核蛋白编码基因标记为例，读者需要提前从Genbank数据库中下载人相对应分子标记的目标片段序列 (命名为“reference.fasta”)，以作为BLAST的数据库，然后运行以下命令:
   $ python 4_tblastx.py (进行TBLASTX)
   注: 读者需要根据目标基因序列的性质 (编码序列或非编码序列)，决定使用哪种BLAST比对方法，如果鉴定非编码的直系同源序列，推荐使用BLASTN。
2. 优化BLAST结果，删除冗余序列。按顺序运行以下命令:
   $ python 5_filter_blast_records.py (优化BLAST结果)
   $ python 6_get_sequences.py (以物种为单位，提取目标基因序列)
   $ python 7_sort_sequences_by_gene.py (以基因为单位，生成目标基因序列文件)
3. 切除引物序列:
   $ python 8_trim_primers.py (切除引物序列)
4. 确保每一条目标编码序列能够翻译成氨基酸。按顺序运行以下命令:
   $ python 9_check_sequence_translatable.py (检查目标基因序列是否可以翻译成氨基酸)
   $ python 10_combine_seqences.py
   $ python 11_check_again.py (再次检查目标基因序列是否可以翻译成氨基酸)

三、序列比对
对于进化速度较慢的基因标记，如核蛋白编码基因，推荐使用PRANK，MUSCLE等软件进行序列比对，对于进化速度较快的基因标记，比如含有内含子序列的EPIC (Exon-primed intron-crossing) 分子标记，推荐使用SATé进行序列比对。

1. 若比对核编码基因等保守序列，以处理一个编码基因A为例:
   $ prank -d=GeneA.fasta -o=GeneA_align.fasta -f=fasta +F -noxml -notree -codon –once
   注: 具体使用方法请参阅http://wasabiapp.org/download/prank/
2. 若比对内含子等高变异序列，以处理一个基因B为例:
   $ python run_sate.py -d DNA --move-to-blind-on-worse-score --start-tree-search-from-current --merger=opal -i GeneB.fasta
   注: 具体使用方法请参阅http://phylo.bio.ku.edu/software/sate/sate.html

四、数据质量控制
为了进一步提高数据质量，确保鉴定的直系同源基因数据集的准确性，本方法对每个alignment 构建单基因树，通过分析基因树中各分枝的长度和判断基因树的拓扑结构，将枝长异常的序列，旁系同源序列以及污染序列从数据集中剔除(如图6)。

1. 单基因树构建:
   $ FastTree -gtr -nt GeneA_align.fasta > GeneA.tree
2. 借助Phylopypruner去除长枝 (含有长枝的基因树拓扑结构可参照图6中的拓扑结构A)。例如，删除某一枝长大于树中所有分枝标准差3倍的序列，运行以下命令:
   $ Phylopypruner --trim-lb 3 GeneA_align.fasta GeneA.tree
   注: 具体使用方法可参阅https://gitlab.com/fethalen/phylopypruner/-/wikis/Introduction
3. 去除旁系同源序列或污染序列，其标准是观察单基因树拓扑结构中是否存在聚类关系怪异的分支，可参照图6中基因树拓扑结构B和C同正确的基因树拓扑结构之间的差异。
4. 对数据质量控制后的alignment再次进行比对，否则，可直接进行下一步。
5. 使用Gblocks软件剔除比对模糊的区域，得到可以直接进行系统发育分析的比对矩阵:
   $ Gblock GeneA_align.fasta -t=c -b5=a
   注: Gblock软件的详细参数介绍可参阅http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_documentation.html

1. 串联树构建: 
   1)

   寻找最佳分区模型:
   $ python PartitionFinder2.py Input_Folder
   注: Input_Folder中包含2个文件: (1) phlipy格式的串联数据集; (2) 控制文件“partition_finder.cfg”，里面包含模型选择，基因分区信息等。具体使用方法可参阅http://www.robertlanfear.com/partitionfinder/tutorial/
   2)

   使用RAxML构建ML树: 
   $ raxmlHPC-PTHREADS-AVX -T 4 -f a -x 12347 -p 12845 -# 500 -m
   GTRGAMMA -s CombineGenes.fasta -n
   CombineGenes_MLtree_500boostrap -q Partion_File
   
2. 物种树构建: 
   1)

   使用RAxML构建单基因树， 以处理一个基因为例:
   $ raxmlHPC-PTHREADS-AVX -T 4 -f a -x 12347 -p 12845 -# 100 -m
   GTRGAMMA -s GeneA _align_cut.fasta -n GeneA _MLtree_100boostrap
   
   2)

   使用ASTRAL构建物种树:
   $ java -jar astral.5.7.4.jar -i Gene_best_trees.tre -b Gene_bs_files -o Gene_astral_result
   注: astral参数的详细使用说明请参阅 https://github.com/smirarab/ASTRAL/blob/master/astral-tutorial.md

溶液配方

1. 特异性标签制备方法
   在0.2 ml离心管中， 建立特异性标签退火体系:
   
   
   充分振荡混匀，利用下述条件进行退火反应，形成局部互补的双链特异性标签: 95°C 5 min; 95°C降到12°C，每秒0.1°C; 保持在12°C，将退火后的双链标签放置于-20°C保存。
   

致谢

实验方案摘自发表的文章Feng, Y. J., Liu, Q. F., Chen, M. Y., Liang, D. and Zhang, P. (2016). Parallel tagged amplicon sequencing of relatively long PCR products using the Illumina HiSeq platform and transcriptome assembly. Molecular Ecology Resources 16(1): 91-102. http://doi:10.1111/1755-0998.12429.

竞争性利益声明

作者声明没有利益冲突

参考文献

1. Che, L. H., Zhang, S. Q., Li,Y., Liang,D., Pang,H., Ślipiński, A. and Zhang, P. (2017). Genome-wide survey of nuclear protein-coding markers for beetle phylogenetics and their application in resolving both deep and shallow-level divergences. Molecular Ecology Resources 17(6): 1342-1358. 
2. Feng, Y. J., Liu, Q. F., Chen, M. Y., Liang, D. and Zhang, P. (2016). Parallel tagged amplicon sequencing of relatively long PCR products using the Illumina HiSeq platform and transcriptome assembly. Molecular Ecology Resources 16(1): 91-102.
3. Kapli, P., Yang, Z. H. and Telford, M. J. (2020). Phylogenetic tree building in the genomic age. Nature Review Genetics 21(7): 428-444.
4. Li, J. N., He, C., Guo, P., Zhang, P. and Liang, D. (2017). A workflow of massive identification and application of intron markers using snakes as a model. Ecology & Evolution 7(24): 10042-10055.
5. Shen, X. X., Liang, D., Feng, Y. J., Chen, M. Y. and Zhang, P. (2013). A versatile and highly efficient toolkit including 102 nuclear markers for vertebrate phylogenomics, tested by resolving the higher level relationships of the Caudata. Molecular Biology and Evolution 30(10): 2235-2248. 
6. Zhang, Y., Deng, S. H., Liang, D. and Zhang, P. (2019). Sequence capture across large phylogenetic scales by using pooled PCR-generated baits: A case study of Lepidoptera. Molecular Ecology Resources 19(4): 1037-1051.