Introduction

Modelling the nucleotide substitution processes that underly the diversification of virus genome sequences lies at the heart of many viral evolutionary analyses. The most widely used nucleotide substitution models belong to the general time reversible (GTR) family (Tavare, 1986) and assume that the Markov process of evolution is time-reversible (Hoff, Orf, Riehm, Darriba, & Stamatakis, 2016), (Lio & Goldman, 1998), (Tavare, 1986).

The GTR model is defined by its instantaneous rate matrix Q_ij equation (1), where Q_ij defines the instantaneous rate of change from nucleotide i ∈ {A, C, G, T} to nucleotide j; subject to the detailed balance condition: q_jiπ_i = q_ijπ_j, with rates q and equilibrium frequencies π (Squartini & Arndt, 2008) (Posada, 2003). The instantaneous rate matrix of the GTR model includes six rate parameters (a, b, c, d, e, and f). Because only products of substitution rates and evolutionary times can be estimated, one of the rate parameters is set to 1 (e.g. b), or the entire matrix is normalized to yield one expected substitution per unit time.

The rate matrix in equation (1) is symmetrical in that, for example, the relative rate at which A changes to G is the same as the relative rate at which G changes to A.

Time reversible nucleotide substitution models such as GTR form the basis of almost all nucleotide sequence-focused evolutionary analyses (including those involving eukaryotes, prokaryotes, and viruses) (Lefort, Longueville, & Gascuel, 2017), (Posada & Crandall, 2001), (Minin, Abdo, Joyce, & Sullivan, 2003).

The reliability of a phylogenetic tree constructed using a particular nucleotide sequence dataset should be maximized when the evolutionary models used to construct the tree accurately reflect the important aspects of the evolutionary process (Buckley & Cunningham, 2002), (Ripplinger & Sullivan, 2008), (Hoff, Orf, Riehm, Darriba, & Stamatakis, 2016). The suitability of different models for describing the evolution of DNA or RNA sequences is, therefore, expected to depend to some degree on the biological and environmental contexts of the sequences being analyzed.

Mutations in viral genomes arise due to diverse biotic (such as replication enzyme infidelities, RNA/DNA editing enzymes) and abiotic (such as ionizing radiation, inorganic oxidizers and chemical mutagens) factors (Sanjuan & Domingo-Calap, 2016). Mutagenic chemical reactions or types of radiation that, for example, cause G to A or C to U mutations in DNA or RNA, are not the same as those that cause A to G or U to C mutations (Cheng, Cahill, Kasai, Nishimura, & Loeb, 1992), (Nguyen, et al., 1992), (Chelico, Pham, Calabrese, & Goodman, 2006), (Sharma, Patnaik, Taggart, & Baysal, 2016). It should not be expected, therefore, that the relative rates of G to A substitution will equal the relative rates of A to G substitution. Instead, in evolving double-stranded (ds) DNA and dsRNA molecules where both strands of the genome are in existence for similar amounts of time, both G to A and C to T substitutions should occur at relatively similar rates. Therefore, for nucleotide sequence datasets derived from any organisms with dsDNA or dsRNA genomes, a non-reversible nucleotide substitution model with a different relative substitution rate category for each of the six possible pairs of complementary nucleotide substitutions (e.g. NREV6 in equation (2)), with q_AC = q_TG, q_AG = q_TC, q_AT = q_TA, q_CG = q_GC, q_CT = q_GA, q_GT = q_CA), might plausibly provide a better description of mutational processes than GTR (Baele, Van de Peer, & Vansteelandt, 2010), (Wickner, 1993).

In the case of single-stranded (ss) RNA viruses, ssDNA viruses, retroviruses, and dsRNA/dsDNA viruses where the two complementary genome strands do not exist for equal amounts of time (Yu, et al., 2004), a model where all twelve different substitutions occur at different rates might be best. Specifically, with ssRNA viruses, ssDNA viruses, and retroviruses, only one of the genome strands (called the virion strand) is packaged into viral particles for transmission and, in many dsRNA viruses, the genome strand that is translated into proteins (called the + strand) exists for longer during the life cycle than does the complementary (or −) strand (Bruslind, 2020), (Onwubiko, et al., 2020). In all these viruses, some degree of strand-specific substitution bias is expected to occur (Van Der Walt, Martin, Varsani, Polston, & Rybicki, 2008), (Polak & Arndt, 2008) such that NREV6 might be anticipated to provide a poorer description of mutational processes than a model such as NREV12 (equation (3)) where each of the twelve different types of substitution has a separate rate (Baele, Van de Peer, & Vansteelandt, 2010).

