Genomic Epidemiology of Imported Cases of COVID-19 in Guangdong Province, China, October 2020 – May 2021

Objective The pandemic of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been engendering enormous hazards to the world. We obtained the complete genome sequences of SARSCoV-2 from imported cases admitted to the Guangzhou Eighth People's Hospital, which was appointed by the Guangdong provincial government to treat coronavirus disease 2019 (COVID-19). The SARS-CoV-2 diversity was analyzed, and the mutation characteristics, time, and regional trend of variant emergence were evaluated. Methods In total, 177 throat swab samples were obtained from COVID-19 patients (from October 2020 to May 2021). High-throughput sequencing technology was used to detect the viral sequences of patients infected with SARS-CoV-2. Phylogenetic and molecular evolutionary analyses were used to evaluate the mutation characteristics and the time and regional trends of variants. Results We observed that the imported cases mainly occurred after January 2021, peaking in May 2021, with the highest proportion observed from cases originating from the United States. The main lineages were found in Europe, Africa, and North America, and B.1.1.7 and B.1.351 were the two major sublineages. Sublineage B.1.618 was the Asian lineage (Indian) found in this study, and B.1.1.228 was not included in the lineage list of the Pangolin web. A reasonably high homology was observed among all samples. The total frequency of mutations showed that the open reading frame 1a (ORF1a) protein had the highest mutation density at the nucleotide level, and the D614G mutation in the spike protein was the commonest at the amino acid level. Most importantly, we identified some amino acid mutations in positions S, ORF7b, and ORF9b, and they have neither been reported on the Global Initiative of Sharing All Influenza Data nor published in PubMed among all missense mutations. Conclusion These results suggested the diversity of lineages and sublineages and the high homology at the amino acid level among imported cases infected with SARS-CoV-2 in Guangdong Province, China.


INTRODUCTION
S evere acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the pathogenic agent of coronavirus disease 2019 (COVID- 19), has caused a global pandemic and conferred a profound effect on all aspects of society. As of June 15, 2021, 175 million infections have been confirmed, and 3.8 million deaths have been reported all over the world [1] . The Chinese government managed to contain the virus through a series of draconian control measures. However, the number of newly confirmed cases outside of China has been rising rapidly. Thus, China prioritized protecting itself from imported cases from other countries or territories.
Guangdong Province is known as the "South Gate of China," with a large number of foreign passengers; hence, the most number of imported cases was detected here compared with other regions of Mainland China [2] . The World Health Organization announced several variants of concern (VOCs) and variants of interest (VOIs) of SARS-CoV-2 on May 31, 2021 [3] . With the increase in the number of genome sequences available in the public database, analysis using more sequences by sampling from a wider area over a longer period was imperative. In the present study, we employed second-generation sequencing to obtain virus sequence information and performed advanced comprehensive methods to analyze the genetic characteristics of SARS-CoV-2 in the imported cases based on the demographic and baseline characteristics of these patients. Our studies enriched the knowledge of the existing epidemic strains of SARS-CoV-2 in the world and provided valuable clues for defense against this infectious disease imported to China.

Ethics Approval Statement
The protocol of the study was reviewed and approved by Guangdong Provincial Center for Disease Control and Prevention (approval number: W96-027E-202121).

Sample Sources
This study included human throat swab samples that were collected and processed in the Guangdong Provincial Center for Disease Control and Prevention, China, and throat swab samples, which were deidentified prior to analysis, collected from the Guangzhou Eighth People's Hospital from October 2020 to May 2021. All the samples were inactivated at 56 °C for 30 min. Demographic characteristics and clinical information were collected from electronic medical records.

Nucleic acid Test (NAT) and Next-generation Sequencing (NGS)
All the throat swab samples were screened by quantitative real-time polymerase chain reaction method using the SARS-CoV-2 NAT kits (BioGerm, Shanghai, China). The positive specimens were processed for NGS using the MGISEQ-2000 platform of the Beijing Genomics Institute (BGI, Shenzhen, China). DNA fragment assembly was performed using the iGenome® VirusDetector system (BGI, Shenzhen, China), with the Wuhan-hu-01 (GenBank access no. MN908947.3) genome as the reference sequence. All ambiguous nucleotides in the alignment analysis were set to N. The ratio of non-N greater than 85% indicated a good fragment assembly in accordance with the manufacturer's instruction.

