Preliminary Identification of Bacterial Species Based on 16S rRNA Gene and Genomic Sequence   

材料与试剂

1. Chelex树脂 (BIO-RAD，catalog number: 143-2832)
2. PCR反应Mix Taq：Premix Taq™ (Ex Taq™ Version 2.0)〔宝生物工程 (大连) 有限公司，Takara，catalog number: RR003A〕
   

仪器设备

1. PCR仪 (德国SENSO，Sensoquest LabCycler PCR仪)
2. 离心机 (日本TOMY, MicroOne离心机)
   

软件和数据库

1. blast 2.8.1+：https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.8.1/；操作系统：windows或linux。
2. FastANI v1.1：https://github.com/ParBLiSS/FastANI；操作系统：linux。
3. R语言bactaxR包：https://github.com/lmc297/bactaxR；操作系统：windows或linux。
4. CompareM：https://github.com/dparks1134/CompareM；操作系统：linux。
5. MAFFT V7：https://mafft.cbrc.jp/alignment/software/；操作系统：windows或linux。
6. MEGA V5.2：https://www.megasoftware.net/downloadform；操作系统：windows。
7. UBCG：https://www.ezbiocloud.net/tools/ubcg；操作系统：linux。
8. 原核生物全基因组数据库：https://www.ncbi.nlm.nih.gov/genome/browse#!/pro karyotes/ 
9. EzBioCloud数据库：https://www.ezbiocloud.net/
   

实验步骤

1. 染色体DNA提取和纯化
   1.1

   用无菌接种环或竹签挑取1个菌落，加入50 µl 5%的Chelex溶液，充分振荡混匀。
   1.2

   沸水浴10-15 min。
   1.3

   2,000 × g离心1 min。上清液即为提取的基因组DNA溶液，-20 °C保存备 用。
2. 16S rRNA基因的PCR扩增
   50 μl反应总体系包含：
   2 × Mix Taq                         25 μl
   Primers (10 μM)
        27F                                  2.0 μl
        1492R                             2.0 μl
   DNA模板 (100 ng/μl)         2.0 μl
   ddH₂O                                 19 μl
   其中扩增引物：27F (5'- AGAGTTTGATCMTGGCTCAG-3') 和1492R (5'- GGTTACCTTGTTACGACTT-3')，为生工生物工程 (上海) 股份有限公司合成。
   PCR扩增参数：94 °C预变性，4 min；94 °C变性，1 min；55 °C复性，1 min；72 °C延伸，1.5 min，共30个循环；72 °C延伸10 min。
3. 16S rRNA基因序列测定
   将16S rRNA基因扩增产物送至测序公司测序。
