李才武， 張 蘭， 張志忠,3， 張貴權(quán)， 魏榮平， 吳代福， 瞿春茂，李 鳳， 金森燕， 董 超， 程建斌， 楊長江， 何勝山， 劉小強(qiáng)， 黃 炎

(1.大熊貓國家公園珍稀動(dòng)物保護(hù)生物學(xué)國家林業(yè)和草原局重點(diǎn)實(shí)驗(yàn)室(中國大熊貓保護(hù)研究中心)邛崍山瀕危野生動(dòng)植物保護(hù)生物學(xué)國家長期科研基地, 都江堰 611830；2.四川大學(xué)生命科學(xué)學(xué)院 動(dòng)物疫病防控與食品安全四川省重點(diǎn)實(shí)驗(yàn)室, 成都 610065；3.國家林業(yè)和草原局野生動(dòng)植物保護(hù)司, 北京 100714； 4.四川臥龍自然保護(hù)區(qū)管理局, 臥龍 623006)

1 引 言

大熊貓是我國國寶級(jí)動(dòng)物,其繁殖率較低，是易危動(dòng)物之一[1].因此，在人工飼養(yǎng)環(huán)境下，提高大熊貓繁殖率十分重要.有研究報(bào)道顯示，某些生殖道微生物可直接或間接影響動(dòng)物繁殖率[2].另外，對(duì)人和動(dòng)物的陰道微生態(tài)的研究表明，陰道微生態(tài)的失衡及微生物感染與不孕存在密切關(guān)系[3]，如人類的外陰陰道炎[4]、宮頸炎[5]、奶牛的臨床陰道炎[6]等.近年來，已有對(duì)豬[7-8]、牛[9]、馬[10]等動(dòng)物的陰道菌群的研究報(bào)道，除此之外，對(duì)大熊貓生殖道微生物菌群也開展了研究.本研究團(tuán)隊(duì)在前期研究中利用微生物分離鑒定及高通量測序[11]揭示了正常大熊貓的陰道和子宮的細(xì)菌微生物組以及大熊貓陰道真菌微生物的多樣性與豐度[12-13]，但是使用細(xì)菌分離培養(yǎng)方法能培養(yǎng)出的細(xì)菌非常有限，而高通量測序方法僅能針調(diào)查部分細(xì)菌和真菌，不能了解其他微生物的組成情況.所以利用這兩個(gè)方法不能夠全面的反映大熊貓生殖道中的微生物菌群[14].

宏基因組測序技術(shù)能對(duì)樣本中所有的微生物多樣性、種群結(jié)構(gòu)、進(jìn)化關(guān)系、功能活性等進(jìn)行全面分析.在本研究中，我們采用宏基因組測序方法對(duì)生活在不同地區(qū)、不同年齡的大熊貓陰道樣本進(jìn)行了微生物全面分析，有助于大熊貓生殖道感染性疾病的防治，為科學(xué)的保護(hù)、繁育大熊貓?zhí)峁?shù)據(jù)支持.

2 材料和方法

2.1 材 料

2018年，根據(jù)大熊貓生活的不同區(qū)域和年齡，采集14個(gè)大熊貓陰道樣本.

2.2 方 法

2.2.1 樣本采集方法 使用長無菌拭子采集妊娠期的大熊貓生殖道(陰道)分泌物[12]，然后將每個(gè)拭子置于無菌的5 mL螺旋蓋管中，運(yùn)送樣品至四川省動(dòng)物疫病防控與食品安全重點(diǎn)實(shí)驗(yàn)室，-80 ℃冷凍保存，直至提取DNA時(shí)取出.

根據(jù)生活地區(qū)和年齡分為兩個(gè)大組：地區(qū)組(WLQ、BFX、DJY)和年齡組(A1、A2、A3).如表1 所示.根據(jù)不同的年齡將這14個(gè)樣品分為3組：其中5～10歲的為A1組(SZD1、SZD2、SZD3、SZD4、SZD5、SZD6、SZD7)；11～16歲的為A2組(SZD8、SZD9、SZD10、SZD11、SZD12)；17～23歲的為A3組(SZD13、SZD14).根據(jù)不同的采樣點(diǎn)將其分成三組.將來自臥龍的5個(gè)陰道樣品分成WLQ組(SZD1、SZD2、SZD8、SZD13、SZD14)，將來自雅安碧峰峽的6個(gè)陰道樣本分成BFX組(SZD4、SZD5、SZD9、 SZD10 、SZD11、 SZD12).同時(shí)將來自都江堰的3個(gè)陰道樣本分成DJY組(SZD3、SZD6、SZD7).

