周江林，孔娜，張琪，胡明達(dá)，周靜，岳俊杰，任洪廣，靳遠(yuǎn)，梁龍

軍事科學(xué)院 軍事醫(yī)學(xué)研究院 生物工程研究所，北京100071

洋蔥伯克霍爾德復(fù)合群（Burkholderia cepa?ciancomplex，Bcc）是以康奈爾大學(xué)植物病理學(xué)家Burkholder 命名的一個細(xì)菌復(fù)合群，它的模式菌就是洋蔥伯克霍爾德菌（B.cepacia），由Burkholder 于1949年首次發(fā)現(xiàn)，因其引起洋蔥表皮腐爛而命名為“洋蔥（cepacia）”[1]，最初被稱作洋蔥假單胞菌（Pseudomonas cepacia），隨后在1992 被重新分類為β變形桿菌綱、伯克霍爾德目、伯克霍爾德科、伯克霍爾德屬而稱為洋蔥伯克霍爾德菌[2]。Bcc是一群非葡萄糖發(fā)酵劑的好氧革蘭陰性桿菌，包括至少22種不同的菌種[3]。這些菌種系統(tǒng)發(fā)生關(guān)系密切相關(guān)，通常使用洋蔥伯克霍爾德菌的recA基因序列進(jìn)行種鑒定[4]。這些種當(dāng)中有9個先前被命名為以下伯克霍爾德基因型：B.cepacia（ge?nomovar Ⅰ）、B.multivorans（genomovar Ⅱ）、B.ceno?cepacia（genomovar Ⅲ）、B.stabilis（genomovar Ⅳ）、B.vietnamiensis（genomovarⅤ）、B.dolosa（genomovarⅥ）、B.ambifaria（genomovar Ⅶ）、B.anthina（genom?ovar Ⅷ）、B.pyrrocinia（genomovar Ⅸ），隨后不斷被確認(rèn)并命名為新的物種。Bcc 廣泛存在于自然水源、土壤和其他潮濕的環(huán)境中，它們可以在這些環(huán)境存活達(dá)數(shù)月之久[5-6]，現(xiàn)有證據(jù)顯示環(huán)境可能是新型Bcc 感染的源頭[7]。自20世紀(jì)50年代首次報告人感染洋蔥伯克霍爾德菌患病以來[1]，現(xiàn)在已是人群中流行程度僅次于ESKAPE病原體（腸球菌、金黃色葡萄球菌、克雷伯菌、不動桿菌、銅綠假單胞菌、腸桿菌）的致病菌。Bcc 頻繁地從囊性纖維化患者的痰液當(dāng)中分離出來，由于其抗生素耐藥性和治療困難，往往會加速囊性纖維化患者肺功能衰竭并增大死亡率[8-9]。同時，Bcc 也會導(dǎo)致慢性肉芽腫病患者的吞噬細(xì)胞無法產(chǎn)生活性氧，從而造成致命性疾病[10-11]。近年來不斷增加的Bcc 導(dǎo)致的院內(nèi)感染報告使得這些細(xì)菌被認(rèn)為是緊急醫(yī)院病原體[12]。

由于Bcc細(xì)菌群體高度的表型和基因型相似性，準(zhǔn)確地識別鑒定這些菌種比較困難，而錯誤的病原菌鑒定可能會導(dǎo)致選擇不恰當(dāng)?shù)闹委煼椒ㄅc用藥。常規(guī)生化反應(yīng)等表型檢測方法無法有效鑒定Bcc菌種，特別是對于很多Bcc菌株甚至?xí)蟾娌煌目茖賉3]。16S rRNA 被廣泛用于細(xì)菌鑒定，但Bcc 中不同種的16S rRNA基因序列相似性高達(dá)98%～100%，無法有效鑒別Bcc菌種[13]。Bcc的recA基因被認(rèn)為具有較高的分辨能力，能夠有效區(qū)分Bcc 群內(nèi)菌種[4]，并且特異性的recA基因引物對通常被臨床上用于recA基因的擴(kuò)增與Bcc的鑒定[14]。為了克服單基因分辨力較低等缺陷，基于7個看家基因的多位點序列分析（multilo?cus sequence analysis，MLSA）被開發(fā)出來并作為最有效的分類工具用于Bcc菌群[15]。此外，基于細(xì)菌全基因組信息的基因組平均核酸一致性（av?erage nucleotide identity，ANI）比較為物種的全面劃分提供了新的解決途徑，并且ANI值已成為現(xiàn)代細(xì)菌分類學(xué)中物種分界的金標(biāo)準(zhǔn)，95%～96%的ANI值被提出且普遍接受為細(xì)菌物種分界線[16-17]。

