SARS-CoV-2病毒全基因組序列比對及進(jìn)化分析

周 翔，趙仁生，崔藝璇，許詩嘉，宋鵬飛，溫 敏，袁 燕，萬春平

(1.云南民族大學(xué) 民族醫(yī)藥學(xué)院 民族藥資源化學(xué)國家民委-教育部重點(diǎn)實(shí)驗(yàn)室,云南 昆明 650504;2.云南中醫(yī)藥大學(xué) 第一附屬醫(yī)院,云南 昆明 650021)

2019年11月新型冠狀病毒在中國武漢首次被發(fā)現(xiàn)，在隨后短短半年，新型冠狀病毒給人類帶來了前所未有的災(zāi)難，無數(shù)的生命的流逝，難以言計(jì)的經(jīng)濟(jì)損失，時(shí)至今日，疫情的影響還在延續(xù)，人們的生命，國家的安全正在受到威脅.SARS-CoV-2(Severe Acute Respiratory Syndrome Coronavirus 2)出現(xiàn)以來，病毒溯源和中間宿主等問題得到了民眾的廣泛關(guān)注，雖然只有約2%的感染者有野生動物接觸史，但研究人員普遍認(rèn)為SARS-CoV-2是1種來自野生動物(如蝙蝠、穿山甲等)的病毒.病毒的跨種屬傳播十分困難，而且往往會引起新宿主的高病死率，但SARS-CoV-2被發(fā)現(xiàn)以來，表現(xiàn)出對人類機(jī)體已經(jīng)有了很好的適應(yīng)性，具有傳播能力強(qiáng)和病死率較低的特性，這與之前2003年在中國香港發(fā)現(xiàn)的SARS病毒和2013年中東地區(qū)爆發(fā)的MERS-CoV病毒所引起新發(fā)、突發(fā)傳染病的大多數(shù)病毒有所不同[1].國際病毒分類學(xué)委員會(International Committee on Taxonomy of Viruses，即 ICTV)在2018年將冠狀病毒(CoV)科分為病毒亞科和正冠狀病毒亞科.病毒亞科僅包括Alphaletovirus 1個(gè)屬，但正冠狀病毒包含有α屬、β屬、γ屬和δ屬，總共4個(gè)屬[2].目前發(fā)現(xiàn)，冠狀病毒感染脊椎動物，主要以哺乳動物，尤其是人類和鳥類以及蝙蝠等為主.感染人類和畜類的冠狀病毒多屬于正冠狀病毒科中α屬和β屬，感染禽類的多屬正冠狀病毒中的γ屬，在豬和禽類的感染病例中也發(fā)現(xiàn)了部分正冠狀病毒科的δ屬.根據(jù)目前冠狀病毒侵?jǐn)_人類歷史還發(fā)現(xiàn)在感染人群和蝙蝠的冠狀病毒中，既有正冠狀病毒α屬也有β屬，造成此次疫情的病毒源新型冠狀病毒(SARS-CoV-2)正是正冠狀病毒科β屬[3-7].新型冠狀病毒(SARS-CoV-2)為直徑60～140 nm 的正鏈單股RNA病毒[8].其主要會感染機(jī)體上皮細(xì)胞，導(dǎo)致免疫功能低下的人群，尤其是老年人和小孩表現(xiàn)出呼吸道，消化道等疾病的臨床癥狀.有些感染癥狀輕微，有些則致命.新型冠狀病毒(SARS-CoV-2)主要傳播途徑為經(jīng)呼吸道飛沫和接觸傳播，相對封閉條件下，高濃度氣溶膠傳播及糞便和尿液對環(huán)境污染造成氣溶膠或接觸傳播[9].鑒于目前對該病毒進(jìn)行的全序列分析相對較少，本論文對該病毒進(jìn)行了全段位的研究分析，旨在闡述病毒的傳播特點(diǎn)，為病毒的追根溯源提供相關(guān)依據(jù)，以及探究新型冠狀病毒的易變異位點(diǎn)，為盡快研究出相關(guān)疫苗及可行性藥物提供參考.

1 SARS-CoV-2生物信息學(xué)分析

1.1 材料篩選