表1 14個(gè)大熊貓陰道樣品分組信息

2.2.2 樣品DNA提取 將拭子從冰箱中拿出來，放在冰上解凍.待解凍完全后，每個(gè)樣品加入2 mL的PBS磷酸緩沖液，然后將拭子渦旋震蕩.隨后使用康為世紀(jì)生物科技有限公司的DNA提取試劑盒提取樣本全基因組DNA，提取方法按照試劑盒實(shí)驗(yàn)操作說明書執(zhí)行.使用1%瓊脂糖凝膠電泳(AGE)分析 DNA 的純度和完整性，同時(shí)使用Qubit 對(duì) DNA 濃度進(jìn)行精確定量.

2.2.3 宏基因組測序及分析 將合格的DNA樣品送至北京諾禾致源科技股份有限公司進(jìn)行宏基因組測序，分析樣品中含有的各種微生物的相對(duì)豐度、功能基因和抗性基因.組裝的基因組或重疊群序列主要使用的是MetaGeneMark、CD-HIT、Bowtie2軟件進(jìn)行的分析.物種注釋主要使用genes 與各功能數(shù)據(jù)庫進(jìn)行比對(duì):DIAMOND軟件將 Unigenes 與從 NCBI 的 NR(Version: 2018.01) 數(shù)據(jù)庫中抽提出的細(xì)菌(Bacteria)、真菌(Fungi)、古菌(Archaea)和病毒(Viruses)序列進(jìn)行比對(duì)(blastp，evalue ≤ 1e-5).功能注釋主要使用DIAMOND軟件將 Unigenes 與各功能數(shù)據(jù)庫進(jìn)行比對(duì)(blastp，evalue ≤ 1e-5).抗性基因注釋主要使用CARD數(shù)據(jù)庫(v2.0.1)提供的Resistance Gene Identifier (RGI)軟件將 Unigenes 與CARD數(shù)據(jù)庫進(jìn)行比對(duì)(RGI內(nèi)置blastp，采用bitscore值對(duì)比對(duì)結(jié)果進(jìn)行評(píng)分).

技術(shù)參數(shù)是利用Illumina技術(shù)測序平臺(tái)和雙末端測序(Paired-End)，建立的文庫類型為350b小片段文庫，測序策略是Illumina PE150，每個(gè)樣品的測序數(shù)據(jù)量(raw data)為6G.

3 結(jié) 果

3.1 測序數(shù)據(jù)概述

本研究采用 Illumina HiSeq 測序平臺(tái)測序，共獲得 94 662.28 Mbp 的原始數(shù)據(jù)(Raw Data)(平均數(shù)據(jù)量為 6 761.59 Mbp)，經(jīng)過質(zhì)控得到 94 455.81 Mbp 的有效數(shù)據(jù)(Clean Data)(平均數(shù)據(jù)量為 6 746.84 Mbp)，經(jīng)過單樣品組裝及混合組裝后，共得到 477 921 089 bp 的 Scaftigs.對(duì)各樣品及混合組裝的結(jié)果，采用 MetaGeneMark 軟件進(jìn)行基因預(yù)測，共得到 583 524 個(gè)開放閱讀框(ORFs)(平均為 38 902)，經(jīng)過去冗余后，共獲得 412 312 個(gè)ORFs，總長為 233.7 Mbp，其中完整基因的個(gè)數(shù)為 74 805，所占比例為 18.14%.非冗余基因集與 MicroNR 庫進(jìn)行 blastp 比對(duì)，運(yùn)用 LCA 算法進(jìn)行物種注釋，注釋到屬和門的比例分別為 84.21% , 96.43%.使用 DIAMOND 軟件對(duì)非冗余基因集進(jìn)行常用功能數(shù)據(jù)庫注釋(e-value ≤10～5)，有 9 978(2.42%) 個(gè) ORFs 比對(duì)到 CAZy 數(shù)據(jù)庫，259 596(62.96%) 個(gè) ORFs 比對(duì)到 KEGG 數(shù)據(jù)庫，240 050(58.22%) 個(gè) ORFs 比對(duì)到 eggNOG 數(shù)據(jù)庫.非冗余基因集與抗性基因數(shù)據(jù)庫(CARD)進(jìn)行注釋(e-value ≤10～30)，有299 個(gè)基因比對(duì)到 CARD 數(shù)據(jù)庫.

