The study protocols were approved by the Ethics Committee of the First Affiliated Hospital of Anhui Medical University (Quick -PJ 2021-12-17). A total of 97 oropharyngeal swab samples were collected from 29 patients who presented with COVID-19 at the First Affiliated Hospital of Anhui Medical University between January 28 and March 8, 2020. The First Affiliated Hospital of Anhui Medical University is comprehensive (first-rate of Level three) hospital and provides complete medical care and health services for patients across Hefei as well as surrounding areas in the Anhui province. It is one of the eight designated hospitals for treating COVID-19 patients in Hefei. At present, the hospital is equipped with 4,990 beds and the annual outpatient volume is about 5 million. Confirmation of COVID-19 and clinical classifications were based on the protocols outlined in the New Coronavirus Pneumonia Prevention and Control Program (4^th Edition), which was published by the National Health Commission of China (37). This program specifies that to be considered as a confirmed COVID-19 case, patients must have detectable SARS-CoV-2 RNA in at least one respiratory sample since illness onset and exhibit acute respiratory infection syndrome and/or abnormalities on computed tomography (CT) scan of the chest. All clinical data on epidemiology (including exposure history), symptoms, underlying comorbidities, and laboratory results were retrospectively extracted from electronic medical records. An illness was considered to be serious when a critical illness notice was present in the medical record.

Oropharyngeal specimens were obtained with flocked swabs and placed in universal transport medium (Beijing Youkang Technology, Beijing, China) at 4°C until processed. All stored samples were processed within 6 h. Nucleic acid was extracted using a viral RNA extraction kit (Da An Gene Co., Ltd. affiliated with Sun Yat-sen University, Guangzhou, China). The presence of SARS-CoV-2 (N and ORF1ab genes) was detected using a SARS-CoV-2 RT-PCR kit (Da An Gene Co., Ltd.) on an Applied Biosystems 7500 system (Thermo Fisher Scientific, Waltham, MA, USA), according to the manufacturer's instructions.

Genome sequences were determined for the 97 SARS-CoV-2 RNA samples isolated as described above. The whole-genome sequence of SARS-CoV-2 was amplified using the Ion AmpliSeq™ DNA custom Panel WG00428_Coronavirus (Thermo Fisher Scientific), containing two pools with 121 primer pairs each. Sequencing libraries were prepared with the AmpliSeq Library Kit 2.0 (Thermo Fisher Scientific) according to the manufacturer's instructions, with barcoding of each sample. PCR amplification was performed as follows: 99°C for 2 min, 26 times (99°C for 15 s and 60°C for 4 min), followed by a hold at 10°C. The PCR amplicons were treated with 2 μL FuPa reagent to partially digest the primer sequences and were phosphorylated at 50°C for 10 min, followed by 55°C for 10 min and 60°C for 20 min. The amplicons were then ligated to adapters with diluted barcodes using the Ion Xpress Barcode Adapters kit (Thermo Fisher Scientific) for 30 min at 22°C and then 72°C for 10 min, followed by purification of the adapter-ligated amplicons (library) using the Agencourt AMPure XP reagent (Beckman Coulter, Brea, CA, USA). Library concentration was evaluated using Real-Time PCR Systems. Each diluted library (100 pM) was amplified through emulsion PCR using the OneTouch™ Instrument (Thermo Fisher Scientific) and enriched with the OneTouch™ ES Instrument (Thermo Fisher Scientific) using the Ion PI Hi-Q OT2 200 kit following the manufacturer's instructions. Finally, sequencing was performed on an Ion Proton instrument (Thermo Fisher Scientific) using an Ion PI Hi-Q Sequencing 200 Kit (Thermo Fisher Scientific). For genome assembly that had reference sequences available, sequencing reads were mapped to the reference using Burrows–Wheeler aligner (Bwa, version 0.7.12-r1039). Reads with excessive variations, which suggest artifacts, were removed from the dataset. Reads with mapped lengths shorter than 30 bp were also removed, and the soft-clipped bases were trimmed from both ends. For ampliseq data, depth was sufficient; redundant/duplicate reads were removed accordingly. Finally, genome assembly was performed with trinity (v1.2.9) with default parameters MEGAHIT (v1.2.9) or with parameter –k-min 15 (38–40).