Phylogenetic Analysis
The complete reference genomes were retrieved from the National Center for Biotechnology Information database (https://www.ncbi.nlm.nih. gov/). Pangolin web tool (https://cov-lineages.org/) was used for the origin of clades. Mattf v7.037 (Kazutaka Katoh, Japan, 2002) and BioEdit v7.0.9.0 (Borland, USA, 2009) were applied to sequence alignment and trimming, respectively. The maximum-likelihood (ML) phylogenetic tree was constructed using IQ-Tree v1. 6 USA, 1987). In addition, the Nextclade web tool (https://clades.nextstrain.org) was used for putting up information about nucleotide mutations and amino acid mutations. R v3.6.1 (Ross Ihaka, Robert Gentleman, New Zealand, 1995) was used to calculate the genetic distance of sequences, obtain the homology, and cleanse the mutation data at the nucleotide and amino acid levels.
Age variables were expressed using "median {Interquartile range [(IQR): 25%-75%]}," and categorical variables including age group, gender, and NAT results were presented as percentages. A two-tailed P -value of < 0.05 was considered statistically significant throughout the study.

Demographic Information of the Patients
A total of 177 patients were included in this study, and their demographic characteristics are shown in Table 1. More than half of the respondents were male (136, 76.84%). A total of 78 individuals (44.07%) were aged 31-40 years, and 6 people were over 61. The minimum age was 16, and the maximum was 85 (median age, 34 years). A total of 123 out of 177 (69.49%) individuals showed positive results for SARS-CoV-2 RNA by the NATs, and among them, 117 cases (95.12%) were sequenced successfully. Exactly 43 patients in our studies (24.29%) were asymptomatic cases (Table 1).

Phylogenetic Tree and Variation Analyses of SARS-CoV-2 at the Nucleotide and Amino Acids Levels
For this study, a GTR + F + I + G4 nucleotide substitution model, the best-fit model based on the Akaike Information Criterion according to the Bayesian information criterion, was used for the phylogenetic analyses. The virus species in this study clustered in two large groups. Most viral sequences from patients with SARS-CoV-2 infection clustered in the African, North American, and European lineages. The rest of the sublineages were scattered in diverse branches of the phylogenetic tree. The diversity of molecular evolution of SARS-CoV-2 cases imported to Guangdong Province was revealed, which was consistent with the case distribution (left part of Figure 3). In the homology analyses, a comparison of the reference sequences with the sequences of respective samples indicated that the homology among genomic sequences was 99.68%-100%. In addition, the viral genomes from infected patients were consistent with the animal reference sequences, including the Canis lupus family, feline, mink, ferret, Panthera leo, tiger, and Neovison vison (99.76%-99.98%). Moreover, we analyzed the correlations between age, gender, clinical classification, and VOCs, with a statistically significant difference between males and females (P < 0.05, not shown in the table).

besjournal.com) lists all substitutions in the nucleotides and amino acids.
Supplementary Tables S2 and S3 (available in www.besjournal.com) show the nonsense and missense mutation distributions in the nonsynonymous mutation. Additionally, Supplementary Tables S4 and S5 (available in www.besjournal.com) list the reported amino acid mutation loci on the Global Initiative of Sharing All Influenza Data (GISAID) website (https://www.gisaid.org/) and publications in PubMed (https://pubmed.ncbi.nlm. nih.gov/), respectively. In total, 1,076 nonsense mutations were observed on 476 nucleotide sites in all the analyzed sequences. A total of 2,010 missense mutations were present in 441 loci, and they were largely disseminated in positions ORF 1a and S. The missense mutations were demonstrated in all 12 gene fragments (right part of Figure 3). We also identified some missense mutations that had neither been reported on the GISAID database nor published in PubMed, including position non-S (replacements at five amino acid residue lineages, including those from the North American (sublineage B.1.

DISCUSSION
We narrated the complete sequences of 117 patients with SARS-CoV-2 infection from the Guangzhou Eighth People's Hospital, Guangdong, China, during the whole study period (June 2020 to May 2021). Based on our basic description and phylogenetic analysis of imported cases, the median age of patients in this study was younger than that of patients in earlier studies [5][6][7] , and we noticed a close identity between the genome of the respondents and animal reference sequences [including Canis lupus family, feline, mink, ferret, Panthera leo, tiger, and Neovison vison (99.76%-99.98%)]. Similar findings were described previously in other studies (homology between rhinolophine, bat, pangolin, and human) [8][9][10][11] . New variants with several critical mutations, such as B. 1.1.7, B.1.351, B.1.617, and P.1, to which the available vaccines or potential drugs may be ineffective, emerged in 2021 [12][13][14] .    [15] . This variant was first isolated and sequenced in Guangzhou in May 2021 [16] . Sublineages B.1.17 and B.1.351 were considered to be the main variants. They had been affiliated with the increased transmissibility of COVID-19 epidemiology and virulence, changes in clinical disease presentation, and reduced effectiveness of available vaccines and therapeutics [17,18] . In addition, we discovered the B.1.427/B.1.429, P.2, and B.1.525 variants in our studied cases. They caused community transmission or multiple clusters of COVID-19, which should arouse our attention [19][20][21][22] . Previous studies showed that D614G endowed a replication advantage to SARS-CoV-2 by increasing the feasibility of human-to-human transmission [23][24][25] .
Given the limitations of epidemiological and clinical data, with a conservative interpretation method, we can only discover the D614G in the spike protein that presented the highest mutation density, presenting another limitation of our research. According to current studies, the ORF7b gene of SARS-CoV-2 with a procaspase-3 cleavage activity has a prerequisite role in viral replication [17,26] . It potentially impedes certain cellular processes that dominate several common symptoms of SARS-CoV-2 infection, such as leucine zipper formation and epithelial cell-cell adhesion, which may contribute to heart-rate disorders [26][27][28] . The ORF9b of SARS-CoV-2, encoding a 97-amino acid protein on the mitochondrial membrane, was located within the nucleocapsid gene (N gene) [27] . Previous studies revealed that SARS-CoV-2 ORF9b, as an accessory protein, inhibits innate immunity by targeting mitochondria. It can potentially suppress the antiviral response of infected cells and regulate the host immune response by suppressing type I interferon synthesis [29][30][31] . Among the various missense mutations in the 117 sequences, we disclosed some missense mutations in several sublineages (such as B.1, B.1.1.25, B.1.1.7, B.1.2,  B.1.214.2, B.1.402, B.1.595, B.1.618, and L3) of partial lineages, including those from North America, Europe, Asia (excluding China), and Africa. They had neither been reported in the GISAID database nor published on PubMed and included five amino acid substitutions of the non-S positions and eight substitution sites in the S position. Given our lack of access to the information on clinical symptoms and dissemination chain in patients infected with the corresponding sublineages, we cannot prove whether these novel mutations contribute to its pathological functions (like serious symptoms: loss of speech or mobility, or confusion), spread, or adaptation to particular surroundings, including the age, pre-existing medical conditions, and living environment of the host. Our findings nevertheless provide important resources for tracking SARS-CoV-2 single-nucleotide variations, lineages, and clades in the imported cases of Southern China.

CONCLUSION
These imported cases in Guangdong, China, were identified, sequenced, and characterized epidemiologically and genetically, portraying the current studies on the molecular evolution of SARS-CoV-2. The descriptive epidemiology approaches  demonstrated the distribution of this virus. Our phylogenetic and homology analyses increased the awareness of the detailed description of this virus. Based on the results of our viral sequence analysis, sublineage diversity and high homology at the amino acid level were present in the imported cases in Guangdong. This finding confirmed the result of a previous study reporting that SARS-CoV-2 exhibited genetic diversity over time [32] . All our results provide an important step in the search for viral variants with reduced pathogenicity or without potential pathogenicity.

AUTHOR CONTRIBUTIONS
No potential conflicts of interest were disclosed. LIANG Dan performed the experiment, drew figures, wrote the manuscript, analyzed phylogenetic relationships, completed statistical tables, and interpreted the results. WANG Tao experimented and undertook fragment assembly. ZHANG Guan Ting and LIANG Yu Feng were in charge of data analysis and statistics and performed the literature search and the experiment. HONG Wen Shan and LI An An performed the experiment and recorded the results. LI Jia Jiao, GUAN Da Wei, DENG Xiao Ling, CHEN Feng Juan, PAN Xing Fei, and LEI Chun Liang were responsible for collecting samples and collating the basic information about patients. WANG Li and CHEN Meng Lin gave advice on the diagrams in this manuscript. JIA Hong Ling was in charge of editing the English text of a draft of this manuscript. KE Chang Wen designed the whole study framework and took responsibility for the integrity and accuracy of the data.