3.1.1 測序原始數(shù)據(jù)處理 采用 Illumina測序平臺(tái)測序獲得的原始數(shù)據(jù)(Raw Data)存在一定比例低質(zhì)量數(shù)據(jù)，為了保證后續(xù)分析的結(jié)果準(zhǔn)確可靠，需要對(duì)原始的測序數(shù)據(jù)進(jìn)行預(yù)處理，獲取用于后續(xù)分析的有效數(shù)據(jù)(Clean Data).具體處理步驟如下：

1) 去除所含低質(zhì)量堿基(質(zhì)量值≤38)超過一定比例(默認(rèn)設(shè)為 40 bp)的 reads；

2) 去除 N 堿基達(dá)到一定比例的 reads(默認(rèn)設(shè)為10 bp)；

3) 去除與 Adapter 之間 overlap 超過一定閾值(默認(rèn)設(shè)為 15 bp)的 reads；

4) 如果樣品存在宿主污染，需與宿主序列進(jìn)行比對(duì)，過濾掉可能來源于宿主(默認(rèn)采用 Bowtie2 軟件，參數(shù)設(shè)置: --end-to-end, --sensitive, -I 200, -X 400)的 reads；

測序數(shù)據(jù)預(yù)處理統(tǒng)計(jì)結(jié)果見表2.

3.1.2 Metagenome組裝 使用SOAP denovo組裝軟件進(jìn)行組裝分析，具體組裝步驟如下：

1) 經(jīng)過預(yù)處理后得到 Clean Data，使用 SOAP denovo 組裝軟件進(jìn)行組裝分析(Assembly Analysis)；

2) 對(duì)于單個(gè)樣品，組裝時(shí)選取 K-mer=55進(jìn)行組裝，得到該樣品的組裝結(jié)果；組裝參數(shù)：-d 1, -M 3, -R, -u, -F；

3) 將組裝得到的 Scaffolds 從 N 連接處打斷，得到不含 N 的序列片段，稱為 Scaftigs (i.e., continuous sequences within scaffolds)；

4) 將各樣品質(zhì)控后的 CleanData 采用 Bowtie2 軟件比對(duì)至各樣品組裝后的 Scaftigs 上，獲取未被利用上的 PE reads：比對(duì)參數(shù)：--end-to-end, --sensitive, -I 200, -X 400；

5) 將各樣品未被利用上的 reads 放在一起，選取 K-mer=55進(jìn)行混合組裝，其他組裝參數(shù)與單樣品組裝參數(shù)相同；

6) 將混合組裝的 Scaffolds 從 N 連接處打斷，得到不含 N 的 Scaftigs 序列；

7) 對(duì)于單樣品和混合組裝生成的 Scaftigs，過濾掉 500 bp以下的片段，并進(jìn)行統(tǒng)計(jì)分析和后續(xù)的分析(表3).

表2 14個(gè)樣品的測序數(shù)據(jù)處理結(jié)果

表3 各樣品組裝結(jié)果Scaftigs 基本信息統(tǒng)計(jì)(>=500 bp)

3.2 物種注釋