A total of 231 genome sequences of betacoronaviruses isolated from mammals were downloaded from NCBI (Supplementary Data 1). The sequences were aligned with SARS-CoV-2 sequences obtained in this study using MAFFT. Phylogenetic trees were reconstructed using MEGA by Maximum Likelihood method. A bat coronavirus sequence (MN996532) was the most closely related sequence to SARS-CoV-2 and was thus used as an outgroup in subsequent analyses (44, 45). Since network reconstruction is extremely sensitive to missing data, only 48 samples (those missing <20% of the full-length genome sequence data) were included in the network analysis (Supplementary Data 1). To filter out rare mutations or alignment errors, network reconstructions used only parsimony-informative sites with binary polymorphisms. Two networks of haplotypes were generated in the NETWORK program using two data sets (46, 47). The first network used all the downloaded sequences from GenBank, as well as the sequences obtained in this study (Supplementary Data 1). Based on these results, sequences that did not distribute in the same branches as the study samples were removed, and the remaining sequences were used to create the second data set (Supplementary Data 1). In the NETWORK program, the median-joining network algorithm with the MP calculation option was used to reconstruct the most parsimonious network. Transversions were given a weight of three, while transitions were given a weight of one. In the network visualizations, nodes were proportional to frequencies of haplotypes and different colors indicated different regions. The network for the first data set was simplified using the star contraction option, as it was initially too complex to read. The data used in the first network was also used to reconstruct the phylogenetic tree of all SARS-CoV-2 sequences using MEGA by Maximum Likelihood method.

Levels of IL-6, IL-8, IL-10, IL-2R, and tumor necrosis factor-α (TNF-α) were assessed using a chemiluminescent immunoassay (Siemens, Munich, Germany); serum ferritin (SF) level was assessed using an electrochemiluminescence immunoassay (Roche, Mannheim, Germany) according to the manufacturer's instructions. The IMMULITE Cytokine Control Module, IMMULITE IL-10 Control Module (Siemens, Germany), and Lyphochek Immunoassay Plus Control (Bio-Rad, CA, USA) were used as internal controls, according to the manufacturer's instructions. The IMMULITE Cytokine Control Module is an assayed, bi-level control (containing different concentrations of selected lyophilized cytokines in a human serum matrix) intended for use with the IMMULITE1000 and IMMULITE 2000 IL-6, IMMULITE1000 IL-8, IMMULITE1000 IL-2R, and IMMULITE1000 TNFa assay. The IMMULITE IL-10 Control Module is a bi-level, synthetic matrix control intended for use with the IMMULITE1000 IL-10 assay. Lyphochek Immunoassay Plus Control (Bio-Rad) is a three-level, human serum matrix control intended for use with the Robas 6000 SF assay kit. The patients underwent multiple inflammatory factors tests during the hospitalization period. Based on the patient's conditions combined with the results of CT scan and nucleic acid test, the inflammatory factor data from samples obtained at the severest disease condition were selected for further analysis.

For each variable site, we divided patients into two groups: “reference” and “variation.” The “reference” group included patients carrying SARS-CoV-2 which was the same as reference genome at the site, whereas the “variation” group included patients carrying SARS-CoV-2 which mutated at the site. For 12 sites, both groups contained at least three patients, so totally 12 sites were used for subsequent analysis. The patients were also separated into “severe case” and “mild case” based on whether the medical record was present. Fisher's exact test was performed using SPSS software version 20.0 (IBM, Armonk, NY, USA) to test whether severe case percentage differed significantly (p < 0.05) between two mutation status. For each of the 12 sites, the expression levels of the six inflammatory factors between groups were compared using parametric t-tests using the SPSS software version 20.0 (IBM, Armonk, NY, USA).

Our analysis of 97 SARS-CoV-2 RNA samples, derived from 29 COVID-19 patients at various time points, identified 263 SNPs, while no indels were found (summarized in Table 1, Supplementary Table 1). Thirty-five of the identified SNPs were parsimony-informative sites and 228 were singleton SNPs. This ratio of parsimony-informative sites to singleton SNPs is comparable to that found in human genomes (48). Notably, there were 188 novel sites found in the study samples that were absent from all SARS-CoV-2 genomes accessed from the NCBI database for this study. All but six of the discovered SNPs were assigned to protein-coding regions. Specifically, non-synonymous, synonymous, and nonsense mutations encompassed 167, 80, and 10 sites, respectively. This suggests a high tolerance for sequence variation at the function - primarily protein-level. Most strikingly, sites 10380, 18060, and 28144 exhibited the highest mutation frequencies of all sites and appeared to be hotspots for mutation. SNPs at sites 10380, 18060, and 28144 were found in 24, 29, and 43 isolates, respectively. Previous studies have also identified SNPs at sites 18060 and 28144 in SARS-CoV-2 isolates from other geographical locations (33, 49, 50). Moreover, a non-synonymous mutation at site 28144 can give rise to a Leu to Ser substitution in the ORF8 protein, whereas an SNP at site 18060 in ORF1b is functionally silent. However, the SNP at site 10380 (harboring a G → T mutation, Gly → Val) has not yet been reported in any other isolates and is likely to account for a novel variation to Hefei. Thus, it may serve as a useful marker in further large-scale molecular epidemiological studies.