隨著全基因組測序技術(shù)的發(fā)展和測序成本的降低，越來越多的Bcc基因組被測序并完整拼接，為生物信息學(xué)分析提供了良好的數(shù)據(jù)基礎(chǔ)。本研究使用了所有目前可用的Bcc細(xì)菌全基因組序列數(shù)據(jù)，基于多種方法進(jìn)行了系統(tǒng)發(fā)生和分類現(xiàn)狀分析，初步糾正了一些Bcc菌株當(dāng)前錯誤的種鑒定。我們的研究結(jié)果提示許多GenBank 公共數(shù)據(jù)庫中的Bcc菌株可能存在錯誤的種鑒定，需要進(jìn)行深入分析和糾正。

1 材料與方法

1.1 菌株基因組序列數(shù)據(jù)篩選與下載

從GenBank基因組數(shù)據(jù)庫（https://www.ncbi.nlm.nih.gov）選取72株注釋為洋蔥伯克霍爾德復(fù)合群內(nèi)菌種的全基因組序列數(shù)據(jù)（截至2019-04-25）進(jìn)行分析，其中模式株或代表株基因組9個。這些基因組按照Parks 等所述方法進(jìn)行去重[18]，惟一例外是本研究中去重是基于基因組序列ANI 而不是原文中bac120 比對的兩兩之間的平均氨基酸一致性（average amino acid identity，AAI）。用checkM 工具評估下載的基因組數(shù)據(jù)的完整度（completeness）和污染程度（contamination）[19]，定義基因組質(zhì)量（quality）為完整度-5×污染程度[18]，去除所有基因組數(shù)據(jù)質(zhì)量值低于50的菌株，去重之后共得到62株細(xì)菌基因組進(jìn)行后續(xù)分析。選取復(fù)合群外的菌株B.pseudomalleiK96243、B.glumaeLMG2196和B.oklahomensisC6786 作 為 外群。

1.2 系統(tǒng)發(fā)生分析

對于recA單基因系統(tǒng)發(fā)生分析，選取新洋蔥伯克霍爾德菌的模式株J2315（=LMG1665）的recA基因全長序列（NCBI 核酸序列訪問號為ALK16523.1）作為搜索請求，用blastn 獲取所有樣本Bcc菌株基因組中的全長recA基因，所有序列使用muscle 工具進(jìn)行多序列對齊[20]，比對結(jié)果用trimAl 工具進(jìn)行剪裁[21]，移除所有g(shù)aps 超過50%的位點，比對結(jié)果共有1071個核苷酸位點。使用MEGA-X 進(jìn)行1000 次bootstrap 構(gòu)建最大似然系統(tǒng)發(fā)生樹[22]，模型方法為General Time Reversible model，位點替換速率模型為Gamma Distributed With Invariant Sites（G+I）。

對于多基因位點系統(tǒng)發(fā)生分析，使用7個管家基因（atpD、gltB、gyrB、recA、lepA、phaC和trpB）進(jìn)行Bcc菌株系統(tǒng)發(fā)生分析。選取新洋蔥伯克霍爾德菌的模式株J2315（=LMG1665）的序列作為搜索請求進(jìn)行基因位點序列獲取，從PubMLST 數(shù)據(jù)庫（https://pubmlst.org/bcc/）下載J2315株的7個位點序列[23]，分別為atpD（443 bp）、gltB（400 bp）、gyrB（454 bp）、recA（393 bp）、lepA（395 bp）、phaC（385 bp）、trpB（301 bp）。分別以下載的7個基因位點序列作為請求，用blastn 從每一個樣本基因組中獲取其對應(yīng)的基因序列，并用muscle 工具做多序列對齊[20]，比對結(jié)果alignment 用trimAl 工具進(jìn)行剪裁[21]，移除所有g(shù)aps 超過50%的位點。之后將7個基因位點的多序列比對結(jié)果用AMAS軟件進(jìn)行拼接[24]，拼接好的超矩陣按上述方法使用MEGA-X 構(gòu)建最大似然系統(tǒng)發(fā)生樹[22]。構(gòu)建好的系統(tǒng)發(fā)生樹用在線工具iTOL（https://itol.embl.de）進(jìn)行標(biāo)注和展示[25]。