3.2.1 樣品的相對(duì)豐度分析 各組以及各樣品中微生物在門水平上最大相對(duì)豐度排名前10的物種分別見表4、5、6.如表所示，在這些樣品中占主導(dǎo)地位的門主要為變形菌門(Proteobacteria), 擔(dān)子菌門(Basidiomyta), 厚壁菌門(Firmicutes)和放線菌門(Actinobacteria),分別占樣品的平均值為34.75%、21.53%、5.52%和3.29%.同時(shí)在屬的水平上分析，其結(jié)果為：以假單胞菌(21.90%)最多，其次為鏈球菌(3.47%)、冷桿菌(1.89%)和變形桿菌(1.38%).從表3、4中可以看出，不同的分組，其優(yōu)勢菌門不同(將物種相對(duì)豐度大于1%的規(guī)定為優(yōu)勢門)，在WLQ組中優(yōu)勢門為變形菌門和厚壁菌門；在BFX組中的優(yōu)勢門為變形菌門、擔(dān)子菌門、厚壁菌門、放線菌門和擬桿菌門；在DJY組中的優(yōu)勢門為變形菌門、厚壁菌門、放線菌門、擬桿菌門和衣原體.在A1組中的優(yōu)勢門為變形菌門、厚壁菌門、放線菌門和擬桿菌門；在A2組中的優(yōu)勢門為變形菌門、擔(dān)子菌門、厚壁菌門和放線菌門；在A3組中的優(yōu)勢門為變形菌門和厚壁菌門.擔(dān)子菌門只在BFX組和A2組中為優(yōu)勢菌門，衣原體在DJY組豐度較高.

表4 在門水平上14個(gè)樣品的物種相對(duì)豐度表

表5 在門水平上WLQ、BFX、DJY這個(gè)三個(gè)組的物種相對(duì)豐度表

表6 在門水平上AGE1、AGE2、AGE3這三個(gè)組的物種相對(duì)豐度表

3.2.2 基于物種豐度的聚類分析 為了進(jìn)一步的研究不同樣品的相似性，對(duì)樣品進(jìn)行聚類分析，構(gòu)建樣品的聚類樹.采用Bray-Curtis 距離系統(tǒng)聚類法，從基因在各樣品中的豐度出發(fā)，以Bray-Curtis距離矩陣進(jìn)行樣品間聚類分析，并將聚類結(jié)果與各樣品在屬水平上的物種相對(duì)豐度整合進(jìn)行展示.(如圖1所示)結(jié)果顯示：總體上，WLQ、DJY、和BFX這三個(gè)地區(qū)組和A1、A2和A3這三個(gè)年齡組都可以單獨(dú)聚類，但在年齡組方面的聚類更加明顯.說明生活地區(qū)和年齡的差異對(duì)大熊貓陰道的微生物菌群都有一定的影響，但年齡對(duì)大熊貓陰道菌群的影響更加明顯.

圖1 在屬水平上基于 Bray-Curtis 距離的聚類樹

A圖代表不同地區(qū)差異(WLQ組、BFX組、DJY組)的基于 Bray-Curtis 距離的聚類樹.B圖代表不同年齡差異(A1組、A2組、A3組)的基于 Bray-Curtis 距離的聚類樹.圖中，左側(cè)是 Bray-Curtis 距離聚類樹結(jié)構(gòu)；右側(cè)的是各樣品在屬水平上的物種相對(duì)豐度分布圖.

Fig.1 Clustering tree based on Bray-Curtis distance at genus level

圖2 在基因功能水平上基于Bray-Curtis距離的聚類樹

A圖代表不同地區(qū)差異(WLQ組、BFX組、DJY組)的基于 Bray-Curtis 距離的聚類樹，B圖代表不同年齡差異(A1組、A2組、A3組)的基于 Bray-Curtis 距離的聚類樹.圖中，左側(cè)是 Bray-Curtis 距離聚類樹結(jié)構(gòu)；右側(cè)的是各樣品在門水平上的物種相對(duì)豐度分布圖.

Fig.2 Clustering tree based on bray-curtis distance at gene function level

3.3 功能注釋

功能數(shù)據(jù)庫注釋結(jié)果概述：原始去冗余后的預(yù)測基因共有 412 312個(gè)，有259 596(62.96%)個(gè)基因能夠比對(duì)上 KEGG 數(shù)據(jù)庫，其中，有143 757(34.87%) 個(gè)基因能夠比對(duì)上數(shù)據(jù)庫中的 7 791個(gè) KEGG ortholog group(KO)；有240 050(58.22%)個(gè)基因能夠比對(duì)上 eggNOG 數(shù)據(jù)庫；有 9 978(2.42%) 個(gè)基因能夠比對(duì)上 CAZy 數(shù)據(jù)庫.

