鄭宏源，閆琳，楊超,3，武雅蓉，秦婧靚，郝彤宇，楊大進(jìn)，郭云昌，裴曉燕，趙彤言，崔玉軍

研究報(bào)告

溶藻弧菌群體基因組學(xué)研究

鄭宏源1，閆琳2，楊超1,3，武雅蓉1，秦婧靚1，郝彤宇1，楊大進(jìn)2，郭云昌2，裴曉燕2，趙彤言1，崔玉軍1

1. 軍事科學(xué)院軍事醫(yī)學(xué)研究院微生物流行病研究所，北京 100071 2. 國(guó)家食品安全風(fēng)險(xiǎn)評(píng)估中心，北京 100022 3. 深圳市疾病預(yù)防控制中心，深圳 518055

溶藻弧菌()是一種能夠?qū)θ祟愐约棒~、蝦、貝類等水產(chǎn)品致病的弧菌，給人類健康帶來威脅，也給水產(chǎn)養(yǎng)殖業(yè)造成巨大的經(jīng)濟(jì)損失。目前該物種基于全基因組的遺傳多樣性和重要遺傳元件研究報(bào)道較少。本研究對(duì)采集自全國(guó)4個(gè)省份的68株溶藻弧菌進(jìn)行高通量測(cè)序，獲得全基因組序列，并結(jié)合113株公開發(fā)表的全球序列數(shù)據(jù)，利用fineSTRUCTURE軟件、VFDB毒力因子庫和CARD、ResFinder耐藥數(shù)據(jù)庫，對(duì)溶藻弧菌的種群結(jié)構(gòu)和毒力、耐藥因子分別進(jìn)行解析。結(jié)果表明：溶藻弧菌可分為譜系1和譜系2。兩個(gè)譜系在美洲和亞洲均有分布，但歐洲僅分離到譜系1菌株；共鑒定發(fā)現(xiàn)12個(gè)克隆群，其中一個(gè)克隆群內(nèi)菌株存在跨洋傳播現(xiàn)象。該物種攜帶、、等多種不同功能的毒力因子；毒力因子在兩個(gè)譜系間的分布無特異性，但存在地域間差異：其中歐洲菌株攜帶、、和的比率低于其他地區(qū)，而基因的攜帶率則明顯高于其他地區(qū)，我國(guó)廣西菌株中基因攜帶率低于其他省份，且不攜帶。多個(gè)基因組攜帶、、、、、和等與多種抗生素耐受相關(guān)的基因，其中和基因在譜系2中的出現(xiàn)頻率遠(yuǎn)高于譜系1；此外，基因在亞洲菌株中的攜帶率高于美洲和歐洲地區(qū)，而在我國(guó)四川菌株中的攜帶率則低于其他省份。在5個(gè)基因組中(VA24、VA28、2014V-1011、ZJ-T和Vb1833)觀察到質(zhì)粒或ICE等攜帶多種耐藥基因的大片段。本研究通過群體基因組學(xué)的研究方法，揭示了溶藻弧菌的種群結(jié)構(gòu)組成和毒力、耐藥相關(guān)元件的分布，為進(jìn)一步了解溶藻弧菌的遺傳特征和致病機(jī)制提供必要基礎(chǔ)，為該病原的監(jiān)測(cè)、預(yù)防和控制工作提供科學(xué)支撐。

溶藻弧菌；全基因組測(cè)序；種群結(jié)構(gòu)；毒力因子；耐藥基因

溶藻弧菌()，是一種廣泛存在于海洋和河口，可以從病人體內(nèi)、海水、水產(chǎn)品等多種不同環(huán)境中分離得到的革蘭氏陰性桿菌[1,2]。溶藻弧菌感染人體后可以引起皮膚感染、外耳道感染等多種疾病，但與霍亂弧菌()和副溶血弧菌()等食源性弧菌不同，臨床上報(bào)導(dǎo)的溶藻弧菌感染病例均有海水暴露史，因此其主要感染途徑是與海水的接觸[3,4]。此外，溶藻弧菌可引起魚、蝦、貝類等水產(chǎn)品患病與大量死亡，帶來巨大的經(jīng)濟(jì)損失[2,5,6]。美國(guó)COVIS (Cholera and OtherIllness Surveillance)和FoodNet (Foodborne Diseases Active Surveillance Network)兩個(gè)監(jiān)測(cè)數(shù)據(jù)庫的數(shù)據(jù)表明：自2007年以來，溶藻弧菌已經(jīng)躍升為第二常見的致病弧菌[7,8]，為公共衛(wèi)生及漁業(yè)經(jīng)濟(jì)帶來巨大壓力。深入研究溶藻弧菌遺傳多樣性和病原相關(guān)特征，將有助于該病原的精準(zhǔn)檢測(cè)和致病機(jī)制研究，對(duì)其監(jiān)測(cè)和防控工作具有重要意義。

多年來，已有多種傳統(tǒng)分子生物學(xué)方法應(yīng)用于溶藻弧菌分型研究，例如核糖體分型、隨機(jī)擴(kuò)增多態(tài)性DNA (random amplified polymorphic DNA, RAPD)、腸桿菌基因間重復(fù)一致序列PCR (entero-bacterial repetitive intergenic consensus PCR, ERIC- PCR)、低頻限制性位點(diǎn)PCR (infrequent restriction site PCR, IRS-PCR)和脈沖場(chǎng)凝膠電泳(pulsed field gel electrophoresis, PFGE)等[9~12]。這些方法在溯源性研究等方面發(fā)揮了一定的作用，具有操作簡(jiǎn)便、獲得結(jié)果較快等優(yōu)勢(shì)，但缺點(diǎn)是分辨率較低，實(shí)驗(yàn)結(jié)果不穩(wěn)定，而且對(duì)實(shí)驗(yàn)條件要求較高。PFGE曾被認(rèn)為是暴發(fā)溯源的金標(biāo)準(zhǔn)[13,14]，但是這種分型方法需要嚴(yán)苛的實(shí)驗(yàn)條件，實(shí)驗(yàn)過程中操作導(dǎo)致的細(xì)微差別則可能得到不同的分型結(jié)果，不利于實(shí)驗(yàn)結(jié)果在不同實(shí)驗(yàn)室間的交流。目前尚無多位點(diǎn)可變數(shù)量串聯(lián)重復(fù)序列分析(multilocus variable-number tandem repeat analysis, MLVA)和多位點(diǎn)序列分型(multilocus sequence typing, MLST)等多靶標(biāo)分子分型方法應(yīng)用于溶藻弧菌的研究報(bào)導(dǎo)，而基于全基因組層面的遺傳多樣性研究和種群結(jié)構(gòu)鑒定的研究報(bào)道較少。另一方面，當(dāng)前對(duì)于溶藻弧菌的毒力和耐藥相關(guān)基因的研究，是利用已知其他弧菌物種的毒力和耐藥基因特異性引物，通過PCR方法在部分菌株中檢驗(yàn)相應(yīng)基因是否存在。研究表明溶藻弧菌攜帶了一些副溶血弧菌等物種常見毒力因子，也表現(xiàn)出多種臨床常見抗生素的耐受性，同時(shí)鑒定到部分菌株攜帶有耐藥基因[1,15~21]。然而，溶藻弧菌全基因組層面毒力和耐藥相關(guān)遺傳因子分布的認(rèn)識(shí)仍不夠深入。

