Objectives: The coronavirus disease 2019 (COVID-19) pandemic, caused by a dynamic virus, has had a profound global impact. Despite declining global COVID-19 cases and mortality rates, the emergence of new severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) variants remains a major concern. This study provides a comprehensive analysis of the genomic sequences of SARS-CoV-2 within the Bhutanese population during the pandemic. The primary aim was to elucidate the molecular epidemiology and evolutionary patterns of SARS-CoV-2 in Bhutan, with a particular focus on genetic variations and lineage dynamics. Methods: Whole-genome sequences of SARS-CoV-2 collected from Bhutan between May 2020 and February 2023 (n=135) were retrieved from the Global Initiative on Sharing All Influenza Database. Results: The SARS-CoV-2 variants in Bhutan were predominantly classified within the Nextstrain clade 20A (31.1%), followed by clade 21L (20%) and clade 22D (15.6%). We identified 26 Pangolin lineages with variations in their spatial and temporal distribution. Bayesian time-scaled phylogenetic analysis estimated the time to the most recent common ancestor as February 15, 2020, with a substitution rate of 0.97×10–3 substitutions per site per year. Notably, the spike glycoprotein displayed the highest mutation frequency among major viral proteins, with 116 distinct mutations, including D614G. The Bhutanese isolates also featured mutations such as E484K, K417N, and S477N in the spike protein, which have implications for altered viral properties. Conclusion This is the first study to describe the genetic diversity of SARS-CoV-2 circulating in Bhutan during the pandemic, and this data can inform public health policies and strategies for preventing future outbreaks in Bhutan.