劉海霞，韓大勇，閆 偉，朱愛文，周明夏，王步忠，吳 昊，耿昕之

（1.江蘇農牧科技職業(yè)學院，江蘇泰州 225300； 2.江蘇西來原生態(tài)農業(yè)有限公司，江蘇泰州 225328； 3.江蘇省阜寧縣農業(yè)農村局，江蘇鹽城 224400）

宏基因組學是以高通量測序技術為基礎的組學研究方法，宏基因組測序技術可以全方位測定分析樣品中所含微生物的種群結構、進化關系、功能活性和多樣性等。郭玉梅等（2022）利用宏基因組測序技術研究發(fā)熱病人咽部微生物菌群，鑒定出4079個物種，歸屬于97個門，951個屬，發(fā)現(xiàn)試驗組較對照組存在消化鏈球菌、擬桿菌、支原體等物種差異，β-內酰胺類酶類基因為耐藥基因。宏基因組技術可以通過挖掘微生物基因組信息，分析其功能及代謝通路，進而挖掘菌群的種群結構、功能特性及種群間的相互關系（彭璽等，2022；孫波等，2021；許悅等，2020；馬海霞等，2015；賀紀正等，2008）。本研究采用宏基因組測序技術，對綿羊呼吸道樣本進行微生物菌群多樣性的測定分析，指導防治綿羊呼吸道感染性疾病，為綿羊健康養(yǎng)殖提供數(shù)據(jù)支撐。

1 材料與方法

1.1 材料

1.1.1 樣品采集 采集8只試驗用綿羊鼻涕樣本（健康羊、呼吸道癥狀病羊各4只）。鼻拭子輕輕蘸取鼻腔中間段涕樣，采集完成后，將鼻拭子插入預置2 mL PBS溶液的進口凍存管中，液氮速凍，-20℃保存。

1.1.2 主要儀器、試劑 Illumina測序儀（6000，NovaSeq），高效拭子基因組DNA提取試劑盒（WH0215，離心柱型，百奧萊博），TIANSeq DNA Library Prep Kit文庫構建試劑盒（NG103，天根生物），TruSeq PE Cluster Kit v3-cBot-HS（PE-401- 3001，Illumina）。

1.2 方法

1.2.1 樣品DNA提取 采用基因組DNA提取試劑盒提取樣本中的基因組DNA，然后分析提取的基因組DNA的純度和完整性（瓊脂糖凝膠電泳），再DNA濃度和質量進行檢測（Nanodrop 2000）。

1.2.2 文庫構建與檢測 用超聲波破碎儀將檢測合格的DNA樣品隨機打斷，片段長度約為350 bp，修復補平打斷的DNA片段末端，在3’端添加堿基“A”，將片段末端轉換為黏性末端，利用堿基互補在黏性末端兩側添加包含Index序列的DNA接頭。磁珠篩選收集一定長度范圍內的目的片段，PCR擴增，在目的片段末端添加Index，構建完成測序文庫。使用Aglient 2100生物芯片分析系統(tǒng)檢測文庫片段大小，并使用Step One PlusTM Real-Time PCR System，以P5、P7接 頭作為引物進行QPCR定量。

1.2.3 上機測序 將檢測合格的不同測序文庫按照目標下機數(shù)據(jù)量的需求和有效濃度進行pooling，在Illumina測序平臺上運行雙端測序程序，得到150 bp的雙端測序reads。

1.2.4 數(shù)據(jù)分析 將測序得到的原始數(shù)據(jù)進行質控后得到有效數(shù)據(jù)。使用CD-HIT、MetaGeneMark、Bowtie2等軟件進行所有樣品的reads的Metagenome組裝、物種注釋、功能數(shù)據(jù)庫注釋、抗性基因注釋等。DIAMOND軟件將Unigenes與從NCBI NR數(shù)據(jù)庫中抽提出的微生物序列進行比對，genes與各功能數(shù)據(jù)庫進行比對進行物種注釋，DFAMOND軟件將Unigenes與各功能數(shù)據(jù)庫進行比對，進行功能注釋，CARD數(shù)據(jù)庫提供的Resistance GeneIdentifier（RGI）軟件將Unigenes與抗性基因數(shù)據(jù)庫（CARD）進行比對，完成抗性基因注釋。

2 結果與分析