1.3 基因組比較

用fastANI 工具計算所有樣本菌株任意2個基因組之間的ANI值，以及Bcc菌株基因組與外群菌株基因組直接的ANI值[26]。對于同一個種的某一群基因組，若其中任意2個基因組之間的ANI值大于或等于99.95%，則被認(rèn)為是高度相似而冗余的菌株，從中挑選一個基因組作為該群的代表株，優(yōu)先選擇參考株（reference strain）或代表株（representative strain），若不存在時隨機(jī)挑選一個完整基因組作為該群的代表株。

2 結(jié)果

2.1 recA基因序列分析

62株Bcc的recA基因平均長度為1072.25 bp，最長的recA基因為2株B.ambifaria菌的1080 bp，最短的recA基因為4株B.cenocepacia菌的1070 bp；序列相似度為94.118%～100%，平均95.898%（數(shù)據(jù)未示）。基于recA序列的系統(tǒng)發(fā)生樹如圖1，B.ubonensis（6株）、B.vietnaminensis（6株）、B.con?taminans（2株）的所有菌株都在同一個進(jìn)化分支上，B.multivorans除了1株外其他6株也均在同一進(jìn)化分支上且分支長度較短，說明這幾株的分類和鑒定相對比較清晰和可靠。新洋蔥伯克霍爾德（B.cenocepacia）菌株主要分布在2 大枝上，分別對應(yīng)圖1 中所標(biāo)示的基因型ⅢA和ⅢB，有個別新洋蔥伯克霍爾德菌零星分布在其他進(jìn)化分支上，需要對照多位點序列系統(tǒng)發(fā)生樹進(jìn)一步分析。值得注意的是，我們發(fā)現(xiàn)有5株細(xì)菌的分布位置與其當(dāng)前鑒定種的進(jìn)化分支差異較大，在圖1 中黑體加粗表示，其種屬情況需要進(jìn)一步分析。

2.2 MLSA 進(jìn)化分析

為了進(jìn)一步分析Bcc菌株的系統(tǒng)發(fā)育和分類情況，我們做了MLSA 分析。62株Bcc菌株都成功找到了所有7個看家基因的同源保守位點，多序列對齊后串聯(lián)長度為2771 bp[atpD（443 bp）、gltB（400 bp）、gyrB（454 bp）、recA（393 bp）、lepA（395 bp）、phaC（385 bp）、trpB（301 bp）]，基于串聯(lián)等位基因序列的最大似然系統(tǒng)發(fā)生樹如圖2。首先，recA基因進(jìn)化樹上處于單一分支的幾個種在圖2的串聯(lián)看家基因片段進(jìn)化樹上依然處于單獨的進(jìn)化分支，因此這幾個Bcc菌種可能是單系發(fā)生的（monophytic），并且其菌株的鑒定和分類情況相對較好。其次，圖2 再次確認(rèn)了新洋蔥伯克霍爾德菌至少分2個基因型ⅢA、ⅢB，這與以往的研究結(jié)果一致。第三，對于recA基因進(jìn)化樹上5個分布位置異常的菌株（圖1），我們發(fā)現(xiàn)在基于MLSA 進(jìn)化樹上分布位置也顯示與其注釋菌種所在進(jìn)化分支存在較大差異（圖2），且這5株菌在2個進(jìn)化樹上的分布位置相互吻合，再次提示我們這幾株菌的種屬鑒定可能存在問題。

2.3 基因組平均核酸一致性計算與比較結(jié)果

