陳嘉源，施勁松，丘棟安，劉 暢，李 鑫，趙 強，阮吉壽，高 山*

(1.南開大學 生命科學學院, 天津300071; 2.東部戰(zhàn)區(qū)總醫(yī)院, 南京 210016;3.英國諾丁漢特倫特大學 生物科學系, 諾丁漢 NG11 8NS;4.南開大學 醫(yī)學院, 天津300071;5.南開大學 數(shù)學科學學院, 天津 300071)

冠狀病毒屬(Coronavirus)的病毒是具有包膜的正鏈RNA病毒。由于2003年嚴重急性呼吸綜合征(Severe Acute Respiratory Syndrome, SARS)和2012年中東呼吸綜合征(Middle East Respiratory Syndrome, MERS)的爆發(fā)，冠狀病毒逐漸成為病毒學領域的一個研究熱點。2019年12月，中國武漢報道了冠狀病毒引起的肺炎，其臨床癥狀與2003年爆發(fā)的SARS不同，因此推斷該病毒可能是冠狀病毒的一個新變種。當前，冠狀病毒溯源相對比較清晰：2003年，Enserink等在國際上首次提出了SARS冠狀病毒(SARS Coronavirus, SARS-CoV)來自蝙蝠，而后通過果子貍傳播給人類[1]；再后來的病毒溯源研究都是根據(jù)蝙蝠-果子貍-人途徑回溯，并重點研究各種蝙蝠體內(nèi)的SARS樣(SARS-like)冠狀病毒與已知SARS冠狀病毒的進化關系；特別是2013年，SARS冠狀病毒(GenBank: AY274119.3)被溯源到中華菊頭蝠(Rhinolophussinicus)[2]。盡管已有足夠的冠狀病毒基因組數(shù)據(jù)公開，有關冠狀病毒跨物種傳播和宿主適應性方面的研究卻一直難以突破，特別是無法解釋冠狀病毒變異快、宿主多且具有較強的宿主適應性等特點。其中一個最主要原因就是冠狀病毒基因組較大，多個基因功能未知。例如，SARS冠狀病毒基因組(GenBank: AY274119.3)的全長是29 751 bp，比HIV-1病毒全長三倍還要多。另外，冠狀病毒基因組注釋中基因命名不統(tǒng)一、且大部分注釋的蛋白質序列都來自預測，無實驗驗證，難以保證正確性。高通量測序的迅速發(fā)展，使大量冠狀病毒基因組序列得以公開。然而，當前大部分的冠狀病毒基因組研究都是簡單使用全基因組或某個病毒結構基因的序列，后一種做法主觀性較大。簡單使用全基因組序列進行進化分析和重要功能或致病位點分析，不僅工作量巨大，而且大量與研究目的無關的基因組區(qū)域會對數(shù)據(jù)分析產(chǎn)生干擾。因此，不同于簡單使用全基因組序列的其它研究，我們于2018年在國際上首次提出分子功能與進化分析相結合的研究思想，并應用于冠狀病毒基因組的研究。分子功能與進化分析相結合，即僅使用與特定功能相關的基因組特征對一定變異范圍內(nèi)的序列進行進化分析[3]。這樣不僅可以簡化進化分析、提高信噪比，使分析結果更可靠，而且可以用于指導功能實驗。2018年, 南開大學高山等在國際上首次報道了SARS冠狀病毒中存在互補回文小RNA(complemented palindromic Small RNAs, cpsRNAs)的現(xiàn)象，并確定了Beta冠狀病毒B亞群(BB冠狀病毒)與功能相關的兩個重要的基因組特征(Genomic feature)[3]，即一段22 bp的互補回文序列(命名為Nankai complemented palindrome)與其所在的編碼區(qū)(命名為Nankai CDS)。在上述研究中，對Nankai complemented palindrome進行突變分析的結果與使用Nankai CDS進行進化分析的結果一致支持SARS冠狀病毒的跨物種傳播途徑是蝙蝠-果子貍-人。在分子功能與進化分析相結合的研究思想的指導下，本研究使用Nankai complemented palindrome和Nankai CDS對新發(fā)布的2019新型冠狀病毒基因組(GenBank: MN908947)進行分析以期準確溯源，并對BB冠狀病毒跨物種傳播和宿主適應性進行初步研究。