4. 16S rRNA基因序列本地化比对分析
   
   4.1

   本地化数据库构建
   用于BLAST分析的原核生物模式菌株的16S rRNA基因序列本地化数据库, 可直接在美国国家生物技术信息中心 (National Center for Biotechnology Information, NCBI) 的ftp站点 (ftp://ftp.ncbi.nih.gov/blast/db/16S_ribosomal_RNA.tar.gz) 下载。本文以构建Cryobacterium属模式菌株16S rRNA基因序列本地化数据库为例，在https://lpsn.dsmz.de/genus/cryobacterium 网页下载Cryobacterium属模式菌株的16S rRNA基因序列，并通过检索EzBioCloud数据库查看该属是否有更新，最终确定Cryobacterium属共15个模式菌株。序列命名为Cryo16S.fasta后，利用BLAST+工具构建本地化数据库。在windows系统CMD命令行窗口中，输入以下命令：makeblastdb.exe -in Cryo16S.fasta -parse_seqids -hash_index -dbtype nucl -out D:\blast\db\Cryo16S，如图1，生成以Cryo16S命名的本地化数据库。
   
   
   图1. 本地数据库构建
   
   
   4.2

   16S rRNA基因序列比对
   本文以Liu et al. (2020) 论文中报道的Cryobacterium属分离菌株数据为例，可在文章所列的菌株信息附表中查询到16S rRNA基因序列的GenBank号。根据分离菌株的GenBank号，利用NCBI的Entrez在线工具 (https://www.ncbi.nlm.nih.gov/sites/batchentrez) 下载原始数据。准备好其中65株菌株的16S rRNA基因序列fasta文件 (命名为Cryo_query.fasta) 后，利用BLAST+工具比对。在windows系统CMD命令行窗口中，输入以下命令：blastn.exe -query Cryo_query.fasta -out Cryo_16S_blast.xls -db D:\blast\db\ Cryo16S -outfmt 6 -evalue 1e-5 -max_target_seqs 1 -num_threads 8，如图2。
   
   
   图2. 本地blast命令
   
   生成Cryo_16S_blast.xls命名的结果文件 (表1)，每条序列输出一条与其相似性最高的物种信息。
   Kim et al. (2014) 通过对数千个基因组和16S rRNA基因序列的统计分析发现，当两菌株16S rRNA基因序列之间相似性低于 ≈ 98.65%时，可判断它们归属于不同的种；而高于 ≈ 98.65%时，可能属于同一个种，也可能属于不同的种，需结合全基因组序列分析等其他方法判定。根据比对结果可知，仅有两个菌株与已知物种间序列相似性低于98.65%，为潜在的新物种，其余菌株与已知物种最高相似在98.76-100%之间，因此，需根据其基因组序列计算全基因组平均核苷酸序列一致性 (average nucleotide identity，ANI)，得到精确的物种鉴定结果。
   
   表1. Cryobacterium 菌株的16S rRNA基因序列比对结果

   菌株	最相似的模式菌株	相似性 (%)
   MDB2-33-2	C. breve TMT4-23T	98.42
   MDB2-10	C. breve TMT4-23T	98.48
   MDB2-B,MDB1-5,MDB1-18-1,MDB2-A-2,TMT1-19,MDB2-A-1,TMT1-23-1,TMT3-29-2,Hz7, TMT1-62,TMT2-48-2,Sr3,Hh4,MDB1-18-2,TMT2-16,TMT2-14,TMT2-17-1,TMT2-4,TMT4-10, TMT2-59,TMT1-21,TMT2-42-4,TMT1-22,TMT2-18-3,TMT2-18-2,TMT2-15-1,TMT2-10, TMT2-23,MDT1-3,HLT2-28,TMT1-51,HLT2-23,RHLT2-21,TMT4-31	C. breve TMT4-23T	98.76-99.93
   Hh14	C. melibiosiphilum Hh39T	99.22
   Sr8	C. psychrotolerans 0549T	99.35
   RHLS22-1, Sr59	C. soli GCJ02T	99.57-99.93
   Sr54,Hh38,Hz9,Hz16,Sr39,TMS1-13-1	C. levicorallinum Hh34T	99.65-99.93
   TMB1-8, TMB1-7,TMB3-15,TMB3-10,TMN-39-1,TMB3-12,TMB3-1-2	C. zongtaii TMN-42T	99.658-99.93
   TMT1-1, Hb1	C. luteum Hh15T	99.86
   TMT1-66-1	C. flavum Hh8T	99.93
   Hh7, TMS1-20-1,Hh11	C. levicorallinum Hh34T	100
   TMT1-3	C. luteum Hh15T	100
   TMT1-2-2	C. flavum Hh8T	100
   TMN-39-2	C. zongtaii TMN-42T	100
   TMT3-12	C. breve TMT4-23T	100
   TMT1-2-1,Sr47	C. psychrotolerans 0549T	100

   
   
5. 基因组平均核苷酸和氨基酸序列一致性分析
   ANI是基于两两基因组之间所有直系同源蛋白编码基因序列比较的平均值，该值可反映基因组之间的进化距离。基因组之间的ANI 值与16S rRNA基因序列相似性分析及DNA-DNA杂交结果相一致，因此，ANI分析方法已经代替繁琐的DNA-DNA杂交技术。两菌株基因组ANI 值为95-96%时，相当于DNA-DNA杂交值70%。因此，当两菌株基因组ANI值大于96%时，鉴定为同一个物种；小于95%时，鉴定为不同物种 (Richter et al., 2009)。
   5.1

   基因组数据准备
   本文所使用全基因组数据来源于Liu et al. (2020) 文章或原核生物全基因组数据库。目前有多个工具可直接对全基因组进行ANI值计算，如，Jspeciesws和ANI Calculator等在线计算工具，Jspecies等图形界面工具，以及FastANI、autoANI等命令行工具 (表2)。本文选择计算速度最快、操作最为方便的FastANI工具介绍，FastANI需在linux系统中运行，可在普通计算机中安装ubuntu系统运行FastANI工具。首先建立Cryobacterium属模式菌株以及待鉴定菌株的全基因组序列文件，准备两个文本文件：将待比对序列文件名称保存至query.list文件 (图3a)，参比序列文件名称保存至ref.list (图3b)，其中每一个基因组文件占用一行，用于FastANI软件计算的输入文件。
   
   表2. ANI计算工具

   名称	软件网页
   FastANI	https://github.com/ParBLiSS/FastANI
   Jspeciesws	http://jspecies.ribohost.com/jspeciesws/#Home
   ANI Calculator	http://www.ezbiocloud.net/tools/ani
   Jspecies	http://imedea.uib-csic.es/jspecies/docu.html
   OAT	http://www.ezbiocloud.net/tools/orthoani
   OrthoANI	http://www.ezbiocloud.net/tools/
   autoANI	https://github.com/osuchanglab/autoANI
   GET_HOMOLOGUES	https://github.com/eead-csic-compbio/get_homologues

   
   
   图3. query.list和ref.list文件部分内容示例
   
   
   5.2

   ANI计算
   打开shell窗口后，进入基因组文件、query.list和ref.list文件所在的目录，输入命令：fastANI --ql query.list --rl ref.list -o ANI.out。运行过程如图4，计算结束后，输出结果到ANI.out文件。
   表3列出了部分ANI计算结果，根据待鉴定菌株与Cryobacterium模式菌株的ANI值可知，其中45株菌ANI值低于95%，无法鉴定为已知的物种，属于潜在的新种，另外20株与已知物种的ANI值高于96%，可鉴定到种，鉴定结果见表3。
   
   
   图4. FastANI运行过程
   
   表3. 基于基因组ANI值对未知的Cryobacterium菌株的鉴定结果

   菌株	ANI值最高的参比菌株	ANI值 (%)	鉴定结果
   MDB2-33-2	C. psychrotolerans 0549T	81.15	新种
   MDB2-10	C. breve TMT4-23T	81.27	新种
   MDB2-B	C. breve TMT4-23T	81.15	新种
   MDB1-5	C. breve TMT4-23T	81.42	新种
   MDB1-18-1	C. psychrotolerans 0549T	81.18	新种
   MDB2-A-2	C. psychrotolerans 0549T	81.20	新种
   TMT1-19	C. breve TMT4-23T	87.00	新种
   MDB2-A-1	C. psychrotolerans 0549T	81.19	新种
   Hh14	C. melibiosiphilum Hh39T	87.43	新种
   TMT1-23-1	C. breve TMT4-23T	87.35	新种
   TmT3-29-2	C. breve TMT4-23T	87.47	新种
   Hz7	C. breve TMT4-23T	86.80	新种
   Sr8	C. psychrotolerans 0549T	98.15	C. psychrotolerans
   TMT1-62	C. breve TMT4-23T	86.76	新种
   TmT2-48-2	C. breve TMT4-23T	87.25	新种
   Sr3	C. breve TMT4-23T	86.86	新种
   Hh4	C. breve TMT4-23T	82.71	新种
   MDB1-18-2	C. psychrotolerans 0549T	81.22	新种
   TMT2-16	C. breve TMT4-23T	86.86	新种
   TMT2-14	C. breve TMT4-23T	86.82	新种
   TMT2-17-1	C. breve TMT4-23T	86.87	新种
   TMT2-4	C. breve TMT4-23T	86.87	新种
   TMT4-10	C. breve TMT4-23T	84.45	新种
   TmT2-59	C. breve TMT4-23T	84.46	新种
   TMT1-21	C. breve TMT4-23T	84.43	新种
   TMT2-42-4	C. breve TMT4-23T	86.88	新种
   TMT1-22	C. breve TMT4-23T	84.34	新种
   TMT2-18-3	C. breve TMT4-23T	86.85	新种
   TMT2-18-2	C. breve TMT4-23T	86.80	新种
   TMT2-15-1	C. breve TMT4-23T	87.00	新种
   TMT2-10	C. breve TMT4-23T	84.47	新种
   TMT2-23	C. breve TMT4-23T	84.21	新种
   RHLS22-1	C. zongtaii TMN-42T	90.73	新种
   MDT1-3	C. breve TMT4-23T	81.64	新种
   Sr54	C. aureum Hh31T	86.47	新种
   TMB1-8	C. zongtaii TMN-42T	96.31	C. zongtaii
   TMB1-7	C. zongtaii TMN-42T	97.87	C. zongtaii
   HLT2-28	C. breve TMT4-23T	82.57	新种
   TMT1-51	C. breve TMT4-23T	89.38	新种
   HLT2-23	C. luteum Hh15T	90.01	新种
   Hh38	C. levicorallinum Hh34T	99.11	C. levicorallinum
   Hz9	C. ruanii Sr36T	86.61	新种
   RHLT2-21	C. breve TMT4-23T	82.56	新种
   TMT1-1	C. aureum Hh31T	86.13	新种
   Hb1	C. aureum Hh31T	85.92	新种
   Hz16	C. ruanii Sr36T	86.83	新种
   TMB3-15	C. zongtaii TMN-42T	97.89	C. zongtaii
   TMB3-10	C. zongtaii TMN-42T	97.86	C. zongtaii
   Sr39	C. aureum Hh31T	87.46	新种
   TMT1-66-1	C. flavum Hh8T	97.33	C. flavum
   TMT4-31	C. zongtaii TMN-42T	97.86	C. zongtaii
   TMN-39-1	C. zongtaii TMN-42T	99.10	C. zongtaii
   Sr59	C. zongtaii TMN-42T	88.16	新种
   TMS1-13-1	C. levicorallinum Hh34T	96.49	C. levicorallinum
   TMB3-12	C. zongtaii TMN-42T	97.90	C. zongtaii
   TMB3-1-2	C. zongtaii TMN-42T	97.87	C. zongtaii
   Hh7	C. levicorallinum Hh34T	97.02	C. levicorallinum
   TMS1-20-1	C. levicorallinum Hh34T	97.09	C. levicorallinum
   Hh11	C. levicorallinum Hh34T	96.94	C. levicorallinum
   TMT1-3	C. luteum Hh15T	99.49	C. luteum
   TMT1-2-2	C. flavum Hh8T	97.46	C. flavum
   TMN-39-2	C. zongtaii TMN-42T	99.07	C. zongtaii
   TmT3-12	C. breve TMT4-23T	99.99	C. breve
   TMT1-2-1	C. psychrotolerans 0549T	99.58	C. psychrotolerans
   Sr47	C. psychrotolerans 0549T	90.92	新种

   
   
   5.3

   平均氨基酸一致性分析
   ANI值可作为种水平的鉴定方法之一，而在更高分类单元的鉴定可通过计算氨基酸一致性 (amino acid identity, AAI) 进行初步判定，在属水平的AAI阈值大约为65-72% (Konstantinidis and Tiedje, 2007)。以本文中的基因组数据为例，AAI值可通过CompareM软件计算，运行操作系统为ubuntu。以~/data/my_genomes文件夹中的基因组序列为输入文件，输入命令：comparem aai_wf ~/data/my_genomes aai_output。
   运行结束后，结果保存至aai_output文件夹，AAI值保存在aai_summary.tsv文件，结果显示Cryobacterium mesophilum DSM 19267T与其他79株菌 (包括14个模式菌株和65个分离菌株) 为64.3-65.86%，其余菌株两两间AAI值为70.86-100%。
   5.4

   ANI聚类分析
   以上过程初步得到了每株菌的鉴定结果，为了分析菌株间亲缘关系，可对两两菌株间ANI值聚类分析。将所有待聚类分析的序列文件路径保存至pairwise.list，打开shell窗口后，进入工作目录，输入命令：
   fastANI --ql pairwise.list --rl pairwise.list -o ANI_pairwise.out
   得到ANI_pairwise.out结果文件，结果包含两两菌株间ANI值。以此文件为输入文件，在R语言bactaxR包中进行聚类分析并生成带有95% 和96%阈值的聚类树状图，操作系统windows或linux均可。代码如下：
   
   
   图5. ANI聚类分析代码
   
   代码逻辑为：

   1. 加载bactaxR程序包；
   2. 读取ANI_pairwise.out文件，命名为ani；
   3. 运行bactaxR程序包中的ANI.dendrogram函数，该函数将两两间ANI值转换为距离矩阵后，利用hclust函数中的非加权组平均法 (UPGMA) 聚类。
      运行结束后，生成ANI聚类图，如图6所示。以95%为阈值，本文中65株待分析菌株和15个已知物种，可划分为34个种，其中45株菌分别归属于19个新种。
      
      
      图6. ANI聚类分析结果
      
      
6. 系统发育分析
   系统发育分析是细菌分类、鉴定和命名的基础，建立在物种间的进化关系上，而不是普遍相似性基础上。本文利用上述Cryobacterium数据为例，构建基于16S rRNA基因序列和全基因组序列的系统发育树。
   6.1

   16S rRNA基因系统发育树构建
   首先，将待分析的16S rRNA基因序列保存为一个fasta格式文件 (16S.fasta)，采用mafft软件对序列进行alignment，命令如下：
   mafft 16S.fasta > 16S_algined.fas
   输出文件为16S_algined.fas。打开MEGA软件，依次点击Phylogeny—Consruct/Test Neighbor-joining Tree，找到16S_algined.fas文件并选中后，选择Nucleotide sequences，点击OK，进入参数选择界面，调整Bootstrap method等参数后 (图7)，点击Compute，生成系统发育树 (图8)。
   
   
   图7. MEGA软件Neighbor-joining系统发育分析参数示例
   
   
   图8. Cryobacterium spp. 基于Neighbor-joining方法的系统发育树 C. mesophilum MSL-15^T设置为外群。图中显示了大于50%的Bootstrap值，标尺代表0.2%核苷酸差异。
   
   
   6.2

   基于基因组的系统发育树构建
   本文以UBCG工具为例，构建基于基因组序列的系统发育树。UBCG可快速抽取细菌92个 (或以下) 共有的保守单拷贝基因序列，然后对每个基因序列和串列序列分别构建基于最大似然法的系统发育树。运行命令如图9所示：
   
   
   图9. 利用UBCG构建基因组系统树命令
   
   代码解释：首先进入UBCG软件所在目录，将~/data/my_genomes文件夹中所有fna后缀的基因组序列转换为bcg格式文件，此处应注意基因组序列文件名称不能含有空格。bcg文件保存至UBCG软件目录下的bcg文件夹，全部转换完成后，执行系统树构建命令，结果保存至~/soft/UBCG/output/Cryo_genome_tree文件夹。基于核心基因组的系统发育树如图10所示。
   
   
   图10. Cryobacterium spp. 基于91个单拷贝持家基因的系统发育树 C. mesophilum MSL-15^T设置为外群。节点处显示出基因支持指数 (GSI) 值，标尺代表2%核苷酸差异。

致谢

本文所使用数据内容由国家自然科学基金青年基金项目 (31600007) 资助，数据来源于Liu et al. (2020)。

参考文献

1. Kim, M., Oh, H. S., Park, S. C. and Chun, J. (2014). Towards a taxonomic coherence between average nucleotide identity and 16S rRNA gene sequence similarity for species demarcation of prokaryotes.Int J Syst Evol Microbiol 64: 346-351.
2. Liu, Q., Song, W. Z., Zhou, Y. G., Dong, X. Z. and Xin, Y. H. (2020). Phenotypic divergence of thermotolerance: Molecular basis and cold adaptive evolution related to intrinsic DNA flexibility of glacier-inhabiting Cryobacterium strains. Environ Microbiol 22: 1409-1420.
3. Richter, M. and Rosselló-Móra, R. (2009). Shifting the genomic gold standard for the prokaryotic species definition. Proc Natl Acad Sci USA 106: 19126-19131.
4. Konstantinidis, K. T. and Tiedje, J. M. (2007). Prokaryotic taxonomy and phylogeny in the genomic era: advancements and challenges ahead. Curr Opin Microbiol 10: 504-9.