全基因組測(cè)序(whole genome sequencing, WGS)在細(xì)菌遺傳多樣性研究中有著最高的分辨率，而且已廣泛運(yùn)用在微生物溯源、傳播及種群結(jié)構(gòu)鑒定等各個(gè)方面[22~25]。本研究通過全基因組測(cè)序，利用fineSTRUCTURE方法[26]鑒定了溶藻弧菌的種群結(jié)構(gòu)，并在此基礎(chǔ)上研究了毒力、耐藥相關(guān)因子及其在物種內(nèi)不同菌株間的分布。

1 材料與方法

1.1 菌株來源

本研究共收集181株溶藻弧菌的全基因組序列，包括68株新測(cè)序的序列和113株公共基因組序列。新測(cè)序的菌株采樣于國(guó)內(nèi)四個(gè)省份(山東、湖北、四川、廣西)的水產(chǎn)品、海水、水底沉積物等不同環(huán)境，采樣時(shí)間均為2014年。公共序列下載于NCBI數(shù)據(jù)庫(National Center for Biotechnology Information Search database，https://www.ncbi.nlm.nih.gov/)，包括截止到2020年10月15日的所有組裝結(jié)果和SRA數(shù)據(jù)(46+67)，這些菌株來自美洲、歐洲和亞洲，分離時(shí)間為1969?2020年間。

1.2 全基因組測(cè)序

使用細(xì)菌基因組DNA提取試劑盒(QIAGEN UltraClean?Microbial DNA Isolation kit)，參照試劑盒說明書提取全基因組DNA。經(jīng)檢測(cè)質(zhì)量合格后送測(cè)序公司(華大基因，深圳)進(jìn)行細(xì)菌全基因組測(cè)序，測(cè)序平臺(tái)為Illumina NovaSeq 6000，測(cè)序文庫大小為350 bp，測(cè)序reads為雙端150 bp，平均測(cè)序深度134X (94~188X)。

1.3 質(zhì)控及組裝

使用軟件Trimmomatic v0.38[27]將測(cè)序結(jié)果中低質(zhì)量的數(shù)據(jù)進(jìn)行過濾，然后利用軟件shovill (https:// github.com/tseemann/shovill)對(duì)過濾后的reads 進(jìn)行組裝。組裝后的基因組草圖序列已上傳至NCBI，BioProject編號(hào)為PRJNA714277。

1.4 變異檢測(cè)和基因預(yù)測(cè)

單核苷酸多態(tài)性(single nucleotide polymorphism, SNP)位點(diǎn)的鑒定參照我們已發(fā)表的方法[28]，以溶藻弧菌FA2(GCA_011801435.1)基因組完成圖為參考序列，使用MUMmer軟件包[29]將181株溶藻弧菌基因組序列分別比對(duì)到參考菌株基因組，獲得核心基因組(所有菌株中均能找到對(duì)應(yīng)堿基狀態(tài))SNP矩陣。為保證SNP鑒定結(jié)果的可靠性，位于重復(fù)區(qū)和測(cè)序質(zhì)量值低的SNP 未納入后續(xù)分析。采用Prokka軟件[30]對(duì)各個(gè)菌株基因組序列進(jìn)行基因預(yù)測(cè)，使用Roary軟件[31]判斷基因獲得缺失。

1.5 系統(tǒng)發(fā)育分析

使用TreeBest軟件(http://treesoft.sourceforge. net/treebest.shtml)對(duì)SNP數(shù)據(jù)集構(gòu)建鄰接(Neighbour- Joining, NJ)樹，bootstrap值設(shè)為100，使用在線工具iTOL[32]進(jìn)行可視化展示。

1.6 種群結(jié)構(gòu)鑒定

首先將SNP矩陣轉(zhuǎn)換為FASTA格式，然后計(jì)算兩兩菌株間的SNP距離，為減小克隆結(jié)構(gòu)的干擾，以10,000個(gè)SNP為閾值去除克隆信號(hào)，并從每個(gè)克隆群中隨機(jī)選取一株菌作為該群的代表株，經(jīng)過篩選最終得到129個(gè)代表株。將129株菌的SNP信息輸入種群結(jié)構(gòu)分析軟件fineSTRUCTURE[26]，設(shè)定默認(rèn)參數(shù)，計(jì)算菌株間共祖率(coancestry)，得到共祖率矩陣，再利用馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo, MCMC)算法實(shí)施矩陣的聚類分析，進(jìn)而推斷種群結(jié)構(gòu)。

1.7 毒力因子鑒定和耐藥基因預(yù)測(cè)

使用BLASTn將所有菌株基因組數(shù)據(jù)與VFDB (Virulence Factor Database)[33]數(shù)據(jù)庫核心庫(setA，僅包含經(jīng)試驗(yàn)驗(yàn)證過的毒力基因)進(jìn)行比對(duì)，鑒定本研究中數(shù)據(jù)集的毒力因子；利用RGI v5.1.1[34]和ResFinder[35]兩種軟件分別基于CARD (Comprehen-sive Antibiotic Resistance Database, https://card. mcmaster.ca/)耐藥基因庫和ResFinder (https://bitbucket. org/genomicepidemiology/resfinder_db/src/master/)耐藥庫鑒定整個(gè)數(shù)據(jù)集的耐藥基因。毒力因子和耐藥基因鑒定結(jié)果中，參考我們已發(fā)表的相關(guān)工作[36]，coverage和identity值均大于70%認(rèn)為該基因存在，否則為缺失。

2 結(jié)果與分析

2.1 溶藻弧菌基因組特點(diǎn)

本研究在我國(guó)4個(gè)省份(山東、湖北、四川、廣西)共分離了68株溶藻弧菌，加上公共數(shù)據(jù)113株，共181株菌基因組數(shù)據(jù)納入分析。新測(cè)序的68株溶藻弧菌基因組序列，經(jīng)過濾后組裝，平均基因組大小為5.1 Mb (4.9~5.4 Mb)，平均GC含量為44.6% (44.4%~44.9%)。整個(gè)數(shù)據(jù)集首先進(jìn)行了泛基因組分析，共得到22,439個(gè)基因，其中核心基因2995個(gè)，附屬基因19,444個(gè)。根據(jù)圖1可以推測(cè)，隨著菌株數(shù)量的增加，整個(gè)物種的泛基因數(shù)量呈上升趨勢(shì)，而核心基因組則趨于穩(wěn)定，說明整個(gè)物種的附屬基因在逐漸增加，即溶藻弧菌的基因組是“開放”的，可以通過水平基因轉(zhuǎn)移等方式不斷獲得新的基因。

圖1 溶藻弧菌基因數(shù)量與基因組數(shù)量的關(guān)系

橫軸為納入分析的基因組數(shù)量，縱軸為基于數(shù)據(jù)集計(jì)算出的基因總數(shù)與核心基因數(shù)量。

2.2 溶藻弧菌種群結(jié)構(gòu)

通過變異檢測(cè)，181株溶藻弧菌共鑒定出335,733個(gè)SNP。根據(jù)鑒定到的SNP，我們構(gòu)建了整個(gè)數(shù)據(jù)集的系統(tǒng)發(fā)育樹(圖2A)，結(jié)果顯示溶藻弧菌的系統(tǒng)發(fā)育樹和副溶血弧菌[24,37]等物種類似，呈現(xiàn)“輪盤狀”結(jié)構(gòu)(圖2A)，說明菌株間水平轉(zhuǎn)移信號(hào)已將垂直遺傳信號(hào)打亂，該物種屬于高頻重組弧菌。使用fineSTRUCTURE軟件對(duì)溶藻弧菌的種群結(jié)構(gòu)做進(jìn)一步解析，最終將溶藻弧菌分成兩個(gè)主要譜系(Lineage)：譜系1和2 (圖2B)。從圖中可發(fā)現(xiàn)，在對(duì)角線上的少量種群中仍呈現(xiàn)明顯的克隆信號(hào)，在我們已發(fā)表的相關(guān)研究中[24]提到，可通過進(jìn)一步去除克隆群菌株以移除克隆信號(hào)干擾。本研究由于菌株數(shù)量相對(duì)較少，且進(jìn)一步去除克隆信號(hào)后得到的結(jié)果仍分為同樣兩個(gè)譜系，故只展示第一次去除克隆信號(hào)的結(jié)果。

分析結(jié)果(表1)顯示：譜系1中的美洲和中國(guó)菌株占比相當(dāng)，而譜系2的中國(guó)菌株占比較高；譜系1菌株與譜系2菌株相比分布更為廣泛，歐洲來源的17株菌僅出現(xiàn)在譜系1，提示譜系1在全球的擴(kuò)散程度更高。

圖2 溶藻弧菌系統(tǒng)發(fā)育樹與種群結(jié)構(gòu)

A：181株菌基于核心基因組SNP構(gòu)建的系統(tǒng)發(fā)育樹。黑色、紅色樹枝分別代表譜系1和譜系2，樹枝頂部紫色圓點(diǎn)代表鑒定到的12個(gè)克隆群(clonal groups, CGs)，右下角紅色虛線框代表的是發(fā)生跨洋傳播的兩株菌，內(nèi)環(huán)代表不同的采樣地點(diǎn)，外環(huán)代表不同采樣來源。B：129株代表菌株基于fineSTRUCTURE軟件的種群鑒定結(jié)果。列代表供體菌株(donor)，行代表受體菌株(recipient)，方格顏色表示供體向受體導(dǎo)入的序列片段數(shù)量，黑色線代表種群間的界線。

表1 不同譜系中菌株分離地點(diǎn)統(tǒng)計(jì)

*代表兩個(gè)譜系中不同采樣地點(diǎn)的菌株數(shù)量，#代表相應(yīng)地點(diǎn)采樣菌株在譜系內(nèi)的占比。

參考我們已發(fā)表的研究方法[25]，以SNP距離小于10為界鑒定出了近期傳播或暴發(fā)的12個(gè)克隆群(CGs，圖2A)，共涉及34株菌。僅在其中一個(gè)克隆群中(圖2A中虛線框)，發(fā)現(xiàn)來自不同國(guó)家的菌株。相關(guān)菌株分別來自英國(guó)和日本(分離時(shí)間不詳)，這兩株菌基因組間 SNP遺傳距離為0，表明溶藻弧菌在上述兩個(gè)國(guó)家之間發(fā)生了跨地域傳播。背景信息顯示這兩株菌分別分離自蛤類和魚類，提示跨國(guó)傳播可能是由于水產(chǎn)品國(guó)際貿(mào)易導(dǎo)致的。本研究中溶藻弧菌跨洋傳播現(xiàn)象的確定對(duì)我們將來的防治思路和策略有指導(dǎo)意義。

2.3 毒力因子的鑒定及分布

181個(gè)基因組中共鑒定出9類37個(gè)毒力因子(表2)，其中鞭毛(flagella)與三型分泌系統(tǒng)(T3SS)相關(guān)基因最多，分別為14個(gè)和16個(gè)，占毒力因子總數(shù)的81%。鞭毛作為一種附著因子(adhesin factor)，可介導(dǎo)溶藻弧菌在宿主表面的粘附；T3SS最初是在耶爾森氏菌屬(spp)中發(fā)現(xiàn)的[38]，而弧菌屬(spp.)中最早發(fā)現(xiàn)T3SS是在副溶血弧菌(血清型O3:K6, strain RIMD2210633)的兩條染色體上，包含兩套T3SS系統(tǒng)，分別命名為T3SS1和T3SS2[39],其中T3SS1在所有副溶血弧菌基因組中都存在，主要功能是誘導(dǎo)細(xì)胞自噬、細(xì)胞圓化和細(xì)胞裂解[40,41]。此外，我們還發(fā)現(xiàn)所有的溶藻弧菌(100%)都攜帶溶血素基因。是5類典型溶血素家族(TDH、HlyA、TLH、δ-VPH、HLX)[42]中的一種，實(shí)驗(yàn)表明這種基因在宿主中介導(dǎo)了凋亡、膜溶解、壞死等多種毒性表現(xiàn)[42,43]。98%的菌株攜帶膜外蛋白基因，這是溶藻弧菌的主要孔蛋白，也是一種重要的毒力因子，有研究已經(jīng)證實(shí)在溶藻弧菌的鐵代謝平衡中起到了重要作用，而且也是溶藻弧菌重要的疫苗研發(fā)靶標(biāo)[44,45]。我們還分別在87%、26%和9%的溶藻弧菌中發(fā)現(xiàn)攜帶多價(jià)粘附分子基因、免疫原性脂蛋白基因、和熱休克蛋白基因。另有個(gè)別菌株表達(dá)莢膜基因和T6SS亞單位。

表2 溶藻弧菌毒力因子及相關(guān)功能

括號(hào)內(nèi)百分比代表整個(gè)數(shù)據(jù)集中攜帶此毒力基因的菌株的百分比。

在鑒定出的37個(gè)毒力因子中，其中23個(gè)在大于90%的菌株中攜帶, 8個(gè)在10%~90%的菌株中攜帶，其余6個(gè)僅在少數(shù)(<10%)菌株中攜帶。統(tǒng)計(jì)這8個(gè)分布差異性較大的毒力因子在不同譜系中、不同國(guó)家和地區(qū)間、中國(guó)不同省份的分布情況(表3)，結(jié)合毒力因子在系統(tǒng)發(fā)育樹上的分布(圖3)，未發(fā)現(xiàn)譜系1或譜系2的群特異性毒力因子；歐洲菌株中、、和的攜帶率低于其他國(guó)家和地區(qū)，而基因的攜帶率較高；除中國(guó)外的亞洲其他地區(qū)的菌株均不攜帶基因；此外，我國(guó)廣西的菌株基因的攜帶率低于其他省份，且均不攜帶基因。因本研究中的菌株采樣量較小，相關(guān)遺傳因子的地域分布特征仍需將來通過擴(kuò)大菌株集進(jìn)行驗(yàn)證。

2.4 耐藥基因的預(yù)測(cè)及分布

RGI和ResFinder均鑒定到32個(gè)耐藥基因，經(jīng)同源比對(duì)、去除冗余后匯總為31個(gè)耐藥基因，分別介導(dǎo)溶藻弧菌對(duì)β內(nèi)酰胺酶類、喹諾酮類、四環(huán)素類等多種常見抗生素的耐受(表4)。我們發(fā)現(xiàn)幾乎所有菌株(>95%)普遍攜帶的耐藥基因有：介導(dǎo)β內(nèi)酰胺酶類耐受的，介導(dǎo)四環(huán)素耐受的和，介導(dǎo)喹諾酮類耐受的，以及調(diào)控耐藥的基因、。

此外，我們發(fā)現(xiàn)介導(dǎo)四環(huán)素耐受的基因在兩個(gè)譜系中的分布存在顯著差異(圖3)，譜系2和譜系1的菌株攜帶這種基因的比率分別為100%和12%。是一種轉(zhuǎn)錄調(diào)節(jié)因子，與四環(huán)素耐藥相關(guān)，而且另一個(gè)四環(huán)素耐藥基因耐藥功能的正常發(fā)揮，依賴于。我們還在9株溶藻弧菌中首次發(fā)現(xiàn)了一種和磷霉素耐受相關(guān)的基因，有趣的是其中8株菌(88.9%)屬于譜系2。此外，經(jīng)統(tǒng)計(jì)發(fā)現(xiàn)，基因在中國(guó)和亞洲其他地區(qū)的菌株中的攜帶率(分別為45%和44%)高于美洲和歐洲菌株的攜帶率(分別為24%和6%)；而我國(guó)四川菌株中該基因的攜帶率(14%)低于其他省(廣西：60%，山東：49%，湖北：47%)。與毒力因子研究類似，對(duì)相關(guān)耐藥因子的分布研究需要后期優(yōu)化采樣策略、擴(kuò)大采樣范圍和樣本量后進(jìn)一步驗(yàn)證。

圖3中，虛線框1代表的是菌株2014V-1011 (SAMN12648291)里的少數(shù)耐藥基因，其中4個(gè)()位于該菌株一個(gè)尚未命名的質(zhì)粒(CP046773.1)上。虛線框2和6中的基因分別代表新測(cè)序菌株VA28和VA24上的少數(shù)耐藥基因，這些基因落在了2~3個(gè)組裝片段上(contig)，因此基于現(xiàn)有數(shù)據(jù)無法確認(rèn)其在基因組上的位置是否相鄰，但其中的與兩個(gè)基因是Ⅰ類整合子(class 1 integron)3’端保守區(qū)域(conserved segments, 3’CS)的兩個(gè)固有基因[46,47]，另外虛線框1中的質(zhì)粒片段上也有這兩個(gè)基因，因此我們推斷虛線框2和6中的少數(shù)耐藥基因可能位于通過水平基因轉(zhuǎn)移(horizontal gene transfer, HGT)而獲得的Ⅰ類整合子片段或質(zhì)粒片段上，這需要后期通過基因組完成圖測(cè)序和分析進(jìn)行確定。虛線框3中的基因位于菌株ZJ-T(SAMN05271497)染色體1上的相鄰位置，是一個(gè)連續(xù)的大片段，我們截取這一段序列并在線BLAST比對(duì)后發(fā)現(xiàn)，該片段與溶藻弧菌的整合共軛元件(integrative conjugative elements, ICEs) ICEVal056-1 (accession number: KR231688.1)，coverage和identity值分別為70%和99.98%，表明該菌株可能通過水平轉(zhuǎn)移獲得了一個(gè)攜帶多種耐藥基因的大片段。虛線框4和5中的基因是位于菌株Vb1833 (SAMN15829766)染色體1上的兩個(gè)大片段,我們分別截取了這兩個(gè)片段在線BLAST后發(fā)現(xiàn)，比對(duì)結(jié)果中多為ICE片段或者質(zhì)粒片段，這也提示菌株Vb1833 (SAMN15829766)在進(jìn)化過程中，通過水平轉(zhuǎn)移的方式獲得了攜帶多種耐藥基因的大片段。本研究在溶藻弧菌中發(fā)現(xiàn)了攜帶多種耐藥基因的大片段，為今后溶藻弧菌所致疾病的治療和防控提供了重要信息。

表3 部分毒力因子在不同譜系、不同國(guó)家和地區(qū)以及中國(guó)不同地點(diǎn)的分布統(tǒng)計(jì)

圖3 181株溶藻弧菌毒力與耐藥因子分布熱圖

紫色代表攜帶毒力因子，藍(lán)色代表攜帶耐藥因子，灰色代表缺失。虛線框代表5株菌中的少數(shù)耐藥基因，紅色箭頭所指的兩個(gè)基因?yàn)樽V系2中攜帶率更高的兩個(gè)耐藥基因和。

3 討論

過去研究發(fā)現(xiàn)，副溶血弧菌的種群結(jié)構(gòu)與不同海域密切相關(guān)[24]。但本研究現(xiàn)有數(shù)據(jù)集下，溶藻弧菌不同采樣地點(diǎn)的菌株在整個(gè)系統(tǒng)發(fā)育樹和不同種群間未呈現(xiàn)地域性聚集，而是呈現(xiàn)出“分散分布”，即：同一采樣地點(diǎn)的菌株(如中國(guó)、美洲等)分散于兩個(gè)不同的譜系中。根據(jù)群體遺傳學(xué)理論，種群分化只在種群間遷徙率很低的情況下發(fā)生(每世代中只有一個(gè)遷徙或者更少)[24]，當(dāng)少量遷徙者(migrant)與當(dāng)?shù)胤N群共存時(shí)，會(huì)逐漸與當(dāng)?shù)胤N群交換遺傳物質(zhì)，進(jìn)而變得與當(dāng)?shù)胤N群越來越相似。本研究中種群已經(jīng)形成，但兩個(gè)種群(譜系1和譜系2)中還是包含了不同地域來源的菌株，可能是遷徙者到達(dá)新環(huán)境后，尚未充分與當(dāng)?shù)胤N群交換遺傳物質(zhì)，依然保留著大量“家鄉(xiāng)”種群的遺傳特征。造成該現(xiàn)象的原因包括：跨地域遷徙發(fā)生在近期，進(jìn)化時(shí)間相對(duì)較短；或者自身遺傳特征限制了遺傳物質(zhì)在不同種群的菌株間進(jìn)行交換；或者是不同種群溶藻弧菌生活在同一地區(qū)的不同微生境下，導(dǎo)致事實(shí)上的生態(tài)條件隔離。將來可通過開展生態(tài)調(diào)查和實(shí)驗(yàn)室研究以檢驗(yàn)上述假說。

