key: cord-0816002-3255edr3
authors: Suzuki, Yoshiyuki
title: Methods for making multiple alignment of genomic sequences for severe acute respiratory syndrome coronavirus 2
date: 2020-08-19
journal: Meta Gene
DOI: 10.1016/j.mgene.2020.100785
sha: 0f43ca8aed22a17efd86a015fadf12c814f64016
doc_id: 816002
cord_uid: 3255edr3

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) emerged in December 2019 and caused a pandemic. To monitor the global transmission pattern of SARS-CoV-2, it is required to constantly update the phylogenetic tree of genomic sequences with 29.9 kb, which may be time consuming. Phylogenetic analysis of SARS-CoV-2 may be accelerated by making a multiple alignment of nucleotide sequences using the CPA (combining pairwise alignments) method, in which a pairwise alignment is made for a reference and each of other sequences, and the pairwise alignments are combined into a multiple alignment. Here it is shown from the analysis of 3729 genomic sequences for SARS-CoV-2 and outgroup strains that the CPA method can produce a multiple alignment with an elevated or a reduced number of variable sites depending on the reference compared to the OMA (ordinary multiple alignment) method, which was considered to be the most reliable. In particular, the topology of the phylogenetic tree constructed from the multiple alignment made using the CPA method adopting the outgroup sequence as the reference was considerably different from that using the OMA method, suggesting that the outgroup sequence may not be suitable as the reference in the CPA method.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) emerged as the etiological agent of pneumonia called coronavirus disease 2019 in Wuhan, China, in December, 2019 (Zhu et al., 2020 . Since then, SARS-CoV-2 caused a pandemic, and as of July 24, 2020, 15,012,731 cases and 619,150 deaths have been confirmed worldwide (World Health Organization, 2020) . SARS-CoV-2 is classified in the genus Betacoronavirus of the family Coronaviridae (Gorbalenya et al., 2020; Lu et al., 2020) . The virion of SARS-CoV-2 is enveloped and 50-200 nm in diameter (Chen et al., 2020) , and is considered to utilize angiotensin-converting enzyme 2 (ACE2) as the cellular receptor (Lan et al., 2020; Shang et al., 2020; Wrapp et al., 2020; Yan et al., 2020) . The genome of SARS-CoV-2 is a non-segmented, linear, and single-stranded RNA of positive polarity with the length of 29.9 kb (Wu et al., 2020; Zhou et al., 2020) .

Phylogenetic analysis of genomic sequences for the strains of SARS-CoV-2 isolated over the world has been conducted to clarify the global transmission pattern in the course of the pandemic (Global Initiative on Sharing All Influenza Data, 2020).

However, hundreds of genomes are daily sequenced, and thus the phylogenetic tree for thousands of genomic sequences is required to be constantly updated to reflect the real time situation of the pandemic. For the construction of phylogenetic tree, it is necessary to make multiple alignment of genomic sequences, which may be time consuming when the number of sequences is large. Therefore, the alignment process has often been modified to accelerate the phylogenetic analysis of SARS-CoV-2; e.g., genomic sequences were partitioned into subsets, which were aligned separately with a reference and combined into a multiple alignment (Hadfield et al., 2018) . J o u r n a l P r e -p r o o f

As the number of genomic sequences increases in the phylogenetic analysis of SARS-CoV-2, it may become critical to devise efficient ways to make multiple alignment. For this purpose, the approach using a reference may be extended such that a pairwise alignment is made for a reference and each of other sequences, and the pairwise alignments obtained are combined into a multiple alignment (Suzuki and Gojobori, 2001) . The aim of the present study was to examine the property of this method in the phylogenetic analysis of SARS-CoV-2.

Genomic sequences for 15,419 and 1,211 strains of SARS-CoV-2 were retrieved from Global Initiative on Sharing All Influenza Data (GISAID) (Shu and McCauley, 2017) (supplementary table S1) and Virus Pathogen Database and Analysis Resource (ViPR) (Pickett et al., 2012) (supplementary table S2 ), respectively, on May 1, 2020. After excluding the sequences containing < 29,500 nucleotide sites and ambiguous nucleotides, 3,728 sequences were retained for the following analysis.

Additionally, the genomic sequence for the bat coronavirus strain isolated from Rhinolophus affinis, RaTG13, which is known to be the most closely related to SARS-CoV-2, was retrieved from International Nucleotide Sequence Database (INSD) (accession number: MN996532) (Zhou et al., 2020) and used as the outgroup sequence for the phylogenetic analysis of SARS-CoV-2. In the phylogenetic analysis of nucleotide sequences, it is customary to use only the sites that are shared by all sequences in multiple alignment, to avoid problems caused by the heterogeneity in evolutionary rates among sites (Nei and Kumar, 2000) . Therefore, multiple alignment without gaps was made for 3,729 nucleotide sequences of SARS-CoV-2 and outgroup strains by each method as follows.

In the OMA (ordinary multiple alignment) method, multiple alignment for J o u r n a l P r e -p r o o f Journal Pre-proof 3,729 nucleotide sequences was made by including all sequences in the input file for the computer program MAFFT (version 7.305b) (Katoh et al., 2002) , in which a multiple alignment was made progressively from closely related sequences to distantly related sequences along with a guide tree. The sites containing gaps were eliminated from the result. Although the OMA method was time consuming, this method was supposed to be the most reliable.

In the CPA (combining pairwise alignments) method, one of 3,729 nucleotide sequences was selected as a reference. Pairwise alignment was made for the reference and each of other sequences with MAFFT (version 7.305b) (Katoh et al., 2002) . The pairwise alignments were investigated to identify the nucleotide positions in the reference that were aligned without gaps to all other sequences. Multiple alignment was generated by aligning the nucleotides at the identified positions in the reference with the nucleotides at the corresponding positions in other sequences. The genomic sequence for the prototype strain of SARS-CoV-2, WH-Human 1 coronavirus (WHCV), which represented the ingroup strain in the phylogenetic analysis of SARS-CoV-2, was selected as the reference (Wu et al., 2020). The genomic sequence for RaTG13, which represented the outgroup strain, was also selected as the reference (Zhou et al., 2020) .

The CPA methods adopting WHCV and RaTG13 as the reference were called the CPA WHCV and CPA RaTG13 methods, respectively, in the present study. Computations for making multiple alignments by the OMA, CPA WHCV , and CPA RaTG13 methods were conducted on Mac OS X (version 10.10.5; 3.1 GHz Intel Core i7; 16 GB 1867 MHz DDR3).

The numbers of conserved and variable sites were counted for the multiple alignments made for 3,729 nucleotide sequences of SARS-CoV-2 and outgroup strains J o u r n a l P r e -p r o o f Journal Pre-proof using the OMA, CPA WHCV , and CPA RaTG13 methods. Since the genetic variation among SARS-CoV-2 strains was effective for resolving the phylogenetic relationships among SARS-CoV-2 strains, the numbers of conserved and variable sites were also counted for 3,728 sequences of SARS-CoV-2 strains.

Each of the multiple alignments made by the OMA, CPA WHCV , and CPA RaTG13 methods was used for constructing the phylogenetic tree by the p distance-based neighbor-joining (NJp) method (Saitou and Nei, 1987) with MEGA (version 7.0.26) , which has been reported to perform better than other methods generally for constructing the phylogenetic tree (Nei and Kumar, 2000; Yoshida and Nei, 2016) . In each phylogenetic tree, the numbers of branches with the length of 0 (0-branches) and > 0 (non-0-branches) were counted separately. In addition, since the topology of the phylogenetic tree is constituted by interior branches (Nei and Kumar, 2000) , the numbers of 0-branches and non-0-branches were divided into those of interior branches (interior 0-branches and non-0-branches, respectively) and exterior branches (exterior 0-branches and non-0-branches, respectively). The total branch length as well as the interior and exterior branch lengths was also computed for each phylogenetic tree. Furthermore, the topologies of the phylogenetic trees constructed from the multiple alignments made using the CPA WHCV and CPA RaTG13 methods were compared with that made using the OMA method by examining the compatibility of partitions supported by interior non-0-branches.

Using the OMA method, multiple alignment for 3,729 nucleotide sequences of SARS-CoV-2 and outgroup strains was made in 28,396 seconds. The multiple alignment contained 28,160 sites, among which 25,059 sites were conserved and 3,101 sites were variable (table 1) Multiple alignments were also made using the CPA WHCV and CPA RaTG13 methods in 1,132 and 1,196 seconds, respectively. Both the numbers of conserved and variable sites were reduced to 25,950 and 2,198, respectively, in the CPA WHCV method (table 1) .

On the other hand, the number of conserved sites was reduced to 25,948, but the number of variable sites was elevated to 2,210 in the CPA RaTG13 method (table 1).

In the phylogenetic tree constructed from the multiple alignment made using the OMA method, the number of 0-branches was 5,088, consisting of 2,498 interior and 2,590 exterior 0-branches, whereas the number of non-0-branches was 2,367, consisting of 1,228 interior and 1,139 exterior non-0-branches (table 2; supplementary figure S1) (Felsenstein, 1986) . The total branch length of the phylogenetic tree was 0.13701, which was divided into 0.03404 and 0.10297 for the interior and exterior branch lengths, respectively. In the CPA WHCV method, the numbers of interior and exterior non-0-branches were decreased to 1,225 and 1,135, respectively (table 2; supplementary   figure S2 ). The interior and exterior branch lengths were also decreased to 0.03402 and 0.10256, respectively. In contrast, in the CPA RaTG13 method, the numbers of interior and exterior non-0-branches were increased to 1,267 and 1,156, respectively (table 2; supplementary figure S3). The interior and exterior branch lengths were also increased to 0.03459 and 0.10370, respectively.

Among 1,225 interior non-0-branches constituting the topology of the phylogenetic tree constructed using the CPA WHCV method, only 23 (2%) were incompatible to the topology of the phylogenetic tree constructed using the OMA method (figure 1). In contrast, 189 (15%) of 1,276 interior non-0-branches constituting the topology of the phylogenetic tree constructed using the CPA RaTG13 J o u r n a l P r e -p r o o f Journal Pre-proof method were incompatible to the topology of the phylogenetic tree constructed using the OMA method (figure 1). Therefore, the incompatibility in the topology of the phylogenetic tree constructed using the CPA RaTG13 method to that using the OMA method was significantly greater than that of the topology of the phylogenetic tree constructed using the CPA WHCV method to that using the OMA method ( = 1.17 × 10 −28 by χ 2 test).

In the present study, the property of the CPA method for making multiple alignment of nucleotide sequences was examined in comparison to the OMA method through the analysis of 3,729 genomic sequences for SARS-CoV-2 and outgroup strains.

Compared to the multiple alignment made using the OMA method, which was supposed to be the most reliable, the multiple alignment made using the CPA WHCV method contained fewer variable sites, resulting in a decrease in the number of interior non-0-branches in the phylogenetic tree. In contrast, the multiple alignment made using the CPA RaTG13 method contained more variable sites, resulting in an increase in the number of interior non-0-branches in the phylogenetic tree. However, many of the interior non-0-branches in the phylogenetic tree constructed using the CPA RaTG13 method were incompatible with those in the phylogenetic trees constructed using the OMA and CPA WHCV methods, which were mostly compatible with each other. selected as the reference, 3,727 pairwise alignments of ingroup sequences and one pairwise alignment of ingroup and outgroup sequences were conducted. On the other hand, in the CPA RaTG13 method, in which an outgroup sequence was selected as the reference, 3,728 pairwise alignments of ingroup and outgroup sequences were conducted. Therefore, the rate of mis-alignment may be higher in the CPA RaTG13 method than in the CPA WHCV method, as observed in the present study, suggesting that the outgroup sequence may not be suitable as the reference in the CPA method. J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f

Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study

The Newick tree format

ViPR: an open bioinformatics database and analysis resource for virology research

Benchmarking tools for the alignment of functional noncoding DNA

The neighbor-joining method: a new method for reconstructing phylogenetic trees

Structural basis of receptor recognition by SARS-CoV-2

GISAID: Global initiative on sharing all influenza datafrom vision to reality

Positively selected amino acid sites in the entire coding region of hepatitis C virus subtype 1b

World Health Organization. 2020. Coronavirus disease

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Efficiencies of the NJp, maximum Likelihood, and Bayesian methods of phylogenetic construction for compositional and noncompositional genes

A pneumonia outbreak associated with a new coronavirus of probable bat origin

A Novel Coronavirus from Patients with Pneumonia in China

Supplementary table S1: Information on the genomic sequences for 15,419 strains of SARS-CoV-2 retrieved from Global Initiative on Sharing All Influenza Data (GISAID

Supplementary table S2: Information on the genomic sequences for 1,211 strains of SARS-CoV-2 retrieved from Virus Pathogen Database and Analysis Resource (ViPR

Supplementary figure S1: Phylogenetic tree for 3,729 genomic sequences of SARS-CoV-2 and outgroup strains constructed from the multiple alignment made using the OMA method. The phylogenetic tree was constructed by the NJp method

Supplementary figure S2: Phylogenetic tree for 3,729 genomic sequences of SARS-CoV-2 and outgroup strains constructed from the multiple alignment made using the CPA WHCV method. The phylogenetic tree was constructed by the NJp method

729 genomic sequences of SARS-CoV-2 and outgroup strains constructed from the multiple alignment made using J o u r n a l P r e -p r o o f the CPA RaTG13 method. The phylogenetic tree was constructed by the NJp method

Conceptualization, Methodology, Software, Investigation, Visualization, Writing

 In CPA method, pairwise alignments are made for reference and other sequences

Genomic sequences of 3,729 SARS-CoV-2 and outgroup strains were analyzed

Tree topology was unreliable using outgroup sequence as reference in CPA method

Outgroup sequence may not be suitable as reference in CPA method

The author thanks two anonymous reviewers for valuable comments. This work was supported by JSPS KAKENHI Grant Number JP19K12221 and AMED Grant Number JP19fk0108033h0003 to Y.S.

The author declares no conflict of interest.J o u r n a l P r e -p r o o f