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

刘庆; 辛玉华; 东秀珠

doi:10.21769/BioProtoc.2003932

Improve Research Reproducibility A Bio-protocol resource

Home
Protocols

Peer-reviewed

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

刘刘庆

辛辛玉华

东东秀珠

Published: Apr 22, 2021 DOI: 10.21769/BioProtoc.2003932 Views: 7600

PDF

Ask a question

How to cite

Favorite

Cited by

材料与试剂

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

仪器设备

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

软件和数据库

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

实验步骤

染色体DNA提取和纯化
1.1
用无菌接种环或竹签挑取1个菌落，加入50 µl 5%的Chelex溶液，充分振荡混匀。
1.2
沸水浴10-15 min。
1.3
2,000 × g离心1 min。上清液即为提取的基因组DNA溶液，-20 °C保存备用。
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。
16S rRNA基因序列测定
将16S rRNA基因扩增产物送至测序公司测序。

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-23^T	98.42
MDB2-10	C. breve TMT4-23^T	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-23^T	98.76-99.93
Hh14	C. melibiosiphilum Hh39^T	99.22
Sr8	C. psychrotolerans 0549^T	99.35
RHLS22-1, Sr59	C. soli GCJ02^T	99.57-99.93
Sr54,Hh38,Hz9,Hz16,Sr39,TMS1-13-1	C. levicorallinum Hh34^T	99.65-99.93
TMB1-8, TMB1-7,TMB3-15,TMB3-10,TMN-39-1,TMB3-12,TMB3-1-2	C. zongtaii TMN-42^T	99.658-99.93
TMT1-1, Hb1	C. luteum Hh15^T	99.86
TMT1-66-1	C. flavum Hh8^T	99.93
Hh7, TMS1-20-1,Hh11	C. levicorallinum Hh34^T	100
TMT1-3	C. luteum Hh15^T	100
TMT1-2-2	C. flavum Hh8^T	100
TMN-39-2	C. zongtaii TMN-42^T	100
TMT3-12	C. breve TMT4-23^T	100
TMT1-2-1,Sr47	C. psychrotolerans 0549^T	100

基因组平均核苷酸和氨基酸序列一致性分析
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 0549^T	81.15	新种
MDB2-10	C. breve TMT4-23^T	81.27	新种
MDB2-B	C. breve TMT4-23^T	81.15	新种
MDB1-5	C. breve TMT4-23^T	81.42	新种
MDB1-18-1	C. psychrotolerans 0549^T	81.18	新种
MDB2-A-2	C. psychrotolerans 0549^T	81.20	新种
TMT1-19	C. breve TMT4-23^T	87.00	新种
MDB2-A-1	C. psychrotolerans 0549^T	81.19	新种
Hh14	C. melibiosiphilum Hh39^T	87.43	新种
TMT1-23-1	C. breve TMT4-23^T	87.35	新种
TmT3-29-2	C. breve TMT4-23^T	87.47	新种
Hz7	C. breve TMT4-23^T	86.80	新种
Sr8	C. psychrotolerans 0549^T	98.15	C. psychrotolerans
TMT1-62	C. breve TMT4-23^T	86.76	新种
TmT2-48-2	C. breve TMT4-23^T	87.25	新种
Sr3	C. breve TMT4-23^T	86.86	新种
Hh4	C. breve TMT4-23^T	82.71	新种
MDB1-18-2	C. psychrotolerans 0549^T	81.22	新种
TMT2-16	C. breve TMT4-23^T	86.86	新种
TMT2-14	C. breve TMT4-23^T	86.82	新种
TMT2-17-1	C. breve TMT4-23^T	86.87	新种
TMT2-4	C. breve TMT4-23^T	86.87	新种
TMT4-10	C. breve TMT4-23^T	84.45	新种
TmT2-59	C. breve TMT4-23^T	84.46	新种
TMT1-21	C. breve TMT4-23^T	84.43	新种
TMT2-42-4	C. breve TMT4-23^T	86.88	新种
TMT1-22	C. breve TMT4-23^T	84.34	新种
TMT2-18-3	C. breve TMT4-23^T	86.85	新种
TMT2-18-2	C. breve TMT4-23^T	86.80	新种
TMT2-15-1	C. breve TMT4-23^T	87.00	新种
TMT2-10	C. breve TMT4-23^T	84.47	新种
TMT2-23	C. breve TMT4-23^T	84.21	新种
RHLS22-1	C. zongtaii TMN-42^T	90.73	新种
MDT1-3	C. breve TMT4-23^T	81.64	新种
Sr54	C. aureum Hh31^T	86.47	新种
TMB1-8	C. zongtaii TMN-42^T	96.31	C. zongtaii
TMB1-7	C. zongtaii TMN-42^T	97.87	C. zongtaii
HLT2-28	C. breve TMT4-23^T	82.57	新种
TMT1-51	C. breve TMT4-23^T	89.38	新种
HLT2-23	C. luteum Hh15^T	90.01	新种
Hh38	C. levicorallinum Hh34^T	99.11	C. levicorallinum
Hz9	C. ruanii Sr36^T	86.61	新种
RHLT2-21	C. breve TMT4-23^T	82.56	新种
TMT1-1	C. aureum Hh31^T	86.13	新种
Hb1	C. aureum Hh31^T	85.92	新种
Hz16	C. ruanii Sr36^T	86.83	新种
TMB3-15	C. zongtaii TMN-42^T	97.89	C. zongtaii
TMB3-10	C. zongtaii TMN-42^T	97.86	C. zongtaii
Sr39	C. aureum Hh31^T	87.46	新种
TMT1-66-1	C. flavum Hh8^T	97.33	C. flavum
TMT4-31	C. zongtaii TMN-42^T	97.86	C. zongtaii
TMN-39-1	C. zongtaii TMN-42^T	99.10	C. zongtaii
Sr59	C. zongtaii TMN-42^T	88.16	新种
TMS1-13-1	C. levicorallinum Hh34^T	96.49	C. levicorallinum
TMB3-12	C. zongtaii TMN-42^T	97.90	C. zongtaii
TMB3-1-2	C. zongtaii TMN-42^T	97.87	C. zongtaii
Hh7	C. levicorallinum Hh34^T	97.02	C. levicorallinum
TMS1-20-1	C. levicorallinum Hh34^T	97.09	C. levicorallinum
Hh11	C. levicorallinum Hh34^T	96.94	C. levicorallinum
TMT1-3	C. luteum Hh15^T	99.49	C. luteum
TMT1-2-2	C. flavum Hh8^T	97.46	C. flavum
TMN-39-2	C. zongtaii TMN-42^T	99.07	C. zongtaii
TmT3-12	C. breve TMT4-23^T	99.99	C. breve
TMT1-2-1	C. psychrotolerans 0549^T	99.58	C. psychrotolerans
Sr47	C. psychrotolerans 0549^T	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聚类分析代码

代码逻辑为：

加载bactaxR程序包；
读取ANI_pairwise.out文件，命名为ani；
运行bactaxR程序包中的ANI.dendrogram函数，该函数将两两间ANI值转换为距离矩阵后，利用hclust函数中的非加权组平均法 (UPGMA) 聚类。
运行结束后，生成ANI聚类图，如图6所示。以95%为阈值，本文中65株待分析菌株和15个已知物种，可划分为34个种，其中45株菌分别归属于19个新种。

图6. ANI聚类分析结果

系统发育分析
系统发育分析是细菌分类、鉴定和命名的基础，建立在物种间的进化关系上，而不是普遍相似性基础上。本文利用上述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)。

参考文献

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.
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.
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.
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.