從美國國立生物技術(shù)信息中心(NCBI)的GenBank核酸序列數(shù)據(jù)庫中篩選出36條來自不同國家和地區(qū)的全基因序列，這些序列分別來自中國、澳大利亞、意大利、波蘭、日本、俄羅斯、美國等不同國家(見表1).這些序列中最長的來自于波蘭的PL_P14、PL_P13和PL_P15，全長 29 903 bp,最短的來自于俄羅斯的SCPM-0-Cdna-07，全長為 29 457 bp.另外，NCBI在參考了各國提交的序列數(shù)據(jù)后，建立了SARS-CoV-2標(biāo)準(zhǔn)全基因參考序列(NC_045512)(RefSeq),長度為 29 903 bp,即總共 29 903 個(gè)堿基，并將其作為標(biāo)準(zhǔn)序列用于與其他全基因組序列比對,從而進(jìn)行SARS-CoV-2的一些特征性分析.

表1 SARS-CoV-2全基因序列及序列來源

1.2 分析方法

鑒于基因組序列的兩端存在非編碼區(qū)，在進(jìn)行36條來自不同國家的全基因組序列比對時(shí)，采用NCBI已確定的SARS-CoV-2標(biāo)準(zhǔn)序列(NC_045512)為標(biāo)準(zhǔn)，運(yùn)用其核實(shí)過的編碼基因區(qū)域，即該序列的266-29 674 位堿基，總共 29 408 個(gè)堿基對，使用序列對比分析軟件snapgene進(jìn)行序列比較.

通過snapgene軟件進(jìn)行36條序列和標(biāo)準(zhǔn)序列的對比分析，找出SARS-CoV-2各基因組序列的變異位點(diǎn)及其在序列上的分布，以及序列上堿基對的缺失、插入和替換.從而分析推測出SARS-CoV-2的易變區(qū)域和保守區(qū)域所對應(yīng)的基因序列和對應(yīng)的編碼蛋白質(zhì).統(tǒng)計(jì)所有的SARS-CoV-2序列上堿基替換數(shù)量和及其變異分布位點(diǎn)，找出其序列位點(diǎn)的變異方式，從而分析出SARS-CoV-2中堿基的易變性.

將36條來自不同國家的新型冠狀病毒(SARS-CoV-2)全基因組序列進(jìn)行序列對比，使用軟件ClustalX1.83,采用NJ(NeighbourJoiningTrees)方法構(gòu)建生物系統(tǒng)進(jìn)化樹，結(jié)合個(gè)地區(qū)分離出的病毒株變異分布，從而分析出SARS-CoV-2地區(qū)進(jìn)化和病毒的流行特點(diǎn).

運(yùn)用軟件Lasergene的程序MegAlign對36個(gè)病毒株進(jìn)行同源性分析，將參考序列NC_045512作為參考對象，對應(yīng)與其進(jìn)行同源性比對分析.

2 分析結(jié)果

2.1 自不同地區(qū)的SARS-CoV-2全基因組序列變異數(shù)量分布

統(tǒng)計(jì)分析來自不同國家地區(qū)的36條SARS-CoV-2全基因組序列中266-29 674 位堿基，總計(jì) 29 408 bp.發(fā)現(xiàn)大部分基因序列都存在變異，且變異的類型主要以堿基替換為主.36條全基因組序列經(jīng)比對分析，共存在128處變異位點(diǎn)，其中包括2個(gè)來自俄羅斯的病毒株在64個(gè)位點(diǎn)上的堿基缺失，33個(gè)病毒株在64個(gè)位點(diǎn)上發(fā)生的堿基替換.統(tǒng)計(jì)發(fā)現(xiàn)，36條全基因組序列，共計(jì)發(fā)生堿基缺失128次，堿基替換次數(shù)達(dá)190次，總計(jì)發(fā)生變異發(fā)生的次數(shù)達(dá)到318次，平均每條序列變異8.8次，變異率在0.3‰(318/29 408/36).36個(gè)病毒株的變異數(shù)量分布見表2和圖1.

表2 SARS-CoV-2全基因組變異數(shù)量分布

圖1 36個(gè)SARS-CoV-2全基因組序列變異數(shù)量比較

其中沒有發(fā)生變異的是來自中國武漢的3個(gè)病毒株，分別為WH-02、WH-04和WH-09；變異數(shù)量較少(≤3)的病毒株有5個(gè)，它們是來自中國武漢WH-01、WH-03、WH-05和來自印度的INMI3、INMI6；變異數(shù)量較多(≥10)的病毒株有3個(gè)，這3個(gè)序列是來自俄羅斯的SCPM-0-Cdna-01、SCPM-0-Cdna-07和來自日本的TKY10256.

