Full-genome sequencing and mutation analysis of SARS-CoV-2 isolated from Makassar, South Sulawesi, Indonesia

Introduction

The first cases of coronavirus disease 2019 (COVID-19) were detected in Wuhan, China, on 31st December 2019 (Zhou et al., 2020; Oley et al., 2021). It was subsequently confirmed that these cases were due to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) (Coronaviridae Study Group of the International Committee on Taxonomy of Viruses, 2020), a respiratory virus that is highly transmissible between humans, and spread rapidly globally. The outbreak was declared a pandemic by the World Health Organization on 30th January 2020 (WHO, 2020b).

Phylogenic studies indicate that the virus first emerged in late November 2019, and has an evolutionary rate of approximately 1.1 × 10⁻⁵ subs/site/year (Duchene et al., 2020). Compared to other RNA viruses, coronaviruses (including SARS-CoV-2) have proofreading activity encoded by non-structural protein (nsp) 14, which reduces RNA replication errors (Sevajol et al., 2014; Smith et al., 2013; Wu et al., 2020). However, the highly transmissible nature of SARS-CoV-2 increases the likelihood of viral mutations and the rise of new variants in each infected individual.

Furthermore, variants that confer advantageous traits for viral survival (such as higher receptor binding capacity and immune evasion) are more likely to dominate the viral population, as observed in the case of the D614G variant (Korber et al., 2020; Lauring & Hodcroft, 2021). Although this variant was already present in China in sequences collected in late January to early February 2020, the mutation was first reported in the UK in early March 2020 (Korber et al., 2020). By June/July 2020, the D614G mutation was present in 63% of the total submitted sequences worldwide (Koyama, Platt & Parida, 2020). Viruses with the D614G variant have a higher viral load, but are not associated with disease severity (Korber et al., 2020). Over time, D614G diverged, leading to new variants in different parts of the world, including the B.1.1.7, B.1.351, P.1, B.1.617.2, and B.1.1.529 variants from the UK, South Africa, Brazil, India, and South Africa respectively. These variants are currently classified as variants of concern (VOC) by the World Health Organization. VOC are those associated with an observed increase in viral load, moderate immune evasion, and disease severity (WHO, 2021b; Robert Koch Institut, 2021; ECDC, 2022). The emergence of new impactful variants highlights the importance of the containment viral spread, genomic surveillance, and the development of broadly cross-protective vaccines (Walensky, Walke & Fauci, 2021; Abdool Karim & de Oliveira, 2021).

Indonesia announced the first two confirmed cases of COVID-19 on 2nd March 2020 (Alwan, 2007). By 1st December 2021, there were 4,256,687 (278 new within December 2021 alone) confirmed cases, 143,480 (10 new within December 2021 alone) deaths, and 4,104,964 recoveries, from 510 districts across all 34 Indonesian provinces (WHO, 2020a). However, the incidence of disease varies between areas, due to the archipelagic nature of the country and the spread of the virus. Genomic surveillance plays an important role in the containment of the virus, and in mitigating the emergence of new variants. From March 2020 to April 2021, 1,210 SARS-CoV-2 sequences were identified, accounting for 0.073% of the 1,668,360 infections recorded in Indonesia, with a monthly average of 140 days between specimen collection and sequence upload to Global Initiative on Sharing All Influenza Data (GISAID) (Pradipta et al., 2022). By 10th December 2021, a total of 10,094 SARS-CoV-2 genome sequences had been submitted to GISAID from Indonesia (GISAID, 2021).

Makassar is the capital city of the South Sulawesi province, which lies at the center of Indonesia. It has the largest economy in eastern Indonesia, due in particular to the Makassar Megapolitan region at its center, which acts as a busy hub for transportation throughout the country. As of 19th April 2021, a total of 50 SARS-CoV-2 whole-genome sequencing data collected from both the Wahidin Sudirohusodo and Universitas Hasanuddin hospitals had been submitted to GISAID (Ministry of Health of Republic Indonesia, 2021).

