楊李玲，郭迎香，位紅穎，王 萌，方藝菲，朱 鵬，姜敬哲

1. 天津農(nóng)學(xué)院，天津 300384

2. 中國水產(chǎn)科學(xué)研究院南海水產(chǎn)研究所/農(nóng)業(yè)農(nóng)村部南海漁業(yè)資源開發(fā)利用重點實驗室，廣東 廣州 510300

3. 廣州美格基因科技有限公司，廣東 廣州 510005

4. 從化區(qū)農(nóng)業(yè)農(nóng)村局，廣東 廣州 510925

5. 上海美吉生物醫(yī)藥科技有限公司，上海 201321

6. 上海海洋大學(xué)，上海 201306

牡蠣屬于軟體動物門、雙殼綱、珍珠貝目，廣泛分布于世界各地，其營養(yǎng)價值高，是全球養(yǎng)殖范圍最廣、產(chǎn)量最高的海水養(yǎng)殖品種。中國作為最大的牡蠣生產(chǎn)國，2019 年牡蠣產(chǎn)量達(dá)到 8 259×104t，占世界總產(chǎn)量的85.3% (數(shù)據(jù)來源：http://www.fao.org)。然而，受環(huán)境變化、陸源污染及種質(zhì)退化等因素影響，牡蠣病害的發(fā)生日趨頻繁，對其養(yǎng)殖業(yè)的健康可持續(xù)發(fā)展造成負(fù)面影響[1]。

1972年Farley等[2]首次報告了牡蠣病毒性疾病病例，證明在高溫條件下牡蠣死亡事件與皰疹病毒存在關(guān)聯(lián)，并將這種病毒命名為1型牡蠣皰疹病毒 (Ostreid Herpesvirus-1, OsHV-1)。該病毒在長牡蠣 (Crassostrea gigas) 中具有高致病性，同時對其他多種雙殼類也造成嚴(yán)重影響[3-4]。除OsHV-1外，其他報道過的與牡蠣相關(guān)的病毒還包括：乳多空病毒科 (Papovaviridae) 病毒，可能導(dǎo)致牡蠣“卵囊炎”，引起卵子和配子細(xì)胞肥大；鰓壞死病毒，屬虹彩病毒科 (Iridoviridae)，可能是導(dǎo)致20世紀(jì)60 年代后期葡萄牙牡蠣 (C. angulata) 種群大規(guī)模死亡的主要原因；另外，披膜病毒 (Togaviridae)、呼腸孤病毒 (Reoviridae) 和小 RNA 病毒 (Picornavridae) 在貝類宿主中也有多次報道[5]。牡蠣等濾食性貝類在過濾大量海水獲取有機顆粒物營養(yǎng)的同時，也會將水體中的其他物質(zhì)如病毒、污染物等帶入體內(nèi)，并通過外套膜進(jìn)行濃縮。因此，如果牡蠣養(yǎng)殖水域受到人類活動影響，牡蠣體內(nèi)可能還會富集諸如諾如病毒 (Norovirus, NV)、甲型肝炎病毒 (Hepatitis A Virus, HAV) 和星狀病毒 (Astrovirus, AV) 等引起人類消化道疾病的食源性病毒，但這些病毒并不是牡蠣的病原[6]。病毒基因組是疾病診斷、流行病檢測和病原生物學(xué)研究的基礎(chǔ)。由于無脊椎動物沒有能夠穩(wěn)定傳代的細(xì)胞系，除了OsHV-1和一些食源性病毒外，其他牡蠣相關(guān)病毒的鑒定證據(jù)主要基于病理學(xué)和電子顯微鏡觀察，很少能夠拿到基因組序列[7]。因此，基因組數(shù)據(jù)的匱乏是牡蠣病害研究進(jìn)展緩慢的一個主要原因。高通量測序和宏病毒組技術(shù)的運用，克服了傳統(tǒng)病毒學(xué)研究對宿主細(xì)胞培養(yǎng)的依賴，極大提高了新病毒的發(fā)現(xiàn)和鑒定效率。近些年，已在海洋環(huán)境、無脊椎動物、低等脊椎動物及人類等相關(guān)研究中得到了廣泛應(yīng)用，鑒定到大量的新型病毒[8-11]，極大地擴展了人們對病毒世界的認(rèn)知。

本研究基于前期工作所獲得的華南沿海多地養(yǎng)殖牡蠣宏病毒組數(shù)據(jù)，進(jìn)一步運用生物信息學(xué)工具對其進(jìn)行深入挖掘，發(fā)現(xiàn)并鑒定到5條新型牡蠣相關(guān)圓環(huán)病毒基因組，為促進(jìn)牡蠣病毒及貝類病害相關(guān)研究提供參考。

1 材料與方法

1.1 數(shù)據(jù)來源與篩選

課題組前期工作通過對中國華南沿海多地養(yǎng)殖牡蠣進(jìn)行宏病毒組測序，得到約25億條測序原始序列 (reads)[12]。對測序數(shù)據(jù)進(jìn)行 Fastp (V0.20.0)[13]質(zhì)控去除低質(zhì)量和接頭序列，用Megahit (V1.2.9)[14-15]將reads組裝成重疊群 (contig) 片段，以NCBI非冗余蛋白質(zhì)數(shù)據(jù)庫為參考，用Diamond[16](V0.9.24.125) 進(jìn)行比對注釋，并用 Megan 6[17]進(jìn)一步對注釋結(jié)果進(jìn)行分類，篩選其中5條基因組完整的、疑似為圓環(huán)病毒科 (Circoviridae) 的病毒序列 (Contig ID：ML2-k141_11943、ML2-k141_27933、QZd1-k141_2132、T5S3-k141_267932、ZHd1-k141_220676)進(jìn)一步深入分析。

1.2 開放閱讀框預(yù)測與比對

基于 NCBI ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/) 進(jìn)行開放閱讀框 (Open reading frame, ORF) 預(yù)測，其中 Minimal ORF length (nt) 選取 75，Genetic code選取 standard，ORF start codon to use 選擇 Any sense codon，忽略 nested ORFs。用預(yù)測到的ORF進(jìn)行Blastp比對，其中數(shù)據(jù)庫選擇nr蛋白數(shù)據(jù)庫，e為0.000 001。對于有比對結(jié)果的ORF，取一致性最高的蛋白用tblastn方法反向比對本研究中病毒的基因組序列，驗證5條ORF序列是否完整，并用DNAstar進(jìn)行翻譯，得到最終的蛋白序列。

