12例海外華人SARS-CoV-2感染者病毒基因組測序及變異分析

曾 永，陳建威，崔恩海，胡國海，王曉光，季巧英，叢 迎，楊 陽，楊偉斌

（1. 青田縣疾控中心，青田 323900；2. 青島華大基因研究院，青島 266555；3. 湖州市中心醫(yī)院，湖州 313000；4. 深圳國家基因庫，深圳 518083；5. 麗水市疾控中心，麗水 323000；6. 麗水市中心醫(yī)院，麗水 323000）

2019年12 月，湖北省武漢市出現(xiàn)嚴(yán)重急性呼吸綜合征冠狀病毒2（severe acute respiratory syndrome coronavirus 2，SARS-CoV-2）引起的新型冠狀病毒肺炎（corona virus disease 2019，COVID-19）［1-2］。病毒感染隨后蔓延到全國各地，導(dǎo)致嚴(yán)重的呼吸道疾病，嚴(yán)重危害人們的生命健康［3］。除中國外，COVID-19疫情在全球多個(gè)國家也出現(xiàn)了，成為國際關(guān)注的突發(fā)公共衛(wèi)生事件［4-5］。同時(shí)，中國開始出現(xiàn)海外輸入性病例，防控病毒的戰(zhàn)斗還在繼續(xù)。

SARS-CoV-2在人群中的持續(xù)傳播可能會(huì)產(chǎn)生適應(yīng)性突變，因此極有可能會(huì)累計(jì)基因突變，進(jìn)而抵抗人類的免疫系統(tǒng)［6］。已有研究表明，基因突變可能會(huì)導(dǎo)致新型冠狀病毒的毒力、傳染性和傳播性的改變［7-8］。因此，對不同來源的新冠肺炎病例進(jìn)行基因組分析有助于了解SARS-CoV-2的變異與感染病程的相關(guān)性，從而為疾病預(yù)防和控制提供建議。

鑒于患者咽拭子和痰液中RNA含量的差異，我們采用最新研發(fā)的捕獲測序［9-10］和宏轉(zhuǎn)錄組測序方法［10］，對意大利歸國和國內(nèi)的新冠肺炎確診病例進(jìn)行基因組突變分析。與宏轉(zhuǎn)錄組測序相比，捕獲測序所需數(shù)據(jù)量僅為宏轉(zhuǎn)錄組測序的1/10或更少，且測序結(jié)果準(zhǔn)確。對病毒基因組的突變分析為深入研究意大利歸國患者SARS-CoV-2的致病力、傳染性以及國內(nèi)毒株的差異打下了基礎(chǔ)。

1 材料與方法

1.1 患者標(biāo)本收集、RNA提取和實(shí)時(shí)逆轉(zhuǎn)錄聚合酶鏈反應(yīng)（qRT-PCR）定量

本項(xiàng)目涉及的所有臨床標(biāo)本（包括咽拭子和痰液）均來自浙江省新冠肺炎確診的19例患者，其中7例為國內(nèi)病例，12例為歐洲輸入性病例（11例為意大利輸入性病例，1例為西班牙輸入性病例）。在11名意大利患者中有7名是意大利貝加莫同一家餐廳的工作人員（詳細(xì)臨床資料見表1）。該項(xiàng)倫理審查已獲麗水市中心醫(yī)院和深圳華大生命科學(xué)研究院倫理委員會(huì)（institutional review board，IRB）批準(zhǔn)。所有參與者在納入這項(xiàng)研究之前都簽署了書面知情同意書。

通過Maxwell?16 Total RNA純化試劑盒進(jìn)行SARS-CoV-2 RNA提取。隨后針對SARS-CoV-2的RdRp基因和N基因進(jìn)行實(shí)時(shí)逆轉(zhuǎn)錄聚合酶鏈反應(yīng)（quantitative real-time PCR，qRT-PCR）。采用HCoV-19核酸檢測試劑盒檢測和定量檢測臨床標(biāo)本中的病毒RNA的載量（捷諾生物，上海）。