1 數(shù)據(jù)與方法

在前期研究中，我們通過互補回文小RNA的發(fā)現(xiàn)從SARS冠狀病毒MA-15株(GenBank: DQ497008)中確定了BB冠狀病毒中與功能相關的兩個重要的基因組特征，即一段22 bp的互補回文序列與其所在的編碼區(qū)。這個編碼區(qū)(DQ497008: 25689-26153)是BB冠狀病毒基因組中一段高度保守的序列，在已知SARS冠狀病毒基因組(如DQ497008)中只對應一個讀碼框ORF 3b；但在其它BB冠狀病毒基因組中普遍對應多個讀碼框，因此ORF 3b無法代表其它BB冠狀病毒基因組中的這段序列。本研究將22 bp的互補回文序列與其所在的編碼區(qū)(可能是一個或多個讀碼框)分別命名為Nankai complemented palindrome和Nankai CDS。DQ497008和2003年爆發(fā)的SARS冠狀病毒基因組(GenBank: AY274119.3)含有相同的Nankai complemented palindrome和Nankai CDS，因此在本研究中唯一代表SARS冠狀病毒基因組。在前期研究中[3]，我們的研究對象僅限于BB冠狀病毒(本研究亦是如此)，并進行如下操作[3]：提取NCBI GenBank數(shù)據(jù)庫中全部BB冠狀病毒的完整基因組(Complete genome)序列(以下簡稱病毒序列或序列)；去除與Nankai CDS相似度在0.999以上的序列后，留下11條病毒序列(GenBank: JX993987、JX993988、GQ153539、GQ153540、GQ153542、DQ071615、DQ412042、DQ412043、AY515512、AY572034和DQ497008)用于進一步研究。本研究增加一條新發(fā)布的2019新型冠狀病毒基因組的序列(GenBank: MN908947)和一條來自蝙蝠的序列(GenBank: MG772934)，其Nankai CDS與2019新型冠狀病毒Nankai CDS相似度最高(92.52%)。另外，由于本論文接收時2019新型冠狀病毒基因組序列(GenBank: MN908947)尚未公開，實際使用的是來自互聯(lián)網(wǎng)公開的序列(見致謝)。在本研究中，13條序列根據(jù)其宿主分為五組用于進一步研究，這五組命名為SARS(DQ497008)、果子貍Civet(AY515512和AY572034)、2019新型冠狀病毒2019-nCoV(MN908947)、蝙蝠群體bat2(MG772934)和蝙蝠群體bat1(MG772934之外8條來自蝙蝠的序列)。在本研究中，進化分析使用軟件PHYLIP v3.69，統(tǒng)計與作圖使用軟件R v2.15.3[4]。其它數(shù)據(jù)處理包括：病毒蛋白質預測和核苷酸三連子分析只考慮一個起始密碼子ATG、三個終止密碼子(TAA、TAG和TGA)和正向三個讀碼框；蛋白質預測不考慮小于20個氨基酸的蛋白質；最后，考慮四個起始密碼子(ATG、GTG、CTG和TTG)和三個終止密碼子重復蛋白質預測和核苷酸三連子分析以檢驗結果的可靠性(見結果)。

2 結 果

