多重耐藥胸膜肺炎放線桿菌GD2107株全基因組測序及生物信息學分析

徐民生，柯海意，楊冬霞，施科達，孫澤儀，牛佳偉，常 鑫，翟少倫，臧瑩安，李春玲

(1.仲愷農業工程學院動物科技學院，廣州 510305；2.廣東省農業科學院動物衛生研究所，廣東省畜禽疫病 防治研究重點實驗室，農業農村部獸用藥物與診斷技術廣東科學觀測實驗站，廣州 510640)

胸膜肺炎放線桿菌(Actinobacilluspleuropneumoniae，APP)是一種引起豬傳染性胸膜肺炎(porcine contagious pleuropneumia，PCP)的革蘭氏陰性桿狀細菌，其定植于豬上呼吸道且具有高度傳染性，可感染各年齡段豬并引發急性或慢性纖維出血性壞死性肺炎癥狀，嚴重情況下會導致死亡[1]。根據對煙酰胺腺嘌呤二核苷酸(nicotinamide adenine dinucleotide,NAD)的生長依賴性可將APP分為生物Ⅰ型和生物Ⅱ型[2]。Apx毒素在APP致病過程中起著重要作用，其中ApxⅣ在所有血清型的APP中均可分泌產生，可用于鑒別APP[3]。目前，基于衣殼多糖(CPS)的組成有19個公認的血清型，各血清型流行率具有一定的地理差異，治療和防控此病仍依賴于抗生素和疫苗[4]。APP引起的急性或慢性病癥導致較高的死亡率或治療費用，給全球養殖業造成了重大損失[5]。

自1995年第一個完整的細菌基因組序列公布以來，全基因組測序技術(whole genome sequencing,WGS)得到了蓬勃發展[6]。從傳統的Sanger測序到基于Illumina平臺的第二代測序和PacBio平臺的第三代測序，不僅提高了完整基因組的組裝速度和質量，所需的時間和費用也逐步縮短和降低。全基因組測序技術作為一種有力測序工具在病原體的表征和流行病學調查中的運用也越來越多[7]。Li等[8]研究表明，全基因組測序技術可用來預測抗生素的耐藥性，識別新的血清型以便用于研究藥物和疫苗的作用靶點。全基因組序列正迅速成為預測APP生物信息的有利工具，尤其是對于了解相關菌株基因特征和進化關系具有重要意義[9]。

本研究以先前從廣東省某發病豬場分離到的1株胸膜肺炎放線桿菌GD2107株為研究對象，藥敏結果顯示其具有多重耐藥性。為進一步探索其致病和耐藥機理，使用Illumina NovaSeq和PacBio SequeI測序平臺進行全基因組測序和組裝，通過生物信息學分析軟件進行基因功能注釋，了解其攜帶的耐藥基因、毒力基因情況，分析菌株間親緣關系，揭示菌株耐藥與進化關系的相關性，以期為豬傳染性胸膜肺炎的基礎研究提供一個新的視角與借鑒。

1 材料與方法

1.1 材料

1.1.1 菌株來源 APP GD2107于2021年7月分離自廣東省某規模化發病豬場感染傳染性胸膜肺炎病死豬的肺部，由廣東省農業科學院動物衛生研究所豬病研究室分離、鑒定和保存。

1.1.2 主要試劑與儀器 胰蛋白大豆肉湯培養基(TSB培養基)、煙酰胺腺嘌呤二核苷酸(NAD)均購自北京索萊寶科技有限公司；小牛血清購自廣州蕊特生物科技有限公司；藥敏紙片購自杭州微生物試劑有限公司；細菌基因組DNA提取試劑盒、2×FastTaqPremix酶、DL2000 DNA Marker均購自TaKaRa公司。電熱恒溫培養箱購自上海精宏實驗設備有限公司；PCR儀購自杭州晶格科學儀器有限公司；電泳儀購自北京市六一儀器廠；顯微鏡購自Olympus公司。

1.2 方法