1.2 宏轉(zhuǎn)錄組文庫的制備和測序

宏轉(zhuǎn)錄組測序技術(shù)可及時(shí)快速地對傳染病進(jìn)行檢測、防控和變異分析［11-12］。首先，我們對每個(gè)RNA樣本用DNase I去除DNA，并對純化的產(chǎn)物進(jìn)行定量（Qubit RNA HS Assay Kit, Thermo Fisher Science, Waltham, MA, USA）；然后將純化的RNA在87℃緩沖液中孵育6 min進(jìn)行打斷，隨后將RNA片段利用隨機(jī)引物合成雙鏈（ds）cDNA并進(jìn)行末端修復(fù)、接頭連接；最后對（ds）cDNA進(jìn)行18個(gè)循環(huán)的PCR擴(kuò)增并加上雙端測序標(biāo)簽（unique dual-indexed adapters，UDI），使用MGIEasy RNA（華大智造，深圳）的試劑盒經(jīng)過環(huán)化、滾環(huán)復(fù)制后構(gòu)建DNA納米球（DNA nanoball，DNB）文庫［10］。使用Qubit 3.0 熒光定量儀 （Invitrogen, Carlsbad， CA, USA）檢測文庫濃度，并使用Agilent 2100 生物分析儀檢測文庫片段大小，對合格文庫采用DNBSEQ T1測序平臺(tái)（華大智造，深圳）或MGISEQ-2000測序平臺(tái)（華大智造，深圳）進(jìn)行雙端100 bp 讀長模式的高通量測序［4］。

1.3 病毒雜交捕獲富集和測序

病毒全基因組雜交捕獲技術(shù)是病毒分離鑒定的常用方法之一，可提升檢測靈敏度和病毒基因組覆蓋度［13-14］。利用MGIEasy外顯子捕獲輔助試劑盒（華大智造，深圳）替代了SARS-CoV-2 DNA/RNA捕獲試劑盒（伯科生物，北京）中的阻滯劑寡核苷酸（blocker oligos）和PCR引物寡核苷酸（PCR primer oligos）后，采用基于混合捕獲的富集方法從宏轉(zhuǎn)錄組逆轉(zhuǎn)錄雙鏈cDNA片段中富集SARS-CoV-2特異性片段，并根據(jù)試劑盒說明對SARS-CoV-2病毒全基因組進(jìn)行擴(kuò)增［10］。給富集擴(kuò)增得到的病毒cDNA片段加上雙端測序標(biāo)簽，經(jīng)過環(huán)化、滾環(huán)復(fù)制后構(gòu)建DNB文庫，對每個(gè)庫檢合格的文庫使用MGISEQ-2000測序平臺(tái)（華大智造，深圳）進(jìn)行雙端100 bp 讀長模式的高通量測序［4］。

1.4 SARS-CoV-2病毒序列鑒定

為了獲得高質(zhì)量的SARS-CoV-2病毒測序序列，我們采取了4個(gè)嚴(yán)格的質(zhì)量控制步驟，從每個(gè)宏轉(zhuǎn)錄組文庫和雜交捕獲文庫測序結(jié)果中獲得高質(zhì)量的SARS-CoV-2片段序列。首先用Kraken（v0.10.5）軟件［15］將原始測序數(shù)據(jù)同SARS-CoV-2參考基因組（GISAID獲取號(hào)：EPI_ISL_402119）［16］數(shù)據(jù)進(jìn)行對比，篩選得到候選的SARS-CoV-2原始測序數(shù)據(jù)；然后依次用fastp（v0.19.5，參數(shù)：-q20-u20-n1-150）、SOAPnuke（v1.5.6，參數(shù)：-l20-q0.2-E50-d）和PRINSEQ（v0.20.4，參數(shù)：LC_Method Dust-LC_Threshold 7）對測序數(shù)據(jù)的低質(zhì)量數(shù)據(jù)、接頭污染、PCR重復(fù)測序序列和低復(fù)雜度序列進(jìn)行過濾。如果該文庫最終高質(zhì)量有效測序序列數(shù)小于1 000條，則舍棄其測序數(shù)據(jù)，不用于下游分析。針對每個(gè)樣本，本文將兩個(gè)不同建庫方法獲得的大規(guī)模并行高通量測序（massively parallel sequencing，MPS）數(shù)據(jù)合并用于變異檢測、組裝和其他下游分析。使用Samtools （v1.9）計(jì)算SARS-CoV-2參考基因組在所有文庫測序數(shù)據(jù)中每個(gè)核苷酸位點(diǎn)的覆蓋深度。對每個(gè)文庫的測序數(shù)據(jù)以200 bp的滑動(dòng)窗口用R（v3.6.1）軟件中的ggplot2包繪制SARS-CoV-2參考基因組的測序覆蓋率和測序深度。