在經(jīng)典定義中，DNA回文(DNA palindrome)和DNA互補回文(DNA complemented palindrome)序列都叫做DNA回文序列，而且大多數(shù)情況下，DNA回文序列的經(jīng)典定義與我們新定義中的互補回文序列相同(例如ATCGCGAT)。根據(jù)我們的新定義[3]，回文序列(DNA或RNA)要求序列(一條鏈)從5’到3’方向與從3’到5’方向的讀取結果相同(例如ATCGGCTA)?；匚暮突パa回文序列應該加以區(qū)別，主要基于三個原因：(1)我們的前期研究發(fā)現(xiàn)了回文和互補回文序列可能具有不同的生物學意義；(2)在動物病毒基因組中，回文和互補回文序列具有不同的分布特征；(3)經(jīng)典定義是從DNA水平定義的，我們的定義考慮了單鏈DNA和RNA(特別針對病毒)情況。2004年，Chew等在國際上首次報道了SARS冠狀病毒基因組中較短互補回文序列(新定義)的分布特征[5]，其主要發(fā)現(xiàn)包括兩點：(1)在所有分析的冠狀病毒基因組中，長度為4 bp的互補回文序列出現(xiàn)頻率顯著較低；(2)長度為6 bp的互補回文序列僅在SARS冠狀病毒(而不是所有冠狀病毒)基因組中出現(xiàn)頻率顯著較低。因此，該研究的結論是SARS冠狀病毒基因組含有較少的6 bp互補回文序列有可能利于該病毒有效規(guī)避宿主細胞內(nèi)的某些防御機制而有利于病毒存活；該研究同時還發(fā)現(xiàn)了兩個較長(14 bp以上)的互補回文序列分別是TCTTTAACAAGCTTGTTAAAGA和TAAAATTAATTTTA。但是，僅僅根據(jù)概率模型計算得到的此類回文或互補回文序列數(shù)量巨大，絕大部分無法與分子功能關聯(lián)，因此沒有進一步研究的價值。直到2018年，我們偶然發(fā)現(xiàn)了互補回文小RNA恰好來自這個22 bp的互補回文序列(命名為Nankai complemented palindrome)，特別是發(fā)現(xiàn)了根據(jù)這個序列預測的RNA二級結構在BB冠狀病毒基因組中高度保守(見圖1A)。RNA水平的證據(jù)是將這個序列與分子功能建立關聯(lián)的關鍵[3]。在前期研究中，我們對11個BB冠狀病毒(見數(shù)據(jù)與方法)的Nankai complemented palindrome進行突變分析的結果與使用其所在的編碼區(qū)(命名為Nankai CDS)進行進化分析的結果一致支持SARS冠狀病毒的跨物種傳播途徑是蝙蝠-果子貍-人(見圖1B紅色方框內(nèi))[3]。

圖1 2019新型冠狀病毒的起源與可變翻譯Fig.1 Origin and alternative translation of the 2019 novel coronavirus genome

