云南新現(xiàn)蝙蝠SARS樣冠狀病毒密碼子偏性及其聚類分析

，， ，，，

嚴(yán)重急性呼吸道綜合癥(Severe acute respiratory syndrome coronavirus, SARS)是由SARS冠狀病毒(SARS-CoV)引起的一種嚴(yán)重的呼吸道疾病[1]，是21世紀(jì)第1個(gè)在全球暴發(fā)流行的傳染病。該病以其起病急、發(fā)展迅猛、死亡率高的特點(diǎn)，引起了全世界民眾的恐慌[2]。SARS-CoV初步認(rèn)為起源于蝙蝠SARS樣冠狀病毒(SARSr-CoV)[3],其中從菊頭蝠體內(nèi)分離出的冠狀病毒基因組結(jié)構(gòu)與SARS-CoV相似，且核苷酸一致性在88%～92%之間[4]。雖然到目前為止SARSr-CoV還未出現(xiàn)暴發(fā)流行，但其核苷酸與SARS-CoV的高度一致性，仍應(yīng)引起相關(guān)部門的高度重視。

密碼子(codon)是生物體內(nèi)遺傳信息傳遞不可缺少的物質(zhì)，在編碼氨基酸過程中，多種密碼子可以編碼同一種氨基酸，即同義密碼子[5]。在不同的生物體內(nèi)，甚至同一種生物不同的蛋白質(zhì)基因?qū)γ艽a子的使用頻率不盡相同，具有一定的偏性，即同義密碼子使用的偏性[6]。密碼子偏性在單細(xì)胞物種和多細(xì)胞物種的基因組和基因形成中起著重要作用[7-8]。同時(shí)，病毒密碼子偏性的使用情況也可揭示不同毒株間的進(jìn)化關(guān)系[5]。

2017年12月，中科院武漢病毒研究所石正麗課題組再次報(bào)道了11株云南新現(xiàn)蝙蝠SARSr-CoV，并與SARS-CoV的基因序列具有高度相似性，且其S蛋白上具有SARS-CoV的受體結(jié)合區(qū)(RBD)[3]。本研究以這11株云南新現(xiàn)的蝙蝠SARSr-CoV為主要對(duì)象，分析了它們與SARS-CoV密碼子偏性的異同，并以密碼子偏性為基礎(chǔ)進(jìn)行聚類分析，探索不同時(shí)期和地點(diǎn)發(fā)現(xiàn)的SARSr-CoV與SARS-CoV密碼子偏性之間的關(guān)系。

1 材料與方法

1.1 材 料

1.1.1目的序列來源 本研究涉及的冠狀病毒基因編碼序列(CDS)均來自NCBI數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/)。