表4 溶藻弧菌耐藥基因與相關(guān)抗生素和耐藥機(jī)制

括號(hào)內(nèi)百分比代表整個(gè)數(shù)據(jù)集中攜帶此耐藥基因的菌株的百分比。

本研究中，無論是從NCBI下載的公共數(shù)據(jù)，還是新測(cè)序的中國(guó)菌株，在采樣地點(diǎn)上都存在一定的偏差。例如公共數(shù)據(jù)中，美國(guó)菌株占60%，而其他國(guó)家或大洲的菌株則相對(duì)較少，新測(cè)序的中國(guó)菌株也僅采集自中國(guó)的山東、湖北、四川和廣西4個(gè)省份，其他省份數(shù)據(jù)缺失。采樣偏差可能影響本研究中的某些結(jié)果。例如研究結(jié)果顯示，同一譜系中(譜系1或譜系2)包含了不同大洲的菌株，提示歷史上曾存在多次遠(yuǎn)距離遷徙事件，但本研究中我們僅在一個(gè)克隆群中觀察到跨洋傳播，可能是由于短時(shí)期遷徙事件頻率低，或者是研究數(shù)據(jù)集不夠完整所致。此外，種群結(jié)構(gòu)研究中，發(fā)現(xiàn)歐洲菌株只存在于譜系1中；毒力、耐藥因子分布研究中，發(fā)現(xiàn)某些毒力、耐藥因子(、、、/、等)在不同的大洲以及中國(guó)不同的省份之間的分布存在差異。這些結(jié)果需要在將來通過完善采樣策略、擴(kuò)大采樣范圍、增加數(shù)據(jù)量等方式消除采樣偏差造成的影響并進(jìn)一步驗(yàn)證。