1.5 病毒基因組組裝及單核苷酸變異（SNVs）檢測

對于參考覆蓋深度大于15 X（X表示全基因組的測序覆蓋深度的層數(shù)）的樣本，截取最多100 X的SARS-CoV-2測序數(shù)據(jù)，使用SPAdes（v3.14.0）軟件［17］默認(rèn)參數(shù)進(jìn)行組裝得到完整的病毒基因序列。對于其他低深度樣本，使用IDBA_hyrid（v1.1.3，參數(shù)為“-mink 21--maxk 91--step 10”）進(jìn)行病毒基因組組裝［18］。

為了檢測每個(gè)病毒樣品的單核苷酸變異（single nucleotide variations，SNVs）信息，首先使用BWA aln（v0.7.16）將每個(gè)文庫的有效測序數(shù)據(jù)與SARSCoV2參考基因組（EPI_ISL_402119）進(jìn)行比對，然后使用Picard （v2.10.10）去除重復(fù)的比對結(jié)果。基于去重的比對結(jié)果使用FreeBayes （v1.3.1，參數(shù)：-p 1-q 20-m 60-V）［19］分析得到每個(gè)測序數(shù)據(jù)的SNVs，最后使用snippy-vcf_filter（參數(shù)：--minqual 100--minfrac 0.8；https://github.com/tseemann/snippy）過濾低置信度SNVs。剩余的高置信度SNVs在SARS-CoV-2參考基因組中用SNVeff（v4.3）［20］的默認(rèn)參數(shù)進(jìn)行注釋。如果測序數(shù)據(jù)的參考基因組覆蓋深度大于15 X，則所有SNVs變異檢測調(diào)用軟件的最小覆蓋參數(shù)設(shè)置為10 X；如果測序覆蓋深度小于15 X，則設(shè)置為4 X。為了評(píng)估每個(gè)樣本的SNVs準(zhǔn)確性，使用pysamstats （v1.1.2）軟件（https://github. com/alimanfoo/pysamstats）和Pilon（v1.23）軟件［21］進(jìn)行變異檢測，并保留至少兩種方法檢測到的SNVs。對于高深度測序（＞40 X）的樣本，對BWA軟件比對結(jié)果使用與Ni等［22］相同方法及參數(shù)進(jìn)行宿主內(nèi)單核苷酸變異（intradividual single nucleotide variations，iSNVs）檢測。

2 結(jié)果與分析

2.1 病毒基因組測序

我們通過qRT-PCR檢測19個(gè)SARS-CoV-2樣本的Ct值，結(jié)果表明所有樣本的Ct值均在22.43～38.84之間（表1）。其中僅3例樣本（C06、F03和F12）檢測到高病毒拷貝（Ct＜24.5），6例樣本（C05、F05、F04、C01、C07和C02）檢測到低病毒拷貝（24.5＜Ct＜28.7），其余10例樣本檢測到極低病毒拷貝（Ct＞28.7）。

表1 19個(gè)新冠病例樣本臨床信息Tab. 1 Detailed clinical information of 19 cases