1.1.2使用軟件 本研究使用EMBOSS(http://emboss.toulouse.inra.fr/)子程序CUSP計(jì)算密碼子Frequency值[9]，CodonW(https://sourceforge.net/projects/codonw/)計(jì)算密碼子ENC、GC、GC3S、RSCU值[10]，Lasergene子程序EditSeq和MegAlign用于比對(duì)和截取蛋白編碼序列。使用Sigmaplot繪圖、SPSS22.0進(jìn)行聚類分析。

1.2 方 法

1.2.1有效密碼子數(shù) 有效密碼子數(shù)(Effective number of codon, ENC)是由 Wright[11]于1990 年提出的一種描述密碼子使用偏離隨機(jī)選擇程度的方法。它是一個(gè)基因的密碼子使用頻率與同義密碼子平均使用頻率的量化值。同時(shí)該值可單獨(dú)由密碼子使用數(shù)據(jù)計(jì)算得出，與基因的長度及氨基酸組成無關(guān)[12]，因此可以對(duì)基因的密碼子偏性程度提供一個(gè)比較客觀的標(biāo)準(zhǔn)。本研究采用CodonW計(jì)算出特定序列中有效密碼子的ENC值。其取值范圍為20(每個(gè)氨基酸只使用一個(gè)密碼子的極端情況)～61(各個(gè)密碼子均被平均使用)。因此，當(dāng)ENC值越接近20時(shí)其偏性越強(qiáng)，而ENC值越接近61，其偏性就越弱甚至沒有偏性[6]。

1.2.2相對(duì)密碼子使用度 相對(duì)密碼子使用度(Relative Synonymous Codon Usage，RSCU)是指對(duì)于某一特定的密碼子在編碼對(duì)應(yīng)氨基酸的同義密碼子間的相對(duì)頻率，它去除了氨基酸組成對(duì)密碼子使用的影響[13]，當(dāng)密碼子不存在偏好時(shí)，該密碼子的RSCU值等于1。當(dāng)某一密碼子的RSCU值大于1時(shí)，代表該密碼子為使用相對(duì)較多的密碼子，即偏性密碼子，反之亦然[14]。

1.2.3ENC-Plot分析 ENC-Plot分析是用于確定密碼子使用偏愛性的影響因素(尤其是突變偏倚/突變壓力)。GC3S值代表同義密碼子的第3個(gè)密碼子位置中的鳥嘌呤或胞嘧啶的頻率，不包括Met、Trp和終止密碼子。ENC-Plot分析方法是用GC3S為橫坐標(biāo)，ENC為縱坐標(biāo)，來預(yù)測ENC和GC3S之間的功能關(guān)系。即針對(duì)GC3S值繪制ENC值[11]。其中，標(biāo)準(zhǔn)曲線是指在不存在選擇壓力，密碼子的偏性完全取決于突變壓力的情況[15]。故當(dāng)基點(diǎn)位于標(biāo)準(zhǔn)曲線上或散布在標(biāo)準(zhǔn)曲線附近時(shí)，說明密碼子的偏愛性受突變影響較大。反之，說明密碼子使用模式受到自然選擇的影響比較大。

1.2.4中性分析 中性分析(Neutrality Plot)是另外一種確定密碼子使用偏性的影響因素(尤其是自然選擇偏倚/自然選擇壓力)，是以每一個(gè)編碼基因的密碼子的GC1、GC2平均值為縱坐標(biāo)，GC3為橫坐標(biāo)繪制的二維散點(diǎn)圖。如果散點(diǎn)的趨勢線分布于對(duì)角線(斜率=1)，則表明基因僅受突變的影響。反之，斜率為0說明該基因在進(jìn)化過程中受到選擇壓力的影響[14]。

1.2.5奇偶規(guī)則分析 奇偶規(guī)則分析(Parity Rule 2 plot analyses)也被應(yīng)用于密碼子使用偏性的分析。此規(guī)則是一種核苷酸內(nèi)鏈規(guī)則，如果在2個(gè)互補(bǔ)鏈DNA之間不存在任何突變或選擇效應(yīng)上的偏倚，則預(yù)測A=T和G=C[16]。如果由4個(gè)密碼子編碼的氨基酸的第3位密碼子中的奇偶規(guī)則存在顯著偏倚，則在進(jìn)化過程中則以自然選擇壓力為主[17]。因此規(guī)則規(guī)定僅選擇由4個(gè)密碼子編碼的氨基酸。以G3/(G3+C3)|4為橫坐標(biāo)，以A3/(A3+T3)|4為縱坐標(biāo)繪制二維散點(diǎn)圖[15]。圖的中心處遵循PR2原則，A=T且G=C，即橫縱坐標(biāo)都為0.5。從這個(gè)中心出發(fā)的矢量表示了PR2偏倚的程度和方向。

1.2.6基于密碼子偏性的聚類分析 本研究采用CondonW計(jì)算出所需的SARS-CoV及SARSr-CoV基因序列的密碼子使用頻數(shù)，使用SPSS 22.0進(jìn)行聚類分析。

2 結(jié) 果

2.1有效密碼子數(shù)目分析 SARSr-CoV的全基因組由29 727個(gè)核苷酸組成，由11個(gè)開放閱讀框(ORF)。共編碼S、E、M和N共4種結(jié)構(gòu)蛋白和ORF1a、ORF1b、ORF3a、ORF3b、ORF6、ORF7a、ORF7b、ORF8等多種功能未知的蛋白[18]。研究表明，新發(fā)現(xiàn)的11株云南蝙蝠SARSr-CoV各蛋白有效密碼子數(shù)目(ENC)值均接近61(如表1所示)，初步表明密碼子偏性總體較低。并且，云南新現(xiàn)蝙蝠SARSr-CoV與SARS-CoV相同蛋白的ENC值非常接近。

表1 云南新現(xiàn)蝙蝠SARSr-CoV與SARS-CoV的ENC值(均值)Tab.1 ENC value of emerging bat SARSr-CoV in Yunnan province and SARS-CoV (Average)

2.2相對(duì)密碼子使用度 運(yùn)用CodonW軟件分別計(jì)算新現(xiàn)11株云南蝙蝠SARSr-CoV與SARS-CoV(SZ3、GZ02、Tor2、BJ01、PC4-227)12條蛋白各同義密碼子的RSCU平均值(結(jié)果見表2)。結(jié)果顯示，S、M、N、E、ORF1a、ORF1ab、ORF3a、ORF3b、ORF6、ORF7a、ORF7b、ORF8分別具有27、25、30、28、26、27、25、27、25、27、18、30個(gè)偏性密碼子(RSCU>1)，并以A、U結(jié)尾的密碼子為主；其中，ACU為新現(xiàn)云南蝙蝠SARSr-CoV的12種蛋白共有的偏愛密碼子，而GGG、CGG、AGC、UGG、UAG、AUC、AUG為非偏性密碼子。

表2 云南新現(xiàn)蝙蝠SARSr-CoV的同義密碼子平均使用度(RSCU)Tab.2 Relative synonymous codon usage (RSCU) of emerging bat SARSr-CoV in Yunnan province

AAC0.60 1.10 0.79 0.80 0.83 0.75 1.05 2.00 0.00 0.84 1.00 0.61 LysAAA1.22 1.33 1.37 2.00 0.99 1.02 1.28 0.91 0.30 1.33 2.00 1.45 AAG0.78 0.67 0.63 0.00 1.01 0.98 0.72 1.09 1.70 0.67 0.00 0.55 AspGAU1.18 0.38 0.96 2.00 1.23 1.25 0.82 0.00 1.53 0.97 1.00 0.81 GAC0.82 1.62 1.04 0.00 0.77 0.75 0.99 0.00 0.47 1.03 1.00 1.19 GluGAA1.07 1.04 1.01 2.00 1.11 1.10 0.86 0.00 0.80 0.79 1.50 1.84 GAG0.93 0.96 0.99 0.00 0.89 0.90 1.14 2.00 1.20 1.21 0.50 0.16 CysUGU1.18 1.27 0.00 0.67 1.29 1.28 0.61 0.09 0.00 0.68 1.00 1.48 UGC0.82 0.73 0.00 1.33 0.71 0.72 1.39 1.91 0.00 1.32 1.00 0.52 TERUGA0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.55 2.77 0.00 0.00 TrpUGG1.00 1.00 1.00 0.00 1.00 1.00 0.91 0.09 1.00 0.09 1.00 1.00 ArgCGU1.29 1.42 0.97 3.00 2.11 1.85 0.79 0.26 0.00 0.90 0.00 2.04 CGC0.56 1.32 1.55 0.00 0.64 0.76 0.00 0.00 0.00 0.08 0.00 2.18 CGA0.36 0.90 1.18 3.00 0.30 0.21 0.00 2.68 0.00 0.18 0.00 0.14 表2(續(xù))AACondonSMNEORF1aORF1abORF3aORF3bORF6ORF7aORF7bORF8CGG0.02 0.45 0.00 0.00 0.08 0.10 0.00 0.08 0.00 0.06 0.00 0.00 SerAGU0.87 0.47 1.23 0.87 1.25 1.27 1.02 0.51 2.09 0.11 0.00 1.65 AGC0.37 0.94 0.90 0.87 0.50 0.49 0.83 0.99 0.00 0.78 0.00 0.27 ArgAGA2.06 1.04 1.83 0.00 2.09 2.25 4.42 1.62 2.00 4.72 0.00 1.37 AGG1.70 0.87 0.48 0.00 0.78 0.82 0.79 1.36 4.00 0.06 0.00 0.27 GlyGGU1.67 1.31 0.90 2.00 2.30 2.17 1.35 0.00 0.00 1.82 0.00 1.15 GGC1.19 0.90 1.41 0.00 0.80 0.83 1.85 4.00 0.00 0.91 0.00 0.55 GGA1.00 0.80 1.33 1.64 0.75 0.86 0.81 0.00 0.00 1.15 0.00 2.31 GGG0.14 1.00 0.36 0.36 0.16 0.14 0.00 0.00 0.00 0.12 0.00 0.00

終止密碼子及Trp和Met(RSCU值為1)未列入表中；AA表示氨基酸

2.3ENC-Plot關(guān)聯(lián)分析 本研究進(jìn)一步繪制了SARS-CoV與新現(xiàn)的11株云南蝙蝠SARSr-CoV各蛋白的ENC和GC3S的散點(diǎn)圖(圖1)，可以看出兩者幾乎落在了同一區(qū)域，這進(jìn)一步說明兩者的密碼子偏好受到了相似因素的影響。其中大部分點(diǎn)落在了遠(yuǎn)離標(biāo)準(zhǔn)曲線的位置，這說明云南新現(xiàn)蝙蝠SARSr-CoV與SARS-CoV受自身突變和自然選擇的雙重影響，但以自然選擇為主。

2.4中性分析 為進(jìn)一步分析11株云南新現(xiàn)蝙蝠SARSr-CoV與SARS-CoV密碼子偏愛性的影響因素，我們以S蛋白、ORF1a兩蛋白作為結(jié)構(gòu)蛋白與非結(jié)構(gòu)蛋白的代表，分析突變及自然選擇壓力對(duì)兩者密碼子偏性的影響(圖2)。

如圖2所示，云南新現(xiàn)蝙蝠SARSr-CoV與SARS-CoV的S、ORF1a兩蛋白的斜率(b)分別是0.208 4和0.257 9，更接近0，表明兩者主要受到自然選擇的影響。

圖1 云南新現(xiàn)蝙蝠SARSr-CoV與SARS-CoV相關(guān)蛋白編碼基因的ENC-Plot分析Fig.1 ENC-Plot analysis about the protein coding sequences of emerging bat SARSr-CoV in Yunnan province and SARS-CoV

圖2 云南新現(xiàn)蝙蝠SARSr-CoV及SARS-CoV相關(guān)蛋白編碼基因的中性繪圖分析Fig.2 Neutrality Plot analysis of the protein coding gene of emerging bat SARSr-CoV in Yunnan province and SARS-CoV

2.5奇偶規(guī)則分析 為進(jìn)一步明確新現(xiàn)的11株云南蝙蝠SARSr-CoV與SARS-CoV密碼子偏性形成中自然選擇壓力的作用，我們分析了由4個(gè)密碼子編碼的氨基酸的第3位密碼子中的奇偶規(guī)則。如圖3所示，兩者大部分基因落在了PR2圖的左下方，這意味著在密碼子第3位C和T(嘧啶)的使用頻率高于G和A(嘌呤)。即其頻率存在顯著偏倚，提示兩者在進(jìn)化過程中主要受自然選擇壓力為主[17]。

2.6基于密碼子偏性的聚類分析 我們進(jìn)一步從密碼子偏愛性角度分析了新發(fā)現(xiàn)的11株蝙蝠SARSr-CoV毒株與以往發(fā)現(xiàn)的SARSr-CoV毒株以及SARS-CoV毒株之間的進(jìn)化關(guān)系。為了能夠和石正麗課題組基于基因序列構(gòu)建系統(tǒng)進(jìn)化樹的結(jié)果進(jìn)行比較，我們同樣選擇了包膜蛋白S和非結(jié)構(gòu)蛋白基因ORF1a 2個(gè)蛋白的密碼子使用頻率進(jìn)行聚類分析，如圖4。

圖3 云南新現(xiàn)蝙蝠SARSr-CoV與 SARS-CoV相關(guān)蛋白編碼基因的奇偶規(guī)則分析(PR2)Fig.3 Parity Rule Analysis (PR2) of the protein coding gene of emerging bat SARSr-CoV in Yunnan province and SARS-CoV

SARS-CoV為紫色字體；云南新現(xiàn)的蝙蝠SARSr-CoV為紅色加粗字體；其余在云南省發(fā)現(xiàn)的SARSr-CoV為黑色加粗字體。圖4 基于S蛋白、ORF1a蛋白密碼子偏愛性的聚類分析Fig.4 Cluster analysis of S and ORF1a proteins based on codon bias

對(duì)S蛋白的聚類分析發(fā)現(xiàn)，云南新報(bào)道的11株蝙蝠SARSr-CoV分散聚類于先前報(bào)道的SARSr-CoV中。其中，新報(bào)道的Rs7327、Rs9401、Rs4084、Rs4231、Rs4874這5株云南蝙蝠SARSr-CoV與既往報(bào)道的3株云南SARSr-CoV (Rs/YN2013、Rs/WIVI1、Rs/WIV16)一起與SARS-CoV聚為一類。其中， Rs/YN2013、Rs/WIVI1、Rs/WIV16與SARS-CoV進(jìn)化關(guān)系最為密切。新報(bào)道的6株云南蝙蝠SARSr-CoV(Rs4081、Rs4255、As6526、Rs4237、Rs4247、Rf4092)則與香港(HKU3-1～ HKU3-13)、廣西(Rs/GX2013)發(fā)現(xiàn)的SARSr-CoV聚為一類。既往報(bào)道的云南株YNLF_31C、YNLF_34C、RS3367與貴州(Rs/Rs672)、廣西(Rs/Rp3)、湖北(Rf1、Rm1、HuB2013)、河北(HeB2013)、吉林(RfJL2012)發(fā)現(xiàn)的病毒株聚在了一起；此外，山西毒株Rf/SX2013和陜西毒株Rp/SAX2011則單獨(dú)聚為一類。

而非結(jié)構(gòu)蛋白ORF1a的聚類分析發(fā)現(xiàn)，云南蝙蝠SARSr-CoV(包括11株新報(bào)道SARSr-CoV及以往報(bào)道的YNLF_31C和YNLF_34 C)均與SARS-CoV聚在了一起，提示這些毒株與SARS-CoV的密碼子偏愛性相似度較高，該分支同時(shí)包括了來自于貴州(Rs/Rs672)、廣西(Rs/Rp3、Rs/GX2013)、湖北(Rf1、Rm1)的毒株。而香港報(bào)道的SARSr-CoV毒株與SARS-CoV差異較大，與一株來自湖北的毒株(HuB2013)共同聚為一類。山西毒株Rf/SX2013、陜西毒株Rp/SAX2011、河北毒株HeB2013、湖北毒株Rf1、吉林毒株RfJL2012等SARSr-CoVs則共同聚為一類。

值得注意的是，與SARS-CoV密碼子偏愛性最為密切的SARSr-CoV毒株均來自云南。其余省份檢測到的SARSr-CoVs與SARS-CoV的距離較遠(yuǎn)。

3 討 論

密碼子是蛋白質(zhì)編碼基因中的基本單位和進(jìn)化單位，生物信息學(xué)的發(fā)展使得分析與計(jì)算大量基因序列數(shù)據(jù)成為可能。本研究中，我們使用一些常用的生物信息軟件圍繞密碼子使用模式，對(duì)新報(bào)道的11株云南蝙蝠SARSr-CoVs密碼子偏性及其與SARS-CoVs以及以往報(bào)道的SARSr-CoVs的進(jìn)化關(guān)系進(jìn)行分析。

本研究發(fā)現(xiàn)，這11株SARSr-CoVs各蛋白間有效密碼子數(shù)目(ENC)值均接近61，密碼子偏性總體較低。此外云南蝙蝠SARSr-CoV與SARS-CoV進(jìn)化過程中，在偏愛密碼子的選擇上非常接近。通過RSCU分析我們發(fā)現(xiàn)云南蝙蝠SARSr-CoV各蛋白偏愛性密碼子以A、U結(jié)尾的密碼子居多，這與中東呼吸綜合征冠狀病毒(MERS-CoV)偏愛以A、U結(jié)尾的密碼子的結(jié)果相似[19]。ENC-plot、中性繪圖分析、PR2分析均顯示，云南蝙蝠SARSr-CoV與SARS-CoV的密碼子偏性主要受到自然選擇等其他因素的影響，該結(jié)果也和MERS-CoV一致[19]。

通過基于密碼子偏性的聚類分析發(fā)現(xiàn)，無論是S蛋白還是ORF1a蛋白，與SARS-CoV密碼子偏性相似度較高的SARSr-CoV毒株均來自云南。另發(fā)現(xiàn)這11株新報(bào)道的云南蝙蝠SARSr-CoV在非結(jié)構(gòu)蛋白基因ORF1a上均彼此接近，且與SARS-CoV緊密地聚在了一起。這些結(jié)果提示，不論是結(jié)構(gòu)蛋白和非結(jié)構(gòu)蛋白，云南蝙蝠SARSr-CoV的密碼子偏性與SARS-CoV具有較高的相似性，該結(jié)果也從密碼子偏性的角度，佐證了云南蝙蝠SARSr-CoV可能是SARS-CoV進(jìn)化起源的觀點(diǎn)[3]。

此外，不同地區(qū)來源的SARSr-CoV其密碼子偏愛性存在差異，并與地理位置之間可能存在一定的關(guān)聯(lián)。該發(fā)現(xiàn)也與石正麗課題組在ORF1蛋白系統(tǒng)進(jìn)化的結(jié)果一致[3]。這可能與鄰近地區(qū)之間蝙蝠的SARSr-CoV相互傳播過程有關(guān)。此外，中國其他省份發(fā)現(xiàn)的SARSr-CoV幾乎均可在云南找到相似密碼子偏性的毒株，這進(jìn)一步驗(yàn)證了石正麗課題組關(guān)于云南蝙蝠SARSr-CoV可能是SARS-CoV、及其他地區(qū)SARSr-CoV重要的天然基因庫的觀點(diǎn)[3]。