從變異數(shù)量的分布來看，來自俄羅斯的2個(gè)病毒株SCPM-0-Cdna-01在29544-29608位點(diǎn)和SCPM-0-Cdna-07在29457～29521位點(diǎn)(相當(dāng)于SARS-CoV-2參考序列的29610～29674位點(diǎn))堿基對缺失較多.由于缺失的問題也導(dǎo)致了這2個(gè)病毒株的序列編碼蛋白和其他序列有所區(qū)別.

2.2 變異位點(diǎn)在SARS-CoV-2全基因組序列上的分布

對36條來自不同國家和地區(qū)的SARS-CoV-2病毒株全基因組序列的全長 29 408 bp 的堿基替換位點(diǎn)數(shù)量和堿基缺失位點(diǎn)數(shù)量以及堿基插入數(shù)量進(jìn)行了統(tǒng)計(jì)比對(見表3)，其分布趨勢(見圖2).結(jié)合了NCBI網(wǎng)站上提供的各編碼蛋白質(zhì)類別及分布位點(diǎn)(見圖3)，對各編碼蛋白質(zhì)基因區(qū)域的變異位點(diǎn)及類型進(jìn)行了統(tǒng)計(jì)(見表4).分析對比SARS-CoV-2各蛋白編碼區(qū)間類別(見表4和圖2)，可以觀察到，在全基因組序列中，堿基替換發(fā)生在3k～4k、14k～15k、23k～24k和28k～29k區(qū)間相對比較多，這些區(qū)間主要是ORF1ab蛋白區(qū)域和S蛋白及N蛋白編碼區(qū).堿基缺失主要發(fā)生在29k～30k區(qū)域，主要是ORF10蛋白區(qū)域.在28k～30k區(qū)域，即ORF10蛋白區(qū)域和N蛋白區(qū)域內(nèi)堿基的缺失和堿基替換的位點(diǎn)數(shù)量達(dá)到最高.

表3 36條SARS-CoV-2全基因組上堿基變異位點(diǎn)和數(shù)量分布

表4 SARS-CoV-2全基因組序列上蛋白編碼區(qū)域及變異位點(diǎn)分布

圖2 36條SARS-CoV-2全基因組序列上變異位點(diǎn)及數(shù)量對比

圖3 各蛋白在全基因組序列中的編碼區(qū)域分布

在SARS-CoV-2全基因序列中，堿基變異多數(shù)集中在序列的后半?yún)^(qū)域，也就是發(fā)生在3′一端的23k后面的區(qū)域，即23k～30k之間.后半段的序列編碼包括S蛋白、M蛋白、E蛋白和N蛋白等重要的結(jié)構(gòu)蛋白.除去ORF10蛋白區(qū)域的堿基缺失，變異的集中區(qū)域在S蛋白區(qū)域(21 563～25 384)和N蛋白區(qū)域(28 273～29 533).在變異的位點(diǎn)和數(shù)量上來看，S蛋白編碼區(qū)域變異位點(diǎn)8個(gè)，總變異數(shù)量達(dá)35個(gè)，N蛋白編碼區(qū)域變異位點(diǎn)3個(gè)，變異數(shù)量達(dá)29個(gè)，N蛋白重復(fù)變異率較高.除此之外，E蛋白和M蛋白的變異位點(diǎn)及變異數(shù)量相對較低，值得注意的是ORF6蛋白區(qū)域和ORF7b蛋白區(qū)域沒有發(fā)生變異.

SARS-CoV-2全基因序列中，變異的類型以堿基替換和堿基缺失為主.堿基的缺失發(fā)生在來自俄羅斯的2個(gè)病毒株(SCPM-0-Cdna-07和SCPM-0-Cdna-01)，其缺失區(qū)域?yàn)镺RF10蛋白編碼區(qū)域，幾乎缺失了整個(gè)編碼區(qū)域.堿基的替換發(fā)生的位置較多，S蛋白區(qū)域的堿基替換以A→G類型為主，N蛋白區(qū)域的堿基替換以G→A為主，位點(diǎn)23k前段區(qū)域的堿基替換以C→T為主.

2.3 SARS-CoV-2生物進(jìn)化分析

使用Cluastal X軟件對36個(gè)SARS-CoV-2病毒株的全基因序列進(jìn)行多序列比對(Do Complete Alignment),再利用其進(jìn)化樹功能繪制NJ進(jìn)化樹圖(Neighbour Joining Trees),繪制出來后再運(yùn)用MEGA7.0進(jìn)行進(jìn)化樹調(diào)整及分析，見圖4.