Preliminary data evaluation showed that the N439K variant with a D416G background dominated the SARS-CoV-2 viral population in Indonesia (Cahyani et al., 2021). This variant contains an N439K mutation in the receptor-binding motif (RBM), leading to enhanced binding affinity for the human angiotensin I-converting enzyme 2 (hACE2) receptor and resistance to several neutralizing monoclonal and polyclonal antibodies. These properties allow it to infect cells more easily and evade the immune system (Thomson et al., 2021). This study aims to: (1) report the full genome sequences of SARS-CoV-2 isolated from COVID-19 patients in Makassar, South Sulawesi, Indonesia, throughout the period of January to April 2021; and (2) perform mutation and phylogenetic analysis of SARS-CoV-2 from Makassar, South Sulawesi, Indonesia.

Materials & Methods

Results

Full genome characterization of SARS-CoV-2 isolated from Makassar, South Sulawesi

The phylogenetic tree of our original isolates (n = 50) and GISAID sequences (n = 31) is shown in Fig. 1. All isolates originating from Makassar were clustered in the main lineage of B.1, where 48 samples were identified from the B.1.466.2 lineage, and two (hCoV-19/Indonesia/SN-MKS_17/2021 and hCoV-19/Indonesia/SN-MKS_36/2021) from the B.1.319 lineage. However, 38 sequences of our isolates formed new clades which were different from other samples from the GISAID reference lineage of B.1.466.2 (EPI_ISL_877456; EPI_ISL_768689; EPI_ISL_728184). All isolates from Makassar were distantly related to other VOC lineages (Alpha/B.1.1.7, Beta/B.1.351, Gamma/P.1, Delta/B.1.617.2, Omicron/B.1.1.529).

Table 1:

The clinical characteristics of patients with COVID-19 in Makassar, South Sulawesi, Indonesia, from January to April 2021.

Variable	N (%)
Gender
Male	28 (56)
Female	22 (44)
Age (year)
0–4	4 (8)
5–17	1 (2)
18–29	9 (8)
30–39	15 (30)
40–49	5 (19)
50–64	10 (20)
65–74	4 (8)
75–84	2 (4)
Clinical severity
Asymptomatic	11 (22)
Mild	17 (34)
Moderate	17 (34)
Severe	5 (10)
Comorbidities*
Hypertension	1 (10)
Coronary Artery Disease	1 (10)
AV Block	1 (10)
Coagulopathy	2 (10)
Malnutrition	1 (10)
Elevated Liver Enzyme	1 (10)
Anxiety Disorder	1 (10)

DOI: 10.7717/peerj.13522/table-1

Notes:

*Data available in 10 samples.

Mutation analysis of the S protein

Mutation analysis of the S protein region showed that the most frequently occurring mutations in our study were D614G (50/50, 100%), followed by N439K (49/50, 98%) and P681R (38/50, 76%) (Fig. 2A and Table 2). Other mutations were found at lower frequencies, including L5F, V143F, Y144del, D253G, S255F, A348S, Q677H, A1070S, and P1263L. Mapping of the identified mutations showed they were predominantly in the S1 subunit (10 mutations), while two were located in the S2 subunit, and two others (A348S and N439K) were found in the receptor-binding domain (RBD). D641G, P681R, and Q677H mutations were located close to the furin cleavage site (FCS) region (Fig. 2B). Four mutations (V143F, Y144del, D253G, and S255F) were identified in the N-terminal domain (NTD) of the S1 subunit (Fig. 2A). The distribution of these specific mutations in our collected samples are provided in Table S2.

Table 2:

The dominant mutations in SARS-CoV-2 proteins identified in this study.

Mutations with a frequency of > 10% in our collected samples are included.

No	Region	Mutation	Frequency [n (%)]
1	nsp3-ORF1ab	S126L	38 (76)
		T350I	50 (100)
		P822L	49 (98)
2	nsp6-ORF1ab	L75F	38 (76)
		V149F	6 (12)
3	nsp12-ORF1ab	P323L	50 (100)
4	nsp13-ORF1ab	S259L	38 (76)
5	Spike	N439K	49 (98)
		D614G	50 (100)
		P681R	38 (76)
6	NS3-ORF3a	Q57H	49 (98)
7	Nucleoprotein	T205I	50 (100)

DOI: 10.7717/peerj.13522/table-2

Discussion