Haplotype network and phylogenetic analysis were carried out to infer the evolutionary relationship between regional samples. Only sequencing data with > 80% total sequence coverage was included in the haplotype network reconstruction, so only 37 haplotypes were identified from the 48 samples used to reconstruct the network (Table 2). All the haplotypes were novel and had not been previously described; five of them were associated with multiple samples.

The whole haplotype network (Figure 1A) could be separated into two big lineages, which were labeled Lineage A and Lineage B. The core of Lineage B and that of Lineage A were distinguished by two mutations: the synonymous mutation T8782C and the non-synonymous mutation C28144T changing a leucine to a serine. These two lineages were well-known as L-type and S-type, respectively, as also reported in other network studies (43, 51). Both core haplotypes were super-spreaders and were distributed in almost all areas included in this study.

Lineage A (S-type) showed a closer relationship with the root (an animal virus sequence, MN996532) that represents the ancestral origin. Sequence MT079847 was used as the representative core haplotype for Lineage A. MT079847 differed from the reference sequence by only one synonymous mutation, C8782T, and one non-synonymous mutation, T28144C (Leu → Ser). Although Lineage A had transmitted to all the areas we included, it was mainly distributed in China. Within Lineage A, only one sublineage with more than five haplotypes was discovered, which was labeled a1 (Figure 1). All these haplotypes had a common synonymous mutation C18060T and were geographically exclusive to Hefei and Yangtze River Delta, thereby forming an endemic cohort.

Lineage B was derived from Lineage A, which consisted of more haplotypes and longer branches. The core haplotype of Lineage B was identical to the reference sequence NC_045512. Within Lineage B, two sublineage s, each with more than five haplotypes, were discovered. Sublineage B1, which had common mutations 26144 or 11083, was mainly distributed in Hong Kong and Southeast Asia. Sublineage B2, which had a common mutations 241, was mainly distributed in Hong Kong and Japan. Based on this analysis and data records, Sublineage B2 emerged more recently than Lineage A or the other Lineage B sublineages. Moreover, the overlap in the geographical distribution of Lineage A and Lineage B indicates frequent traveling as an important accelerating factor for the spread of infection.

Although Lineage A was more ancestral than Lineage B, their transmission to human and global circulation occurred synchronously (Figure 1B) and the genome of Lineage B was sequenced first (2). In Chinese mainland, both lineage s accounted for almost one half of the samples, however, Lineage B was represented at much higher frequency than Lineage A in Hong Kong and Japan. Especially, a newly emerging Sublineage B2, which was mainly isolated in most samples in Hong Kong and Japan, was rare in Chinese mainland. By contrast, Lineage a1 was local concentrated distribution in Hefei and the Yangtze River Delta, indicating that this cohort is an endemic S-type variation of the virus.

The results of phylogenetic tree (Supplementary Figure 1) showed the same results as network. All the sequence could be split into two Lineage A and B. Within Lineage A, a distinguished Sublineage labeled “a1” was limited in Hefei and Yangtze River Delta. Within Lineage B, two distinguished Sublineage labeled “B1” and “B2” were found. While Sublineage B1 was mainly distributed in Hong Kong and Southeast Asia. With a particular long branch, Sublineage B2 was transmitted in Japan quickly and extensively. The root (MN996532) had a very long branch, but it is nested in Sublineage a1. As for the samples collected in Hefei, they are distributed in nine different evolutionary branches in both phylogenetic network (Figure 2) and phylogenetic tree, indicting the introduction of SAR-CoV-2 to Hefei happened at least nine times.

One practical application of phylogenetic network is to reconstruct potential pathways of human-to-human transmission of SARS-CoV-2. We investigated the transmission paths specifically associated with the infection in the Hefei region by focusing on locally obtained samples and comparing with other sequences distributed in same branches with them (Figure 2). The phylogenetic network analysis indicated that the initial infections likely occurred independently on at least nine separate occasions. Different viral types appeared to be introduced sequentially to the region, beginning with the L-type and followed by the S-type. A more detailed picture of potential transmission paths began to emerge when epidemiological data was incorporated into the analysis.

A minimum of six independent initial infections of the L-type virus likely occurred within a short period in January 2020 (labeled L1–L6 in Figure 2).

