JC69 model is the simplest model of nucleotide substitution (Liò and Goldman, 1998). It is based on two simple assumptions. The first assumption is that each residue of DNA is equally likely to change to any of the other three nucleotide bases. The second assumption is that all four bases have the same frequency. Hence, the rate of transition is equal to the rate of transversions (Pevsner, 2015). Because of its simple assumptions, the JC69 model is unlikely to be applicable in most of the data sets and works reasonably only in closely related sequences. Though it can be utilized in distantly related sequences, the correction made can sometimes be too significant to be reliable. As this model involves a single parameter, α , for both the rate of substitution for each nucleotide (3 α per unit time) and the rate of substitution in each of the three possible directions of change ( α , it is called a one-parameter model. Kimura 2 Parameter (K80) is an extension of the JC69 model (Jukes and Cantor, 1969). As in real data, transversions generally occur at lower rates than transitions; (Kimura, 1980) proposed a model, which assumes that the transitions rate is different from the transversions rate. However, like the JC69 model, Kimura also assumes that all four bases have the same frequency.

Later several more robust nucleotide substitution models, like, F81 (Felsenstein, 1981), HKY85 (Hasegawa et al., 1985) and TN93 (Tamura and Nei, 1993) models, were developed. The F81 model was developed by American scientist (Felsenstein, 1981). Unlike Jukes & Cantor or Kimura-2 parameters, the F81 model assumes that the base frequency of all bases is different. However, like Jukes & Cantor, the F81 model assumes that the base substitution occurs with equal probability (Felsenstein, 1981). Hasegawa-Kishino-Yano 85 (HKY85) model assumes unequal base frequencies as well as different substitution rates between transversions and transitions (Hasegawa et al., 1985). Tamura and Nei’s 1993 (TN93) model assumes unequal base frequencies, but all transversions are assumed to take place at an equal rate, but the transition rate between purine differs from that of pyrimidine (Tamura and Nei, 1993).

For the first time in 1986, Simon Tavaré described a general independent, finite-sites, neutral, and time-reversible model called the general time-reversible (GTR) model or the general reversible (REV) model, which assumes different substitution rates for each pair of nucleotide and unequal base frequencies (Tavaré, 1986). Additionally, the rate of variation across sites (+G) (Yang, 1994a) and/or a proportion of invariable sites (+I) (Shoemaker and Fitch, 1989) can also be included in any model. Recently several other DNA substitution models comprising of non-stationary (nucleotide composition can change over time) and non-reversible (asymmetric) matrices (Boussau and Gouy, 2006; Jayaswal et al., 2011) or even involving neighbor interactions (Lunter and Hein, 2004) were developed for inferring phylogenetic trees more accurately.

The two commonly used amino acid substitution matrices are the PAM matrices (Dayhoff et al., 1978) and the Blocks amino acid substitution matrices (BLOSUM) matrices. Margaret Dayhoff and the team aligned closely related protein sequences of seventy-one groups (Dayhoff et al., 1978). As all the sequences were closely related homologs, mutation detected in them were less likely to change the function of the protein, and hence the matrix designed was named PAM, which is an abbreviated form of Percent Accepted Mutations, where “accepted” designate the mutation favored via natural selection in the sequence (Xiong, 2006). The PAM matrices were generated on the basis of the evolutionary divergence amongst sequences of the same group. For instance, one PAM unit is described as 1% of amino acids have been modified (Xiong, 2006) and PAM60 is generated when the PAM1 matrix is multiplied by itself sixty times. Thus, PAM with a lower serial number is suitable for aligning closely related sequences, and PAM with a higher serial number is suitable for divergent sequences (Xiong, 2006). Later, Jones and the team utilized PAM matrices and developed a more advanced replacement matrix, namely the JTT model based on a large sequences dataset. After constructing a phylogenetic tree of each protein family, this method identified sequence pairs that are >85% identical and nearest-neighbors. Further, it also calculated the evolutionary distance among them. This pair of sequences were subsequently removed for avoiding recounting modifications on any given branch of a phylogeny. Likewise, this complete process was repeated for all such pairs of sequences in all protein families until the JTT matrix was finally developed (Jones et al., 1992).

BLOSUM is constructed on the basis of >2,000 conserved amino acid arrangements representing 500 groups of diverse protein sequences (Henikoff and Henikoff, 1992). Unlike PAM matrices, the BLOSUM matrices indicate the actual identity percentage amongst sequences selected for constructing the matrices (Henikoff and Henikoff, 1992). For instance, BLOSUM52 represents that sequences nominated for generating matrix share an average identity value of 52%. Hence, a higher BLOSUM number represents less divergent sequences. As PAM matrices, except PAM1, are generated from an evolutionary model and BLOSUM matrices are generated from direct observations, PAM matrices have more evolutionary meaning as compared to the BLOSUM matrices. Hence, PAM matrices are generally utilized for reconstructing phylogenetic trees. Nevertheless, due to the mathematical extrapolation technique utilized, the PAM matrices are less realistic for divergent sequences. The BLOSUM matrices are generated from local sequence alignments of conserved sequence blocks, while the PAM1 matrix is generated based on the global alignment of full-length sequences comprising both variable and conserved regions. Hence, BLOSUM matrices are more advantageous during database searching as well as finding conserved domains within proteins (Henikoff and Henikoff, 1992).

Later, (Adachi and Hasegawa, 1996), Yanga & team (Yang et al., 1998), and Adachi & team (Adachi et al., 2000) utilized the Maximum Likelihood method for developing vertebrate mitochondrial, mammalian mitochondrial (mtMAM), and chloroplast sequences (cpREV) specific amino acid replacement models, respectively. As Adachi & team, Yanga & team, and Adachi & team utilized only 20, 23 and 10 sequences, respectively, for constructing an amino acid replacement model, the accuracy of their matrices is always under question (Whelan and Goldman, 2001). Later Whelan and Goldman combined the best attributes of both Maximum Likelihood (ML) and counting methods for developing a more powerful amino acid replacement model from an extensive database of different globular protein families (Whelan and Goldman, 2001). Recently, Le and the team have also developed a more robust amino acid substitution model for metazoan mitochondrial (mtMet), vertebrate mitochondrial (mtVer) and invertebrate mitochondrial (mtInv) (Le et al., 2017). Amino acid substitution models have also been developed for Influenza virus (FLU) (Dang et al., 2010), HIV between-patient matrix HIV-Bm (HIVb) (Nickle et al., 2007), HIV within-patient matrix HIV-Wm (HIVw) (Nickle et al., 2007), arthropod mitochondrial (mtART) (Abascal et al., 2007), retrovirus (rtREV) (Dimmic et al., 2002) and general ‘Variable Time’ matrix (VT) (Müller and Vingron, 2000).

A codon is a continuous three DNA/RNA bases sequences, which encodes a specific amino acid or stop signal during protein synthesis. As there are four different nucleotides, there are only 64 possible codons. Out of these 64, only 61 code for specific amino acids, while rest three codes act as a stop codon. Since there are only 20 amino acids, more than one codon encodes one amino acid. This degeneracy property of genetic code enables us to distinguish between synonymous (do not alter encoded amino acid) and non-synonymous (alter encoded amino acid) substitution at the nucleotide level (Yang et al., 2000). Codon models are generally utilized to estimate evolutionary pressures on proteins across divergent lineages via comparing the ratio of substitution rates at non-synonymous (dN) and synonymous sites (dS) in the protein-coding regions (ω = dN/dS). Employing synonymous polymorphisms as a proxy of neutral diversity, one can estimate if non-synonymous polymorphisms are hindered or favored by natural selection. In the neutral evolving genes, the fixation rate of non-synonymous and synonymous mutation will be the same (ω = 1). During negative (purifying) selection, the non-synonymous mutation is not favored by natural selection and thus is eliminated, causing the fixation rate of non-synonymous mutation to be lower than the synonymous rate (ω< 1). During positive (adaptive) selection, the non-synonymous mutation is favored via Darwinian selection, thereby causing the fixation rate of non-synonymous mutation to be higher than the synonymous rate (ω> 1) (Gupta and Vadde, 2020).

A study reported that ancient proteins are under strong purifying selection, while newly developed proteins are under positive selection (Vishnoi et al., 2010). As newly developed young genes perform either highly specialized (if generated de novo or via horizontal transfer) or redundancy (if generated via duplication) functions, they are more at risk of either losing their function or gaining novel functions in succeeding lineages (Domazet-Loso and Tautz, 2003; Daubin and Ochman, 2004; Wolf et al., 2009; Vishnoi et al., 2010). Though initially, young genes experience a large number of adaptive mutations, the substitution of some of the mutations will slowly optimize the function of the gene in due course of time and hence the supply of new adaptive mutations will also reduce; hence, ω value of a young gene will decline over time (Vishnoi et al., 2010; Moutinho et al., 2022). On the contrary, functions of old genes, like diabetic genes, are highly optimized and they are likely to have already exhausted all beneficial mutations by recent times and; thus, they are expected to evolve under negative selection and fix only neutral and/or nearly neutral mutations (Vishnoi et al., 2010).

Though ω was originally designed for detecting selective pressure acting on a protein across divergent lineage, ω can also be utilized for detecting selective pressure acting on a protein in a single population (Kryazhimskiy and Plotkin, 2008). However, selective pressure estimated via ω on sequences sampled from a single population differs from that of the divergent lineages. For instance, though ω< 1 is a clear signature of negative selection across divergent lineages, weak negative or strong positive selection between population samples is also expected to produce ω< 1 (Roumagnac et al., 2006; Holt et al., 2008; Kryazhimskiy and Plotkin, 2008). Strong positive selection in a population will generate speedy sweeps at selected sites (but not at neutral sites, which are presumed to be independent). Thus, two individuals from the same population under strong positive selection are likely to contain identical alleles at each selected site, generating ω< 1 (Kryazhimskiy and Plotkin, 2008).

Mechanistic codon models detect selective pressure on the coding region of a gene utilizing a finite set of parameters, for instance, synonymous/non-synonymous rate ratio, transversion/transition rate ratio, and codon frequencies at equilibrium. As the mechanistic codon model utilizes a finite set of parameters, it is also known as the parametric codon substitution model. Mechanistic codon models focus mainly on silent-transversion, silent-transition, replacement-transversion rates, and replacement-transitions amongst sense codons and codon frequencies. Considering all parameters in a single codon model will be computationally more demanding. Thus to avoid this problem, several mechanistic codon models have been developed to date. Each mechanistic codon substitution models have distinctive parameters that differentiate the substitution rate at the nucleotide level and selective pressure at the protein level. Thus, each mechanistic model has the capacity to estimate selective forces acting on any protein in their unique way (Whelan et al., 2001; Delport et al., 2009). If selective pressure at the protein level is not considered, codon models will be equivalent to nucleotide substitution models. If the substitution rate at nucleotide is not considered, the codon model will be equivalent to amino acid substitution models (Miyazawa, 2011a). Several studies utilizing a large set of protein-coding sequences reported that codon substitution models are statistically more powerful than nucleotide and amino acid models (Seo and Kishino, 2009; Miyazawa, 2011b). However, the codon model having a larger substitution rate was reported to be equivalent to the amino acid substitution model (Seo and Kishino, 2008).

The first two mechanistic codon substitution models (Goldman and Yang, 1994; Muse and Gaut, 1994) were capable of estimating only the global ω of the coding region of a gene. These two models considered transition/transversion ratio and codon frequencies for estimating ω (Goldman and Yang, 1994; Muse and Gaut, 1994). Besides, Goldman and Yang (Goldman and Yang, 1994) also considered replacement probabilities amongst amino acids on the basis of the Grantham physicochemical distance matrix (Grantham, 1974). Later Nielsen and Yang (Nielsen and Yang, 1998) & Yang and the team (Yang et al., 2000) developed more robust mechanistic Bayesian models individually. It is pertinent to note that in the models developed by Goldman and Yang (Goldman and Yang, 1994) and Nielsen and Yang (Nielsen and Yang, 1998), the rate of substitution is proportional to the frequency of the target codon (which is not very mechanistic), and later many models employes these models to “explain” the stationary distribution in codons, whereas in the model developed by Muse and Gaut (Muse and Gaut, 1994), it is proportional to the target nucleotide, which is much more mechanistic considering the mutation process, and based on which later different type of mutation-selection (MutSel) model was developed [described below].

In 1994, Yang (Yang, 1994b) developed two approximation approaches for Maximum Likelihood phylogenetic estimation, which allow for varying substitution rates across nucleotide sites. The first, known as the “discrete gamma model,” approximates the gamma distribution by using many rate categories with equal probability for each category. Each category’s mean is employed to depict all of the rates within that category. This method’s performance has been shown to be rather acceptable, with four such categories seeming to be adequate to achieve both an optimum or near-optimal fit by the model to the data, as well as an acceptable approximation to the continuous distribution. The second strategy, dubbed the “fixed-rates model,” divides sites into multiple groups based on the rates anticipated by the star tree. When evaluating alternative tree topologies, sites in various classes are considered to evolve at these constant rates. Analyses of the data sets indicated that this approach might yield good results; however, it seems to have certain aspects with a least-squares pairwise comparison (Yang, 1994b). These models, however, overlook the fact that substitution rates of each amino acid differ distinctly. For instance, as only one transversion is required to convert the phenylalanine codon (UUU) into a leucine codon (UUG) as well as the tryptophan codon (UGG) into a leucine codon (UUG), they consider their substitution rate to be same (Doron-Faigenboim and Pupko, 2007). But in reality, the probability of occurring a former event is approximately 5 times higher than a later event (Doron-Faigenboim and Pupko, 2007).

Considering this lacuna, for the first time, in 2004, Whelan and Goldman developed a complete parametric model that considers several instantaneous substitutions (Whelan and Goldman, 2004). This model estimated substitution rate matrice for single-, double- and triple nucleotide mutation individually utilizing transition to transversion ratio and equilibrated frequency of mutated nucleotides. Later, these three matrices were joined together to estimate the general codon rate matrix. This method is reported to estimate the likelihood of parameters more accurately in comparison to other mechanistic models (Whelan and Goldman, 2004). Double and triple nucleotide substitutions are reported to occur through the mechanistic process, like during repairing DNA break (Sakofsky et al., 2014) or error-prone polymerase activity (Harris and Nielsen, 2014). Although double and triple substitutions rates are predicted to be two to three orders of magnitude lower than single substitutions (Smith et al., 2003; Whelan and Goldman, 2004; Tamuri et al., 2012), the model which included double and triple substitutions were reported to fit better in real data. Later several different models were developed by Doron-Faigenboim and Pupko (Doron-Faigenboim and Pupko, 2007), Kosiol, Holmes and Goldman (Kosiol et al., 2007), De Maio & team (De Maio et al., 2013), Miyazawa (Miyazawa, 2011b), Zoller and Schneider (Zoller and Schneider, 2012), Zaheri, Dib and Salamin (Zaheri et al., 2014), Venkat & team (Venkat et al., 2018) and Jones & team (Jones et al., 2018), which included double and triple substitution between codon.

Later, the model developed via Goldman and Yang (Goldman and Yang, 1994) was modified to include various nucleotide models (Pond et al., 2005; Pond and Frost, 2005; Arenas and Posada, 2014), estimate ω variation across sites (Yang, 2007) and branches (Yang, 2007; Dutheil et al., 2012). In the model developed by Pond and Muse, they consider the possibility of site-to-site variation in synonymous and non-synonymous substitution rates in protein-coding DNA sequences and observed that within-gene variability in synonymous substitution rates is common (Pond and Muse, 2005). Another model developed by Mayrose et al. (2007) that uses two hidden Markov models and function on the spatial dimension. First and the second model depicts the dependency between adjacent non-synonymous and synonymous rates rates, respectively. They demonstrate that taking into consideration synonymous rate variability and dependence substantially improves the accuracy of ω estimate, in particular for positively selected sites. In some models, codons were partitioned based on the physicochemical properties of the encoded amino acids (e.g., polarity or charge) (Sainudiin et al., 2005; Wong et al., 2006), codon bias (Yang and Nielsen, 2008) or the effects of GC contents (Misawa, 2011). Models in which codons were partitioned based on the physicochemical properties of the encoded amino acids explicitly parameterized physiochemical constraints due to non-synonymous substitution. The model developed by Goldman & Yang (Goldman and Yang, 1994) and Yang, Nielsen & Hasegawa (Yang et al., 1998) applied mathematical functions for modeling association amongst physiochemical properties and ω parameter. Yang (Yang, 2000) permitted the effect of the physicochemical property to fluctuate among sites. Sainudiin & team (Sainudiin et al., 2005) and Wong & team (Wong et al., 2006) developed two separate models that at first divided non-synonymous substitutions into small groups in accordance with the pre-defined physiochemical property. As the main objective of these two models is to examine the impact of certain physicochemical properties of amino acids on the structure and function of a protein, their parameterization is focused on comparing the property-modifying substitutions rate with the property-conserving substitutions rate. Conant and Stadler (Conant and Stadler, 2009) estimated multiple amino acid properties via modeling exchangeability amongst non-synonymous codons as a linear combination of five pre-specified measures of physiochemical properties. This model enabled us to investigate the association between selection pressure and physicochemical properties while avoiding over parameterization of the codon model.

In 2008, Yang and Nielsen (Yang and Nielsen, 2008) developed the FMutSel model in which the amino acid frequencies are determined by the functional requirements of the protein (Rodrigue et al., 2008; Beaulieu et al., 2019). In the FMutSel model, each codon was allocated a fitness parameter. Dissimilarities in fitness parameters amongst two codons are utilized for specifying substitution rates in the Markov matrix via altering the rates specified by the standard mutation models (Yang and Nielsen, 2008). The FMutSel/FMutSel0 model combination has only been implemented in PAML4 with the M0 and M3 models so far. Model M0 implies that ω across all branches and sites is constant, whereas Model M3 allows ω to vary between sites (Du et al., 2014). Likewise, in 2010, Rodrigue and the team developed a complex extension of this model in which site-specific amino acid propensity scores are utilized for estimating scaled selection coefficients, which in turn was utilized for identifying substitution rates (Rodrigue et al., 2010). In 2013, De Maio and team (De Maio et al., 2013) reported that when some models were employed to compute ω heterogeneity on data, where both multiple-nonsynonymous rates and double & triple codon modification occur, they yielded high false-positive rates. Recently Venkat & team (Venkat et al., 2018) reported that when branch-site codon models are employed in branch-specific tests to detect positive selection, double modification may cause high false-positive rates. To avoid this problem, recently, Dunn and the team developed a statistically more powerful general-purpose parametric modeling framework for codons (Dunn et al., 2019). By including information about all possible instantaneous codon substitutions, along with instantaneous double and triple nucleotide substitution and multiple non-synonymous rates, both accuracy, as well as statistical power was highly improved (Dunn et al., 2019).

Though empirical codon models are highly useful in understanding protein evolution as well as in phylogenetic applications, only a few models have been developed to date (Kosiol et al., 2007). Substitution rates amongst codons were empirically determined to utilize a large set of protein-coding sequences (Schneider et al., 2005; Kosiol et al., 2007). Unlike mechanistic models, empirical codon substitution cannot distinguish between the substitution rate at the nucleotide level and selective pressure at the protein level (Miyazawa, 2011a). Thus, there is no parameter except codon frequencies for tailoring of each protein family (Miyazawa, 2011a). Delport and team reported that empirical substitution matrices represent average propensities of substitutions across several protein families via sacrificing gene-level resolution (Delport et al., 2010).

For the first time in 1990, Schöniger and the team constructed counted codon-codon substitutions matrix on the basis of ∼800 pairwise alignments of 41 actin genes (Schöniger et al., 1990). However, due to a lack of a better electronic facility, this matrix lost its fame in a short interval of time (Cannarozzi and Schneider, 2012). Additionally, as this matrix was developed based on a small number of sequences, it was less reliable. Later in 2005, Schneider and the team (Schneider et al., 2005) developed another empirical codon model utilizing a somewhat similar approach employed by Gonnet and the team (Gonnet et al., 1992) for constructing an amino acid substitution matrix. Since the development of the first amino acid substitution matrix (Dayhoff et al., 1978), it was for the first time Goonet & team (Gonnet et al., 1992) and Jones & team (Jones et al., 1992) developed two individual models based on the sufficiently large amount of sequences and thus are more reliable. In 2007, Kosiol and the team (Kosiol et al., 2007) developed an empirical codon substitution matrix utilizing an extensive database of protein-coding DNA sequences. They reported that the accuracy of the model gets significantly improved by considering instantaneous double and triple substitution. Additionally, the amino acid encoded by each codon, associations amongst codons, and physicochemical properties of amino acids is key factors for driving the process of codon evolution (Kosiol et al., 2007). Empirical codon substitution matrix is reported to outperform mechanistic codon substitution matrix when utilized in likelihood-based phylogenetic analysis (Kosiol et al., 2007). Empirical codon models can also be utilized to detect different lineages sampled in a single phylogenomic dataset (De Maio et al., 2013) rather than depending on a general sequence database for instance Pandit database (Kosiol et al., 2007). In 2014, Bloom developed another novel model that depicts the experimental determination of a parameter-free evolutionary model via deep sequencing, mutagenesis, and functional selection (Bloom, 2014). Employing this, Bloom build a model of influenza nucleoprotein evolution that represents the gene phylogeny in a far better way as compared to earlier existing models with nearly hundreds of free parameters. He also emphasized that the data provided by these types of high-throughput experiments can significantly improve the accuracy of both phylogenetic as well as genetic studies.

The semi-empirical codon substitution matrix is often called a mixed empirical codon substitution matrix because it combines empirical substitution rate with mechanistic parameters of codon evolution (Dunn et al., 2019). For the first time, Doron-Faigenboim & Pupko combined existing empirical amino acid substitution matrices with mechanistic parameters (Doron-Faigenboim and Pupko, 2007). They assumed that the substitution rate between non-synonymous substitution amongst codons was equal to the pre-estimated amino-substitution rate, which was obtained by utilizing 189 parameters and a huge amount of amino acid sequences (Doron-Faigenboim and Pupko, 2007). Later Kosiol and team, utilizing 1830 codon substitution parameters and large datasets, developed the first fully empirical codon model and then appended those models with mechanistic parameters for codon evolution (Kosiol et al., 2007). Subsequently, De Maio and team (De Maio et al., 2013), developed another model with almost the same accuracy but was less complex than the model developed by Doron-Faigenboim & Pupko (Doron-Faigenboim and Pupko, 2007). The empirical matrices in those studies denote a wide range of amino acid change propensity (Dunn et al., 2019). Later, Zoller & Schneider (Zoller and Schneider, 2012) and Miyazawa (Miyazawa, 2011a) developed different methods for tailoring information contained in an empirical substitution matrix to a specific dataset, and the benefit of these two approaches is that they can easily distinguish between substitution rate at the nucleotide level and selective pressure at the protein level.

Recent advancements in high throughput sequencing technologies have enabled the generation of a huge amount of sequence data (Gupta et al., 2017; Gupta and Vadde, 2019b; 2020). This enormous amount of sequence data provides an opportunity to detect a direct association between selective pressure and the function of any protein in a more comprehensive way (Anisimova and Kosiol, 2009). Codon models are commonly utilized for identifying candidate genes and their variants under positive selection (Ouyang and Liang, 2007; Parto and Lartillot, 2018; Dunn et al., 2019). Initially, codon models presume that synonymous and non-synonymous substitution rates among sites as well as throughout the phylogenetic history, are constant (Anisimova and Kosiol, 2009). Though the majority of proteins are evolved under purifying selection, the positive selection may affect a few lineages. During this adaptive evolution, only a few protein sites have the capability to increase protein fitness during amino acid substitution (Pupko and Galtier, 2002). Hence, these codon model approaches presuming constant selective pressure over time as well as across sites lack the power to detect genes evolving under positive selection (Anisimova and Kosiol, 2009). Subsequently, several situations of variation in the selective pressure was included with the model developed by Muse & Gaut (Muse and Gaut, 1994) and Goldman & Yang (Goldman and Yang, 1994). These models were later utilised extensively to detect positive selection by likelihood ratio test comparing two nested models. One model (null hypothesis) do not permit positive selection while other model (alternative hypothesis) permit positive positive selection. Positive selection is identified when model permitting sites/lineages under positive selection fits data significantly better than the model restricting the site/lineages under positive selection. Nevertheless, if few parameters become invaluable or because of boundary problems, the asymptomatic null distribution may differ from the standard (Anisimova and Kosiol, 2009).

Additionally, the codon substitution model can be utilized for detecting site-specific positive selection in proteins. Later, this information can be used for testing the biological hypothesis through laboratory experiments (Anisimova and Kosiol, 2009). For instance, in 2005, Swayer and the team reported that a small portion of TRIM5α, an immune defense protein, was recognized to be under positive selection. Later, through functional analysis, they confirmed the significance of the peptide segment in species-specific viral inhibition (Sawyer et al., 2005). The conditional selection model developed by Chen and the team may be utilized particularly for detecting interaction amongst sites during drug resistance (Chen and Lee, 2006). Considering this earlier, we also employed a phylogenetic approach implemented in the PAML’s CODEML modeling tool to identify the kind of selection operating on T2D genes in the Drosophila genus. The data showed that the gene sequences encoding T2D are evolving under purifying selection. However, few membrane protein sites, including those encoded by CG8051, ZnT35C, and kar, are substantially evolving under positive selection. This may be due to adaptive evolution in response to changes in the niche, food, or other environmental conditions (Gupta and Vadde, 2020). Thus, the identification of selective pressure via codon substitution models may provide detailed insight into disease progression, pathogenic drug resistance, and epidemic dynamics (Anisimova and Kosiol, 2009).

Gene expression is modulated through transcription (DNA to mRNA) and translation (mRNA to protein) mechanisms (Zhou et al., 2016). Promoter strength & RNA stability are mainly responsible for mRNA concentrations and transcript levels & protein stability is responsible for protein concentrations in any cell (Ikemura, 1985; Sharp et al., 1986). During translation, the information is transmitted as codons. This genetic code is degenerate in nature, i.e., except for tryptophan and methionine, more than one codon (synonymous codons) can encode a single amino acid (Wang et al., 2018). In coding sequences of many organisms, these synonymous codons are utilized at unequal frequencies (Chakraborty et al., 2017). This phenomenon is called codon usage bias. Preferred codons are more frequently utilized in highly expressed genes (Zhou et al., 2016). The degree of codon usage bias differs amongst genes & species and is mainly affected via neutral selection, directional mutation, tRNA abundance (Olejniczak and Uhlenbeck, 2006), selection for efficient translation initiation (Zalucki et al., 2007), gene length (Sun et al., 2009), an expression level (Hiraoka et al., 2009), DNA replication initiation site (Huang et al., 2009), etc. Codon usage bias can also be utilized in detecting phylogenetic trees amongst species (Wu et al., 2007; Zhao et al., 2016). In 2016, SENCA (site evolution of nucleotides, codons, and amino acids), a codon substitution model, was developed that distinctly describes (a) preferences amongst synonymous codons, (b) amino acids, and (c) nucleotide processes that apply on all sequence sites such as the mutational bias (Pouyet et al., 2016). This model assumes that the vast majority of synonymous substitutions are not neutral and can predict more accurate estimates of selection in comparison to more traditional codon sequence models (Pouyet et al., 2016).

Codon substitution models are also utilized for reconstructing ancestral coding sequences through parsimony and Maximum Likelihood approaches (Anisimova and Kosiol, 2009). These ancestral sequences can further be utilized to detect alterations that have been experienced in every branch of phylogeny and at each individual site of the gene sequence. Several studies have utilized ancestral state information to understand protein evolution and episodic or lineage-specific base composition (Long and Langley, 1993; Akashi, 1996; Eanes et al., 1996; Fitch et al., 1997; Takano-Shimizu, 2001). For instance, the evolution of steroid receptors (Thornton et al., 2003) and ancestral archosaur visual pigment rhodopsin (Chang et al., 2002). Ancestral sequence reconstruction is also employed in studying HIV evolution (Gaschen et al., 2002), protein engineering (Cole and Gaucher, 2011), and understanding variation in DNA turnover because of indels and substitutions amongst eutherian mammalian lineages (Blanchette et al., 2004). Additionally, numerous population genetic tests depend on this ancestral reconstruction to understand the impact of natural selection on the functional classes of mutations or genetic regions (Akashi, 1995; Templeton, 1996; Akashi, 1999; Suzuki and Gojobori, 1999) and also identify coevolving nucleotides/amino acids (Osada and Akashi, 2012; Liao et al., 2013).

Codon models reconstruct phylogenetic trees by considering genetic code and the rate of non-synonymous & synonymous base substitutions. In almost every protein-coding gene, the incidence of non-synonymous substitution is less and is mainly involved in early divergence. Synonymous substitutions are higher and are responsible for recent divergence. By considering this information, the codon models may be utilized in reconstructing phylogenetic trees more accurately (Ren et al., 2005). Earlier studies have reported that though nucleotide substitution models are modified to accommodating differences in the evolutionary dynamics at three codon positions (Yang, 1996), the accuracy of this model is lower as compared to codon models. Nevertheless, due to the lack of efficient codon-based tree search methods, tree inference from coding sequence data is generally performed under DNA and AA (amino acid) models. Because of the 61 × 61 matrix, tree generation utilizing codon models is computationally more demanding. To date, no efficient methods have been developed for phylogeny reconstruction utilizing the codon substitution model on a large dataset. For the small dataset, phylogeny can be constructed using CODEML from the PAML package (Yang, 2007). However, Yang has implemented a heuristic algorithm in PAML, which is not the most efficient approach. One possible way to reconstruct an efficient phylogenetic tree is by initially generating numerous phylogenetic trees utilizing both DNA and amino acid substitution models. Later, these trees can be utilized for constructing more accurate trees under efficient Maximum Likelihood (ML) heuristics under codon models (Anisimova and Kosiol, 2009).

Another significant approach in the reconstruction of the phylogenetic tree is by implementing codon models with a Bayesian framework and sampling topological space with an efficient Markov chain Monte Carlo (Anisimova and Kosiol, 2009). By using the Bayesian framework, we can either get similar or even better tree topology in comparison to ML approaches. The main reason for achieving better tree topology via the Bayesian framework is that ML approaches search for a single best tree while the Bayesian framework scan cluster of best trees. The benefits of the Bayesian framework can also be explained via the matter of probability. The best tree generated via ML may ∼90% probability of demonstrating the real information. On the contrary, the Bayesian framework generates hundreds/thousands of near-optimal/optimal having ∼90% probability representing the real information. Hence, the phylogenetic tree generated via the Bayesian framework is more realistic than the phylogenetic tree generated via ML methods. However, in this Bayesian framework, the rate of substitution will differ for each three codon sites because it considers different data partitions. Thus, codon model usage may serve as an important asset while comparing several candidate trees inferred under either DNA or amino acid models (Anisimova and Kosiol, 2009).