2.1 測序結果概述本研究測序在Illumina NovaSeq測序平臺完成，共獲得625347.73 Mbp的原始數(shù)據(jù)，經過質控得到536874.34 Mbp的有效數(shù)據(jù)，混合組裝單樣品、多樣品，共得到989926098 bp的Scaftigs。使 用MetaGeneMark軟件進行基因預測，共得到561403個開放閱讀框（ORFs），去冗余后共獲得418423個ORFs，完整基因的個數(shù)為96776，所占比例為20.13%。非冗余基因集注釋到屬的比例為72.21%，非冗余基因集注釋到門的比例83.36%。使用DIAMOND軟件對非冗余基因集進行KEGG、eggNOG、CAZy、CARD等功能數(shù)據(jù)庫進行注釋，有572318個ORFs比 對到KEGG數(shù) 據(jù) 庫，560255個ORFs比對 到eggNOG數(shù) 據(jù) 庫，10112個ORFs比 對 到CAZy數(shù)據(jù)庫。通過CARD注釋，比對到CARD數(shù)據(jù)庫有299個基因。

2.2 物種注釋通過物種豐度分析，試驗組共得到103個 門，159個 綱，413個 目，840個 科，1789個屬，4790個種；對照組共得到99個門，159個綱，409個 目，836個 科，1780個 屬，4783個 種。試驗組與對照組比較發(fā)現(xiàn)，變形菌綱豐度對照組顯著高于試驗組；乳桿菌目豐度對照組顯著高于試驗組，而桿菌目豐度試驗組顯著高于對照組；葡萄球菌科豐度試驗組顯著高于對照組，鏈球菌科豐度對照組顯著高于試驗組；乳球菌屬豐度試驗組顯著高于對照組，葡萄球菌屬豐度對照組顯著高于試驗組；金黃色葡萄球菌豐度試驗組顯著高于對照組，格氏乳球菌豐度對照組顯著高于試驗組。綜合分析，試驗組致病菌豐度相對較高，對照組益生菌豐度相對較高。

2.3 功能注釋

2.3.1 樣品KEGG注釋結果 本研究發(fā)現(xiàn)，樣品中636981個不同基因被注釋，從屬于KEGG數(shù)據(jù)庫6大類46個KEGG通路中，分別是新陳代謝（187129，29.37%）、人類疾?。?73291，27.20%）、生命系統(tǒng)（169859，26.66%）、細胞進程（48558，7.62%）、環(huán)境信息處理（29505，4.63%）、遺傳信息處理（28639，4.52%）。富集差異表達基因最多的6個KEGG通路，分別是信號傳導（52127）、氨基酸代謝（42336）、免疫系統(tǒng)（27816）、轉運與分解代謝（20982）、傳染性疾?。?0029）、感知系統(tǒng)（17920）。環(huán)境信息處理類通路中的信號傳導通路、新陳代謝類通路中氨基酸代謝通路所含基因數(shù)最多。（如圖1）

圖1 基因KEGG代謝通路統(tǒng)計

2.3.2 樣品eggNOG數(shù)據(jù)庫注釋結果 eggNOG數(shù)據(jù)庫注釋發(fā)現(xiàn)，本試驗相對含量在前5位的功能大類是細胞內轉運、分泌和囊泡轉運，翻譯后修飾、蛋白質翻轉，復制、重組與修復，信號轉導機制，轉錄，基 因數(shù)分別為50128、24833、12680、12635、11206個。eggNOG功能分類統(tǒng)計結果（如圖2）

圖2 基因eggNOG功能分類

2.3.3 樣品CAZy數(shù)據(jù)庫注釋結果 經過與CAZy數(shù)據(jù)庫比對，本試驗樣品中共鑒定出6905個碳水化合物活性酶，包含3238個糖基轉移酶，2217個糖苷水解酶類，數(shù)量最多，占據(jù)樣品檢測數(shù)據(jù)總數(shù)的77%。1002個碳水化合物結合模塊，258個碳水化合物酯酶，112個輔助氧化還原酶，78個多糖裂合酶最少，占據(jù)樣品檢測數(shù)據(jù)總數(shù)的0.63%。