注：A：Nankai CDS是BB冠狀病毒基因組中一段高度保守的序列，包括一段22 bp的互補回文序列(用*表示)。Nankai CDS還包括一個8或11 bp的可變區(qū)，分別屬于兩種不同長度(是465或468 bp，用#表示)的Nankai CDS；B：進化樹構建使用13條去除可變區(qū)的Nankai CDS。使用多種方法進行進化分析的結果一致，這里只展示鄰接(Neighbor Joining，NJ)法的結果：2019新型冠狀病毒源自蝙蝠某個群體(藍色方框內(nèi))，但與SARS冠狀病毒差異巨大(紅色方框內(nèi))；用P1到P17(見圖1C)可以對13個BB冠狀病毒進行基因分型(小括號內(nèi))，其結果與進化分析的結果一致。C：可變翻譯導致從Nankai CDS可以預測出至少17個蛋白質，分別命名為P1到P17。這里展示的是使用一個起始密碼子ATG和三個終止密碼子進行蛋白質預測的結果，使用四個起始密碼子(ATG、GTG、CTG和TTG)和三個終止密碼子會預測更多的蛋白質，但不改變可變翻譯相關結論。

Notes:A: The 22-bp complemented palindrome (named Nankai complemented palindrome) and the CDS (named Nankai CDS) have a high degree of evolutionary conservation in genomes of betacoronavirus subgroup B(BB coronavirus). An 8- or 11-bp long variable region in Nankai CDS resides in two types of Nankai CDS. The first type has a length of 465 bp, while the second type has a length of 468 bp.B: Consistent results were obtained by different phylogenetic analysis methods. Here, we only show the result using the NJ (Neighbor Joining) method. The 2019 novel coronavirus (in blue box) with large differences from the SARS coronavirus (in red box), may originate from BB coronaviruses in Chinese horseshoe bats. Phylogenetic analysis and genotyping (in parentheses) of 13 viruses using 17 putative proteins (Fig.1C) achieved consistent results. C: The putative proteins (named as P1 to P17) were predicted from the Nankai CDS using ATG as the start codon, and TAA, TAG and TGA as stop codons. The prediction using ATG, GTG, CTG and TTG as start codons, and TAA, TAG and TGA as stop codons did not change all conclusions in the present study.

在本研究中，所有13條Nankai CDS(見數(shù)據(jù)與方法)都含有一個完整的Nankai complemented palindrome和一個緊鄰其3'端的長度為8或11 bp的可變區(qū)(Variable region)。使用MN908947(2019新型冠狀病毒)中Nankai complemented palindrome預測的RNA二級結構仍然能形成穩(wěn)定的莖環(huán)(Stem-loop)結構，再次驗證了Nankai complemented palindrome在BB冠狀病毒基因組中的保守性[3]。Nankai CDS中的可變區(qū)導致了長度變異(見圖1A)：MG772934(來自蝙蝠)和MN908947(2019新型冠狀病毒)的Nankai CDS全長為468 bp(含11 bp可變區(qū))；其余11條序列的Nankai CDS全長為465 bp(含8 bp可變區(qū))。由于無法確定每條Nankai CDS的可變區(qū)的變異方式，進行進化分析前，13條Nankai CDS去除可變區(qū)后即可對齊，長度統(tǒng)一到457 bp。使用多種方法(UPGMA、鄰接法、最大簡約法、最大似然法和貝葉斯法)對13條Nankai CDS進行進化分析的結果(見圖1B)一致支持：(1)MN908947(2019新型冠狀病毒)與MG772934(來自蝙蝠)在進化樹的一個分支內(nèi)(見圖1B藍色方框內(nèi))，說明2019新型冠狀病毒源自蝙蝠；(2)2019新型冠狀病毒(見圖1B藍色方框內(nèi))與SARS冠狀病毒(見圖1B紅色方框內(nèi))在不同分支，而且距離很大，這一結果與兩者臨床癥狀差異一致；(3)SARS冠狀病毒的跨物種傳播途徑“蝙蝠-果子貍-人”(見圖1B紅色方框內(nèi))與2018年分析結果完全相同[3]，證明了我們方法的可靠性。

Nankai CDS是BB冠狀病毒基因組中一段高度保守的序列(見數(shù)據(jù)與方法)，靠近其3'端的99個堿基(編碼33個氨基酸)極端保守，不僅沒有起始或終止密碼子突變，也極少有其它突變，具有重要的研究和應用(如引物設計)價值。Nankai CDS在來自SARS和果子貍的beta冠狀病毒基因組中是一個完整的讀碼框，編碼蛋白P16(見圖1C)，但在來自其它宿主的冠狀病毒基因組中(不包括3'端的99個堿基)出現(xiàn)了大量的起始或終止密碼子突變(見表1)，并導致Nankai CDS出現(xiàn)讀碼框的多樣性，即可變翻譯。由于去冗余后的完整基因組數(shù)據(jù)較少，需要從兩個方面排除這些突變是來自測序錯誤等導致的假陽性。一方面，用13條Nankai CDS比對到NCBI NT數(shù)據(jù)庫，看是否有其它冠狀病毒序列支持這些突變。比對結果顯示，大部分序列上的突變都能得到NT數(shù)據(jù)庫中其它數(shù)據(jù)支持，只有三個例外：JX993988在Nankai CDS上的三個突變(C80A、T80A和T254A)，MG772934(來自蝙蝠)和MN908947(2019新型冠狀病毒)在Nankai CDS上共有的四個突變(T182A、C182A、T237A和T236G)沒有其它數(shù)據(jù)支持。JX993988代表的病毒類型比較特殊，與其它來自蝙蝠的病毒差異較大(見圖1B)，這里不做進一步研究。另外一方面，統(tǒng)計13條Nankai CDS中核苷酸三連子出現(xiàn)次數(shù)，次數(shù)最高的前14位依次是：ATG(123次)、ACT(106次)、AAG(98次)、TTG(83次)、CAG(77次)、CAA(71次)、GTG(68次)、AAA(66次)、CTA(65次)、ACA(60次)、TTA(59次)、ATT(58次)、CGA(57次)和CTG(55次)。核苷酸三連子統(tǒng)計結果顯示：前14位三連子出現(xiàn)次數(shù)之和超過三連子總數(shù)一半以上(1 046次/2 016次)；6個三連子與起始密碼子ATG只有0或1個核苷酸的差異；8個三連子與終止密碼子(TAA、TAG和TGA)只有1個核苷酸的差異。統(tǒng)計結果說明了通過1個核苷酸突變(點突變)導致起始或終止密碼子的丟失和獲得的概率非常高。根據(jù)以上兩個方面的分析結果，基本可以排除Nankai CDS中的大量起始或終止密碼子突變是假陽性(見表1)。