雖然本研究中鑒定到多種不同功能的毒力因子，但是除3株菌采集自病人外，其余均為環(huán)境菌株，因而難以將該病原的毒力因子與菌株致病性表型進(jìn)行關(guān)聯(lián)。此外，雖然本研究中鑒定到多種耐藥相關(guān)基因，但部分鑒定結(jié)果可能為耐藥基因的同源體(homology)，而是否攜帶耐藥基因同源體與菌株的臨床耐藥表型并不一定具有相關(guān)性，因?yàn)椴糠植≡谂R床上表現(xiàn)出對(duì)抗生素的耐受，是由染色體上的具體的基因突變導(dǎo)致的，包括點(diǎn)突變、基因片段獲得缺失(Indel)等。另一方面，一個(gè)基因的表達(dá)與否還受到多種因素的調(diào)控，因此攜帶特定耐藥基因(或同源體)的菌株是否表現(xiàn)出耐藥表型，還要結(jié)合實(shí)驗(yàn)測(cè)定結(jié)果進(jìn)一步研究。

綜上所述，本研究利用全基因組測(cè)序技術(shù)和群體遺傳學(xué)方法，對(duì)181株溶藻弧菌的遺傳多樣性和毒力、耐藥特征進(jìn)行了研究。研究結(jié)果表明溶藻弧菌是一種高頻重組致病弧菌，由兩個(gè)主要譜系組成，其中譜系1菌株的地域分布更為廣泛。兩個(gè)譜系菌株的地理分布特點(diǎn)和克隆群分析均表明該物種存在近期跨洋傳播，但種群跨地域遷徙后尚未觀察到種群融合。溶藻弧菌基因組包含多種不同功能的毒力因子，主要與溶藻弧菌自身代謝、對(duì)宿主細(xì)胞的粘附、裂解、免疫反應(yīng)等功能相關(guān)，兩個(gè)譜系中均未發(fā)現(xiàn)種群特異性毒力因子。此外，溶藻弧菌基因組既有多種核心耐藥基因，也存在攜帶多種耐藥基因的質(zhì)粒或者ICE等移動(dòng)遺傳元件。本研究存在采樣偏差以及缺乏表型實(shí)驗(yàn)結(jié)果等不足，但上述結(jié)果為更深入研究溶藻弧菌遺傳特征和致病機(jī)制提供了基礎(chǔ)數(shù)據(jù)，將加深對(duì)該病原的認(rèn)識(shí)并促進(jìn)其防控工作。此外，本研究中的研究思路和方法，也可供霍亂弧菌、河流弧菌()等其他高頻重組致病菌的群體基因組學(xué)研究參考，通過揭示相關(guān)病原菌的遺傳多樣性、致病因子、傳播模式等特征，在溯源、防控和治療等方面提供科學(xué)支撐并發(fā)揮積極作用。

