洋蔥伯克霍爾德菌復合群部分菌株基于基因組的物種水平重鑒定

周江林，孔娜，張琪，胡明達，周靜，岳俊杰，任洪廣，靳遠，梁龍

軍事科學院 軍事醫學研究院 生物工程研究所，北京100071

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

由于Bcc細菌群體高度的表型和基因型相似性，準確地識別鑒定這些菌種比較困難，而錯誤的病原菌鑒定可能會導致選擇不恰當的治療方法與用藥。常規生化反應等表型檢測方法無法有效鑒定Bcc菌種，特別是對于很多Bcc菌株甚至會報告不同的科屬[3]。16S rRNA 被廣泛用于細菌鑒定，但Bcc 中不同種的16S rRNA基因序列相似性高達98%～100%，無法有效鑒別Bcc菌種[13]。Bcc的recA基因被認為具有較高的分辨能力，能夠有效區分Bcc 群內菌種[4]，并且特異性的recA基因引物對通常被臨床上用于recA基因的擴增與Bcc的鑒定[14]。為了克服單基因分辨力較低等缺陷，基于7個看家基因的多位點序列分析（multilo?cus sequence analysis，MLSA）被開發出來并作為最有效的分類工具用于Bcc菌群[15]。此外，基于細菌全基因組信息的基因組平均核酸一致性（av?erage nucleotide identity，ANI）比較為物種的全面劃分提供了新的解決途徑，并且ANI值已成為現代細菌分類學中物種分界的金標準，95%～96%的ANI值被提出且普遍接受為細菌物種分界線[16-17]。

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

1 材料與方法

1.1 菌株基因組序列數據篩選與下載

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

1.2 系統發生分析

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

對于多基因位點系統發生分析，使用7個管家基因（atpD、gltB、gyrB、recA、lepA、phaC和trpB）進行Bcc菌株系統發生分析。選取新洋蔥伯克霍爾德菌的模式株J2315（=LMG1665）的序列作為搜索請求進行基因位點序列獲取，從PubMLST 數據庫（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 從每一個樣本基因組中獲取其對應的基因序列，并用muscle 工具做多序列對齊[20]，比對結果alignment 用trimAl 工具進行剪裁[21]，移除所有gaps 超過50%的位點。之后將7個基因位點的多序列比對結果用AMAS軟件進行拼接[24]，拼接好的超矩陣按上述方法使用MEGA-X 構建最大似然系統發生樹[22]。構建好的系統發生樹用在線工具iTOL（https://itol.embl.de）進行標注和展示[25]。

1.3 基因組比較

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

2 結果

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%（數據未示）。基于recA序列的系統發生樹如圖1，B.ubonensis（6株）、B.vietnaminensis（6株）、B.con?taminans（2株）的所有菌株都在同一個進化分支上，B.multivorans除了1株外其他6株也均在同一進化分支上且分支長度較短，說明這幾株的分類和鑒定相對比較清晰和可靠。新洋蔥伯克霍爾德（B.cenocepacia）菌株主要分布在2 大枝上，分別對應圖1 中所標示的基因型ⅢA和ⅢB，有個別新洋蔥伯克霍爾德菌零星分布在其他進化分支上，需要對照多位點序列系統發生樹進一步分析。值得注意的是，我們發現有5株細菌的分布位置與其當前鑒定種的進化分支差異較大，在圖1 中黑體加粗表示，其種屬情況需要進一步分析。

2.2 MLSA 進化分析

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

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

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

圖1 62株Bcc菌株基于recA基因序列的系統發生樹

3 討論

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

圖2 62株Bcc菌株基于7個看家基因片段串聯序列的系統發生樹

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

表2 當前種鑒定存在疑問的菌株信息與種鑒定糾正

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