1.2.1 細菌培養及ApxⅣ基因擴增 將菌株GD2107凍存液以1∶100的比例加入含0.1% NAD和5%牛血清的TSB培養液中，37 ℃、220 r/min培養12 h。使用細菌DNA提取試劑盒提取基因組DNA。參考Zhang等[10]合成ApxⅣ基因特異性引物，引物序列為：F：5′-TTATCCGAA-CTTTGGTTTAGCC-3′；R：5′-CATATTTGATA-AAACCATCCGTC-3′，預期擴增片段長度為418 bp。引物由生工生物工程(上海)股份有限公司合成。以提取的基因組DNA為模板對菌株進行鑒定，PCR反應體系25 μL：2×rTaqMix 酶12.5 μL，上、下游引物(10 pmol/μL)各0.5 μL，DNA 模板1.0 μL，ddH2O 10.5 μL。PCR反應條件：94 ℃ 預變性5 min；94 ℃變性30 s，58 ℃退火30 s，72 ℃延伸1 min，共35個循環；72 ℃延伸10 min。PCR產物經1.0%瓊脂糖凝膠電泳后，將ApxⅣ基因陽性PCR產物送至生工生物工程(上海)股份有限公司測序，菌株GD2107基因組DNA送至上海派森諾生物科技股份有限公司進行全基因組測序。

1.2.2 藥敏試驗 采用K-B藥敏紙片法對分離菌株進行藥敏試驗。將分離菌株接種至含有5%牛血清及0.1% NAD的TSB培養基中，制成含菌量約為0.5麥氏濁度的菌懸液，吸取100 μL并用無菌棉簽均勻涂布于TSA平板(含0.1% NAD和5%牛血清)表面，待液體吸干后用無菌鑷子夾取抗生素藥敏紙片貼在平板表面，置于37 ℃溫箱中培養30 min，后倒置培養12～18 h，記錄抑菌圈直徑并參照杭州微生物試劑有限公司提供的判斷標準進行結果判定。

1.2.3 全基因組測序 采用全基因組鳥槍法(whole genome shotgun，WGS)策略構建PE和S10K2個不同片段的文庫，分別利用Illumina NovaSeq和PacBio SequeI測序平臺對GD2107株進行全基因組測序。采用FastQc(v 0.11.7)[11]軟件和AdapterRemoval(v 2.1.7)[12]和SOAPec(v 2.0)軟件[13]對測序數據進行質量控制和過濾。使用HGAP(v 4.0)[14]和CANU(v 1.7.1)[15]軟件進行拼裝，得到contig序列后用pilon(v 1.22)[16]軟件對結果進行校正，最終拼接得到完整序列。

1.2.4 基因組圈圖繪制 首先將基因組序列、基因組預測及非編碼RNA的預測信息整合成標準的GBK(GenBank)格式文件；然后采用cgview[17](http:∥stothard.afns.ualberta.ca/cgview_server)繪制該基因組的圈圖；最后采用Photoshop CS對該圖進行編輯。

1.2.5 功能元件分析 采用不同生物信息學分析軟件對組裝后的基因組、非編碼RNA(ncRNA，主要類型包括sRNA、rRNA、tRNA、snRNA、snoRNA及miRNA等)、CRISPRs功能元件、原噬菌體(prophage)、基因島(genomics islands，GIs)進行預測。用GeneMarkS(v 4.32)[18]軟件預測GD2107全基因組基因序列；用tRNAscan-SE(v 1.3.1)[19]預測全基因組中的tRNA基因；用Barrnap(v 0.9)預測rRNA基因，其余ncRNA主要通過與Rfam(v 14.1)[20]數據庫比對進行預測；用CRISPR finder(v 20170509)[21]預測全基因組中的正向重復序列和間隔區；用PhiSpy[22](http:∥phaster.ca)預測基因組中存在的原噬菌體；用IslandViewer4[23](http:∥www.pathogenomics.sfu.ca/islandviewer)預測基因組中存在的GIs。

1.2.6 蛋白編碼基因功能注釋 將所有預測得到的蛋白編碼基因與數據庫中包含的蛋白編碼基因進行比對，若兩者具有顯著的序列相似性，則可初步確定其編碼合成的蛋白具有相似甚至是相同的功能；采用diamond(v 0.9.10.111)軟件完成序列比對[24]。序列比對的臨界值選取為1e-6。序列的功能判別規則為：E-value<1e-6。

1.2.7 亞系統分析 對組裝后的GD2107基因組進行碳水化合物活性酶(CAZy)、毒力因子(VFDB)和耐藥因子(CARD)相關亞系統的信息學分析：采用hmmscan[25]軟件與CAZy數據庫比對來預測全基因組中存在的CAZy酶類基因；通過與VFDB[26]和CARD數據庫[27]進行BLAST比對，預測存在的毒力基因和耐藥基因。

1.2.8 系統進化樹構建 將測序結果在GenBank數據庫中進行BLAST比對。選取NCBI數據庫中與GD2107株相似性相近的14條序列(表1)，利用Mega 7.0軟件構建系統進化樹。

表1 參考菌株信息Table 1 Reference strain information

2 結 果

2.1 ApxⅣ基因擴增結果