[1] Lajnef R, Snoussi M, Romalde JL, Nozha C, Hassen A. Comparative study on the antibiotic susceptibility and plasmid profiles of Vibrio alginolyticus strains isolated from four Tunisian marine biotopes., 2012, 28(12): 3345–3363.

[2] Xie JS, Bu LF, Jin S, Wang XY, Zhao QS, Zhou SM, Xu YJ. Outbreak of vibriosis caused by Vibrio harveyi and Vibrio alginolyticus in farmed seahorse Hippocampus kuda in China., 2020, 523: 735168.

[3] Reilly GD, Reilly CA, Smith EG, Baker-Austin C. Vibrio alginolyticus-associated wound infection acquired in British waters, Guernsey, July 2011., 2011, 16(42): 19994.

[4] Zhou K, Sun F, Zha DJ, Liu JY, Ma YY, Hao XK, Xu XL. Biological characteristics of a Vibrio alginolyticus strain isolated from a patient with chronic otitis externa., 2020, 158(2): 146–151.周柯, 孫菲, 查定軍, 劉家云, 馬越云, 郝曉柯, 徐修禮. 慢性外耳道炎患者分離溶藻弧菌的鑒定及生物學(xué)特性分析. 中國(guó)病原生物學(xué)雜志, 2020, 158(2): 146–151.