Because non-reversible models consider the directionality of evolution, they could, in some cases, be used to identify root nodes of phylogenetic trees (Yap & Speed, 2005), (Boussau & Gouy, 2006). It is, however, unclear whether non-reversible models might, in certain situations at least, perform better than reversible models in the context of phylogenetic inference. Although it is possible to use non-reversible nucleotide substitution models such as NREV6 and NREV12 during maximum likelihood-based phylogenetic inference with computer programs such as IQ-TREE (Nguyen, Schmidt, Von Haeseler, & Minh, 2015), these models are not routinely used for phylogenetic inference. This is in part because non-reversible models render several commonly used algorithmic techniques for efficient likelihood computation inapplicable making inference slower. It is also in part because it remains undetermined whether, under conditions where strandspecific substitution biases are evident, non-reversible models consistently yield substantially more accurate phylogenetic trees than reversible models.

Here we present evidence that strand-specific nucleotide substitution biases are common within virus genomic sequence datasets such that NREV12 generally provides a significantly better fit than both GTR and NREV6 for such datasets. We then use simulations to demonstrate that whereas strand-specific nucleotide substitution biases reduce the accuracy of phylogenetic inference under both GTR or NREV12, when these biases become extreme use of NREV12 can yield significantly more accurate phylogenetic trees than GTR.

Results and Discussion

Non-reversible nucleotide substitution models generally provide a better fit than reversible models to virus sequence datasets