L1. Patient HF1 marked the first confirmed COVID-9 case admitted to our hospital. The patient traveled back from Wuhan, China on January 19, 2020, developed a fever on January 23, and was hospitalized on January 26. Patient HF3 had also returned from Wuhan and was hospitalized at approximately the same time as HF1. There was only one mutation that distinguished the HF1 isolate from the HF3, suggesting that HF1 and HF3 should be counted as one introduction event, which represents the earliest infections of the L-type SARS-CoV-2 in Hefei.

L2. Patient HF2 accounted for the second introduction event. The patient traveled back from Wuhan and was hospitalized on January 28, 2020.

L3. The third introduction occurred with patients HF4 and HF5. Both patients reported attending a social event with several individuals from Wuhan, soon developed COVID-19 symptoms, and were hospitalized on January 28. After 9 days, HF9 was admitted to our hospital. Although no direct contact was found between patient HF9 and patients HF4 and HF5, they shared similar haplotypes.

L4. The fourth introduction involved six cases (HF6, HF7, HF10, HF13, HF14, and HF15). Patient HF6 returned to Hefei from Wuhan on January 21 and was admitted to the hospital on January 28 with serious conditions. Soon after, patients HF7, HF10, HF13, HF14, and HF15 were hospitalized. Notably, patients HF10, HF13, and HF14 were from the same family and had no historical direct contact with anybody in Wuhan. This suggests the possibility of human-human transmission, a theory supported by evidence that one of the two viral haplotypes identified in the HF13 isolate was shared by HF14 and that the sequences of the HF10 and HF14 isolates differed by a single nucleotide.

L5 and L6. HF8 and HF12 account for two separate introduction events.

During the same period in early 2020, the S-type SARS-CoV-2 was introduced into Hefei on at least three independent occasions (labeled S1–3 in Figure 2).

S1. HF16 represented the first introduction, which occurred on January 21 when the patient-already exhibiting serious symptoms-traveled back to Hefei from Wuhan. The viral haplotype of this isolate has no direct relationship to any other case of S-type infections under study and was thus defined as a separate introduction.

S2. The second introduction involved nine individual cases (HF17, HF18, HF19, HF20, HF21, HF22, HF23, HF25, and HF26). With no clear historical direct contact with infection cases in Wuhan, patients HF17 and HF18 were admitted to the hospital on February 2 and 3, respectively, and represent the earliest patients in S2 (Figure 2) within our healthcare system. The viral haplotypes of the HF17 and HF18 isolates were differed by only one mutation and were closely associated with all other isolates of this cohort. Notably, HF21 and HF22 cases occurred in the same household, which also indicate a pattern of human-to-human viral transmission within this cohort. This second introduction of an S-type SARS-CoV-2 was responsible for nine individual cases involving viral haplotypes that shared a synonymous mutation, C18060T. These haplotypes constitute the majority of Sublineage a1, as highlighted in the haplotype network analysis (Figure 1A), which appears to be concentrated in Hefei and the Yangtze River Delta region, indicating that this cohort is an endemic S-type variant of the virus.

S3. HF24 accounted for the third S-type introduction event and was distinguished from other S-type isolates.

As previously mentioned, an illness was only considered to be serious if a critical illness notice was present in the medical record. The level of inflammatory factors and the grouping information for each patient are summarized in Supplementary Table 2. In total, 12 sites (8782, 10380, 11083, 13394, 14418, 16954, 17614, 18060, 26144, 26885, 28144, and 28253) were selected for a Fisher's exact test for comparing the rate of severe cases among different mutation status (Table 3). In this regard, no significant differences were observed among these sites, suggesting that the differences in mutation status might contribute to other viral functions and properties, such as transmissivity, or the sample size was insufficient to infer significant results.

Recent reports indicate that the initiation and progress of the sickness caused by COVID-19 are driven by cytokine (e.g., IL-6 and IL-8) responses. Thus, therapy that targets cytokines may improve the health of COVID-19 patients (9). Similarly, numerous studies have shown that high SF level plays a key role in inflammation and is significantly correlated to the severity of the disease (52). Therefore, we investigated whether mutation status at genomic sites (Table 3) could affect serum concentrations of inflammatory factors (Figure 3). We monitored levels of several inflammatory factors, including IL-6, IL-8, IL-10, IL-2R, TNF-α, and SF, in the patient blood samples and found that inflammatory responses were significantly reduced (p < 0.05) in patients infected with the isolates bearing mutations at sites 18060(harboring a C → T synonymous mutation, Leu) and 28253 (harboring a C → T synonymous mutation, Phe). Basing on these findings, we speculated that these mutations, especially those at 18060 (i.e., the local endemic variation), may facilitate viral transmission and contribute to global public health initiatives by exposing populations to asymptomatic or mildly symptomatic versions of this deadly virus.