對所有SARS-CoV-2 RNA樣本構(gòu)建逆轉(zhuǎn)錄文庫，并對其全基因組進(jìn)行了測序。這些樣品獲得的原始測序數(shù)據(jù)雙端序列條數(shù)在8 995萬條到140 824萬條之間，平均約為64 968萬條（130 Gb）。用Kraken和BWA對SARS-CoV-2病毒測序數(shù)據(jù)進(jìn)行分類，并計(jì)算每個(gè)文庫數(shù)據(jù)的覆蓋深度。雖然每個(gè)文庫均產(chǎn)出了大量的測序數(shù)據(jù)，但由于病毒拷貝數(shù)較低，在Ct＜28.7的6個(gè)樣本中，僅有0.004 5%～0.055 8%的測序數(shù)據(jù)被歸類為SARSCoV-2病毒序列（比例最低的樣本C01和最高的樣本C05除外），而在Ct＞28.7的11個(gè)樣本中，僅有0.001 4%～0.006 2%的測序數(shù)據(jù)被歸類為SARSCoV-2。在19個(gè)測序結(jié)果中，只有10個(gè)樣本的有效參考基因組覆蓋率超過95%，且測序覆蓋深度大于5 X（表2）。

然后，我們構(gòu)建了14個(gè)基于雜交捕獲的文庫并進(jìn)行測序以獲得病毒的基因組數(shù)據(jù)，其中包括9個(gè)宏轉(zhuǎn)錄測序數(shù)據(jù)基因組覆蓋度不足的樣本（F05、F10、F01、C04、F02、F06、C03、F09和F11）和5個(gè)宏轉(zhuǎn)錄測序數(shù)據(jù)基因組高覆蓋度的樣本（C06、F12、C05、F04和F07，包含不同的Ct值）。雜交捕獲測序的每個(gè)樣本平均得到了約7 000萬條病毒測序序列條數(shù)，分類鑒定得到SARS-CoV-2病毒測序數(shù)據(jù)的平均比率為8.4%（表2），其比例顯著高于宏轉(zhuǎn)錄組測序所得的結(jié)果，表明捕獲測序法更適合于高Ct（低病毒載量）樣本的測序。該方法測序的所有樣品覆蓋深度均大于5 X，基因組覆蓋率為95%（表2），可滿足用于突變和組裝分析。

2.2 宏轉(zhuǎn)錄組測序和雜交捕獲測序的SNVs檢測比較

為了比較兩種不同建庫方法大規(guī)模并行MPS結(jié)果的準(zhǔn)確性，我們使用宏轉(zhuǎn)錄組和雜交捕獲方法分別對5個(gè)樣本（C06、F12、C05、F04和F07）進(jìn)行了文庫構(gòu)建及測序。分析表明，每個(gè)文庫測序數(shù)據(jù)的參考基因組深度都超過了5 X，且參考基因組覆蓋率均大于96%（圖1a、表2）。使用Pilon和FreeBayes方法對每個(gè)測序數(shù)據(jù)進(jìn)行SNVs檢測，取兩個(gè)方法都認(rèn)為是可靠的SNVs作為該測序數(shù)據(jù)檢測到的SNVs。結(jié)果表明，不論在高病毒拷貝數(shù)樣本（C06、F12）中，還是在低病毒拷貝數(shù)樣本（C05、F04）和極低病毒拷貝數(shù)樣本（F07）中，兩個(gè)不同建庫方法測序數(shù)據(jù)檢測到的SNVs完全一樣，具有很高的一致性（表3）。上述結(jié)果顯示，相對直接進(jìn)行宏轉(zhuǎn)錄組測序，雜交捕獲測序?qū)Σ煌珻t值的樣本均具有準(zhǔn)確的突變檢測結(jié)果，且所需樣本量僅為宏轉(zhuǎn)錄組測序的1/10或更少。高Ct（低病毒載量）的樣品即使測序超過100 Gb也不能獲得有效的SARS-CoV-2病毒測序信息，但捕獲測序可獲得有效的信息，表明雜交捕獲測序是一種更適合于低病毒載量樣本的檢測方法。