3.3.1 樣品的聚類分析 根據(jù)Unigenes在KEGG數(shù)據(jù)庫六大代謝通量中的注釋情況并基于Bary-Curtis距離繪制了聚類樹(圖2A、B).從圖中可以看出樣品在基因功能上的Bary-Curtis距離與在物種上的Bary-Curtis距離相似，結(jié)果一致：相比于按地區(qū)分組來說，按照年齡分組的樣品更能體現(xiàn)各組間的功能差異.表明年齡的差異對(duì)大熊貓陰道菌群的影響較大.

3.3.2 基于功能豐度的降維分析 為研究地區(qū)組和年齡組中功能組成的相似度，在KEGG水平上對(duì)功能豐度值進(jìn)行了PCA分析.如圖3所示，在地區(qū)組中，WLQ組和DJY組的功能具有一定的相似性，BFX組則與WLQ組和DJY組存在一定的差異性.在年齡組中，A1組中的樣品在功能上存在一定的相似性，而A2組和A3組中的樣品功能差異較明顯.總的來說，盡管個(gè)體間差異較明顯，但每一個(gè)組都表現(xiàn)為獨(dú)立的分布.

PCA結(jié)果展示，橫坐標(biāo)表示第一主成分，百分比則表示第一主成分對(duì)樣品差異的貢獻(xiàn)值；縱坐標(biāo)表示第二主成分，百分比表示第二主成分對(duì)樣品差異的貢獻(xiàn)值；圖中的每個(gè)點(diǎn)表示一個(gè)樣品，同一個(gè)組的樣品使用同一種顏色表示.A圖代表不同地區(qū)差異(WLQ組、BFX組、DJY組)的PCA 分析，B圖代表不同年齡差異(A1組、A2組、A3組)的PCA 分析.

Fig.3 PCA analysis based on functional abundance

A、B、C圖按照地區(qū)分組的分別為WLQ組、BFX組和DJY組中的抗性基因物種歸屬分析圈圖.圈圖中的內(nèi)圈為ARO 的物種分布情況，外圈為組內(nèi)所有樣品基因的物種分布情況.

Fig.4 Resistance gene species assignment analysis circle map and resistance mechanism and species Overview circle diagram

3.4 抗性基因注釋

3.4.1 抗性基因的物種歸屬關(guān)系 根據(jù)樣品的基因物種注釋結(jié)果，可以分析各耐藥基因所對(duì)應(yīng)的物種信息，通過比較耐藥基因所屬的物種信息可以直觀的反映攜帶耐藥基因的優(yōu)勢菌群.如圖4所示，在WLQ組中(圖4A)，變形菌門、厚壁菌門為占主導(dǎo)地位的兩個(gè)菌種，分別占菌量的30%、6%，而各自攜帶的耐藥基因比例分別為62%、4%；在BFX組中(圖4B)，變形菌門、擔(dān)子菌門，放線菌門為占主導(dǎo)地位的三個(gè)菌種，分別占菌量的58%、6%、4%，而各自攜帶的耐藥基因比例分別為40%、12%、4%；在DJY組中(圖4C)，變形菌門，擬桿菌門為占主導(dǎo)地位的兩個(gè)菌種，分別占菌量的16%、5%，而各自攜帶的耐藥基因比例分別為63%、3%.這種不對(duì)稱的結(jié)果表明，變形菌門比厚壁菌門和擬桿菌門更易攜帶耐藥基因，但各組中占主導(dǎo)地位的菌種與在物種注釋中各組的優(yōu)勢菌門結(jié)果一致.

圖5 ARO分布情況熱圖Fig.5 ARO distribution heat map注：橫軸為樣品名稱，右側(cè)縱軸為抗性基因類型 ARO 名稱，黑色表示樣品中含有該 ARO，白色表示樣品中沒有該 ARO.