The COVID-19 pandemic continues to present a pressing challenge to global health. Despite possessing a proofreading enzyme, SARS-CoV-2 mutates and adapts rapidly, resulting in new variants with enhanced fitness. By 29th March 2022, the World Health Organization had identified a total of five VOCs, including previously circulating VOCs (Alpha, Beta, Gamma) and VOCs in current circulation (Delta and Omicron) (WHO, 2021b). VOCs tend to demonstrate increased transmissibility, and mitigation measures, such as vaccines, therapeutics, diagnostics, and other public or social health measures, are generally less effective for them (Thomson et al., 2021). Although the parental SARS-CoV-2 strains overwhelmingly affected older age groups in particular, our results show that the virus also affects adolescents and young adults in comparable proportions. Global genomic surveillance of SARS-CoV-2 plays an important role in mitigating and monitoring the emergence and spread of new variants, as well as in providing information for public health measures (WHO, 2021; WHO, 2021a).

This study reports SARS-CoV-2 genomic surveillance using whole-genome sequencing in Makassar, the capital city of the South Sulawesi province of Indonesia, from January to April 2021. We performed phylogenetic and mutation analysis for in-depth characterization of the SARS-CoV-2 variants circulating in this region.

Phylogenetic analysis showed that most isolates originated from the B.1.466.2 lineage, while the others (two samples) belonged to the B.1.319 lineage. These results were higher than those of another study carried out in West Java, another province in Indonesia, between September 2020 to June 2021, which found that 27% of isolates were of the B.1.466.2 lineage (Fibriani et al., 2021). B.1.466.2 variants are more locally dispersed in Indonesia, and have been detected at lower rates in other Southeast Asian countries (Cahyani et al., 2021) B.1.466.2 is one of the Indonesian lineages identified by PANGO, and although it has not been considered a cause for concern on a global scale, it has been designated on alert for further monitoring by the World Health Organization (WHO, 2021b).

The results of the mutation analysis show that ORF1ab contains the greatest number of amino acid mutations, followed by the ORF3a, Spike, and N. The ORF1ab combines codes of a polyprotein that are cleaved into 16 NSPs, which play a variety of roles in viral replication (Yoshimoto, 2020). Amino acid mutations are also able to affect the structure of proteins, which can impact their function in viral replication. The ORF1ab consists of two amino acid mutations (T350I and P323L) observed in 100% of the genome sequences in the cohort. The ORF3a gene, which encodes an ion channel protein related to NLRP3 inflammasome reactivation (Yoshimoto, 2020), was found to have a total of 16 amino acid mutations in the cohort, with the Q57H mutation present in 49 (98%) sequences. The nucleocapsid (N) protein encoded by the N gene encapsulates the RNA genome of SARS-CoV-2 (Yoshimoto, 2020), and all 50 sequences collected in this cohort were found to contain the T205I mutation.

The S protein, encoded by the S gene, is a SARS-CoV-2 surface glycoprotein of particular interest, as it plays an essential role in viral entry into the host cell by binding to cell surface angiotensin-converting enzyme 2 (ACE2) proteins. The S protein is a target for antibody neutralization; thus, amino acid mutations in antibody binding sites within the S protein allow for the evasion of antibody-based host defense mechanisms (Walls et al., 2020). In addition, the RNA sequence of the S region is employed as the target of the primer used in the PCR viral detection method. Consequently, nucleic acid mutations in these areas will result in false negative PCR results (Wang et al., 2020; Vogels et al., 2020). All genomic sequences collected in this cohort contain the background D614G (B.1. lineage) mutation, and 98% of them contain the N439K mutation.

The D614G mutation quickly became prominent in SARS-CoV-19 strains (Korber et al., 2020; Zhou et al., 2022; Viedma et al., 2021), and remains present in VOCs circulating currently. D614G was detected in early cases of COVID-19 in several Europe regions, and quickly became dominant (55–85% by the end of April 2020) (Viedma et al., 2021; Tang et al., 2020). In East Asia, the transition from D to G was slower, with almost 90% of viral samples from Chinese patients still displaying D614 at the end of April 2020 (Tang et al., 2020). In Indonesia, D614G was reported to be dominant (94.2%) until June 2021 (Cahyani et al., 2021), and this is consistent with our findings in Makassar city during our study period of January–April 2021. In other Indonesian provinces (namely Yogyakarta and Central Java) the D614G variant circulating since May 2020 finally replaced the D614 variant completely as of September 2020 (Wijayanti et al., 2022). Several studies have attempted to examine the effects of this mutation, and results suggest that it is a crucial to the FCS region. Structural and biochemical research on a full-length G614 spike trimer revealed that there are connections in D614G not seen in D614, which prevent the S1 subunit from prematurely losing its binding to the ACE2 (Zhang et al., 2021). Other reports found that the change of amino acid at position 614 from aspartic acid (D) to glycine (G) results in a greater ability to penetrate ACE2-expressing cells and replicate in human respiratory cells (Plante et al., 2021; Zhang et al., 2020). This ability essentially increases the number of spikes that can allow infection to occur, and indicates that the D614 viral spike is less stable than the G614 variant. The G614 mutation is thought to confer a moderate degree of infectivity and transmissibility, which could account for its persistent circulation worldwide (Korber et al., 2020; Plante et al., 2021; Zhang et al., 2020).