為了進(jìn)一步確認(rèn)這5株菌的種屬信息，我們計算比較了其各自與注釋物種、本文進(jìn)化樹上分支所在菌種的參考株或代表株的ANI值，結(jié)果見表1?？梢钥吹剑闐DS 7H-2的當(dāng)前鑒定菌種是B.cepacia，但是其在系統(tǒng)發(fā)生樹上處于B.ceno?cepacia基因型ⅢA 所在的進(jìn)化分支，與B.cenocepa?cia模式株J2315基因組的ANI值為98.949%，遠(yuǎn)遠(yuǎn)大于其與當(dāng)前鑒定菌種B.cepacia模式株ATCC 25416的ANI值91.9433%，也顯著大于細(xì)菌物種分界的ANI 閾值95%[17]。結(jié)合ANI值比較和進(jìn)化樹分析，我們有足夠理由認(rèn)為菌株DDS 7H-2 實際上是一株新洋蔥伯克霍爾德菌（B.cenocepacia），且其基因型為ⅢA。同樣的分析，我們可以糾正菌株FDAARGOS 496、LO6、DWS 37E-2、DDS 22E-1的菌種鑒定，詳細(xì)信息見表2。

圖1 62株Bcc菌株基于recA基因序列的系統(tǒng)發(fā)生樹

3 討論

準(zhǔn)確鑒定細(xì)菌種屬信息非常重要，尤其是當(dāng)這些細(xì)菌會引起人或動物感染發(fā)病并且需要選擇對應(yīng)治療方案的時候，準(zhǔn)確可靠的鑒定信息對于應(yīng)對傳染病暴發(fā)和疾病防控至關(guān)重要。洋蔥伯克霍爾德菌復(fù)合群由于其廣泛的分布環(huán)境和經(jīng)常性感染人群，特別是在肺纖維化患者間的相互傳播給人們帶來較大的生命財產(chǎn)損失，引起了研究者的關(guān)注。該復(fù)合群已知菌種至少有22種，各個種的外形和生化特性比較相似，16S rRNA 序列相似度也非常高，臨床上常規(guī)生化方法和廣泛使用的16S rRNA基因分類方法不能準(zhǔn)確有效地分辨各個菌種[3]，甚至?xí)a(chǎn)生錯誤的鑒定注釋信息，從而對后續(xù)治療方案選擇帶來誤導(dǎo)；另一方面，現(xiàn)有菌株的測序結(jié)果往往會上傳到公共序列數(shù)據(jù)庫，這在方便研究人員基于大規(guī)模序列數(shù)據(jù)集綜合分析的同時，也可能會由于之前不可靠方法鑒定的菌株注釋信息給研究者的分析結(jié)果帶來偏差甚至錯誤。

圖2 62株Bcc菌株基于7個看家基因片段串聯(lián)序列的系統(tǒng)發(fā)生樹

表1 當(dāng)前種鑒定存在疑問的菌株與Bcc相關(guān)菌種基因組ANI值比較（%）

表2 當(dāng)前種鑒定存在疑問的菌株信息與種鑒定糾正

本研究中我們首先通過構(gòu)建recA單基因系統(tǒng)發(fā)生樹，發(fā)現(xiàn)有5株細(xì)菌在進(jìn)化樹上分布位置與其本身注釋菌種所在分支不同并且進(jìn)化距離較大，而與其他種的細(xì)菌分布在一個進(jìn)化分支且分支長度較小，這提示我們這幾株細(xì)菌在公共數(shù)據(jù)庫里的分類注釋可能存在異常。進(jìn)一步，我們用7個看家基因片段串聯(lián)進(jìn)化分析和細(xì)菌全基因組ANI值比較等多種方法，綜合分析了當(dāng)前Bcc細(xì)菌的分類情況，結(jié)果顯示這5株菌當(dāng)前的菌種注釋存在錯誤，并且根據(jù)多種方法更新糾正了這5株Bcc菌株的種鑒定信息。我們的分析和結(jié)果表明細(xì)菌全基因組信息的ANI值比較與MLSA 合作分析可以良好地鑒定Bcc菌種并避免可能的種鑒定錯誤，臨床上應(yīng)該以這2種方法的鑒定結(jié)果作為標(biāo)準(zhǔn)，避免僅僅依靠單個基因而對Bcc 復(fù)合群的細(xì)菌做出種鑒定；同時提示當(dāng)前公共數(shù)據(jù)庫中可能存在許多種鑒定錯誤的Bcc菌種，這些序列注釋信息需要研究人員進(jìn)一步分析和糾正。