采用針對APP菌種的特異性引物ApxⅣ基因擴增，結果獲得大小為418 bp的特異性條帶(圖1)，與預期結果相符。

2.2 藥敏試驗結果

藥敏試驗結果顯示，菌株GD2107只對氨芐西林、阿莫西林/克拉維酸、頭孢唑啉(先鋒Ⅴ)敏感，對青霉素、頭孢拉定(先鋒Ⅵ)、卡那霉素、鏈霉素等14種抗菌藥耐藥(表2)。

2.3 GD2107株的基因組結構與分析

由圖2可知，GD2107株的全基因組包括1條大小為2 271 987 bp大小的環狀染色體，GC含量為41.21%，另存在2條大小為5 027和3 497 bp的環狀質粒，GC含量分別為44.06%和36.06%。最后將基因組序列提交至NCBI數據庫中，得到對應的染色體GenBank的登錄號為CP097377，2條質粒的GenBank登錄號分別為CP097378和CP097379。

M，DL2000 DNA Marker；-，陰性對照；+，陽性對照；1，分離株 M，DL2000 DNA Marker; -，Negative control; +，Positive control; 1，The isolate圖1 分離株ApxⅣ基因PCR鑒定結果Fig.1 PCR identification results of ApxⅣ gene in the isolates

表2 藥敏試驗結果Table 2 Drug sensitivity test result mm

從內到外，第1圈代表刻度；第2圈代表GC Skew；第3圈代表GC含量；第4、7圈代表每一個CDS所屬的COG；第5、6圈代表CDS、tRNA及rRNA在基因組上的位置 From inside to outside,the first circle represents the scale;The second circle represents GC Skew;The third circle represents GC content;The fourth and seventh circles represent the COG of each CDS;The positions of CDS,tRNA and rRNA in the genome in the fifth and sixth cycles圖2 胸膜肺炎放線桿菌GD2107全基因組染色體圈圖Fig.2 Circular map of Actinobacillus pleuropneumoniae GD2107 chromosome

測序得到的開放閱讀框(ORF)數量有2 290個，總長度為1 972 470 bp，占全基因組的86.82%；在非編碼RNA(ncRNA)中含有tRNA基因21個，rRNA基因19個(其中5S rRNA基因7個，16S rRNA基因6個，23S rRNA基因6個)，其他ncRNA基因20個；使用PhiSpy預測基因組中存在的4個原噬菌體；經CRISPR finder預測發現菌株含有2組CRISPR相關序列；通過IslandViewer4預測到菌株GD2107基因組中有23個基因島；基因組結構預測結果見表3。

表3 基因組序列基本信息Table 3 Basic information of genome sequence

2.4 蛋白編碼基因功能注釋

將GD2107株基因的蛋白序列與表4中的功能數據庫進行比對分析，對于每一條序列選取得分最高的比對結果進行注釋；GD2107株的全基因組共有2 290個蛋白質編碼基因，其中，分別有2 263、1 979、2 112、1 549、1 866、1 735、514、557個基因在NR、Swiss-Prot、COG、KEGG、GO、Pfam、TCDB、PHI數據庫中得到注釋(表4)。

表4 蛋白編碼基因功能注釋匯總Table 4 Summary of protein coding gene function annotations

2.4.1 COG數據庫注釋結果 結果顯示，共注釋到2 112個蛋白，占全部基因的92.23%，分布于25個COG的條目中(圖3)，其中，參與復制、重組和修復的基因有166個，占7.25%；參與翻譯、核糖體結構和生物合成的基因有179個，占7.82%；與細胞壁、細胞膜、胞外被膜生物合成相關的基因有200個，占8.73%；參與氨基酸轉運和代謝的基因有173個，占7.55%；與轉錄相關的基因有84個，占3.67%。

圖3 COG功能分類圖Fig.3 COG function classification diagram

2.4.2 KEGG和GO注釋結果 GD2107菌株主要富集在碳水化合物代謝、膜運輸、輔助因子、維生素代謝和氨基酸這些代謝途徑，也是維持細菌新陳代謝所必需的。經KEGG注釋可分為8類45個功能條目(圖4)，其分別注釋為31.9%新陳代謝、7.1%遺傳信息處理、6.5%環境信息處理、3.1%細胞過程、2.9%人類疾病及1.6%生物體系統，其中每一類型又包含有各自的亞型。

圖4 GD2107基因組KEGG功能注釋結果Fig.4 Results of KEGG functional annotation of GD2107 genome