Similarly to B.1.466.2, the N439K mutation was found in moderately high frequency in this study. The N439K mutation was first detected in Scotland in March 2020 and increased in prevalence until June 2020, before disappearing entirely from the country following stringent public health measures. Despite this, the mutation has reappeared independently several times in other parts of the world. It was first sampled in Romania in May 2020, and as of January 2021 had been detected in 34 countries as the second most observed receptor-binding domain (RBD) mutation worldwide (Thomson et al., 2021). The N439K mutation was most dominant in Indonesia between January and April 2021, and formed its own characteristic lineage (B.1.466.2) in the Indonesian population (Rambaut et al., 2020; O’Toole et al., 2021a; O’Toole et al., 2021b), as reflected in the samples obtained for this study. Experiments involving recombinant N439K, hACE2 proteins, and X-ray structural analysis have shown that the N439K mutation results in a 2-fold increase in hACE2 binding affinity (Thomson et al., 2021). The effect of this on transmissibility and infectivity measures requires further evaluation. The clinical outcomes and viral replication of N439K mutants in cell tissue culture were comparable to those of N439 SARS-CoV-2 (Thomson et al., 2021). The N439K mutants demonstrated more significant 2-fold reduction in response to polyclonal antibody in a small proportion (6.8%) of sera collected from 442 individuals who recovered from COVID-19. A reduction in binding of more than 2-fold was also demonstrated in 16.7% of 140 mAbs isolated from recovered individuals (Thomson et al., 2021). In addition, Zhou et al. (2021) found that the N439K RBD was significantly resistant to the SARS-CoV-2 neutralizing antibody REGN10987, which could potentially result in neutralization failure. This information is important for ensuring the efficacy of vaccines and immunotherapy. Amino acid changes appear to accumulate gradually over time. N439K is the most prominent mutation in the RBD domain, and should be used preferentially to determine the effect of mutations on pathogen evolution.

Another amino acid mutation that should be highlighted is found in the spike protein–P681R. The frequency of P681R was found to be moderately high in Makassar during the study period, compared to data from the rest of Indonesia (42.6%) and the world (2.1%). The potential causes of this, including natural conditions, lifestyle, genetic factors of hosts, and adherence to health protocols, should be further investigated. Studies on viruses with P681R have revealed that the mutation could result in a reduction in neutralizing antibodies. Furthermore, it has been shown to increase cell–cell fusion by promoting S protein cleavage, mediated by furin (Saito et al., 2021). Interestingly, P681R was found frequently in the Delta variant, based on data obtained from GISAID from June to August 2021 (92.7% of all data submitted to GISAID and 99.6% of all Delta variants submitted to GISAID). This mutation may have benefitted from a cumulative fitness advantage over other variants, due to its high frequency in the Indonesian B.1.466.2 lineage.

These findings highlight the necessity of ongoing surveillance to identify novel COVID-19 mutations. The ability to detect new variants and their spread is critical for appropriate public health response and policy making.

Conclusion

We report the full-genome sequences of SARS-CoV-2 collected from Makassar, South Sulawesi, Indonesia, between January and April 2021. Lineage B.1.466.2 represents the most dominant (96%) of our isolates, and none of these variants are classified as VOC. All SARS-CoV-2 isolates were found to contain the D614G mutation, while the N439K mutation was also found in high proportions (98%) compared to the lineage in the total Indonesian SARS-CoV-2 sequences submitted to GISAID. Our study therefore highlights the importance of continuous genomic surveillance of SARS-CoV-2, particularly in the South Sulawesi region, as well as in other parts of Indonesia, to monitor the emergence and spread of novel mutations and variants with possible impacts on global health.

Supplemental Information