圖4 基于SARS-CoV-2全基因組的遺傳進(jìn)化分析

通過調(diào)整之后得出的進(jìn)化樹可以將36個(gè)序列分成2群組，第2群組里面又可以分成2小群組.從圖中可以看出，第1群組(Group 1)8個(gè)病毒株大部分都來自SARS-CoV-2病毒首次被發(fā)現(xiàn)的地方，即中國武漢，而且序列提交的時(shí)間都在2—3月.另外還有2條來自澳大利亞的病毒株，與源自中國武漢的病毒株親緣性較高，序列提交的時(shí)間在5—6月，可以認(rèn)為病毒株在此時(shí)間段內(nèi)并未發(fā)生太大的變化.第2群組(Group 2)涵括有28個(gè)病毒株，而這28個(gè)病毒株又可以分為2小組，即Group 2-1和Group 2-2.Group 2-1該群組包括有9個(gè)病毒株分別來自日本、波蘭、俄羅斯和意大利，SARS-CoV-2病毒株的提交時(shí)間都在5月20日～6月10日之間，其中5條來自日本，另一分支上的4條序列來自其他3個(gè)國家，親緣性較高，可以看出病毒株極有可能在這個(gè)時(shí)間段內(nèi)伴隨著國際交通路線四處傳播.另外1個(gè)群組(Group 2-2),包含了剩下的19個(gè)SARS-CoV-2病毒株，其中6個(gè)來自美國的病毒株和2個(gè)來自澳大利亞的病毒株占據(jù)了1條分支，提交時(shí)間在5—6月中旬，此群組親緣性較高.剩下的一條進(jìn)化分支上包含了3個(gè)來自波蘭的病毒株，1個(gè)來自日本的病毒株，2個(gè)來自澳大利亞的病毒株，以及5個(gè)來自意大利的病毒株，提交時(shí)間都在5—6月中旬.綜合進(jìn)化樹可以看出，流行于美國等國家的SARS-CoV-2病毒株與中國武漢傳播的SARS-CoV-2病毒株有著一定的區(qū)別，其親緣性并不是很高.

2.4 SARS-CoV-2同源性分析

通過Lasergene軟件的Megalign小程序?qū)?6株選自不同地區(qū)國家的病毒株進(jìn)行同源性分析，其基因組同源性分析結(jié)果如下表5.

表4 SARS-CoV-2參考病毒株NC_045512與其他病毒株的同源性

以病毒株NC_045512為參考序列，與其他病毒株進(jìn)行比對，得出同源性差異.從表中可以看出參考序列NC_045512與其他病毒株同源性是相當(dāng)高的，基本都在99%以上.除了來自美國的序列MA_MGH_00278與參考序列的同源性只有94.9%，以及來自澳大利亞的序列VIC1057與參考序列同源性為98.6%.可以說，在同源性上，所有病毒株與參考病毒株并沒有多大差異，具有高度同源性.

3 討論與分析

