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

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

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

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

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

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

1 材料與方法

1.1 患者標本收集、RNA提取和實時逆轉錄聚合酶鏈反應（qRT-PCR）定量

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

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

1.2 宏轉錄組文庫的制備和測序

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

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

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

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

為了獲得高質量的SARS-CoV-2病毒測序序列，我們采取了4個嚴格的質量控制步驟，從每個宏轉錄組文庫和雜交捕獲文庫測序結果中獲得高質量的SARS-CoV-2片段序列。首先用Kraken（v0.10.5）軟件［15］將原始測序數(shù)據(jù)同SARS-CoV-2參考基因組（GISAID獲取號：EPI_ISL_402119）［16］數(shù)據(jù)進行對比，篩選得到候選的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ù)的低質量數(shù)據(jù)、接頭污染、PCR重復測序序列和低復雜度序列進行過濾。如果該文庫最終高質量有效測序序列數(shù)小于1 000條，則舍棄其測序數(shù)據(jù)，不用于下游分析。針對每個樣本，本文將兩個不同建庫方法獲得的大規(guī)模并行高通量測序（massively parallel sequencing，MPS）數(shù)據(jù)合并用于變異檢測、組裝和其他下游分析。使用Samtools （v1.9）計算SARS-CoV-2參考基因組在所有文庫測序數(shù)據(jù)中每個核苷酸位點的覆蓋深度。對每個文庫的測序數(shù)據(jù)以200 bp的滑動窗口用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］默認參數(shù)進行組裝得到完整的病毒基因序列。對于其他低深度樣本，使用IDBA_hyrid（v1.1.3，參數(shù)為“-mink 21--maxk 91--step 10”）進行病毒基因組組裝［18］。

為了檢測每個病毒樣品的單核苷酸變異（single nucleotide variations，SNVs）信息，首先使用BWA aln（v0.7.16）將每個文庫的有效測序數(shù)據(jù)與SARSCoV2參考基因組（EPI_ISL_402119）進行比對，然后使用Picard （v2.10.10）去除重復的比對結果。基于去重的比對結果使用FreeBayes （v1.3.1，參數(shù)：-p 1-q 20-m 60-V）［19］分析得到每個測序數(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］的默認參數(shù)進行注釋。如果測序數(shù)據(jù)的參考基因組覆蓋深度大于15 X，則所有SNVs變異檢測調用軟件的最小覆蓋參數(shù)設置為10 X；如果測序覆蓋深度小于15 X，則設置為4 X。為了評估每個樣本的SNVs準確性，使用pysamstats （v1.1.2）軟件（https://github. com/alimanfoo/pysamstats）和Pilon（v1.23）軟件［21］進行變異檢測，并保留至少兩種方法檢測到的SNVs。對于高深度測序（＞40 X）的樣本，對BWA軟件比對結果使用與Ni等［22］相同方法及參數(shù)進行宿主內單核苷酸變異（intradividual single nucleotide variations，iSNVs）檢測。

2 結果與分析

2.1 病毒基因組測序

我們通過qRT-PCR檢測19個SARS-CoV-2樣本的Ct值，結果表明所有樣本的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個新冠病例樣本臨床信息Tab. 1 Detailed clinical information of 19 cases

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

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

2.2 宏轉錄組測序和雜交捕獲測序的SNVs檢測比較

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

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

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

2.3 基因組組裝及SNVs分析

合并兩個不同建庫方式的測序數(shù)據(jù)后，總共有10個樣品的測序深度大于15 X，其余9個樣品測序深度為5～15 X，全部樣品的參考基因組覆蓋率均大于93%（圖1b）。使用SPAdes或IDBA_hyrid進行組裝，共得到10個完整的SARS-CoV-2病毒基因組序列和9個SARS-CoV-2草圖基因組序列，其中9個草圖基因組序列含有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個不同Ct樣本宏轉錄組測序和雜交捕獲測序參考基因組覆蓋深度；（b）19個樣品全部測序數(shù)據(jù)統(tǒng)計參考基因組覆蓋深度。(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個樣本進行SNVs檢測，取至少2種方法檢測到的高置信SNVs為最終鑒定的SNVs （表5）。在19個樣品中總共檢測到了ORF1a:C409S、ORF1a:P1207H、ORF8:L84S和S:D614G 4種突變型，分布在基因ORF1a、ORF8和S上（表4、5）。突變分支S（Clade S，ORF8:L84S）包括4個中國樣本（C01、C02、C07和C03）和唯一的西班牙樣本F01，這些樣品都是在基因ORF8的第251位堿基由T變成C，導致病毒第84位編碼氨基酸由亮氨酸（Leu，L）突變?yōu)榻z氨酸（Ser，S）。突變分支G（Clade G，S:D614G）包括來自同一家意大利餐廳的4個樣本（F07、F08、F09和F10）以及另外5個意大利樣本（F05、F02、F03、F04和F12）。這9個意大利樣本的S蛋白具有相同的突變，都是SRAS-CoV-2病毒第614位氨基酸由天冬氨酸（Asp，D）變成了甘氨酸（Gly，G）。此外，具有相同突變ORF1a:P1207H的樣本包括2個來自中國慶元市的樣本（C05和C06）和2個來自同一家意大利餐廳的樣本（F06和F12），具有突變ORF1a:C409S的只有1個來自中國的樣本C04。變異分析表明，在中國及意大利的新型冠狀病毒肺炎患者檢測到的SARS-CoV-2病毒基因組均存在不一樣的變異類群。同時樣本間檢測具有差異的SNVs，表明了SARS-CoV-2病毒變異的多態(tài)性。

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

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

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

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

2.4 iSNVs變異檢測

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

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

3 討論

本研究從新冠肺炎病例咽拭子和痰液樣本中提取Total RNA，并采用宏轉錄組和雜交捕獲方法進行高通量全基因組測序，最終獲得了19個樣本有效的高通量測序數(shù)據(jù)，組裝得到了10個完整的新冠病毒基因組序列和9個新冠病毒基因組草圖。通過比較宏轉錄組測序和捕獲測序的參考基因組覆蓋度以及變異檢測情況，證實了捕獲測序對于不同病毒載量的樣品均可獲得有效、準確的組裝和變異檢測結果。

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

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