2.3.4 樣品抗性基因注釋結果 利用CARD數(shù)據(jù)庫進行抗生素抗性基因的預測，繪制ARO分布熱圖，可反映各抗性基因類型ARO在各樣品中的分布情況。如圖3可知，葡萄球菌和埃希氏桿菌屬耐藥基因分布最廣，存在于測定的8個樣品中，耐藥基因APH（3）-Ib只存在于樣 品D2中，耐 藥 基 因novA、antibiotic、msbA、Mycobacterium、macB廣泛存在于除樣品D1外的其他樣品中。耐藥基因smeR、aadA只存于B2和D1樣品中。

圖3 ARO分布熱圖

3 討論

3.1 試驗樣品的物種注釋隨著生物信息學和現(xiàn)代分子生物學的發(fā)展，宏基因組測序技術的應用領域更廣。為探究試驗樣本微生物組成特征，使用Diamond軟件將unigene序列與NR數(shù)據(jù)庫比對進行物種注釋，可獲得樣本在分類水平各物種組成及豐度信息，通過物種豐度分析反映試驗樣本微生物結構組成多樣性變化情況。李娟（2020）等預測到麥洼牦牛瘤胃微生物基因2065638個開放閱讀框（ORFs）完整基因數(shù)為333371個，舍飼組厚壁菌門、擬桿菌門、軟壁菌門、浮霉菌門相對放牧組豐度均增加，纖維桿菌門相對豐度舍飼組顯著低于放牧組。本研究發(fā)現(xiàn)，綿羊呼吸道微生物菌群多樣性特征明顯，在門水平上，試驗組優(yōu)勢菌群主要為厚壁菌門、擬桿菌門、變形菌門、放線菌門、螺旋體門，對照組優(yōu)勢菌群主要為厚壁菌門、擬桿菌門、放線菌門、梭桿菌門、變形菌門。

3.2 試驗樣品的功能注釋通過對基因在KEGG生物通路數(shù)據(jù)庫上進行注釋分析，可以查看基因參與的生物通路、體現(xiàn)的生物功能。對比分析eggNOG數(shù)據(jù)庫與基因翻譯獲得的蛋白質序列，獲得相應的蛋白質功能注釋信息。將基因序列與CAZy數(shù)據(jù)庫進行比較，得到碳水化合物活性酶類的物種來源、酶功能EC分類、基因、蛋白序列及蛋白結構信息。Pitta等（2016）利用宏基因測序研究泌乳期奶牛瘤胃微生物群，發(fā)現(xiàn)基因序列對碳水化合物蛋白質代謝等近50%的代謝有重要作用。本研究從KEGG數(shù)據(jù)庫分別注釋到636981個基因，環(huán)境信息處理類通路中的信號傳導通路所含基因數(shù)最多，其次是新陳代謝類通路中氨基酸代謝通路所含基因量。CAZy數(shù)據(jù)庫比對發(fā)現(xiàn)，樣品中共鑒定出碳水化合物活性酶6905個，其中，糖基轉移酶3238個和糖苷水解酶類2217個，數(shù)量最多，占據(jù)樣品檢測數(shù)據(jù)總數(shù)的77%。

3.3 試驗樣品的抗性基因注釋濫用抗生素會導致動物體內產生抗生素耐藥性細菌，Ma等（2015）在豬、雞和人類糞便中檢測出了含量較高的四環(huán)素和紅霉素等抗生素。Hitch等（2018）在綿羊瘤胃中檢測到達托霉素和粘菌素。利用CARD數(shù)據(jù)庫進行抗生素抗性基因預測，能反映各抗性基因類型ARO物種中的分布情況。本研究發(fā)現(xiàn)，葡萄球菌和埃希氏桿菌屬耐藥基因分布最廣，存在于檢測的8個樣品中，耐藥基因APH3-1b只存在于樣品D2中，耐藥基因novA、antibiotic、msbA、Mycobacterium、macB廣 泛 存在于除了樣品D1以外的其他樣品中，耐藥基因smeR、aadA只存于B2和D1樣品中。

4 結論

本研究通過宏基因組技術對綿羊呼吸道微生物群落分進行分析，發(fā)現(xiàn)厚壁菌門、變形菌門、擬桿菌門、放線菌門為綿羊呼吸道優(yōu)勢菌群。葡萄球菌和埃希氏桿菌屬耐藥基因在試驗動物中分布較廣。病羊鼻腔菌群中的支原體顯著高于健康羊，藥敏試驗表明，羊肺炎支原體對泰妙菌素敏感，為臨床首選藥物。