表2 19個(gè)新冠病例樣本宏轉(zhuǎn)錄組測序數(shù)據(jù)和雜交捕獲測序數(shù)據(jù)與參考基因組比對統(tǒng)計(jì)Tab. 2 Sequencing data and SARS-CoV-2 reference genome mapping statistics

表3 宏轉(zhuǎn)錄組測序和雜交捕獲測序檢測的SNVs變異信息Tab. 3 SNVs detection by metatranscriptomic and hybrid capture

2.3 基因組組裝及SNVs分析

合并兩個(gè)不同建庫方式的測序數(shù)據(jù)后，總共有10個(gè)樣品的測序深度大于15 X，其余9個(gè)樣品測序深度為5～15 X，全部樣品的參考基因組覆蓋率均大于93%（圖1b）。使用SPAdes或IDBA_hyrid進(jìn)行組裝，共得到10個(gè)完整的SARS-CoV-2病毒基因組序列和9個(gè)SARS-CoV-2草圖基因組序列，其中9個(gè)草圖基因組序列含有12～2 174 bp的gaps（表4）。除F09的覆蓋率為93%外，所有組裝的基因組相對于參考基因組的覆蓋率均在95%以上。

圖1 SARS-CoV-2參考基因組（EPI_ISL_402119）覆蓋深度情況Fig. 1 The sequencing coverage depth of SARS-CoV-2 reference genome (EPI_ISL_402119)（a）5個(gè)不同Ct樣本宏轉(zhuǎn)錄組測序和雜交捕獲測序參考基因組覆蓋深度；（b）19個(gè)樣品全部測序數(shù)據(jù)統(tǒng)計(jì)參考基因組覆蓋深度。(a) The sequencing depth of 5 different Ct samples by metatranscriptomic sequencing and hybrid capture sequencing; (b) Statistical sequencing depth for all sequencing data of 19 samples.

使用FreeBayes、Pysamstats和Pilon 3種方法對19個(gè)樣本進(jìn)行SNVs檢測，取至少2種方法檢測到的高置信SNVs為最終鑒定的SNVs （表5）。在19個(gè)樣品中總共檢測到了ORF1a:C409S、ORF1a:P1207H、ORF8:L84S和S:D614G 4種突變型，分布在基因ORF1a、ORF8和S上（表4、5）。突變分支S（Clade S，ORF8:L84S）包括4個(gè)中國樣本（C01、C02、C07和C03）和唯一的西班牙樣本F01，這些樣品都是在基因ORF8的第251位堿基由T變成C，導(dǎo)致病毒第84位編碼氨基酸由亮氨酸（Leu，L）突變?yōu)榻z氨酸（Ser，S）。突變分支G（Clade G，S:D614G）包括來自同一家意大利餐廳的4個(gè)樣本（F07、F08、F09和F10）以及另外5個(gè)意大利樣本（F05、F02、F03、F04和F12）。這9個(gè)意大利樣本的S蛋白具有相同的突變，都是SRAS-CoV-2病毒第614位氨基酸由天冬氨酸（Asp，D）變成了甘氨酸（Gly，G）。此外，具有相同突變ORF1a:P1207H的樣本包括2個(gè)來自中國慶元市的樣本（C05和C06）和2個(gè)來自同一家意大利餐廳的樣本（F06和F12），具有突變ORF1a:C409S的只有1個(gè)來自中國的樣本C04。變異分析表明，在中國及意大利的新型冠狀病毒肺炎患者檢測到的SARS-CoV-2病毒基因組均存在不一樣的變異類群。同時(shí)樣本間檢測具有差異的SNVs，表明了SARS-CoV-2病毒變異的多態(tài)性。

表4 病毒基因組組裝和突變類型統(tǒng)計(jì)Tab. 4 Genome assembly and mutation clade statistics

表5 （續(xù)）Tab. 5 (continued)

表5 19個(gè)新冠病例樣本SNVs變異檢測結(jié)果Tab. 5 SNVs detection of 19 sequencing SARS-CoV-2 samples