3.4.2 抗性基因類型的分布 為了反映各抗性基因類型(ARO)在各樣品中的分布情況，繪制ARO分布的黑白熱圖(圖5)，從圖中可以看出，沙門氏菌耐藥基因APH3-lb分布最廣，存在于所有的14個(gè)樣品中.同時(shí)還發(fā)現(xiàn)耐藥基因mexQ只存在樣品SZD2中，耐藥基因vanWl只存在樣品SZD4中.而與之相反的是，耐藥基因APH6-ld、AAC6-lla、Sul1分別廣泛存在于除了樣品SZD6、SZD1、SZD3的其他樣品中.

4 討 論

4.1 樣品的物種注釋

宿主中微生物種群數(shù)量的不同會(huì)對(duì)宿主造成不同的影響.比如高豐度的厚壁菌門或者低豐度的擬桿菌門會(huì)對(duì)宿主的肥胖有一定的影響[15].同時(shí)微生物多樣性的變化從側(cè)面反映了機(jī)體的正常情況.患者腸道中的微生物多樣性會(huì)出現(xiàn)顯著的降低，并且厚壁菌門和梭菌屬所占的比例也顯著降低[16].

在我們所研究樣本中，占主導(dǎo)地位的門主要為變形菌門、擔(dān)子菌門、厚壁菌門和放線菌門.這與前人關(guān)于正常大熊貓生殖道菌群的結(jié)果一致.但檢測出來衣原體，并且衣原體感染可以影響生物的繁殖能力[17].因此推測衣原體可能對(duì)大熊貓的繁殖能力具有潛在的不利影響.在物種樣品聚類分析的這一方面發(fā)現(xiàn)生活地區(qū)和年齡的差異對(duì)大熊貓陰道的微生物菌群都有一定的影響，但年齡的差異對(duì)大熊貓陰道菌群的影響更加明顯.

4.2 樣品的功能注釋

從分類和功能組成的角度來看，腸道微生物群可能與許多復(fù)雜的疾病有關(guān).有研究發(fā)現(xiàn)泌乳期奶牛中瘤胃微生物群會(huì)隨年齡和生理狀態(tài)的變化而變化，并且通過對(duì)其進(jìn)行宏基因測序發(fā)現(xiàn)其基因序列對(duì)碳水化合物，蛋白質(zhì)代謝等近50%的代謝具有潛在貢獻(xiàn)[18].

同時(shí)我們發(fā)現(xiàn)大熊貓生殖道的基因除了未被注釋的功能序列之外，主要功能是新陳代謝，其次是環(huán)境信息處理.這與其他有關(guān)微生物菌群宏基因組的研究結(jié)果一致[19].說明大多數(shù)生物中的基因主要表現(xiàn)為新陳代謝這一功能，維持生物的正常生命活動(dòng).此外，根據(jù)物種和功能豐度的聚類分析可以得出：樣品按照年齡分組進(jìn)行的聚類分析能更加的明顯的看出各組之間的不同，因此與地區(qū)差異相比，年齡差異對(duì)大熊貓生殖道中微生物菌群的物種和功能豐度的影響較大.這與在物種注釋中關(guān)于樣品聚類分析得到的結(jié)果一致.

4.3 樣品的抗性基因注釋

抗生素能夠在防治動(dòng)物疾病方面發(fā)揮重要作用，然而最近幾年抗生素的濫用導(dǎo)致許多動(dòng)物體內(nèi)產(chǎn)生了具有抗生素耐藥性的細(xì)菌，比如在綿羊瘤胃中發(fā)現(xiàn)了達(dá)托霉素和粘菌素[20].同時(shí)在豬、雞和人類的糞便中也發(fā)現(xiàn)了高水平的四環(huán)素和紅霉素等[21].

在抗性基因的注釋中，雖然并沒有發(fā)現(xiàn)特別的耐藥基因，但我們發(fā)現(xiàn)耐藥基因所屬物種的優(yōu)勢菌門與樣品在物種注釋中所呈現(xiàn)的優(yōu)勢菌門一致，均為變形菌門、厚壁菌門、擔(dān)子菌門和放線菌門.并且變形菌門比厚壁菌門和擬桿菌門更易攜帶耐藥基因.從抗性基因數(shù)目來看，地區(qū)組和年齡組中含量最多的為BFX組和A2組，并且由于BFX組中的樣品年齡均在A2組的范圍中，因此推測年齡的差異對(duì)各樣品中抗性基因的數(shù)目有較大的影響.