通過GO注釋，共有1 866個基因，其基因產物共分為三大類：生物過程(biological process,BP)占45.4%；分子功能(molecular function,MF)占23.1%；細胞組分(cellular component,CC)占31.5%(圖5)。表明該菌株的基因產物主要集中在生物過程方面，其中大量基因與生物過程、細胞氮化合物代謝、生物合成和小分子代謝有關；在細胞組分方面，GD2107主要與細胞、細胞位置、胞外區、細胞質和漿膜等相關；在分子功能方面，主要與分子功能、離子結合等相關。

圖5 GD2107基因組GO功能注釋結果Fig.5 Results of GO functional annotation of GD2107 genome

2.4.3 TCDB和PHI注釋結果 TCDB第一級統計結果見圖6。由圖6可知，共有514個基因注釋到TCDB數據庫的9個功能分類中，分別為：通道/孔隙(channels/pores)80個基因；電化學勢驅動轉運體(electrochemical potential-driven transporters)116個基因；主要活性轉運蛋白(primary active transporters)221個基因；轉錄組(group translocators)21個基因；跨膜電子載體(transmembrane electron carriers)17個基因；運輸中涉及的附屬因素(accessory factors involved in transport)6個基因；不完全表征的運輸系統(incompletely characterized transport systems)53個基因。

圖6 TCDB功能分類圖Fig.6 TCDB functional classification diagram

病原體PHI表型突變類型基因數目的統計情況見圖7。由圖7可知，這557個基因中有3個突變后致病菌株具有化學耐受性、2個具有化學敏感性、3個突變成效應子、42個突變后毒力增加、39個喪失致病性、338個突變后毒力降低、122個突變后不受影響。

圖7 病原體PHI表型突變類型分布圖Fig.7 Distribution of PHI phenotype mutation types of pathogens

2.5 亞系統分析結果

2.5.1 CAZy數據庫注釋結果 由圖8可知，GD2107株中有72個基因注釋到CAZy數據庫中，其中包含與糖基轉移酶(GTs)相關基因35個、多糖裂解酶基因2個、糖類脂解酶基因10個、氧化還原酶基因6個、碳水化合物基因5個、糖苷水解酶基因14個。

圖8 CAZy功能分類圖Fig.8 CAZy function classification diagram

2.5.2 VFDB基因預測分析結果 將GD2107株全基因序列與致病菌毒力因子綜合數據庫VFDB進行比對，共注釋到48個毒力基因，其中毒力因子脂寡糖(lipooligosaccharide，LOS)和莢膜(capsule)的基因占比最大(70.8%)。質粒中并不存在毒力基因。

2.5.3 CARD分析 將測序序列與CARD進行比對分析，共獲得了22個耐藥基因，分別為aph(3″)-Ⅰb、macA、macB、sul3、tet34、tet35、tetO、parC、parE、gyrA、gyrB、floR、rpoB、rpoC、dfrA3、hns、tuf2、golS、folP、crp和cpxR，且分別對應表5中的抗菌藥，其中氟苯尼考藥物的耐藥基因floR位于質粒上。

表5 GD2107株耐藥基因及相對應抗菌藥Table 5 Resistance genes and corresponding antibiotics of GD2107 strain

2.6 系統進化樹分析

將PCR擴增得到的APP菌株經ApxⅣ基因擴增后進行測序，并將測序結果與14株參考菌株進行序列相似性比對分析，系統進化樹結果顯示，GD2107株與中國APP5株(CP063424.1)進化關系最近，處在同一分支；其次是中國分離株JL03株(CP000687.1);與韓國KL16、瑞士CVJ13216、瑞士WF83、愛爾蘭S405等地的分離株進化關系均較近；而與美國NCTC11384、英國MDG2331、加拿大L20、匈牙利A-85/14、德國AP76、中國APP6、瑞士P1875、丹麥NCTC10976等APP分離株進化關系較遠(圖9)。

圖9 GD2107分離株基于ApxⅣ基因序列構建的系統進化樹Fig.9 Phylogenetic tree of GD2107 isolate based on ApxⅣ gene sequence

3 討 論

APP是一種臨床上常見的呼吸道病原菌，已成為目前影響豬場經濟效益的重要細菌性疾病之一[28]。近年來隨著飼料“禁抗”政策的推行和養殖場抗生素的濫用導致耐藥性APP的滋生和擴散，這提示應長期對APP耐藥性進行監測，盡量減少由APP感染引起的潛在危險[29]。

全基因組測序技術的快速發展給研究微生物多樣性、進化及物種間相互作用帶來了新的手段，利用全基因組測序技術進行DNA水平的測序分析已成為近年來的發展趨勢[30]，盡管細菌基因組信息的獲取變得越來越容易，然而對APP完整基因序列的報道和研究相對較少。本研究所用GD2107株是從臨床上1頭發病豬的肺臟中分離得到的多重耐藥菌株。通過全基因組測序和生物信息學分析了解GD2107株的基因組結構和功能，豐富了APP的基因組數據庫，為進一步探明其致病機制和耐藥機制提供重要的生物信息學依據。