表5 （續(xù)）Tab. 5 (continued)

2.4 iSNVs變異檢測

我們對高測序深度（＞40 X）的9個(gè)樣本進(jìn)行了iSNVs分析。結(jié)果表明：在F03、F08和C06 3個(gè)樣本中沒有檢測到任何iSNVs；在C02、C07、C05和F04這4個(gè)樣本中至少檢測到一個(gè)iSNVs，但沒有iSNVs位于基因區(qū)或非翻譯區(qū)（untranslated region，UTR）。此外，在F06和F11這2個(gè)樣本中均檢測到了4個(gè)iSNVs，都具有相同的突變c.1841A＞G（p.Asp614Gly，即S:D614G）（表6）。F06和F11樣本基因區(qū)的iSNVs與來自同一家意大利餐廳的其他5個(gè)病例樣本的SNVs突變結(jié)果意外一致，表明在同一患者中可能存在重復(fù)感染。

表6 高深度測序樣品iSNVs變異信息統(tǒng)計(jì)Tab. 6 The statistics of iSNVs for high depth sequencing samples

3 討論

本研究從新冠肺炎病例咽拭子和痰液樣本中提取Total RNA，并采用宏轉(zhuǎn)錄組和雜交捕獲方法進(jìn)行高通量全基因組測序，最終獲得了19個(gè)樣本有效的高通量測序數(shù)據(jù)，組裝得到了10個(gè)完整的新冠病毒基因組序列和9個(gè)新冠病毒基因組草圖。通過比較宏轉(zhuǎn)錄組測序和捕獲測序的參考基因組覆蓋度以及變異檢測情況，證實(shí)了捕獲測序?qū)τ诓煌《据d量的樣品均可獲得有效、準(zhǔn)確的組裝和變異檢測結(jié)果。

此外，通過變異結(jié)果分析，11個(gè)意大利歸國華人中9個(gè)病例均檢測到了S蛋白D614G突變，且具有D614G突變的樣本大部分都伴有5'UTR中的C到T突變（相對于參考基因組EPI_ISL_402119的241位）、3 037位的C到T突變、14 408位的C到T突變。另外2個(gè)病例雖然直接檢測到的變異型不屬于D614G，但在iSNVs檢測中也均發(fā)現(xiàn)了D614G突變。具有D614G突變的SARS-CoV-2病毒株2020年2月初開始在歐洲傳播［23］。我們3月初采集的來自意大利不同地區(qū)的F02、F03、F04、F07、F08、F10和F12 7個(gè)樣品都能檢測到同時(shí)含有這4種突變的相同突變型（表5），表明在3月初之前具有這4個(gè)遺傳連鎖突變的單倍型已在意大利不同地區(qū)開始傳播流行。根據(jù)目前研究，D614G突變毒株已成為世界范圍內(nèi)的主流毒株，在歐洲和北美的測序樣本中占近70%，且其具有更強(qiáng)的感染性和傳播性［24］，這可能是意大利疫情暴發(fā)的一個(gè)重要因素。

我們對19個(gè)樣本進(jìn)行SNVs檢測，結(jié)果共發(fā)現(xiàn)4種突變，其中意大利患者中分離的病毒株是具有高傳染性的S:D614G突變病毒株，而ORF1a:C409S、ORF1a:P1207H和ORF8:L84S這3個(gè)突變主要在中國人群中被檢測到，可見SARS-CoV-2病毒突變在意大利和中國并不相同。鑒于SARS-CoV-2基因組會(huì)不斷的發(fā)生突變，那么對SARS-CoV-2基因組的測序就顯得極為重要，而我們新研發(fā)的捕獲測序的方法與傳統(tǒng)的宏轉(zhuǎn)錄組測序相比所需的樣本量更少，是一種更適合于低病毒載量樣本的檢測方法。通過捕獲測序獲得SARS-CoV-2病毒的突變信息并研究突變的功能，對于后續(xù)抗體和疫苗設(shè)計(jì)具有重要的指導(dǎo)意義。