[5] Gómez-León J, Villamil L, Lemos ML, Novoa B, Figueras A. Isolation of Vibrio alginolyticus and Vibrio splendidus from aquacultured carpet shell clam (Ruditapes decussatus) larvae associated with mass mortalities., 2005, 71(1): 98–104.

[6] González-Escalona N, Blackstone GM, Depaola A. Characterization of a Vibrio alginolyticus strain, isolated from Alaskan oysters, carrying a hemolysin gene similar to the thermostable direct hemolysin-related hemolysin gene (trh) of Vibrio parahaemolyticus., 2006, 72(12): 7925–7929.

[7] Newton A, Kendall M, Vugia DJ, Henao OL, Mahon BE. Increasing rates of vibriosis in the United States, 1996-2010: review of surveillance data from 2 systems., 2012, 54 Suppl 5(05): S391–S395.

[8] Jacobs Slifka KM, Newton AE, Mahon BE. Vibrio alginolyticus infections in the USA, 1988-2012., 2017, 145(7): 1491–1499.

[9] Hu NX, Ma C, Li YJ. Research progress on molecular typing technology of Vibrio alginolyticus., 2014, 35(13): 1753–1755.胡南霞, 馬聰, 李艷君. 溶藻弧菌分子分型技術(shù)的研究進(jìn)展. 國(guó)際檢驗(yàn)醫(yī)學(xué)雜志, 2014, 35(13): 1753–1755.

[10] Sudheesh PS, Jie K, Xu HS. Random amplified polymorphic DNA-PCR typing of Vibrio parahaemolyticus and V. alginolyticus isolated from cultured shrimps., 2002, 207(1–2): 11–17.

[11] Ouyang JL, Zhou YC, Wu YG, Dai XL, Wang SF, Xie ZY. RAPD analysis of genetic diversity of Vibrio alginolyticus isolated from Hainan province., 2011, 35(7): 30–36.歐陽吉隆, 周永燦, 吳學(xué)貴, 戴小連, 王世鋒, 謝珍玉. 海南地區(qū)溶藻弧菌遺傳多樣性的RAPD分析. 海洋科學(xué), 2011, 35(7): 30–36.

[12] Ren CH, Hu CQ, Luo P, Chen C, Jiang X, Wang QB. Genotyping of Vibrio alginolyticus isolates from Daya Bay by infrequent-restriction-site PCR and pulsed-field gel electrophoresis., 2008, 22(4): 267– 271.

[13] Talon D, Cailleaux V, Thouverez M, Michel-Briand Y. Discriminatory power and usefulness of pulsed-field gel electrophoresis in epidemiological studies of Pseudomonas aeruginosa., 1996, 32(2): 135–145.

[14] Inglis TJ, Garrow SC, Adams C, Henderson M, Mayo M, Currie BJ. Acute melioidosis outbreak in Western Australia., 1999, 123(3): 437–443.

[15] Xie ZY, Hu CQ, Chen C, Zhang LP, Ren CH. Investigation of seven Vibrio virulence genes among Vibrio alginolyticus and Vibrio parahaemolyticus strains from the coastal mariculture systems in Guangdong, China., 2005, 41(2): 202–207.

[16] Hernández-Robles MF, álvarez-Contreras AK, Juárez- García P, Natividad-Bonifacio I, Curiel-Quesada E, Vázquez-Salinas C, Qui?ones-Ramírez EI. Virulence factors and antimicrobial resistance in environmental strains of Vibrio alginolyticus., 2016, 19(4): 191–198.