表1 Nankai CDS內(nèi)起始或終止密碼子突變Table 1 Mutations in start or stop codons in Nankai CDS

注：可變翻譯導致從Nankai CDS可以預測出至少17個蛋白質，分別命名為P1到P17(見圖1C)。所有突變均來自17個預測的蛋白質，PP表示預測的蛋白質(Putative Protein)的名稱；Start一列記錄的是起始密碼子的第一個堿基的位置；Stop一列記錄的是終止密碼子的最后一個堿基位置；Allele一列記錄所在位置的密碼子的全部類型；Mutation一列只記錄導致獲得起始或終止密碼子的突變，突變格式為“突變前堿基-突變位置-突變后堿基”；AA表示PP的長度，單位是氨基酸(Amino Acid)。這里展示的是只使用一個起始密碼子ATG和三個終止密碼子進行蛋白質預測的結果；使用四個起始密碼子(ATG、GTG、CTG和TTG)和三個終止密碼子會預測更多的蛋白質，但不改變可變翻譯相關結論。

Notes:The putative proteins (named as P1 to P17) were predicted from the Nankai CDS (Fig.1C). This table only includes the mutations in start or stop codons in 17 putative proteins. “PP” is Putative Protein; “Start” records the positions of the first nucleotides in start codons; “Stop” records the positions of the last nucleotides in stop codons; “Allele” records all types of codons on the positions; “Mutation” only records mutations which cause the acquisition of start or stop codons in the format “nucleotide before mutaion-position-nucleotide after mutaion”; “AA” is Amino Acid. The putative proteins were predicted from the Nankai CDS using ATG as the start codon, and TAA, TAG and TGA as stop codons. The prediction using ATG, GTG, CTG and TTG as start codons, and TAA, TAG and TGA as stop codons did not change all conclusions in the present study.

可變翻譯導致從Nankai CDS可以預測出至少17個蛋白質(Putative Protein)，分別命名為P1到P17(見圖1C)。用P1到P17可以對13個病毒進行基因分型(基于可變翻譯的基因分型)，例如，MN908947(2019新型冠狀病毒)可以表示為P10P11P12P13P17。我們發(fā)現(xiàn)基因分型相同或相近的病毒都在進化樹上的同一個分支內(nèi)(見圖1B)。例如，bat1分支的病毒的基因分型都由三個蛋白質組成；SARS和果子貍分支都是P16；MN908947(2019新型冠狀病毒)和MG772934(來自蝙蝠)都有P10。根據(jù)以上事實，可以得到以下推論：(1)BB冠狀病毒(不限于Nankai CDS內(nèi))存在大量的可變翻譯；(2)BB冠狀病毒變異快、多樣性高，僅僅13個病毒即可得到多達9種基因分型(見圖1BC)；(3)在SARS冠狀病毒從蝙蝠到人傳播過程中，Nankai CDS內(nèi)多個點突變積累才形成P16；(4)在2019新型冠狀病毒從蝙蝠到人傳播過程中，P15只需一個點突變(C355T)即形成P12和P13，P14只需兩個點突變(T170A和T191C)即形成P11和P17；(5)可變翻譯可用于(但不限于)BB冠狀病毒的基因分型，用P1到P17進行基因分型的結果與用Nankai CDS進行進化分析的結果一致。(6)BB冠狀病毒可能通過可變翻譯以適應不同宿主。

3 結 論

本研究證明了我們提出的分子功能與進化分析相結合的研究思想的可行性與可靠性，對生物大分子的功能或進化研究具有一定的指導意義，并對病毒或細菌等病原體的研究提供了新的路線和方向。因此，本研究除了解決2019新型冠狀病毒溯源等具體問題，更主要貢獻是提供了研究思想和方法。

1)使用Nankai CDS溯源分析的結果支持2019新型冠狀病毒源自蝙蝠，但與SARS冠狀病毒差異巨大；