1.3 基因組進(jìn)化樹的構(gòu)建及比較基因組分析

將5條病毒全基因組序列提交至VIPTree(V1.9)[18](https://www.genome.jp/viptree/)，選擇核酸類型為ssDNA，選擇宿主病毒類型為Eukaryote，選擇基因預(yù)測的遺傳密碼為Eukaryotic，構(gòu)建基因組進(jìn)化樹。利用MAFFT[19]將5條病毒基因組序列與近源的病毒基因組序列進(jìn)行比對，并使用TrimAI[20]進(jìn)一步修剪以去除對齊不明確的區(qū)域，利用 MEGA[21]構(gòu)建 Neighbor-joining Tree，設(shè)置 Bootstrap Replications為 1 000，選擇 substitution model為 p-distance，設(shè)置 Site Coverage Cutoff為50%。

根據(jù)與5條病毒復(fù)制酶蛋白一致性最高的參考蛋白序列 (YP_006281010.1)，計算參考蛋白與5條病毒序列之間的 Amino Acid Identity (AAI)。下載YP_006281010.1對應(yīng)的基因組序列 (NC_017843.2)，利用 Orthologous Average Nucleotide Identity Tool (OAT) 對 NC_017843.2 和 5 條病毒序列計算兩兩序列的 Average Nucleotide Identity (ANI)，得到兩兩序列的ANI，使用R的ggcorrplot包對數(shù)值進(jìn)行可視化。同時在VIPtree中進(jìn)行基因組兩兩比較分析：設(shè)置 positioning of each sequence 為auto，colors為 default，ticks of genomes為 shown，vertical dashed lines為 shown，representation of tBLASTx hits為 normal。

1.4 復(fù)制酶蛋白進(jìn)化樹分析

因5條病毒復(fù)制酶蛋白序列Blatsp比對的結(jié)果存在重疊，因此僅以其中1條序列(ML2-k141_11 943)比對結(jié)果前50的序列和5條序列一起構(gòu)建系統(tǒng)進(jìn)化樹，建樹過程參考1.3。

1.5 復(fù)制酶蛋白結(jié)構(gòu)域分析

在NCBI上下載ARI44308.1[22]，用Jalview對ARI44308.1、YP_006281010.1和這5條病毒復(fù)制酶蛋白序列進(jìn)行結(jié)構(gòu)位點分析，使用ClustalW進(jìn)行比對，下載保存SVG格式文件。用SMART (http://smart.embl-heidelberg.de/) 軟件進(jìn)行結(jié)構(gòu)域預(yù)測，同時采用Swiss model進(jìn)行三維結(jié)構(gòu)預(yù)測。

1.6 病毒豐度分析

為了計算每個病毒的相對豐度值，首先合并5 條病毒基因組序列用 salmon (V0.13.1) index[23]命令生成參考基因組數(shù)據(jù)集，然后再利用salmon quant[23]命令將全部牡蠣病毒組文庫的clean reads逐一映射 (Mapping) 到參考基因組上，統(tǒng)計映射上的reads數(shù)，根據(jù)調(diào)整后的FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) 計算公式計算各病毒的相對豐度值。

式中：G為比對的reads數(shù)；T為測序reads總數(shù)(×106)；L為基因組長度 (kb)。

2 結(jié)果

2.1 病毒注釋結(jié)果

通過對中國華南沿海多地養(yǎng)殖牡蠣樣品進(jìn)行宏病毒組測序，消除低質(zhì)量序列并組裝，共獲得3 375 091 條非冗余 contig序列，進(jìn)一步與 NCBI的NR庫進(jìn)行blastp比對和分類注釋最終獲得731 536條疑似病毒來源的contig，篩選其中5條基因組相對完整的、被注釋為圓環(huán)病毒的基因組序列進(jìn)一步分析，具體信息見表1。

表1 牡蠣相關(guān)圓環(huán)病毒基因組信息Table 1 Genome information of oyster-related circoviruses

2.2 基因組水平進(jìn)化樹分析

為了明確病毒的分類，用5條基因組序列與全部單鏈DNA病毒一起構(gòu)建基因組水平的進(jìn)化樹(VIPTree，圖1)?；蚪M樹中內(nèi)圈顏色代表病毒科(Family)，主要包括圓環(huán)病毒科、雙生病毒科(Geminiviridae)、類雙生病毒科 (Genomovirdae) 和細(xì)項圈病毒科 (Anellovridae)。從科水平分類看，圓環(huán)病毒科在整個單鏈DNA基因組樹中占據(jù)絕大部分。候選的5條序列均聚在了一個小的分支 (紅色五角星標(biāo)識，圖1)，且在圓環(huán)病毒大分支內(nèi)，說明它們屬于圓環(huán)病毒科成員。圖1的最外圈代表病毒的宿主，包括脊索動物門、節(jié)肢動物門、綠藻門和軟體動物門。從宿主來源看，候選病毒分支周圍主要是節(jié)肢動物宿主和其他來源宿主，而節(jié)肢動物與軟體動物同屬無脊椎動物。

圖1 單鏈DNA病毒基因組進(jìn)化樹 (VIPtree)Fig. 1 Phylogenetic tree of single-stranded DNA virus genome (VIPtree)

為了解5條病毒序列的進(jìn)化關(guān)系，將5條病毒基因組序列與近源的病毒基因組序列進(jìn)行比對。5 條病毒序列與 NC 017843.2、NC 027797.1、NC 030466.1均聚在了一支，其中除了NC 030468.1來源于海水樣品外，NC 017843.2、NC 027797.1、NC 030466.1分別來源于橈足動物、珊瑚、大猩猩樣品(圖2)。該聚類結(jié)果與5條候選序列的來源背景基本一致 (海洋、動物)。值得注意的是，QZd1-k141_2132與T5S3-k141_267932在全基因組水平關(guān)系最接近。

圖2 牡蠣相關(guān)圓環(huán)病毒基因組與近源基因組序列構(gòu)建的系統(tǒng)發(fā)育樹Fig. 2 Phylogenetic tree of five oyster-related circoviruses and relative viral genome sequences

2.3 全基因組序列和復(fù)制酶蛋白序列一致性分析

用候選序列中預(yù)測到的復(fù)制酶蛋白序列在NCBI上進(jìn)行Blastp比對，發(fā)現(xiàn)5條病毒序列均與YP_006281010.1的蛋白序列相似性較高。該序列注釋為圓環(huán)病毒復(fù)制酶蛋白，該蛋白來源于一種海洋節(jié)肢動物——夏唇角水蚤 (Labidocera aestiva)。用其基因組序列 (NC_017843.2) 和5條序列進(jìn)行兩兩比較，NC_017843.2基因組序列與5條候選序列的ANI介于0%～72%，其中最高的是QZd1-k141_2132 (63%)，最低的是 ML2-k141_27933 (0%，圖3)。5條病毒序列兩兩間的一致性為0%～94%，其中QZd1-k141_2132與T5S3-k141_267932的最高(94%)，ML2-k141_27933和其余4條病毒序列均為0%。根據(jù)ICTV規(guī)定的圓環(huán)病毒科在種水平劃分閾值80%來看，5條序列可能屬于4個不同的新種，且均與NC_017843.2不同，即QZd1-k141_2132與T5S3-k141_267932可能屬于同一種，而其余3條序列屬于另外3個種。復(fù)制酶蛋白序列一致性分析結(jié)果一致。

圖3 圓環(huán)病毒基因組序列間的ANI (a) 和復(fù)制酶蛋白序列間的AAI (b)Fig. 3 ANI value among genomic sequences (a) and AAI value among replication proteins of circoviruses (b)

2.4 比較基因組分析

比較基因組分析表明，夏唇角水蚤圓環(huán)病毒基因組序列 (NC_017843.2) 與5條候選基因組序列均在復(fù)制酶蛋白區(qū)域 (即YP_006281010.1) 具有相似性 (圖4)，但在另一個ORF，即疑似衣殼蛋白區(qū)域，均未見明顯的相似性。可見候選病毒的衣殼蛋白屬于新型的衣殼蛋白，與牡蠣宿主細(xì)胞受體的特異性識別相關(guān)。

圖4 牡蠣相關(guān)圓環(huán)病毒比較基因組分析圖Fig. 4 Comparative genome analysis of oyster-related circoviruses

2.5 Neighbor-joining 進(jìn)化樹構(gòu)建

因5條病毒復(fù)制酶蛋白序列Blatsp比對的結(jié)果存在重疊，因此僅以其中1條序列(ML2-k141_11943) 比對結(jié)果的前50的序列和5條序列一起構(gòu)建系統(tǒng)進(jìn)化樹。ML2-k141_11943、ZHd1-k141_220676、QZd1-k141_2132、T5S3-k141_267932 這4個病毒和YP_006281010.1、GAC77785.1、QNG41123.1、AXQ66182.1形成了一個分支，其中YP_006281010.1、QNG41123.1和 AXQ66182.1的來源是動物宏基因組 (具體物種不明)，GAC77785.1來源于海洋沉積物 (圖5)。ML2-k141_27933序列與AXH76045.1、AXH74760.1形成了一個分支，該序列來源同樣是動物宏基因組 (具體物種不明)。

圖5 牡蠣相關(guān)圓環(huán)病毒復(fù)制酶蛋白進(jìn)化樹Fig. 5 Phylogenetic tree of replication proteins of oyster-related circoviruses

2.6 功能分析

為查明復(fù)制酶ORF中是否有相關(guān)功能位點和結(jié)構(gòu)域，用ARI44308.1和YP_006281010.1作為參考序列與5條候選序列進(jìn)行全局比對 (圖6)。雖然5條候選序列在典型環(huán)狀病毒的復(fù)制酶結(jié)構(gòu)域(HUH phage replication sub-family) 處與參考蛋白序列差異較大。但SMART分析發(fā)現(xiàn)，除ZHd1-k141_220676外，余下的4條候選序列與2條參考序列都鑒定到明確的HUH結(jié)構(gòu)域。復(fù)制酶蛋白結(jié)構(gòu)域三維結(jié)構(gòu)預(yù)測表明，5條候選序列，尤其是ML2-k141_11943、ML2-k141_27933和T5S3-k141_267932，與YP_006281010.1復(fù)制酶蛋白的空間結(jié)構(gòu)存在較高的相似性 (圖7)。

圖6 牡蠣相關(guān)圓環(huán)病毒復(fù)制酶蛋白序列比對Fig. 6 Sequence alignment of oyster-related circoviruses

圖7 牡蠣相關(guān)圓環(huán)病毒復(fù)制酶蛋白三維結(jié)構(gòu)預(yù)測圖Fig. 7 Three-dimensional structure prediction of replication protein of oyster-related circoviruses

2.7 病毒豐度分析

為了計算每個病毒的相對豐度值，用salmon統(tǒng)計映射上的reads數(shù)。在牡蠣測序文庫中均檢測到了5條備選病毒序列的reads，F(xiàn)PKM豐度各不相同 (表2)。整體來看，T5S3-k141_267932的分布最廣泛，ML2-k141_11943在特定樣品中的豐度最高。具體來看，ML2-k141_11943和ML2-k141_27933 在 ML1 中豐度最高 (148.67 和 2 080.33)；QZd1-k141_2132和T5S3-k141_267932在ChYJ1509Da中豐度最高 (542.87和76.93)；ZHd1-k141_220676在 ChZH1511Da 豐度最高 (55.17)。

表2 牡蠣相關(guān)圓環(huán)病毒豐度分析Table 2 Analysis of abundance of oyster-related circovirus virus

3 討論

近年來，我國沿海地區(qū)養(yǎng)殖貝類大規(guī)模死亡事件頻發(fā)，涉及的貝類品種、養(yǎng)殖面積和區(qū)域范圍不斷擴大，使貝類養(yǎng)殖業(yè)遭受重大經(jīng)濟損失[24]。傳統(tǒng)病毒學(xué)研究嚴(yán)重依賴宿主細(xì)胞的培養(yǎng)，由于缺乏可穩(wěn)定傳代的無脊椎動物細(xì)胞培養(yǎng)技術(shù)，加之相關(guān)研究投入有限，貝類病毒研究進(jìn)展十分緩慢[25]。由于高通量測序和病毒組技術(shù)的應(yīng)用，近年幾項海洋無脊椎動物病毒研究成果的發(fā)表，加深了人們對于海洋動物病毒的了解。2016年，Shi等[26]報道了關(guān)于無脊椎動物RNA病毒組的重要研究成果，對包括多種貝類生物在內(nèi)的約400種無脊椎動物 (包括牡蠣) 組織混合后的約200個文庫進(jìn)行了測序，發(fā)現(xiàn)了1 400多種新型RNA病毒，極大地豐富了RNA病毒的多樣性和種類。Rosani和Gerdol[27]從長牡蠣 (C. gigas)和地中海貽貝 (Mytilus galloprovincialis) 宿主轉(zhuǎn)錄組數(shù)據(jù)中鑒定了7個幾乎完整的RNA病毒基因組，屬于小核糖核酸病毒目 (Picornavirales)，其中有6種是全新的病毒。Ian等[28]對健康的與被感染的海星進(jìn)行了比較研究，鑒別出1種細(xì)小病毒科 (Parvoviridae) 的疑似致病病原，并且明確其在浮游生物和海洋沉積物中也廣泛存在。Cui[29]從不同海域共采集了3門6綱58種的海洋無脊椎動物樣品 (包括牡蠣)，并利用宏轉(zhuǎn)錄組的測序方法對不同物種的RNA病毒組進(jìn)行研究，鑒定出涵蓋9個病毒家族 (Durnavirales、Totiviridae、Bunyavirales、Chuviridae、Picornavirales、Flaviviridae、Hepelivirales、Solemoviridae、Tombusviridae) 的363個RNA病毒，其中的315個與已知的RNA病毒相似度較低，被認(rèn)為是全新的RNA病毒，其中牡蠣來源病毒4個。牡蠣等海洋無脊椎動物體內(nèi)可能存在著非常豐富的病毒種類。但前述研究多基于轉(zhuǎn)錄組測序數(shù)據(jù)，發(fā)現(xiàn)的新病毒大部分為RNA病毒，對DNA病毒的相關(guān)報道則十分有限。

圓環(huán)病毒科病毒是已知基因組最小的一類動物致病性DNA病毒，呈二十面體對稱球形，無囊膜，直徑為17～22 nm，基因組由大約1.7～2.3 kb的單鏈環(huán)狀DNA組成。該家族成員種類繁多，可感染各類動物宿主，包括蝙蝠、嚙齒動物、畜禽動物和人類等，其中了解比較清楚的是感染畜禽動物和人類的種類，如雞傳染性貧血病毒 (Chicken Infectious anemia)、豬圓環(huán)病毒 (Porcine circovirus)、鴿圓環(huán)病毒 (Pigeon circovirus)等[30-31]。截止到本研究前，在牡蠣中尚未有圓環(huán)病毒的報道。經(jīng)查詢Virus-Host數(shù)據(jù)庫確認(rèn) (數(shù)據(jù)截止至2021年5月)，目前已知的所有圓環(huán)病毒宿主都是動物界內(nèi)的兩側(cè)對稱動物，加之基因組和蛋白序列進(jìn)化樹表明，它們與夏唇角水蚤圓環(huán)病毒最為接近，同屬無脊椎動物宿主范疇。這意味著這些圓環(huán)病毒不可能是牡蠣體內(nèi)的某些共生細(xì)菌或原生動物的病毒。另一方面，圓環(huán)病毒的宿主特異性是由其衣殼蛋白序列決定的[32]。而本研究發(fā)現(xiàn)5條病毒基因組序列中，僅在ZHd1-k141_220676序列中比對到了疑似衣殼蛋白序列，且與水蚤圓環(huán)病毒衣殼蛋白無相似性。因此，這5個新型圓環(huán)病毒的宿主是水蚤等節(jié)肢動物的可能性不高。當(dāng)然，現(xiàn)有證據(jù)尚不足以確認(rèn)牡蠣是其真正宿主，還需要結(jié)合水體浮游蚤類篩查、擴大養(yǎng)殖牡蠣調(diào)查范圍及牡蠣人工感染實驗來進(jìn)一步明確其真正的宿主。

根據(jù)ICTV標(biāo)準(zhǔn)，同一圓環(huán)病毒科物種的成員在其整個基因組中應(yīng)該具有大于80%的核苷酸同一性[33]。本研究中的5條候選序列與已知圓環(huán)病毒在整個基因組水平的相似性低于80%，在衣殼蛋白方面未觀察到明顯的相似性。因此，本研究鑒定的5個病毒可能與現(xiàn)有圓環(huán)病毒科家族成員有較大差異，具體的分類界定還有待發(fā)現(xiàn)更多家族成員序列后的分析。

病毒是海洋中最豐富的生命體，以貝類為代表的軟體動物門也是海洋中最大的動物類群。但目前對于兩者的交叉領(lǐng)域——貝類病毒的了解仍十分有限。宏病毒測序方法已經(jīng)在眾多生物和環(huán)境樣品中得到廣泛應(yīng)用，突顯了高通量測序技術(shù)在發(fā)現(xiàn)新病毒方面的巨大潛力[12]。然而，病毒組在未知序列解析和新病毒分類方面還有很多不足尚待解決。本研究發(fā)現(xiàn)的5個牡蠣圓環(huán)病毒的分類地位、家族成員以及宿主范圍、致病性、流行情況等方面還有很多工作要做。后續(xù)研究應(yīng)該在發(fā)揮病毒組技術(shù)優(yōu)勢的基礎(chǔ)上，同時配合以經(jīng)典的病毒學(xué)研究方法，獲得更為全面和深入的數(shù)據(jù)，進(jìn)而為提升為牡蠣等貝類病害的防控水平、支撐相關(guān)產(chǎn)業(yè)發(fā)展發(fā)揮作用。