COVID-19 is a respiratory disease triggered by an RNA virus inclined to mutations. Since December 2020, variants of COVID-19 (especially Delta and Omicron) continuously appeared with different characteristics that influenced death and transmissibility emerged around the world. To address the novel dynamics of the disease, we propose and analyze a dynamical model of two strains, namely native and mutant, transmission dynamics with mutation and imperfect vaccination. It is also assumed that the recuperated individuals from the native strain can be infected with mutant strain through the direct contact with individual or contaminated surfaces or aerosols. We compute the basic reproduction number for each strain independently and take the maximum for R 0 . We prove the nonexistence of backward bifurcation using the center manifold theory, and global stability of disease-free equilibrium when the basic reproduction number R 0 < 1, that is, vaccine is effective to eliminate the native and mutant strains even if it cannot provide full protection. Hopf bifurcation appears when the endemic equilibrium loses its stability. An intermediate mutation rate ν 1 leads to oscillations. When ν 1 increases over a threshold, the system regains its stability and exhibits an interesting dynamics called endemic bubble. An analytical expression for vaccine-induced herd immunity is derived. The epidemiological implication of the herd immunity threshold is that the disease can be effectively eradicated if the minimum herd immunity threshold is attained in the community. Furthermore, the model is parameterized using the Indian data of the cumulative number of confirmed cases and deaths of COVID-19 from March 1 to September 27 in 2021, using MCMC method. The cumulative cases and deaths can be reduced by increasing the vaccine efficacies to both native and mutant strains. We observe that by considering the vaccine efficacy to native strain as 90%, the cumulative cases and deaths would be reduced by 3.27% and 5.2%, respectively; and by considering the vaccine efficacy to mutant strain as 90%, the cumulative cases and deaths would be reduced by 0.9% and 2.5%, respectively. Our study demonstrates that the COVID-19 pandemic may be worse due to the occurrence of oscillations for certain mutation rates (i.e., outbreaks will occur repeatedly) but better due to stability at a lower infection level with a larger mutation rate. We perform sensitivity analysis using the Latin Hypercube Sampling (LHS) methodology and partial rank correlation coefficients (PRCCs) to illustrate the impact of parameters on the basic reproduction number, the number of cumulative cases, and the number of deaths, which ultimately sheds light on disease mitigation.