\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{Sequence analysis of the Spike, nsp12 (RNA-dependent RNA polymerase),
nsp3 (PLpro), and nsp5 (3CLpro) genes reveals a distinct evolutionary
pattern of SARS-CoV-2 variants circulating in Yogyakarta and Central
Java provinces, Indonesia}



\author[1]{Mohamad Hakim}%
\author[1]{* Gunadi}%
\author[1]{Ayu Rahayu}%
\author[2]{Hendra Wibawa}%
\author[1]{Laudria Stella Eryvinka}%
\author[1]{Endah Supriyati}%
\author[1]{Khanza Adzkia Vujira}%
\author[3]{Kristy Iskandar}%
\author[1]{* Afiahayati}%
\author[1]{Edwin Widyanto Daniwijaya}%
\author[1]{Farida Nur Oktoviani}%
\author[1]{Luthvia Annisa}%
\author[1]{Fadila Dyah Trie Utami}%
\author[1]{Verrell Christopher Amadeus}%
\author[1]{Setiani Silvy Nurhidayah}%
\author[1]{Tiara Putri Leksono}%
\author[1]{Fiqih Vidiantoro Halim}%
\author[4]{Eggi Arguni}%
\author[1]{Titik Nuryastuti}%
\author[1]{Tri Wibawa}%
\affil[1]{Universitas Gadjah Mada}%
\affil[2]{Kementerian Pertanian Republik Indonesia}%
\affil[3]{Universitas Gadjah Mada Fakultas Kedokteran Kesehatan Masyarakat dan Keperawatan}%
\affil[4]{Universitas Gadjah Mada Departemen Ilmu Kesehatan Anak}%