綜合來看，篩選的36條SARS-CoV-2全基因組序列大部分都存在變異的現(xiàn)象，除開兩端非編碼區(qū)，全基因組序列總的變異率在0.300‰(318/29 408/36)，ORF1ab蛋白編碼區(qū)的變異率在0.142‰(109/21 289/36)，S蛋白編碼區(qū)域的變異率在0.254‰(35/3 821/36)，N蛋白編碼區(qū)的變異率在0.639‰(29/1 260/36),M蛋白編碼區(qū)的變異率在0.208‰(5/668/36),E蛋白編碼區(qū)的變異率在0.122‰(1/227/36).S蛋白是糖蛋白的1種，由兩個(gè)重要的結(jié)構(gòu)域構(gòu)成，S1亞基纖突蛋白受體結(jié)合域和S2亞基糖蛋白，具有多個(gè)抗原表位，能刺激機(jī)體從而使機(jī)體產(chǎn)生抗體，該蛋白與機(jī)體宿主細(xì)胞ACE2受體的親和力是SARS病毒的10-20倍，因而SARS-CoV-2病毒能夠更快的與人體細(xì)胞結(jié)合使機(jī)體出現(xiàn)病癥，以及在人群中更快的傳播[10].M蛋白是病毒跨膜蛋白，即糖基化的基質(zhì)蛋白的1種，在病毒的組裝和出芽中起作用[11].N蛋白是結(jié)構(gòu)蛋白的1種，主要負(fù)責(zé)SARS-CoV-2病毒的組裝[12].E蛋白控制SARS-CoV-2病毒的組裝，是其包膜的組成成分，E蛋白如果變異缺失，可以導(dǎo)致SARS-CoV-2病毒株喪失使機(jī)體致病的能力[13].結(jié)合上述各蛋白的變異率來看，E蛋白的變異率最低，N蛋白的變異率最高，S蛋白的變異率其次，可以看出病毒在傳播過程中，為了適應(yīng)環(huán)境從而感染不同地區(qū)和國家的人群，在結(jié)構(gòu)上發(fā)生進(jìn)行了相應(yīng)的變異，以達(dá)到和機(jī)體宿主細(xì)胞ACE2受體的結(jié)合.根據(jù)其致病力一直居高不下，可以看出其E蛋白編碼區(qū)的高度穩(wěn)定保守.SARS-CoV-2病毒基因組在ORF1ab蛋白編碼區(qū)具有高度變異性，雖然基本都是沉默變異.在SARS-CoV-2全基因組中還存在另外一處存在高度變異性的區(qū)域，即ORF8蛋白編碼區(qū)，變異率為0.076‰(1/365/36)，由于ORF8蛋白存在多態(tài)性，其變異會導(dǎo)致2個(gè)變異體ORF8-L和ORF8-S，從而導(dǎo)致其蛋白質(zhì)的結(jié)構(gòu)異常[14].除了上述蛋白編碼區(qū)域外，36條SARS-CoV-2病毒全基因組序列值得注意的其他區(qū)域還有ORF6蛋白編碼區(qū)，高度保守沒有發(fā)生任何變異.其次，在36條SARS-CoV-2全基因組序列中，來自俄羅斯的兩條基因組序列尤其值得注意，兩條全基因組序列全部缺少ORF10蛋白編碼區(qū)，是整個(gè)編碼蛋白區(qū)的缺失，可以猜想，在俄羅斯傳播的SARS-CoV-2病毒株大多不含有ORF10蛋白.從病毒株在缺少ORF10蛋白的情況下，仍能感染人群來看，ORF10蛋白在病毒株的重組過程中是可以舍棄的.當(dāng)然也存在另外一種可能，由于其缺失的主要是3′端的基因片段，那就是與測序準(zhǔn)確率和基因組拼接等有關(guān).

雖然現(xiàn)在還未確定SARS-CoV-2病毒株的來源，但根據(jù)已有的研究數(shù)據(jù)及進(jìn)展來看，該病毒不太可能來源于已知家畜、家禽以及犬貓等寵物冠狀病毒[15].在疫情早些時(shí)候，Zhou等發(fā)現(xiàn)SARS-CoV-2基因組序列與蝙蝠身上攜帶的冠狀病毒(Bat CoV RaTG13)的序列一致性達(dá)到了96.3%[16].同時(shí)，中科院等院校通過分析Global Initiative on Sharing All Influenza Data(GISAID)數(shù)據(jù)庫的93個(gè)SARS-CoV-2基因組數(shù)據(jù)發(fā)現(xiàn)：病毒是外源流入該市場內(nèi),并引起大面積傳播與蔓延[17].結(jié)合36株來源于不同國家地區(qū)的SARS-CoV-2病毒株的全基因組序列分析及通過軟件分析的出的生物進(jìn)化樹，可以確定在美國等地區(qū)流行的SARS-CoV-2病毒與中國地區(qū)流行的SARS-CoV-2病毒的親緣性并不是很高.運(yùn)用Lasergene軟件的Megalign程序?qū)λ行蛄羞M(jìn)行了同源性比較，可以得出的結(jié)論就是流行于美國等地區(qū)與流行于中國地區(qū)的SARS-CoV-2病毒是同一種病毒，都是SARS-CoV-2病毒.綜上所述，從美國國立生物技術(shù)信息中心(NCBI)的GenBank核酸序列數(shù)據(jù)庫中篩選出36條來自不同國家和地區(qū)的全基因組序列都是屬于SARS-CoV-2病毒基因組序列，再根據(jù)進(jìn)化樹上的進(jìn)化分支及節(jié)點(diǎn)來看，流行于美國的SARS-CoV-2病毒與中國流行的SARS-CoV-2病毒來源于同一祖先，但是處在不同的分支，其進(jìn)化歷程方向是不同的，因此，傳聞SARS-CoV-2病毒是“以中國為源頭，向全世界傳播”這種說法是沒有依據(jù)的.