[17] Mechri B, Medhioub A, Medhioub MN, Aouni M. Genotypic diversity, antimicrobial resistance and screening of Vibrio cholerae molecular virulence markers in Vibrio alginolyticus strains recovered from a Tunisian Ruditapes decussatus hatchery., 2013, 62(3): 263– 272.

[18] Wei S. Detection of Vibrio alginolyticus and its virulence- related genes[Dissertation]. Jinan University, 2013.魏霜. 溶藻弧菌及其毒力相關(guān)基因的檢測(cè)[學(xué)位論文]. 暨南大學(xué), 2013.

[19] Xian YY. Analysis of the distribution of virulence-related genes and antibiotic resistance of Vibrio parahaemolyticus and Vibrio alginolyticus[Dissertation]. Jinan University, 2015.冼鈺茵. 副溶血弧菌與溶藻弧菌毒力基因分布與抗生素耐藥性分析[學(xué)位論文]. 暨南大學(xué), 2015.

[20] Mechri B, Monastiri A, Medhioub A, Medhioub MN, Aouni M. Molecular characterization and phylogenetic analysis of highly pathogenic Vibrio alginolyticus strains isolated during mortality outbreaks in cultured Ruditapes decussatus juvenile., 2017, 111: 487–496.

[21] Chen HQ. Analysis of drug resistance and drug resistance genes of a pathogenic Vibrio alginolyticus from large yellow croaker., 2020, 35(1): 63–67.陳洪清. 一株大黃魚致病性溶藻弧菌的耐藥性及耐藥基因分析. 漁業(yè)信息與戰(zhàn)略, 2020, 35(1): 63–67.

[22] Limmathurotsakul D, Golding N, Dance DAB, Messina JP, Pigott DM, Moyes CL, Rolim DB, Bertherat E, Day NPJ, Peacock SJ, Hay SI. Predicted global distribution of Burkholderia pseudomallei and burden of melioidosis., 2016, 1: 15008.

[23] Chewapreecha C, Holden MTG, Vehkala M, V?lim?ki N, Yang ZR, Harris SR, Mather AE, Tuanyok A, De Smet B, Le Hello S, Bizet C, Mayo M, Wuthiekanun V, Limmathurotsakul D, Phetsouvanh R, Spratt BG, Corander J, Keim P, Dougan G, Dance DAB, Currie BJ, Parkhill J, Peacock SJ. Global and regional dissemination and evolution of Burkholderia pseudomallei., 2017, 2: 16263.

[24] Yang C, Pei XY, Wu YR, Yan L, Yan YF, Song YQ, Coyle NM, Martinez-Urtaza J, Quince C, Hu QH, Jiang M, Feil E, Yang DJ, Song YJ, Zhou DS, Yang RF, Falush D, Cui YJ. Recent mixing of Vibrio parahaemolyticus populations., 2019, 13(10): 2578–2588.

[25] Yang C, Zhang XL, Fan H, Li YH, Hu QH, Yang RF, Cui YJ. Genetic diversity, virulence factors and farm-to-table spread pattern of Vibrio parahaemolyticus food-associated isolates., 2019, 84: 103270.

[26] Lawson DJ, Hellenthal G, Myers S, Falush D. Inference of population structure using dense haplotype data., 2012, 8(1): e1002453.

[27] Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data., 2014, 30(15): 2114–2120.

[28] Cui YJ, Yu C, Yan YF, Li DF, Li YJ, Jombart T, Weinert LA, Wang ZY, Guo ZB, Xu LZ, Zhang YJ, Zheng HC, Qin N, Xiao X, Wu MS, Wang XY, Zhou DS, Qi ZZ, Du ZM, Wu HL, Yang XW, Cao HZ, Wang H, Wang J, Yao SS, Rakin A, Li YR, Falush D, Balloux F, Achtman M, Song YJ, Wang J, Yang RF. Historical variations in mutation rate in an epidemic pathogen, Yersinia pestis., 2013, 110(2): 577–582.

[29] Delcher AL, Salzberg SL, Phillippy AM. Using MUMmer to identify similar regions in large sequence sets., 2003, 0(1): 10.3.1–10.3.18.

[30] Seemann T. Prokka: rapid prokaryotic genome annotation., 2014, 30(14): 2068–2069.

[31] Page AJ, Cummins CA, Hunt M, Wong VK, Reuter S, Holden MTG, Fookes M, Falush D, Keane JA, Parkhill J. Roary: rapid large-scale prokaryote pan genome analysis., 2015, 31(22): 3691–3693.

[32] Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees., 2016, 44(W1): W242– W245.

[33] Liu B, Zheng DD, Jin Q, Chen LH, Yang J. VFDB 2019: a comparative pathogenomic platform with an interactive web interface., 2019, 47(D1): D687– D692.

[34] Alcock BP, Raphenya AR, Lau TTY, Tsang KK, Bouchard M, Edalatmand A, Huynh W, Nguyen ALV, Cheng AA, Liu SH, Min SY, Miroshnichenko A, Tran HK, Werfalli RE, Nasir JA, Oloni M, Speicher DJ, Florescu A, Singh B, Faltyn M, Hernandez-Koutoucheva A, Sharma AN, Bordeleau E, Pawlowski AC, Zubyk HL, Dooley D, Griffiths E, Maguire F, Winsor GL, Beiko RG, Brinkman FSL, Hsiao WWL, Domselaar GV, Mcarthur AG. CARD 2020: antibiotic resistome surveillance with the com-prehensive antibiotic resistance database., 2020, 48(D1): D517–D525.

[35] Bortolaia V, Kaas RS, Ruppe E, Roberts MC, Schwarz S, Cattoir V, Philippon A, Allesoe RL, Rebelo AR, Florensa AF, Fagelhauer L, Chakraborty T, Neumann B, Werner G, Bender JK, Stingl K, Nguyen M, Coppens J, Xavier BB, Malhotra-Kumar S, Westh H, Pinholt M, Anjum MF, Duggett NA, Kempf I, Nyk?senoja S, Olkkola S, Wieczorek K, Amaro A, Clemente L, Mossong J, Losch S, Ragimbeau C, Lund O, Aarestrup FM. ResFinder 4.0 for predictions of phenotypes from genotypes., 2020, 75(12): 3491–3500.