We tested for evidence of non-reversibility of the nucleotide substitution process in 141 virus sequence datasets (33 ssDNA virus datasets, 30 dsDNA virus datasets, 31 dsRNA virus datasets, and 47 ssRNA virus datasets, all consisting of either full genome sequences (for unsegmented viruses), or complete genome component sequences (for viruses with segmented genomes). Specifically, for each dataset, we compared the goodness-of-fit of the GTR+G, NREV6+G, and NREV12+G models (where G represents gamma-distributed nucleotide substitution rates among sites (Yang, 1994)).

Given that dsDNA viruses such as adenoviruses, papillomaviruses and herpesviruses have both their DNA strands in existence for similar amounts of time before DNA-dependant-DNA polymerase enzymes copy both their + and − DNA strands during replication (Hanson, 2009), we had anticipated that the best fitting substitution model for sequence datasets of these viruses would be NREV6. Using weighted AlC-c scores to reveal trends of model support (Figure 1), it is surprising that NREV12 was overall the best-supported model (illustrated by the redder hues around the top corner of the dsDNA plot in Figure 2. Out of the 30 dsDNA datasets considered, we found that NREV6 provided the best fit to five datasets (HPV18, HPV45, HPV16, BPV, and SV40) and GTR provided the best fit to five (Alphapapilloma virus 6, JC polyomavirus, DPV, RTBV, and DBAV). NREV12 was the best fitting model for the remaining 20 datasets(Table 1). Further, LRTs revealed strong overall support for NREV12, with this model providing a significantly better fit (p < 0.05) than NREV6 for 25/30 of the dsDNA datasets and a significantly better fit than GTR for 24/30 of the datasets.

Figure 1

Figure 2

Table 1

As NREV6 was not the best fitting model for most of the dsDNA virus datasets we infer that, in most dsDNA virus species, strand-specific substitution biases are not ignorable. Further, the datasets where NREV6 was not the best fit are from species in families containing other species where NREV6 was the best fit, indicating that such strand-specific substitution biases are unlikely to be a consequence of some broadly conserved feature of viral life cycles in these families (such as ssDNA replicative intermediates). It is instead plausible that these differences may relate to:

1. differences in replication fidelity and/or proofreading efficiency on the leading and trailing DNA strands in some dsDNA virus species (Grigoriev, 1999). These differences are common in eukaryotes (Youri, Newlon, & KunkelThomas, 2002) (Furusawa, 2012) and prokaryotes (Fijalkowska, Jonczyk, Tkaczyk, Bialoskorska, & Schaaper, 1998) and, considering that the replication processes of dsDNA viruses analyzed here mirror those of their eukaryote hosts, it is perhaps unsurprising that most of these viruses also display some evidence of strand-specific substitution biases.

2. extra exposure to DNA damage of displaced template strands during unidirectional rolling circle replication in some papillomavirus species such as HPV16 could be a contributor to strand-specific nucleotide substitution biases (Kusumoto-Matsuo, Kanda, & Kukimoto, 2011).

3. extra time spent by non-coding strands in single-stranded dissociated states during RNA transcription in some papillomavirus and polyomavirus species (Fernandes & de Medeiros Fernandes, 2012): during transcription processes, the dissociated non-coding strand is transiently more exposed to damage than the coding strand (Wei, et al., 2010) which might also contribute to strand-specific substitution biases.

Similarly, and equally surprising, we found that NREV12 was overall the best supported model for dsRNA viruses (illustrated by the redder hues around the top corner of the dsRNA plot in Figure 1).

NREV6 fit only two of the 31 dsRNA datasets better than both NREV12 and GTR (Human rotavirus A set H and Fiji virus). NREV12 was found to be the best fitting model for 21/31 datasets and GTR was the best fitting of 8/31 (Table 2). In all three Birnaviridae family datasets (which contains virus species with two genome segments) and in 17/22 of Reoviridae family datasets (which contain virus species with 10-12 genome segments) the NREV12 model provided the best fit. Based on the LRTs, strong overall support for NREV12 was found, with this model providing a significantly better fit (p < 0.05) in 27/31 dsRNA virus datasets relative to NREV6, and 23/31 datasets relative to GTR.

Table 2

We anticipated that NREV12 might fit many of these dsRNA datasets better than NREV6 simply because, during their infection cycles, the coding +strand of dsRNA viruses (the one from which protein translation occurs) tends to exist for longer periods within an infected cell than the noncoding -strand. Specifically, there are two main steps during double-stranded RNA virus replication (Wickner, 1993). Firstly, synthesis of the viral +strands from a dsRNA template occurs in the cytoplasm within viral particles. These +strands exist within the cell for prolonged periods in the absence of complementary -strands and are used as templates for translation of viral proteins. In the second step the +strands remaining after translation act as templates for -strand synthesis resulting in the formation of new dsRNA molecules. The +strands of dsRNA viruses are therefore likely more impacted by mutational processes, which in turn could explain the pervasive strand specific substitution biases seen in this group of viruses.

For the ssRNA and ssDNA viruses where one genome strand exists during the virus life cycle for far longer periods of time than the other such that complementary substitutions would not be expected to occur at similar rates, we anticipated that NREV12 should provide a better fit than both NREV6 and GTR. Indeed, for ssRNA viruses, NREV12 was a better fit than NREV6 and GTR for 40/47 of the ssRNA datasets and 24/33 of the ssDNA virus datasets (Figure 1). Of the nine ssDNA virus datasets where NREV12 was not the best fitting model, GTR fit 7/9 better and NREV6 fit 2/9 better (Table 3). Of the seven ssRNA datasets where NREV12 was not the best fitting model, GTR fit 6/7 better and NREV6 fit 1/7 better.

Table 3

Based on the LRTs, strong overall support for NREV12 was found with this model providing a significantly better fit (p < 0.05) than NREV6 for 45/47 of the ssRNA virus datasets (Table 4) and 31/33 of the ssDNA virus datasets (Error! Reference source not found.). Similarly, based on LRTs, NREV12 provided a significantly better fit than GTR for 27/33 of the ssDNA virus datasets (Error! Reference source not found.) and 40/47 of the ssRNA virus datasets (Table 4).

Table 4

We found that DNR estimates alone did not cleanly differentiate between datasets for which NREV12 was or was not best supported (Tables 1, 2, 3 and 4). For the 107 nucleotide sequence datasets with a model preference of NREV12, ten had estimated DNRs that were greater than 0.5, 13 had DNRs between 0.25 and 0.5 and 84 had DNRs between 0.0225 and 0.25. For the ten nucleotide sequence datasets with a model preference of NREV6, one had an estimated DNR greater than 0.5, four had estimated DNRs between 0.25 and 0.5 and five had estimated DNRs between 0.064 and 0.25 (Figure 1). For the 24 nucleotide sequence datasets with a model preference of GTR, none had estimated DNRs greater than 0.5, one had an estimated DNR between 0.25 and 0.5 and the remainder had estimated DNRs between 0.0335 and 0.25.

The dsDNA virus dataset with the highest DNR was Human mastadenovirus D (DNR = 0.646), the dsRNA virus dataset with the highest estimated DNR was Porcine_rotavirus_B (0.351), the ssRNA virus dataset with the highest DNR was SARS-CoV-2 (DNR = 1.536) and the ssDNA virus dataset with the highest DNR was Torque teno sus virus (DNR = 1.56).

Therefore, while NREV12 appears to be generally more appropriate than either NREV6 or GTR for describing mutational processes in ssRNA, ssDNA, dsDNA, and dsRNA viruses, this might only be particularly relevant from a practical perspective when datasets of these viruses yield DNR estimates that are greater than 0.25. For such datasets, NREV12 (and possibly NREV 6 in some instances) might be especially useful for both determining the direction of evolution across phylogenetic trees (i.e., it could potentially be used to root these trees) and for quantifying genomic strand-specific nucleotide substitution biases (Harkins et al., 2009).

Assessing the impacts of model misspecification on phylogenetic tree inference

To determine whether it is worthwhile to use NREV12 rather than GTR for phylogenetic inference when NREV12 is the best fitting nucleotide substitution model, we used simulated datasets to compare the accuracy of phylogenetic trees inferred using these models.

We found that, regardless of dataset diversity and the nucleotide substitution model used, phylogenetic inference tended to become less accurate (i.e. wRF scores increased) as DNR increased (Figure 2). This tendency was, however, more pronounced when using a (mis-specified) GTR model than when using a (correctly specified) NREV12 model with, for any given dataset having DNR > 0, the use of NREV12 tending to yield more accurate phylogenetic trees than when GTR was used. There were, however, only statistically significant improvements (p <0.05 paired t-test) in the accuracy of phylogenetic trees inferred using NREV12 relative to those inferred using GTR in lower diversity datasets (i.e. those with APIs of 85%, 90% and 95%); and then only for DNR > 8. It is noteworthy that the highest estimated DNR in any of the empirical datasets that we analysed was 1.56 more than four-fold lower than the point where statistically significant differences in phylogenetic-inference accuracies became apparent in the simulated datasets.

Materials and Methods

Virus sequence datasets and phylogenetic trees

We obtained viral nucleotide sequences from the National Centre for Biotechnology Information Taxonomy database (http://www.ncbi.nlm.nih.gov/taxonomy) and the Los Alamos National Laboratory HIV sequence database (https://www.hiv.lanl.gov/content/index). These included gene and whole-genome sequences for viruses with ssRNA, ssDNA, dsRNA, and dsDNA genomes (datasets are summarized in Error! Reference source not found.). An outgroup sequence from a closely related virus species was added to each dataset to help root phylogenetic trees. The sequences in each of the datasets were aligned using MUSCLE (Edgar, 2004) implemented in Aliview (Larsson, 2014), and maximum likelihood phylogenetic trees were constructed from each alignment using RAxML v8.2 (Stamatakis, 2016).

Model Testing

We evaluated the fit of NREV12, NREV6, and GTR to the 141 individual sequence datasets using a previously published model test (Harkins, et al., 2009) implemented as a HyPhy package module (Pond, Frost, & Muse, 2005). This script (https://github.com/veg/hyphy-analyses/tree/master/NucleotideNonREV) and is also available as a module in the Datamonkey web server (Weaver, et al., 2018) which takes as input a rooted maximum likelihood phylogenetic tree (minus the rooting sequence) and its corresponding nucleotide sequence alignment. The three models described above: GTR, NREV6, and NREV12 are then fitted to the data using maximum likelihood (ML). The equilibrium frequencies (EF) of the GTR model match those empirically observed in the alignment, while EFs for NREV6 and NREV12 are inducted by the model, satisfying the condition πQ = 0. An additional model NREV12 +F is estimated, where the distribution of nucleotides at the root of the tree is estimated by ML, instead of being set to the empirical frequencies. Nested models were compared by the likelihood ratio test (LRT) with significance evaluated using the distribution with d = difference in degrees of freedom. For all models we also computed the small sample corrected Akaike Information Criterion (AIC-c) score.

Quantification of non-reversibility

We further defined the degree of non-reversibility (DNR) as the absolute difference between the relative rate differences of two nucleotide pairs; i.e. for two nucleotides, x and y, there exists a relative rate of x to y substitutions that we will refer to as m, and a relative rate y to x substitutions that we will refer to as n. Under the NREV12 model, the degree of non-reversibility (DNR) between x and y is defined simply as the absolute difference between m and n:(|m-n|) and will hereby be referred to as ij_DNR were i and j are two nucleotides. We use DNR as a mathematical representation of the degree to which the rates of all pairs of reverse substitutions differ from one another. For each of the 140 individual viral alignments we calculated the average DNR of the six ij_DNR estimates inferred using the NREV12 model.

Simulations for testing the impact of non-reversible evolution on phylogenetic inference

We tested the accuracy of phylogenetic tree inference under reversible and non-reversible models using simulated datasets with varying average pairwise nucleotide sequence identities (APIs) evolved under the NREV12 model with different degrees of non-reversibility (DNR). The goal of these tests was not to exhaustively evaluate model misspecification issues during phylogenetic tree inference, but rather to check, in instances where viral taxa are known to be evolving in a detectably non-reversible manner (i.e. where NREV12 or NREV6 fits the data better than GTR), whether not accounting for this might decrease the accuracy of phylogenetic inference. Using IQTREE, a phylogenetic inference program that has the option to apply an NREV12-like model (referred to in IQTEE as the UNREST model), a phylogenetic tree was inferred from an alignment of real sequences (Avian Leukosis virus) (Figure 3) with an average sequence identity (API) of ~90%. The branch lengths on this tree were then scaled to create four other phylogenetic trees representing sequences with approximately 95%, 85%, 80% and 75% API. These five trees are hereafter referred to as “true” trees and each individual tree was used as the starting point of a different set of simulations.

Figure 3

Phylogenetic trees were inferred from these 5500 simulated datasets and compared to the phylogenetic trees used to simulate the datasets (i.e., the true trees) using wRF distances to assess the impact of varying DNR on the accuracy of phylogenetic inference. We further tested whether the accuracy of phylogenetic inference could be improved for sequences that had evolved under DNR > 0 by using NREV12 instead of GTR. Specifically, for every simulated dataset a phylogenetic tree was inferred using GTR and another using NREV12 and the wRF distances of each of these trees to the true tree was determined. For each of the analysed DNRs a paired t-test was then used to compare the wRF scores of trees inferred using GTR and NREV12. We were particularly interested in determining whether trees inferred using a mis-specified model (i.e. GTR in this case) would be less accurate than trees inferred with a correctly specified model (i.e. NREV12).

To test whether failure to account for non-reversibility might decrease the accuracy of phylogenetic inference we simulated the evolution of 5,500 nucleotide sequence alignments evolved non-reversibly under varying DNR along the five true phylogenetic trees: 100 datasets per true tree per simulated degree of non-reversibility (DNR). Specifically, simulations were done using HyPhy (Pond, Frost, & Muse, 2005) with relative rates ranging from a completely reversible matrix (equation (4))

representing DNR=0 - through matrices with DNR=2, 4, 6, 8, 10, 12, 14, 16, 18 and 20 (Table 5). These baselines simulated substitution rates are reflective of those seen in empirical viral nucleotide sequence datasets.

Table 5

At each DNR, the relative rates used conformed to standard measures of non-reversibility under the Kolmogorov conditions according to which non-reversibly evolving sequence datasets should yield three irreversibility indices (IRI1, IRI,2 and IRI3) that are all non-zero (Squartini & Arndt, 2008). It should be noted that all simulations under NREV12 were performed under the stationarity criterion: πe^Qt = π (where Q is the rate matrix and π is the nucleotide frequency distribution and t≥0).

Quantifying the accuracy of phylogenetic inferences

We used the weighted Robinson-Foulds metric (wRF; implemented in the R phangorn package (Schliep, 2011)) to quantify differences between the true trees used to simulate datasets and the trees inferred from these datasets using the GTR or NREV12 models. wRF considers differences between both the topology and branch lengths of actual and inferred trees (Kuhner, 2015), (Robinson & Foulds, 1981).

Conclusion

The non-reversible nucleotide substitution model, NREV12, provides a substantially better fit to most virus nucleotide sequence datasets than does the widely used reversible substitution model, GTR. NREV12 also provides a better fit to most virus nucleotide sequence datasets than does NREV6; a non-reversible model that would be expected to best describe the evolution of double stranded genome sequences that display no strand-specific nucleotide substitution biases. This suggests that, contrary to our expectations, substantial strand-specific nucleotide substitution biases (i.e. estimated DNRs > 0.25) are common during viral evolution irrespective of genome type. Such biases should be expected for any viruses where one genome strand either is in existence for substantially longer periods of time than the other, or is more exposed to mutagenic processes than the other during transmission, replication, or gene expression.

We had anticipated that, given evidence of sequences evolving both non-reversibly and with strand-specific substitution biases, inferring trees using a model such as NREV12 that appropriately accounts for this might: (1) minimize the impact of increasing DNR on the accuracy of phylogenetic inference (i.e. wRF scores presented in blue in Figure 2 might have been expected to not increase with increasing DNR); and (2) yield significantly more accurate phylogenetic inferences than when using GTR for all datasets where NREV12 was the most appropriate model and DNRs were greater than zero. However, increasing DNR clearly decreased the accuracy of phylogenetic inference even when using NREV12, and, for datasets where DNRs were greater than zero, using GTR did not consistently yield significantly less accurate phylogenetic inferences than those attained using NREV12. From a practical perspective, choosing a non-reversible nucleotide substitution model to construct phylogenetic trees from virus genome sequences that display strand-specific nucleotide substitution biases is not guaranteed to yield more accurate phylogenetic trees. Nevertheless, in instances where strand-specific substitution biases are higher than ~0.5 (such as are found in our SARS-CoV-2, Torque teno sus virus, and Banana bunchy top virus datasets) it may be prudent to select a model such as NREV12 (such as is implemented in programs like IQTREE) over GTR as the better of two suboptimal choices.

The lack of available data regarding the proportions of viral life cycles during which genomes exist in single and double stranded states makes it difficult to rationally predict the situations where the use of models such as GTR, NREV6 and NREV12 might be most justified: particularly in light of the poor over-all performance of NREV6 and GTR relative to NREV12 with respect to describing mutational processes in viral genome sequence datasets. We therefore recommend case-by-case assessments of NREV12 vs NREV6 vs GTR model fit when deciding whether it is appropriate to consider the application of non-reversible models for phylogenetic inference and/or phylogenetic model-based analyses such as those intended to test for evidence of natural section or the existence of molecular clocks.

Data availability

All sequence alignments used during model testing in this study are included in the manuscript supplementary table one. These datasets were obtained from nucleotide sequences from the National Centre for Biotechnology Information Taxonomy database (http://www.ncbi.nlm.nih.gov/taxonomy) and also obtained from the Los Alamos National Laboratory HIV sequence database (https://www.hiv.lanl.gov/content/index). These sequence alignments have been made available as supporting files datasets.

Supplementary Files

Table S1.

Acknowledgements

The authors wish to thank the South African National Research Foundation (NRF) for funding this research under the South African Centre for Epidemiological Modelling and Analysis (SACEMA) bursary. SLKP and SW were supported in part by the U.S. National Institutes of Health (R01 AI134384 and AI140970, GM144468). DPM is supported by the Wellcome Trust (222574/Z/21/Z).

Additional information

Ethics

The University of Cape Town ethics committee declared that this research did not need ethics approval due to the use of freely accessible nucleotide sequences obtained from the National Centre for Biotechnology Information Taxonomy database (http://www.ncbi.nlm.nih.gov/taxonomy) and the Los Alamos National Laboratory HIV sequence database (https://www.hiv.lanl.gov/content/index).

Abbreviations

• API: Average pairwise nucleotide sequence identity

• NREV6: 6 Rate Non-Reversible Model

• NREV12: 12 Rate Non-Reversible Model

• GTR: General time Reversible Model

• dsDNA: Double-stranded deoxyribonucleic acid

• dsRNA: Double-stranded ribonucleic acid

• ssDNA: Single-stranded deoxyribonucleic acid

• ssRNA: Single-stranded ribonucleic acid

• DNR: Degrees of Non-Reversibility

• wRF: Weighted Robinson-Foulds distance

• LRT: likelihood ratio test

• AIC: Akaike information criterion test