通過PacBio及Illumina測序平臺的相互結合，高效并準確地完成了對GD2107株全基因組的測序。測序到的蛋白質編碼基因經GO、KEGG和COG數據庫注釋分析發現，3個數據庫的注釋結果高度一致，基因富集程度最高的都在代謝過程，這說明了代謝過程對于維持細菌的生命至關重要[31]。碳水化合物結合結構域能通過結合CAZy的底物提高CAZy催化底物的活性[32]。從注釋結果可知，糖基轉移酶和糖苷水解酶在GD2107株的碳水化合物活性酶中占比較大(68.1%)。特定的糖苷水解酶可使結構多糖酶解，這是碳水化合物跨環境轉運的基本過程，而糖基轉移酶是參與合成天然糖基結構的CAZy，高效的糖基轉移酶對天然產物的糖基化至關重要，會影響其藥理特性，從而對宿主的臨床治療起指導作用，也為生物技術和制藥開辟了廣闊的可能性[33]。

毒力因子是衡量微生物毒性的重要指標，是引起宿主發病的重要因素[34]。GD2107株有48個毒力因子得到注釋，其中毒力因子中脂寡糖和莢膜的基因最多(70.8%)。由黏性表層組成的莢膜是APP的基本結構成分和致病因子，在感染過程中為保護細菌免受宿主防御提供了一種免疫學機制[35]。脂寡糖是革蘭氏陰性菌細胞壁上的一種結構成分，同脂多糖結構相似，生物功能也相近。其他被注釋的毒力因子還包括溶血素和黏菌素，與APP產生相互依賴的代謝生態系統相關，與蛋白酶相關的毒力因子具有防止吞噬作用和誘導膿腫形成的膠囊的作用，水解酶和脂多糖等毒力因子與降解免疫球蛋白和補體成分有關[36]。在毒力基因中，由某些基因共同調控形成的生物膜構成了抗生素耐藥的第一道防線。

APP的多重耐藥表現通常與其同時攜帶多種耐藥基因有關。Bossé等[33]和Liu等[37]研究顯示，除了大環內酯類藥物外，全基因組測序技術可準確地預測APP對所測試的抗菌藥的耐藥性。經全基因組測序技術分析發現，在GD2107株中存在22種耐藥基因。四環素是最早用于臨床治療PCP的主要藥物，APP對四環素的耐藥現象也已被廣泛發現[38]。本試驗藥敏檢測結果和測序到的耐藥基因(tetO基因等)同樣顯示,GD2107株對四環素類藥物耐藥。Chang等[39]報道位于質粒上的耐藥基因blaROB-1可介導APP菌株間阿莫西林的耐藥性。GD2107株攜的耐藥基因(floR基因)也位于質粒上，推測可能與APP自身的耐藥性和不同菌株間的傳播有關，這增加了臨床防控APP等細菌感染的難度。Kwong等[40]對肺炎克雷伯菌研究顯示，抗生素的耐藥基因在人類攜帶的分離株中也很常見，這些感染通常缺乏與侵襲性疾病相關的基因；COG注釋結果中GD2107株存在大量與基因重組、基因轉移相關的基因。這提示毒力與抗性基因的整合可能會導致變異強毒株的出現。因此，未來在疫苗研發時，除了需要充分了解菌種的生物學信息外，也需敲除其自身所攜帶的耐藥基因和毒力基因等，以防止或減少細菌水平基因轉移和轉化，保證疫苗產品的安全高效[41]。

本研究通過構建GD2107株的基因ApxⅣ系統進化樹發現，GD2107株與來自中國的APP分離株(CP063424.1)處在同一分支，進化關系最近；同國外的APP分離株進化關系較遠，較符合疾病的流行規律，這與Donà等[42]研究APP的流行具有一定的地域流行性這一結論一致，為該病的流行病學調查提供了一定參考。

4 結 論

本研究完成了對菌株GD2107的全基因組測序與生物信息學分析，全面認識了該菌株基因組的結構和功能并探究了與耐藥和致病機制中的相關基因，結果表明菌株GD2107呈多重耐藥性且攜帶多種毒力和耐藥基因，進化關系分析發現該菌株與國內分離株親緣關系最近，具有一定的地域流行性。本研究為預防PCP的疾病流行和探索APP的致病機制提供了參考。