[36] Wang HQ, Yang C, Sun Z, Zheng W, Zhang W, Yu H, Wu YR, Didelot X, Yang RF, Pan JC, Cui YJ. Genomic epidemiology of Vibrio cholerae reveals the regional and global spread of two epidemic non-toxigenic lineages., 2020, 14(2): e0008046.

[37] Cui YJ, Yang C, Qiu HL, Wang H, Yang RF, Falush D. The landscape of coadaptation in Vibrio parahaemolyticus., 2020, 9: e54136.

[38] Zhao Z, Liu JX, Deng YQ, Huang W, Ren CH, Call DR, Hu CQ. The Vibrio alginolyticus T3SS effectors, Val1686 and Val1680, induce cell rounding, apoptosis and lysis of fish epithelial cells., 2018, 9(1): 318–330.

[39] Makino K, Oshima K, Kurokawa K, Yokoyama K, Uda T, Tagomori K, Iijima Y, Najima M, Nakano M, Yamashita A, Kubota Y, Kimura S, Yasunaga T, Honda T, Shinagawa H, Hattori M, Iida T. Genome sequence of Vibrio parahaemolyticus: a pathogenic mechanism distinct from that of V cholerae., 2003, 361(9359): 743–749.

[40] Zhou XH, Konkel ME, Call DR. Type III secretion system 1 of Vibrio parahaemolyticus induces oncosis in both epithelial and monocytic cell lines., 2009, 155(Pt 3): 837–851.

[41] Burdette DL, Seemann J, Orth K. Vibrio VopQ induces PI3-kinase-independent autophagy and antagonizes pha-gocytosis., 2009, 73(4): 639–649.

[42] Jia AR, Woo NYS, Zhang XH. Expression, purification, and characterization of thermolabile hemolysin (TLH) from Vibrio alginolyticus., 2010, 90(2): 121–127.

[43] Wong SK, Zhang XH, Woo NYS. Vibrio alginolyticus thermolabile hemolysin (TLH) induces apoptosis, membrane vesiculation and necrosis in sea bream erythrocytes., 2012, 330–333: 29–36.

[44] Cai SH, Yaol SY, Lu YS, Wu ZH, Jian JC, Wang B. Immune response in Lutjanus erythropterus induced by the major outer membrane protein (OmpU) of Vibrio alginolyticus., 2010, 90(1): 63–68.

[45] Lv TT, Dai F, Zhuang QT, Zhao XL, Shao YN, Guo M, Lv ZM, Li CH, Zhang WW. Outer membrane protein OmpU is related to iron balance in Vibrio alginolyticus., 2020, 230: 126350.

[46] Cambray G, Guerout AM, Mazel D. Integrons., 2010, 44: 141–166.

[47] Rajpara N, Nair M, Bhardwaj AK. A highly promiscuous integron, plasmids, extended spectrum beta lactamases and efflux pumps as factors governing multidrug resistance in a highly drug resistant Vibrio fluvialis isolate BD146 from Kolkata, India., 2018, 58(1): 60–67.

Population genomics study of

Hongyuan Zheng1, Lin Yan2, Chao Yang1,3, Yarong Wu1, Jingliang Qin1, Tongyu Hao1, Dajin Yang2, Yunchang Guo2, Xiaoyan Pei2, Tongyan Zhao1, Yujun Cui1

is a Gram-negative bacillus that causes vibriosis to human and aquatic products, including fish, shrimp and shellfish. It poses a threat to public health and causes enormous economic losses to the aquaculture industry. However, research on genetic diversity and pathogenicity-related genetic elements based on whole genome is still lacking.In this study, sixty-eight strains ofwere collected from four provinces of China and the whole genome sequences were obtained. Combined with 113 publicly available genome sequences downloaded from NCBI, we inferred the population structure ofby using fineSTRUCTURE software, and identified the virulence and antibiotic resistance factors using the VFDB, CARD and ResFinder database. The results indicated thatincluded two main lineages, named Lineage 1 and Lineage 2. Both lineages distributed in America and Asia, but all the European genomes were classified into Lineage 1. A single cross-ocean transmission event was inferred from one of the 12 identified clonal groups in our dataset.genome contains a variety of virulence factors, such as,, andetc. The distribution of virulence factors revealed no lineage-specificity, but some of which revealed differences in their geographical distribution. A lower frequency of,,,and a higher frequency ofwere observed in genomes of Europe than other continents. In China, a lower frequency of,and nowere observed in genomes from Guangxi province. Among the identified antibiotic resistance genes,andare significantly enriched in Lineage 2. In addition,is more common in genomes from Asia, compared with the American and European genomes. But in China, the frequency ofin Sichuan genomes is much lower than in other provinces. We also found that large fragments of plasmids or ICEs that carried multiple drug resistance genes were present in fivegenomes (VA24, VA28, 2014V-1011, ZJ-T and Vb1833). Based on population genomics analysis, our study delineated the population structure, distribution of virulence and antibiotic resistance related factors of, which lays a foundation for future study of genetic characters and pathogenesis mechanism of this pathogen and will improve the works on monitoring, prevention and control of this pathogen.

;whole genome sequencing; population structure; virulence factor; antimicrobial resistance gene

2021-02-09;

2021-03-16

國(guó)家重點(diǎn)研發(fā)計(jì)劃(編號(hào)：2018YFC1603902，2017YFC1601500)和中國(guó)博士后科學(xué)基金項(xiàng)目(編號(hào)：2020M672836)資助[Supported by the National Key R&D Program (Nos. 2018YFC1603902, 2017YFC1601500) of China, and China Postdoctoral Science Foundation (No. 2020M672836)]

鄭宏源，在讀碩士研究生，專業(yè)方向：細(xì)菌基因組流行病學(xué)。E-mail: zhenghongyuan.chn@qq.com

崔玉軍，博士，研究員，研究方向：微生物進(jìn)化與基因組流行病學(xué)。E-mail: cuiyujun.new@gmail.com

10.16288/j.yczz.21-061

2021/4/1 15:08:36

URI: https://kns.cnki.net/kcms/detail/11.1913.R.20210401.1103.004.html

(責(zé)任編委: 高海春)