2)從BB冠狀病毒可變翻譯中獲取的信息可應用于其快速檢測、基因分型、疫苗開發(fā)以及藥物設計；

3)BB冠狀病毒可能通過可變翻譯以適應不同宿主；

4)對BB冠狀病毒可變翻譯的研究將有助于研究該病毒的變異、宿主適應性、感染以及致病機制等問題。

4 討 論

病毒或細菌等原核生物的可變翻譯現(xiàn)象很早就已經(jīng)發(fā)現(xiàn)，但是其生物學意義尚未清楚，特別是缺乏基于較大數(shù)據(jù)的實證研究。我們在前期研究植物病毒[6]、昆蟲病毒[7]、人病毒[8]、哺乳動物病毒(特別是非洲豬瘟病毒[9])和布魯氏菌[10]的過程中，都沒有明顯觀察到可變翻譯現(xiàn)象，特別是BB冠狀病毒Nankai CDS中這樣大量而集中存在的可變翻譯(僅僅13個病毒即可得到多達9種基因分型)極其罕見。此次實證研究結果達到了基于可變翻譯的基因分型與進化分析的一致性，是一次較大突破，主要得益于以下幾點：(1)分子功能與進化分析相結合的研究思想得以應用；(2)有蝙蝠-果子貍-人這一跨物種傳播途徑作為方向性指導；(3)無意中發(fā)現(xiàn)Nankai complemented palindrome和Nankai CDS；(4)來自多種宿主(蝙蝠、果子貍和人)的病毒數(shù)據(jù)，特別是陽性樣本(如來自SARS和武漢肺炎患者)數(shù)據(jù)的積累。經(jīng)典理論認為病毒與宿主共用一套翻譯系統(tǒng)，然而病毒的蛋白質翻譯也有其特殊性，導致本研究從Nankai CDS中預測的17種蛋白質不能提供準確的序列。蛋白質序列的準確度不影響本研究的所有結論和推論，而且為下一步蛋白質、細胞和動物水平的實驗指出方向。

基于可變翻譯的基因分型相同的病毒都在進化樹上的同一個分支(同種宿主)內(nèi)預示了beta冠狀病毒可能通過可變翻譯以適應不同宿主(見圖1B)。但是，由于目前能獲得的病毒數(shù)據(jù)僅僅來自少量宿主(蝙蝠、果子貍和人)，這一推斷的驗證還需要積累來自更多種類宿主的病毒的數(shù)據(jù)。我們推測，病毒通過可變翻譯以適應不同宿主作為一種普遍機制，并不限于冠狀病毒或原核生物；真核生物中某些基因的可變剪接可能也有相似或相近的生物學功能。

BB冠狀病毒通過可變翻譯產(chǎn)生的多種蛋白質的致病性研究是一個更為重要的課題。當前公開的基因組數(shù)據(jù)雖然總量很大，但是，實際有效的陽性樣本(有實驗或臨床記錄證明其致病性)種類不足和缺乏陰性(臨床實驗確認的非SARS病毒)樣本。致病性最明確的2019新型冠狀病毒和SARS冠狀病毒在Nankai CDS上各產(chǎn)生五個小蛋白和一個大蛋白(見圖1C)，并沒有顯示出蛋白質某個屬性(如大小)與致病性的關系。在預測的17個蛋白質中，只有P16(長度為155個氨基酸)和P9(長度為115個氨基酸)是大蛋白質(見表1)；P16(來自SARS冠狀病毒)的致病性已知，而P9(來自蝙蝠的冠狀病毒)未知。因此，來自蝙蝠的其它冠狀病毒產(chǎn)生的蛋白(如P9)未必沒有致病性，或許也具有導致2019肺炎或SARS的能力。本研究預測的17種蛋白質大小適中而且差異明顯，可以很容易通過蛋白質組或外源表達等技術進一步驗證其致病性。

致謝：感謝新型冠狀病毒相關工作人員(省疾控中心等)的基礎性工作，特別復旦大學張永振等公開的2019新型冠狀病毒基因組數(shù)據(jù)，感謝南開大學生命科學學院動物系卜文俊、張濤、黃大衛(wèi)、劉燕強和賀秉軍等各位老師對我們生物信息學研究的長期支持。