\vspace{-1em}



  
  \date{May 21, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
During the Covid-19 pandemic, the resurgence of SARS-CoV-2 was due to
the development of novel variants of concern (VOC). Thus, genomic
surveillance is essential to monitor continuing evolution of SARS-CoV-2
and to track the emergence of novel variants. In this study, we
performed phylogenetic, mutation, and selection pressure analyses of the
\emph{Spike, nsp12, nsp3,} and \emph{nsp5} genes of SARS-CoV-2 isolates
circulating in Yogyakarta and Central Java provinces, Indonesia from May
2021 to February 2022. Various bioinformatics tools were employed to
investigate the evolutionary dynamics of distinct SARS-CoV-2 isolates.
During the study period, 213 and 139 isolates of Omicron and Delta
variants were identified, respectively. Particularly in the \emph{Spike}
gene, mutations were significantly more abundant in Omicron than in
Delta variants. Consistently, in all of four genes studied, the
substitution rates of Omicron were higher than that of Delta variants,
especially in the \emph{Spike} and \emph{nsp12} genes. In addition,
selective pressure analysis revealed several sites that were positively
selected in particular genes, implying that these sites were
functionally essential for virus evolution. In conclusion, our study
demonstrated a distinct evolutionary pattern of SARS-CoV-2 variants
circulating in Yogyakarta and Central Java provinces, Indonesia.%
\end{abstract}%



\sloppy


\textbf{Sequence analysis of the\emph{Spike, nsp12 (RNA-dependent RNA
polymerase),} \emph{nsp3 (PLpro),}and \emph{nsp5 (3CLpro)} genes reveals
a distinct evolutionary pattern of SARS-CoV-2 variants circulating in
Yogyakarta and Central Java provinces, Indonesia}

Mohamad Saifudin Hakim\textsuperscript{1,*}, Gunadi\textsuperscript{2},
Ayu Rahayu\textsuperscript{1}, Hendra Wibawa\textsuperscript{3}, Laudria
Stella Eryvinka\textsuperscript{2}, Endah Supriyati\textsuperscript{4},
Khanza Adzkia Vujira\textsuperscript{2}, Kristy
Iskandar\textsuperscript{5}, Afiahayati\textsuperscript{6}, Edwin
Widyanto Daniwijaya\textsuperscript{1}, Farida Nur
Oktoviani\textsuperscript{1}, Luthvia Annisa\textsuperscript{1}, Fadila
Dyah Trie Utami\textsuperscript{2}, Verrell Christopher
Amadeus\textsuperscript{2}, Setiani Silvy Nurhidayah\textsuperscript{2},
Tiara Putri Leksono\textsuperscript{2}, Fiqih Vidiantoro
Halim\textsuperscript{2}, Eggi Arguni\textsuperscript{7}, Titik
Nuryastuti\textsuperscript{1}, Tri Wibawa\textsuperscript{1}

1 Department of Microbiology, Faculty of Medicine, Public Health and
Nursing, Universitas Gadjah Mada, Yogyakarta, Indonesia

2 Pediatric Surgery Division, Department of Surgery and Genetics Working
Group/Translational Research Unit, Faculty of Medicine, Public Health
and Nursing, Universitas Gadjah Mada, Yogyakarta, Indonesia

3 Disease Investigation Center Wates-Yogyakarta, Directorate General of
Livestok Services, Ministry of Agriculture, Indonesia

4 Centre for Tropical Medicine, Faculty of Medicine, Public Health and
Nursing, Universitas Gadjah Mada, Yogyakarta, Indonesia

5 Department of Child Health and Genetics Working Group, Faculty of
Medicine, Public Health and Nursing, Universitas Gadjah Mada/UGM
Academic Hospital, Yogyakarta, Indonesia

6 Department of Computer Science and Electronics, Faculty of Mathematics
and Natural Sciences, Universitas Gadjah Mada, Yogyakarta, Indonesia

7 Department of Child Health, Faculty of Medicine, Public Health and
Nursing, Universitas Gadjah Mada/Dr. Sardjito Hospital, Yogyakarta,
Indonesia

\textbf{*Corresponding author:}

\textbf{Mohamad S. Hakim, M.D., Ph.D.}

Department of Microbiology, Faculty of Medicine, Public Health and
Nursing

Universitas Gadjah Mada, 55281 Yogyakarta, Indonesia

Email: m.s.hakim@ugm.ac.id

\textbf{Running title:} Sequence analysis of the \emph{Spike, nsp12
(RdRp),} \emph{nsp3 (PLpro),} and \emph{nsp5} \emph{(3CLpro)} genes of
SARS-CoV-2 variants



\textbf{Abstract}

During the Covid-19 pandemic, the resurgence of SARS-CoV-2 was due to
the development of novel variants of concern (VOC). Thus, genomic
surveillance is essential to monitor continuing evolution of SARS-CoV-2
and to track the emergence of novel variants. In this study, we
performed phylogenetic, mutation, and selection pressure analyses of the
\emph{Spike, nsp12, nsp3,} and\emph{nsp5} genes of SARS-CoV-2 isolates
circulating in Yogyakarta and Central Java provinces, Indonesia from May
2021 to February 2022. Various bioinformatics tools were employed to
investigate the evolutionary dynamics of distinct SARS-CoV-2 isolates.
During the study period, 213 and 139 isolates of Omicron and Delta
variants were identified, respectively. Particularly in the \emph{Spike}
gene, mutations were significantly more abundant in Omicron than in
Delta variants. Consistently, in all of four genes studied, the
substitution rates of Omicron were higher than that of Delta variants,
especially in the \emph{Spike} and \emph{nsp12} genes. In addition,
selective pressure analysis revealed several sites that were positively
selected in particular genes, implying that these sites were
functionally essential for virus evolution. In conclusion, our study
demonstrated a distinct evolutionary pattern of SARS-CoV-2 variants
circulating in Yogyakarta and Central Java provinces, Indonesia.

\textbf{Keywords}

3CLpro; Delta variant; nsp12; Omicron variant; PLpro; RdRp; SARS-CoV-2;
Spike

\emph{}

\textbf{Introduction}

During the current pandemic, the resurgence of severe acute respiratory
syndrome coronavirus 2 (SARS-CoV-2) cases was due to the development of
novel variants of concern.\textsuperscript{1} Continuous viral
transmission globally has indeed enabled the virus to fine-tune its more
efficient adaptation in human populations, enabling its sustained
transmission.\textsuperscript{2,3} Genomic surveillance is thus
essential to monitor the continuous evolution of SARS-CoV-2 and to track
the emergence of novel variants or sub-variants. Since January 2020, a
collaborative network of the World Health Organization (WHO), national
authorities, and research institution have been tracking the evolution
of SARS-CoV-2.\textsuperscript{4}

The \emph{Spike} (S) gene of SARS-CoV-2 is under intense selective
pressure and thus, it is highly mutated. These mutations may result in
increased infectivity, transmissibility, and more resistant to
neutralizing antibodies of the novel
variants.\textsuperscript{5}Consistently, the emergence of variants of
interest (VOI) or concern (VOC) was characterized by unique sets of
mutations in the \emph{Spike}gene.\textsuperscript{2} Currently, the
only VOC is the Omicron variant, first identified in South Africa in
early November 2021.\textsuperscript{4} In fact, the emergence of the
Omicron variant that successfully replaced the previously dominant Delta
variant was also marked by a number of unique mutations in the
\emph{Spike}gene.\textsuperscript{1,6,7}

Multiple mutations and recombination events in the \emph{Spike} gene
have led to the diversification of the Omicron variant into various
sub-lineages or sub-variants with possibly different phenotypic
characteristics.\textsuperscript{8,9} Few subvariants, including BQ.1
and XBB, may spread quickly across countries and evolve to be the most
resistant variant against neutralizing
antibodies.\textsuperscript{10,11} Possible recombination events between
Omicron and Delta variants have also been
identified.\textsuperscript{12}

In addition to the \emph{Spike} gene, mutations in the gene encoding the
non-structural proteins also play an essential role in SARS-CoV-2
evolution and adaptation.\textsuperscript{13} As a positive-sense RNA
virus, SARS-CoV-2 encodes its own enzyme to replicate the genome. The
primary RNA polymerase, namely RNA-dependent RNA polymerase (RdRp), is
encoded by \emph{nsp12} gene\textsuperscript{14} which forms the core
polymerase complex with nsp7 and nsp8 proteins.\textsuperscript{15}
The\emph{nsp12} gene of SARS-CoV-2 contains both RdRp and N-terminal
nidovirus RdRp-associated nucleotidyltransferase (NiRAN)
domain.\textsuperscript{16} Several mutations have been identified in
the \emph{nsp12} gene, including S6L, P323L, and P323F, although the
effects of these mutations are currently
unclear.\textsuperscript{13,17}However, it has been shown that the P323L
mutation may increase the mutation rate of
SARS-CoV-2.\textsuperscript{18} RdRp has also been the target of various
SARS-CoV-2 inhibitors, such as remdesivir, although comprehensive
analyses of the global SARS-CoV-2 genomic dataset showed that potential
remdesivir-escape mutations were very rare, indicating its little
selective pressure.\textsuperscript{19} Indeed, the RdRp region
demonstrated strong negative (purifying) selection during the
pandemic.\textsuperscript{20}

Other two important proteins for the SARS-CoV-2 replication cycle are
proteases, encoded by \emph{nsp3} (papain-like protease, PLpro)
and\emph{nsp5} genes (chymotrypsin-like main protease, 3CLpro; also
known as main protease, Mpro).\textsuperscript{21} PLpro and 3CLpro have
been shown to modulate the host innate immune responses by inhibiting
interferon regulatory factor 3 (IRF3) that resulted in attenuated
interferon (IFN) response.\textsuperscript{22} Both proteins have also
been identified as the targets of IgM antibodies and were associated
with the survival of critical Covid-19 patients.\textsuperscript{23} In
addition, both PLpro and 3CLpro are potential targets for antiviral drug
development.\textsuperscript{24} A systematic mutational scanning has
identified several residues within the 3CLpro protein which are critical
for its functionality. These residues are mutation-sensitive and could
serve as promising targets for inhibitors with low likelihood of
resistant development.\textsuperscript{25,26} Nirmatrelvir (PF-07321332)
is an orally available antiviral drug with selective binding to the
active site of 3CLpro\textsuperscript{27} and has shown benefits for
treating Covid-19 patients.\textsuperscript{28} We and others have
identified various mutations in PLpro and 3CLpro
regions\textsuperscript{17,29}, which may pose challenge for antiviral
drug development.

Since the early Covid-19 pandemic, we and others have successfully
sequenced the full-genome of SARS-CoV-2 circulating in several regions
in Indonesia and performed genetic analysis of the isolated
viruses.\textsuperscript{30-40} These studies demonstrated that
B.1.466.2 lineage predominantly circulated in Indonesia during early
pandemic and only few isolates belonged to the B.1.319
lineage.\textsuperscript{35,38} The Delta variant emerged in April 2021,
subsequently replaced B.1.466.2 variant, and posed serious challenge to
public health control.\textsuperscript{35,37} Our previous study showed
that SARS-CoV-2 isolates detected during early pandemic had prominent
mutations in particular gene, including \emph{Spike, nsp3,}
and\emph{nsp12} genes.\textsuperscript{31,38}

Our current study aims to analyse the evolutionary pattern and mutation
rate of the \emph{Spike, nsp12, nsp3,} and \emph{nsp5} genes of various
SARS-CoV-2 variants circulating in Yogyakarta and Central Java
provinces, Indonesia. We also performed selective pressure analysis to
investigate which amino acid changes that are potentially maintained by
positive selection to provide insight of the importance of these sites
in virus evolution. This study provides more insights on genetic
variability of SARS-CoV-2 Delta and Omicron variants, particularly those
circulating in our regions (Yogyakarta and Central Java, Indonesia). To
our limited knowledge, our study is the first comprehensive and
comparative evolutionary analysis of Delta and Omicron variants
circulating in Indonesia.

\textbf{Materials and Methods}

2.1. Study design

The study was part of our previous retrospective study of Covid-19
patients in Yogyakarta and Central Java provinces,
Indonesia.\textsuperscript{36} Briefly, the study was conducted from May
2021 to February 2022. We collected nasopharyngeal swabs of outpatient
or hospitalized patients with Covid-19 for RT-PCR testing in the
Laboratory of Covid-19, Faculty of Medicine, Public Health and Nursing,
Universitas Gadjah Mada, Yogyakarta. During the study period, 352
confirmed Covid-19 patients (based on RT-PCR), consisting of 164 males
and 188 females, were recruited.

2.2. Whole-genome sequencing of SARS-CoV-2

We conducted whole-genome sequencing (WGS) of SARS-CoV-2 for
RT-PCR-positive samples with a Ct value of [?]30 based on our previous
studies.\textsuperscript{30-33,36} Briefly, SuperScript III
First-Strand Synthesis System (Thermo Fisher Scientific, MA, United
States) was used to synthesize single-stranded complementary DNA (cDNA)
from the total RNA. Subsequently, the second strand was synthesized
using Covid-19 ARTIC v3 primer pool design by SARS-CoV-2 ARTIC Network
using Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific,
MA, United States). The library preparations were performed using the
Illumina DNA Prep (Illumina, California, United States). The Illumina
MiSeq next-generation sequencer was used to perform the whole-genome
sequencing of SARS-CoV-2. The genomes of our samples were assembled and
mapped into the reference genome from the first Wuhan (China) isolate
(hCoV-19/Wuhan/Hu-1/2019, GenBank accession number: NC\_045512.2) using
Burrow-Wheeler Aligner (BWA) algorithm embedded in UGENE v. 1.30
(https://doc.ugene.net/wiki/display/UUOUM30/About+UGENE). The lineage of
each sequence was evaluated using Pangolin v4.1.2. The GISAID ID of each
submitted isolate generated from our study was shown in Supplementary
Table 1.

2.3. Phylogenetic analysis

A total of 213 and 139 cleaned sequences of SARS-CoV-2 Delta and Omicron
variants and five reference sequences of each VOC (Alpha, Beta, Delta,
Gamma, and Omicron variants, listed in Supplementary Table 2) were used
to construct four different Maximum Likelihood (ML) trees based on
the\emph{Spike, nsp12, nsp3,} and \emph{nsp5} genes. All reference
sequences were selected based on the earliest submission on GISAID
databases and identified as the respective variants. All ML trees were
generated using\emph{RaxmlGUI v2.0} with 1,000 bootstraps using the
GTR-Gamma model as the substitution model.\textsuperscript{41} The ML
trees were visualized and annotated using \emph{MEGA11}
.\textsuperscript{42}

2.4. Mutation analysis

Mutation analysis was performed using two datasets: 213 Omicron and 139
Delta sequences. Analysis was conducted specifically for each gene
(\emph{nsp3, nsp5, nsp12,} and \emph{Spike} genes). The reference
sequence was retrieved from GenBank with Accession Number NC\_045512.2
(Wuhan-Hu-1). Gene sequence extraction was performed using

\emph{samtools v.1.15.1.}

All sequences were mapped to the genome reference sequence
using\emph{minimap2 v.2.4} . Variant calling was performed
using\emph{Lofreq} . Mutation annotations were performed by \emph{snpEff
v.5.1} . Data processing was performed using the \emph{awk}
and\emph{grep} command.

2.5. Clock rate analysis

The Delta and Omicron sequences were extracted based on the region shown
in the Supplementary Table 4. The datasets were aligned using MAFFT
v.7.5.0 using local pair algorithm. The sequences were filtered by
removing the redundant sequence and dropping all sequences containing
``N''. The clean datasets were used for clock rate analysis.

Clock rate analysis parameter was created using \emph{Beauti
v.1.10.4}and Bayesian analysis was performed using \emph{BEAST v1.10.4}
. HKY-Gamma substitution model and uncorrelated log-normal relaxed
molecular clock were set for the run. Priors were set to classic priors.
The chain length is 4 x 50 million with sampling every 50 thousand. Log
files were combined using LogCombiner \emph{v1.10.4.} Run results were
evaluated using \emph{Tracer v1.7.1} . The ESS value higher than 100
shown acceptable quality run. Molecular Clock was observed from a report
generated by \emph{Tracer} .

2.6. Selective pressure analysis

Selective pressure analysis was performed using MEME in
Hyphy\emph{v.2.5.6.} 213 Omicron and 139 Delta sequences (\emph{nsp3,
nsp5, nsp12,} and \emph{Spike} genes) were subjected to selective
pressure analysis. Each dataset was first aligned using Neighbor-Joining
(NJ) algorithm with 1,000 bootstraps using Maximum Composite Likelihood
Model. The aligned sequences were analyzed using Hyphy and the results
were visualized using HyPhy Vision
(http://vision.hyphy.org/).\emph{p-value} lower than 0.1 was considered
as significant.

2.7. Mutation effect analysis

Mutation effect analysis was performed using
DynaMut.\textsuperscript{43} To compute the stability of SARS-CoV-2
wild-type (WT) and mutated (mutant) proteins, the reference structures
were retrieved from the RCSB database (https://www.rcsb.org/), with the
following IDs: 7CMD (Nsp3); 6M2Q (Nsp5); 6M71 (RdRp); 6VXX (Spike) and
were then uploaded to the DynaMut server
(https://biosig.lab.uq.edu.au/dynamut/). At specific site, a point
mutation was inserted. The effect in the form of total energies and
vibrational entropy energies between wild-type and mutated proteins was
recorded.

2.8. Ethical approval

The study was approved by the Medical and Health Research Ethics
Committee of the Faculty of

Medicine, Public Health, and Nursing, Universitas Gadjah Mada
(KE/FK/0563/EC/2020). Written informed consent was obtained from all
participants and parents or guardians for participating in this study.

\textbf{Results}

3.1. The genome sequencing quality

As mentioned, the GISAID IDs of submitted isolates generated from our
study as well as the date of specimen collections, place, age, and sex
were shown in Supplementary Table 1. The genome sequencing quality of
all sequences was shown in Supplementary Table 3.

3.2. Phylogenetic analyses of SARS-CoV-2 circulating in Yogyakarta and
Central Java provinces based on the Spike, nsp12, nsp3, and nsp5 genes

The Maximum likelihood (ML) phylogenetic tree of the \emph{Spike} gene
of our original isolates (n=352) is shown in Figure 1. Based on
the\emph{Spike} gene, all original viral isolates collected from
Yogyakarta and Central Java provinces were clustered in two main
lineages of the Delta and Omicron variants, while the previous VOCs
(Alpha, Beta, and Gamma) were separated in distinct clusters. Similarly,
the ML tree based on \emph{nsp12} (Figure 2A) and \emph{nsp5} genes
(Figure 2C) also showed that our isolates were segregrated into two
independent clusters of Delta and Omicron variants. While in \emph{nsp3}
gene, the Delta and Omicron variants did not form specific clusters, but
more diverse segregation was observed (Figure 2B).

3.3. Mutation analyses of the Spike, nsp12, nsp3, and nsp5 genes of
SARS-CoV-2 Delta and Omicron variants

The complete mutational analysis was shown in Supplementary Table 5.
High frequency (hotspot) mutations
(\textgreater{}10\%)\textsuperscript{44} identified in either Delta and
Omicron variants were shown in Figure 3, while visualization of these
mutations were depicted in Figure 4. Mutational events highly occurred
in the \emph{Spike} gene. The most frequent amino acid substitutions in
the Spike protein of Delta variant were T19R, G142D, A222V, L452R,
T478K, Q613H, D614G, P681R, D950N, and V1264L. While the Omicron variant
had more mutations than the Delta variant. The amino acid mutations
T478K, D614G, D950N, and N969K were found in Delta and Omicron variants.
However, the frequency of D950N mutation in the Omicron and N969K
mutation in the Delta was low (\textless{}10\%).

The Delta variant had more frequent amino acid substitutions than the
Omicron in the \emph{nsp12} gene. Several mutations found in Delta
variant were V42V, R249M, P323L, G671S, and N600N. While in Omicron
variant, the prevalent mutations were P323L, N600N, and L758L.
Interestingly, Delta and Omicron variants had similar amino acid
substitutions P323L with a frequency of 100\%. The amino acid
substitution in the \emph{nsp5} gene only occurred in the Omicron
variant. There are three amino acid changes with \textgreater{}10\%
frequency, including L27L, R131R, and P132H. In the \emph{nsp5} gene,
the amino acid substitutions occurred in both variants. However, Omicron
had more prevalent mutations than Delta, with nine amino acid
substitutions identified within this gene. F106F mutation occurred in
both variants with a high frequency (\textgreater{}95\%). For P1228L
mutation, the frequency was higher in the Delta (85,9\%) than in Omicron
(2,9\%) variants.

3.4. Substitution rates of the Spike, nsp12, nsp3, and nsp5 genes of
SARS-CoV-2 Delta and Omicron variants

We then determined the susbtitution rate for each gene of SARS-CoV-2
Delta and Omicron variants circulating in Yogyakarta and Central Java
provinces, Indonesia. In four genes studied, the substitution rates of
Omicron were higher than that of Delta variants, particularly in
the\emph{Spike} and \emph{nsp12} genes. The substitution rates for
the\emph{Spike} gene were 2.41 x 10\textsuperscript{-2} subs/site/yr
(95\% HPD: 6.82 x 10\textsuperscript{-4} to 5.67 x
10\textsuperscript{-2}) and 1.05 x 10\textsuperscript{-3} subs/site/yr
(95\% HPD: 3.65 x 10\textsuperscript{-6} to 3.92 x
10\textsuperscript{-3}) for Omicron and Delta variants, respectively.
For \emph{nsp12} gene, the substitution rates were 3.50 x
10\textsuperscript{-2} subs/site/yr (95\% HPD: 3.36 x
10\textsuperscript{-7} to 1.24 x 10\textsuperscript{-1}) and 3.44 x
10\textsuperscript{-3} subs/site/yr (95\% HPD: 1.25 x
10\textsuperscript{-6} to 9.49 x 10\textsuperscript{-3}) for Omicron and
Delta variants, respectively. Similar findings were found in
the\emph{nsp3} and \emph{nsp5} genes \textbf{(Table 1)} .

3.5. Selective pressure analysis

Selective pressure analysis was performed for each gene to comprehend
which sites that evolve under selective pressure, implying that the
sites are functionally essential for virus evolution. Selective pressure
analysis on \emph{Spike} gene revealed that several sites were under
positive selection, including 95 (T), 142 (G), 222 (A), 452 (L), 614
(D), 1264 (V) in Delta and 213 (V), 339 (G), 375 (S), 417 (K), and 440
(N) in Omicron variants. In \emph{nsp12} gene, only one site (249
{[}R{]}) was under positive selection in Delta variant. \emph{Nsp3} gene
revealed 1183 (A), 1206 (S), 1228 (P), and 1733 (E) that were under
positive selection in Delta variant. No sites in \emph{nsp5} gene that
were under positive selection. Thus, in Omicron variant, we only
identified positively selected sites in the \emph{Spike} gene.

3.6. The effect of mutations to the encoded proteins

Some common mutations in Delta and Omicron variants were analyzed to
examine their thermodynamic effect on the Spike protein (Supplementary
Table 6). Mutations in the \emph{Spike} gene of Delta variant, including
T95I, G142D, A222V, Q613H, and D614G, demonstrated a stabilizng effect
on the structure of the Spike protein. We then analyzed common mutations
in the \emph{Spike} gene of Omicron variant. We found that these
mutations have either a stablizing and destabilizing effect on the
structure of the Spike protein. In the \emph{nsp12} gene, amino acid
substitutions P323L found in both Delta and Omicron variants with a
frequency of 100\% has a stabilizing effect on the RdRp protein.

In the \emph{nsp3} gene, we analyzed three mutations (G28R, P77L, and
L120I) in Delta variant, which all have a stabilizng effect on the PLpro
protein. In Omicron variant, T259I mutation has destabilizng effect. For
the \emph{nsp5} gene, we analyzed L75F (Delta) as well as P132H and
P241L (Omicron) mutations, which all showed a stabilizing effect on the
3CLpro protein.

\textbf{Discussion}

Our current study showed that by analyzing the \emph{Spike, nsp12,
nsp3,} and \emph{nsp5} genes, we clearly demonstrated a distinct
evolutionary pattern of SARS-CoV-2 Delta and Omicron variants
circulating in Yogyakarta and Central Java provinces, Indonesia, during
May 2021 to February 2022. In Indonesia, the Omicron upsurge was
occurred from late January until February 2022 and subsequently replaced
the Delta variant.\textsuperscript{45} We have previously shown that
Delta- and Omicron-infected patients had similar hospitalization and
mortality rates.\textsuperscript{36} Our current study emphasized the
need of continued and extensive SARS-CoV-2 sequencing surveillance as
the primary method to quickly detect and respond the emergence of new
variants.

A latest study showed that compared to the previous VOCs (Alpha, Beta,
Gamma, and Delta), Omicron had the highest enrichment of amino acid
substitutions within the \emph{Spike} gene.\textsuperscript{46} Indeed,
new strains with higher transmissibility and infectivity have emerged
due to mutations in the Spike protein.\textsuperscript{2} Alterations of
the amino acid in the Spike protein influenced the binding affinity and
the viral fusion process.\textsuperscript{47} Our results showed that
D614G mutation occurred in all (100\%) Delta and Omicron variants
identified in this study, indicating that it was fixed in the viral
population. Indeed, all VOCs were generated from the G614 variant
lineage. Three mutations always accompany the D614G variant: a C-to-T in
the 5' UTR; a silent C-to-T mutation at position 3.037; and a C-to-T
mutation at position 14,408 that results in an amino acid change P323L
in the \emph{nsp12} gene.\textsuperscript{48} The other mutations at
position 477 (S477G, S477N, and S477R) of the Spike protein were
prominent among monoclonal antibody (mAb)-escape
mutation.\textsuperscript{2} Our results showed that S477N mutation was
found only in the Omicron variant (46.8\%). A computational analysis
showed that N477 had an increased binding affinity to
ACE2.\textsuperscript{49} Other studies reported that T19R, E156G,
L452R, T478K, and P681R variants in the Spike protein exhibited a
stabilizing effect on protein structure, facilitating the binding
affinity for more stable interactions with the human ACE2. This
stability effect may increase the transmission of the virus in human
populations.\textsuperscript{50} Our studies found T19R, L452R, T478K,
and P618R mutations in Delta with a high frequency above 70\% and L452R
in Omicron with 46.8\% frequency.

The nsp12-P323L mutation first appeared in January 2020 and became the
predominant globally (\textgreater{}90\%) by late April. Other mutations
were found in a low frequency (3.5-4.0\%), including E254D, A423V,
A656S, V720I, and V776L.\textsuperscript{51} It has been shown that
P323L mutation may increase the mutation rate of
SARS-CoV-2\textsuperscript{18} and was associated with Covid-19
severity.\textsuperscript{52} Molecular dynamic simulations showed that
P323L mutation led to tighter binding with antiviral drug remdesivir
(RDV).\textsuperscript{53} Interestingly, it has been shown that the
Spike D614G and RdRp P323L (G/L variant) have co-evolved to become more
superior than the original D/P variant.\textsuperscript{54} In our
study, the P323L mutation is similarly conserved between Delta and
Omicron variants, indicating its beneficial effect for viral evolution.

An early study analyzing SARS-CoV-2 isolates circulating in Indonesia
until September 2020 revealed mutations in the PLpro (P77L and V205I)
and 3CLpro (M49I and L50F) genes.\textsuperscript{34} However, we did
not find these mutations in our current study, suggesting that these
mutations were not fixed in the viral population circulating in our
region. Within the \emph{3CLpro} gene, several unique (signature)
mutations are identified, which are different across SARS-CoV-2
lineages. P132H, K90R, and G15S are prevalent mutations found in Omicron
(B.1.1.529), Beta (B.1.351), and Lambda (C.37) variants, respectively.
In contrast, Delta variant had no unique mutations within this
gene.\textsuperscript{55} Consistently, our finding showed that all
Omicron variant circulating in our region harboring P132H mutation,
while Delta variant had no prevalent mutation identified. One mutation
(L75F) found in Delta variant was very low in frequency (1.4\%).
Notably, P132H mutation did not lead to nirmatrelvir
resistance.\textsuperscript{55}

Several mutations have been identified to confer resistance to
nirmatrelvir, including Y54C, L167F, and
E166V.\textsuperscript{26,56}These mutations were not identified in our
isolates. Interestingly,\emph{in vitro} studies combining Y54C and L167F
with P132H resulted in a functional protein.\textsuperscript{56} This
finding provides an insight that additional mutations in the
\emph{3CLpro} gene that may emerge in the Omicron variant may be
clinically relevant in the future. In addition, mutations in
\emph{PLpro} and \emph{3CLpro} genes were shown to be associated with
the clinical course of Covid-19 patients. P108S mutation in the
\emph{3CLpro} gene resulted in reduced activity of the protease protein
and was associated with milder clinical course.\textsuperscript{57} In
the \emph{PLpro} gene, it has been shown that P78L and K233Q mutations
were associated with increased risk of death.\textsuperscript{58}
Interestingly, our current study found the synonymous mutation F106F in
the \emph{PLpro} gene both in Delta (97.7\%) and Omicron (99.3\%)
variants.

It has been estimated that the background nucleotide substitution rate
of SARS-CoV-2 was about 1.1 x
10\textsuperscript{-3}subs/site/year.\textsuperscript{59} The emergence
of VOCs can be ascribed to a 4-fold increase of the substitution rate
above its background level that may have lasted for several weeks or
months.\textsuperscript{60} In addition to their defining (unique)
mutations, each VOC may have a distinct evolutionary
rate.\textsuperscript{60} A previous study analyzing the phylodynamic of
Delta and Omicron variants also showed that a mean rate of nucleotide
substitution was higher in Omicron compared to Delta variants, i.e.
3.898 x 10\textsuperscript{-3} subs/site/year (range: 2.686 x
10\textsuperscript{- 3}~to 5.102 x 10\textsuperscript{- 3})~ and 3.677 \selectlanguage{ngerman}×
10\selectlanguage{english}\textsuperscript{- 4}~subs/site/year (range: 1.311 x
10\selectlanguage{english}\textsuperscript{- 4} to 6.144 x 10\selectlanguage{english}\textsuperscript{- 4}~) for the
Omicron and Delta variants, respectively.\textsuperscript{61}

Alterations in amino acids that decrease viral fitness are often
eliminated by negative selection, while changes that improve viral
fitness are retained by positive selection. On the other hand,
amino-acid alterations are regarded as ``neutral'' when they have no
effects on viral fitness. Because the existence of negative or positive
selection suggests that specific sites are functionally significant, it
is crucial to determine which sites evolve under selective pressure,
especially in the case of novel emerging pathogens such as SARS-CoV-2.
In this study, several positively selected sites in the \emph{Spike,
nsp12,} and \emph{nsp3} genes were identified in Delta variant. However,
in Omicron variant, we identified positively selected sites only in
the\emph{Spike} gene. This finding may imply a distinct genetic
``hot-spot'' in SARS-CoV-2 tropism, replication, and infectivity, and
need to be further investigated.

\textbf{Conclusion}

During the study period, we identified 213 and 139 isolates of Omicron
and Delta variants, respectively, co-circulating in our region.
Particularly in the \emph{Spike} and \emph{nsp5} genes, high frequency
amino acid substitutions were significantly more abundant in Omicron
than in Delta variants. Consistently, in all of four genes studied, the
substitution rates of Omicron were higher than that of Delta variants,
especially in the \emph{Spike} and \emph{nsp12} genes. In addition,
selective pressure analysis revealed several sites in particular genes
that were positively selected, implying that these sites were
functionally essential for virus evolution, including 95 (T), 142 (G),
222 (A), 452 (L), 614 (D), 1264 (V) in Delta and 213 (V), 339 (G), 375
(S), 417 (K), and 440 (N) in the \emph{Spike} gene of Omicron variants.
Our study demonstrated a distinct evolutionary pattern of SARS-CoV-2
Delta and Omicron variants circulating in Yogyakarta and Central Java
provinces, Indonesia. Thus, our study emphasized the need of continued
and extensive SARS-CoV-2 genomic surveillance to quickly detect and
respond the emergence of new variants.

\textbf{Acknowledgments}

We thank Balai Besar Teknik Kesehatan Lingkungan dan Pengendalian
Penyakit, Yogyakarta (BBTKL-PP) and Balai Besar Veteriner (BBVet),
Wates, Yogyakarta for providing the NGS Illumina MiSeq instrument
facility to sequence the virus isolates. We thank the Collaborator
Members of the Yogyakarta-Central Java Covid-19 study group: Fatin
Asfarina, Sri Fatmawati, and Sumartinigsih (Integrated Research
Laboratory, Faculty of Medicine, Public Health and Nursing, Universitas
Gadjah Mada {[}FKKMK UGM{]}), Sumardi, Heni Retnowulan, Eko Budiono,
Bambang S. Riyanto, Satria Maulana, Munawar Gani, Osman Sianipar, Ira
Puspitawati, (FKKMK UGM/Dr. Sardjito General Hospital), Herdiyanto
Mulyawan and Lestari (BBVet Wates, Yogyakarta), Havid Setyawan and
Indaryati (BBTKL-PP, Yogyakarta), Muhammad T. Soekarno and Safitriani
(PT. Pandu Biosains). We gratefully acknowledge the authors, the
originating, and submitting laboratories for their sequences and
metadata shared through GISAID platform. Our study was financially
supported by the Ministry of Health, Indonesia and Dr. Sardjito General
Hospital. The funders had no role in study design, data collection and
analysis, as well as manuscript preparation.

\textbf{Author contributions}

MSH and GUN conceived and designed the study. LSE, ES, KAV, KI, AFI,
FDTU, VCA, SSN, TPL, and FVH, performed RT-PCR, the library preparation,
and NGS. MSH, GUN, AR, ES, EWD, FNO, LA, EA, TN, and TW collected the
data. MSH, GUN, AR, and HW analyzed the data. MSH and AR drafted the
manuscript. HW and TW critically revised the manuscript for important
intellectual content. All authors have read and approved the final
version of the manuscript.

\textbf{Statements and declarations}

The authors declare no conflict of interest.

\textbf{Data availability}

Datasets, output files, and other data generated during this study are
available from the corresponding author upon request.

\textbf{ORCID}

\emph{Mohamad S. Hakim:} http://orcid.org/0000-0001-8341-461X

\emph{Tri Wibawa} : https://orcid.org/0000-0003-0477-4575

\textbf{References}

1. Lambrou AS, Shirk P, Steele MK, et al. Genomic surveillance for
SARS-CoV-2 variants: predominance of the Delta (B.1.617.2) and Omicron
(B.1.1.529) variants - United States, June 2021-January 2022. \emph{MMWR
Morb Mortal Wkly Rep} . 2022;71(6):206-211. doi:10.15585/mmwr.mm7106a4.

2. Harvey WT, Carabelli AM, Jackson B, et al. SARS-CoV-2 variants, spike
mutations and immune escape. \emph{Nat Rev Microbiol} .
2021;19(7):409-424. doi:10.1038/s41579-021-00573-0.

3. Pepin KM, Lass S, Pulliam JR, Read AF, Lloyd-Smith JO. Identifying
genetic markers of adaptation for surveillance of viral host
jumps.\emph{Nat Rev Microbiol} . 2010;8(11):802-813.
doi:10.1038/nrmicro2440.

4. WHO. Tracking SARS-CoV-2 variants. 2023. Accessed January 18, 2023 at
https://www.who.int/activities/tracking-SARS-CoV-2-variants.

5. Li Q, Wu J, Nie J, et al. The impact of mutations in SARS-CoV-2 Spike
on viral infectivity and antigenicity. \emph{Cell} .
2020;182(5):1284-1294 e1289. doi:10.1016/j.cell.2020.07.012.

6. Dhama K, Nainu F, Frediansyah A, et al. Global emerging Omicron
variant of SARS-CoV-2: Impacts, challenges and strategies. \emph{J
Infect Public Health} . 2023;16(1):4-14. doi:10.1016/j.jiph.2022.11.024.

7. Viana R, Moyo S, Amoako DG, et al. Rapid epidemic expansion of the
SARS-CoV-2 Omicron variant in southern Africa. \emph{Nature} .
2022;603(7902):679-686. doi:10.1038/s41586-022-04411-y.

8. Cao Y, Jian F, Wang J, et al. Imprinted SARS-CoV-2 humoral immunity
induces convergent Omicron RBD evolution. \emph{Nature} . 2022.
doi:10.1038/s41586-022-05644-7.

9. Cao Y, Yisimayi A, Jian F, et al. BA.2.12.1, BA.4 and BA.5 escape
antibodies elicited by Omicron infection. \emph{Nature} .
2022;608(7923):593-602. doi:10.1038/s41586-022-04980-y.

10. Wang Q, Iketani S, Li Z, et al. Alarming antibody evasion properties
of rising SARS-CoV-2 BQ and XBB subvariants. \emph{Cell} . 2022.
doi:10.1016/j.cell.2022.12.018.

11. Uraki R, Ito M, Furusawa Y, et al. Humoral immune evasion of the
omicron subvariants BQ.1.1 and XBB. \emph{Lancet Infect Dis} .
2023;23(1):30-32. doi:10.1016/S1473-3099(22)00816-7.

12. Lacek KA, Rambo-Martin BL, Batra D, et al. SARS-CoV-2 Delta-Omicron
recombinant viruses, United States. \emph{Emerg Infect Dis} .
2022;28(7):1442-1445. doi:10.3201/eid2807.220526.

13. Peacock TP, Penrice-Randal R, Hiscox JA, Barclay WS. SARS-CoV-2 one
year on: evidence for ongoing viral adaptation. \emph{J Gen Virol} .
2021;102(4):00158. doi:10.1099/jgv.0.001584.

14. Wu A, Peng Y, Huang B, et al. Genome composition and divergence of
the novel coronavirus (2019-nCoV) originating in China. \emph{Cell Host
Microbe} . 2020;27(3):325-328. doi:10.1016/j.chom.2020.02.001.

15. Peng Q, Peng R, Yuan B, et al. Structural and biochemical
characterization of the nsp12-nsp7-nsp8 core polymerase complex from
SARS-CoV-2. \emph{Cell Rep} . 2020;31(11):107774.
doi:10.1016/j.celrep.2020.107774.

16. Dwivedy A, Mariadasse R, Ahmad M, et al. Characterization of the
NiRAN domain from RNA-dependent RNA polymerase provides insights into a
potential therapeutic target against SARS-CoV-2. \emph{PLoS Comput Biol}
. 2021;17(9):e1009384. doi:10.1371/journal.pcbi.1009384.

17. Baloch Z, Ikram A, Hakim MS, Awan FM. The impact of mutations on the
pathogenic and antigenic activity of SARS-CoV-2 during the first wave of
the COVID-19 pandemic: a comprehensive immunoinformatics
analysis.\emph{Vaccines (Basel)} . 2021;9(12):1410.
doi:10.3390/vaccines9121410.

18. Eskier D, Karakulah G, Suner A, Oktay Y. RdRp mutations are
associated with SARS-CoV-2 genome evolution. \emph{PeerJ} .
2020;8(e9587. doi:10.7717/peerj.9587.

19. Mari A, Roloff T, Stange M, et al. Global genomic analysis of
SARS-CoV-2 RNA dependent RNA polymerase evolution and antiviral drug
resistance. \emph{Microorganisms} . 2021;9(5):1094.
doi:10.3390/microorganisms9051094.

20. Kochan N, Eskier D, Suner A, Karakulah G, Oktay Y. Different
selection dynamics of S and RdRp between SARS-CoV-2 genomes with and
without the dominant mutations. \emph{Infect Genet Evol} .
2021;91(104796. doi:10.1016/j.meegid.2021.104796.

21. Snijder EJ, Decroly E, Ziebuhr J. The nonstructural proteins
directing coronavirus RNA synthesis and processing. \emph{Adv Virus Res}
. 2016;96(59-126. doi:10.1016/bs.aivir.2016.08.008.

22. Moustaqil M, Ollivier E, Chiu HP, et al. SARS-CoV-2 proteases PLpro
and 3CLpro cleave IRF3 and critical modulators of inflammatory pathways
(NLRP12 and TAB1): implications for disease presentation across
species.\emph{Emerg Microbes Infect} . 2021;10(1):178-195.
doi:10.1080/22221751.2020.1870414.

23. Cheng L, Zhang X, Chen Y, et al. Dynamic landscape mapping of
humoral immunity to SARS-CoV-2 identifies non-structural protein
antibodies associated with the survival of critical COVID-19
patients.\emph{Signal Transduct Target Ther} . 2021;6(1):304.
doi:10.1038/s41392-021-00718-w.

24. Lv Z, Cano KE, Jia L, Drag M, Huang TT, Olsen SK. Targeting
SARS-CoV-2 proteases for COVID-19 antiviral development. \emph{Front
Chem} . 2021;9:819165. doi:10.3389/fchem.2021.819165.

25. Flynn JM, Samant N, Schneider-Nachum G, et al. Comprehensive fitness
landscape of SARS-CoV-2 M(pro) reveals insights into viral resistance
mechanisms. \emph{Elife} . 2022;11(e77433. doi:10.7554/eLife.77433.

26. Iketani S, Hong SJ, Sheng J, et al. Functional map of SARS-CoV-2 3CL
protease reveals tolerant and immutable sites. \emph{Cell Host Microbe}
. 2022;30(10):1354-1362.e1356. doi:10.1016/j.chom.2022.08.003.

27. Owen DR, Allerton CMN, Anderson AS, et al. An oral SARS-CoV-2 M(pro)
inhibitor clinical candidate for the treatment of
COVID-19.\emph{Science} . 2021;374(6575):1586-1593.
doi:10.1126/science.abl4784.

28. Arbel R, Wolff Sagy Y, Hoshen M, et al. Nirmatrelvir use and severe
Covid-19 outcomes during the Omicron surge. \emph{N Engl J Med} .
2022;387(9):790-798. doi:10.1056/NEJMoa2204919.

29. Yuan F, Wang L, Fang Y, Wang L. Global SNP analysis of 11,183
SARS-CoV-2 strains reveals high genetic diversity. \emph{Transbound
Emerg Dis} . 2021;68(6):3288-3304. doi:10.1111/tbed.13931.

30. Gunadi, Wibawa H, Hakim MS, et al. Molecular epidemiology of
SARS-CoV-2 isolated from COVID-19 family clusters. \emph{BMC Med
Genomics} . 2021;14(1):144. doi:10.1186/s12920-021-00990-3.

31. Gunadi, Wibawa H, Marcellus, et al. Full-length genome
characterization and phylogenetic analysis of SARS-CoV-2 virus strains
from Yogyakarta and Central Java, Indonesia. \emph{PeerJ} .
2020;8(e10575. doi:10.7717/peerj.10575.

32. Gunadi, Hakim MS, Wibawa H, et al. Is the infection of the
SARS-CoV-2 Delta variant associated with the outcomes of COVID-19
patients? \emph{Front Med (Lausanne)} . 2021;8(780611.
doi:10.3389/fmed.2021.780611.

33. Gunadi, Hakim MS, Wibawa H, et al. Association between prognostic
factors and the outcomes of patients infected with SARS-CoV-2 harboring
multiple spike protein mutations. \emph{Sci Rep} . 2021;11(1):21352.
doi:10.1038/s41598-021-00459-4.

34. Ulfah M, Helianti I. Bioinformatic analysis of the whole genome
sequences of SARS-CoV-2 from Indonesia. \emph{Iran J Microbiol} .
2021;13(2):145-155. doi:10.18502/ijm.v13i2.5973.

35. Fibriani A, Stephanie R, Alfiantie AA, et al. Analysis of SARS-CoV-2
genomes from West Java, Indonesia. \emph{Viruses} . 2021;13(10):2097.
doi:10.3390/v13102097.

36. Gunadi, Hakim MS, Wibawa H, et al. Comparative analysis of the
outcomes of COVID-19 between patients infected with SARS-CoV-2 Omicron
and Delta variants: a retrospective cohort study. \emph{medRxiv} .
2022:2022.2004.2030.22274532. doi:10.1101/2022.04.30.22274532.

37. Prasetyoputri A, Dharmayanthi AB, Iryanto SB, et al. The dynamics of
circulating SARS-CoV-2 lineages in Bogor and surrounding areas reflect
variant shifting during the first and second waves of COVID-19 in
Indonesia. \emph{PeerJ} . 2022;10(e13132. doi:10.7717/peerj.13132.

38. Massi MN, Abidin RS, Farouk AE, et al. Full-genome sequencing and
mutation analysis of SARS-CoV-2 isolated from Makassar, South Sulawesi,
Indonesia. \emph{PeerJ} . 2022;10(e13522. doi:10.7717/peerj.13522.

39. Massi MN, Sjahril R, Halik H, et al. Sequence analysis of SARS-CoV-2
Delta variant isolated from Makassar, South Sulawesi,
Indonesia.\emph{Heliyon} . 2023;9(2):e13382.
doi:10.1016/j.heliyon.2023.e13382.

40. Rantam FA, Prakoeswa CRS, Tinduh D, et al. Characterization of
SARS-CoV-2 East Java isolate, Indonesia. \emph{F1000Res} . 2021;10(480.
doi:10.12688/f1000research.53137.1.

41. Edler D, Klein J, Antonelli A, Silvestro D. raxmlGUI 2.0: A
graphical interface and toolkit for phylogenetic analyses using RAxML.
.\emph{Methods Ecol Evol} 2021;12(373-377. doi:10.1111/2041-210X.13512.

42. Tamura K, Stecher G, Kumar S. MEGA11: Molecular Evolutionary
Genetics Analysis version 11. \emph{Mol Biol Evol} .
2021;38(7):3022-3027. doi:10.1093/molbev/msab120.

43. Rodrigues CH, Pires DE, Ascher DB. DynaMut: predicting the impact of
mutations on protein conformation, flexibility and
stability.\emph{Nucleic Acids Res} . 2018;46(W1):W350-W355.
doi:10.1093/nar/gky300.

44. Alouane T, Laamarti M, Essabbar A, et al. Genomic diversity and
hotspot mutations in 30,983 SARS-CoV-2 genomes: Moving toward a
universal vaccine for the ''confined virus''? \emph{Pathogens} .
2020;9(10):829. doi:10.3390/pathogens9100829.

45. Chen B, Zhao Y, Jin Z, He D, Li H. Twice evasions of Omicron
variants explain the temporal patterns in six Asian and Oceanic
countries. \emph{BMC Infect Dis} . 2023;23(1):25.
doi:10.1186/s12879-023-07984-9.

46. Ou J, Lan W, Wu X, et al. Tracking SARS-CoV-2 Omicron diverse spike
gene mutations identifies multiple inter-variant recombination
events.\emph{Signal Transduct Target Ther} . 2022;7(1):138.
doi:10.1038/s41392-022-00992-2.

47. Yerukala Sathipati S, Shukla SK, Ho SY. Tracking the amino acid
changes of spike proteins across diverse host species of severe acute
respiratory syndrome coronavirus 2. \emph{iScience} . 2022;25(1):103560.
doi:10.1016/j.isci.2021.103560.

48. Korber B, Fischer WM, Gnanakaran S, et al. Tracking changes in
SARS-CoV-2 Spike: Evidence that D614G increases infectivity of the
COVID-19 virus. \emph{Cell} . 2020;182(4):812-827 e819.
doi:10.1016/j.cell.2020.06.043.

49. Mondeali M, Etemadi A, Barkhordari K, et al. The role of S477N
mutation in the molecular behavior of SARS-CoV-2 spike protein: An
in-silico perspective. \emph{J Cell Biochem} . 2023.
doi:10.1002/jcb.30367.

50. Anwar MZ, Lodhi MS, Khan MT, Khan MI, Sharif S. Coronavirus genomes
and unique mutations in structural and non-structural proteins in
Pakistani SARS-CoV-2 Delta variants during the fourth wave of the
pandemic. \emph{Genes (Basel)} . 2022;13(3):552.
doi:10.3390/genes13030552.

51. Showers WM, Leach SM, Kechris K, Strong M. Longitudinal analysis of
SARS-CoV-2 spike and RNA-dependent RNA polymerase protein sequences
reveals the emergence and geographic distribution of diverse
mutations.\emph{Infect Genet Evol} . 2022;97(105153.
doi:10.1016/j.meegid.2021.105153.

52. Biswas SK, Mudi SR. Spike protein D614G and RdRp P323L: the
SARS-CoV-2 mutations associated with severity of COVID-19.\emph{Genomics
Inform} . 2020;18(4):e44. doi:10.5808/GI.2020.18.4.e44.

53. Mohammad A, Al-Mulla F, Wei DQ, Abubaker J. Remdesivir MD
simulations suggest a more favourable binding to SARS-CoV-2 RNA
dependent RNA polymerase mutant P323L than wild-type.\emph{Biomolecules}
. 2021;11(7):919. doi:10.3390/biom11070919.

54. Ilmjarv S, Abdul F, Acosta-Gutierrez S, et al. Concurrent mutations
in RNA-dependent RNA polymerase and spike protein emerged as the
epidemiologically most successful SARS-CoV-2 variant. \emph{Sci Rep} .
2021;11(1):13705. doi:10.1038/s41598-021-91662-w.

55. Ullrich S, Ekanayake KB, Otting G, Nitsche C. Main protease mutants
of SARS-CoV-2 variants remain susceptible to nirmatrelvir. \emph{Bioorg
Med Chem Lett} . 2022;62(128629. doi:10.1016/j.bmcl.2022.128629.

56. Heilmann E, Costacurta F, Moghadasi SA, et al. SARS-CoV-2 3CL(pro)
mutations selected in a VSV-based system confer resistance to
nirmatrelvir, ensitrelvir, and GC376. \emph{Sci Transl Med} .
2022:eabq7360. doi:10.1126/scitranslmed.abq7360.

57. Abe K, Kabe Y, Uchiyama S, et al. Pro108Ser mutation of SARS-CoV-2
3CL(pro) reduces the enzyme activity and ameliorates the clinical
severity of COVID-19. \emph{Sci Rep} . 2022;12(1):1299.
doi:10.1038/s41598-022-05424-3.

58. Tan J, Wu Z, Hu P, Gan L, Wang Y, Zhang D. Association between
mutations in papain-like protease (PLpro) of SARS-CoV-2 with COVID-19
clinical outcomes. \emph{Pathogens} . 2022;11(9):1008.
doi:10.3390/pathogens11091008.

59. Duchene S, Featherstone L, Haritopoulou-Sinanidou M, Rambaut A,
Lemey P, Baele G. Temporal signal and the phylodynamic threshold of
SARS-CoV-2. \emph{Virus Evol} . 2020;6(2):veaa061.
doi:10.1093/ve/veaa061.

60. Tay JH, Porter AF, Wirth W, Duchene S. The emergence of SARS-CoV-2
variants of concern is driven by acceleration of the substitution
rate.\emph{Mol Biol Evol} . 2022;39(2):msac013.
doi:10.1093/molbev/msac013.

61. Benazi N, Bounab S. Comparison of the evolutionary phylodynamic of
Delta and Omicron variants of SARS-CoV-2. \emph{Research Square} . 2022.
doi:10.21203/rs.3.rs-1926171/v1.



\textbf{Figure Legends}

\textbf{Figure 1.} Phylogenetic analysis of SARS-CoV-2 from Yogyakarta
and Central Java, Indonesia based on the \emph{Spike} gene. A
phylogenetic tree was constructed from the \emph{Spike} gene of
SARS-CoV-2 using the maximum likelihood statistical method, with 1,0000
bootstrap replications and the best substitution model for the dataset
(GTR+Gamma).

\textbf{Figure 2.} Phylogenetic analysis of SARS-CoV-2 from Yogyakarta
and Central Java, Indonesia based on the \emph{nsp12}
\textbf{(A),}\emph{nsp3 (PLpro)} \textbf{(B),} and \emph{nsp5 (3CLpro)}
\textbf{(C)}genes.

\textbf{Figure 3.} High frequency (hotspot) mutations
(\textgreater{}10\%) identified in the \emph{Spike}
\textbf{(A),}\emph{nsp12} \textbf{(B),} \emph{nsp3 (PLpro)}
\textbf{(C),} and\emph{nsp5 (3CLpro)} \textbf{(D)} genes of Delta and
Omicron variants.

\textbf{Figure 4.} Visualization of the 3D structure of the
Spike\textbf{(A),} RdRp \textbf{(B),} nsp3 (PLpro) \textbf{(C),} and
nsp5 (3CLpro) \textbf{(D)} containing identified mutations from the top
and side view.

\textbf{Table 1.} Rates of nucleotide substitutions of the \emph{Spike,
nsp12, nsp3,} and \emph{nsp5} genes of SARS-CoV-2 Delta and Omicron
variants found in this study.

\textbf{Table 2.} Selective pressure analysis on the \emph{Spike, nsp12,
nsp3,} and \emph{nsp5} genes of SARS-CoV-2 Delta and Omicron variants.
Positively selected sites were shown.

\textbf{Supplementary Table 1.} The GISAID accession numbers of
SARS-CoV-2 isolates identified in this study.

\textbf{Supplementary Table 2.}The GISAID accession numbers of the
reference sequences of Alpha, Beta, Delta, Gamma, and Omicron variants.

\textbf{Supplementary Table 3.} The genome sequencing quality of all
sequences identified in this study.

\textbf{Supplementary Table 4.} The genome position of the \emph{nsp3,
nsp5, RdRp} and \emph{Spike} genes of SARS-CoV-2 based on the reference
genome from Wuhan, China (hCoV-19/Wuhan/Hu-1/2019, GenBank accession
number: NC\_045512.2)

\textbf{Supplementary Table 5.} Synonymous and nonsynonymous mutations
identified within the \emph{Spike, nsp12, nsp3,} and \emph{nsp5} genes
of SARS-CoV-2 Delta and Omicron variants.

\textbf{Supplementary Table 6.} Mutation effect analysis of the
Spike\textbf{(A)} , and RdRp, nsp3, and nsp5 proteins \textbf{(B)} of
SARS-CoV-2 Delta and Omicron variants.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-1/Figure-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2/Figure-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-3/Figure-3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-4/Figure-4}

\end{center}
\end{figure}

\textbf{Hosted file}

\verb`Table 1.docx`  available at \url{https://authorea.com/users/616696/articles/642544-sequence-analysis-of-the-spike-nsp12-rna-dependent-rna-polymerase-nsp3-plpro-and-nsp5-3clpro-genes-reveals-a-distinct-evolutionary-pattern-of-sars-cov-2-variants-circulating-in-yogyakarta-and-central-java-provinces-indonesia}

\textbf{Hosted file}

\verb`Table 2.docx`  available at \url{https://authorea.com/users/616696/articles/642544-sequence-analysis-of-the-spike-nsp12-rna-dependent-rna-polymerase-nsp3-plpro-and-nsp5-3clpro-genes-reveals-a-distinct-evolutionary-pattern-of-sars-cov-2-variants-circulating-in-yogyakarta-and-central-java-provinces-indonesia}

\selectlanguage{english}
\FloatBarrier
\end{document}

