{'fraction_non_alphanumeric': 0.04025157232704402, 'fraction_numerical': 0.008176100628930818, 'mean_word_length': 4.815356489945155, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Epistasis produces an excess of invariant sites in neutral molecular evolution

Substitution rate variability among sites is a common feature of protein evolution. This variability is frequently modeled by a gamma (Γ) distribution (1), which is proposed to be an emergent property of epistasis by de la Paz et al. (2). However, the conclusion of de la Paz et al. (2) is based on the analysis of highly divergent sequences, as many as 24 substitutions per site. We examined if this conclusion holds for more biologically realistic sequence divergences (0.1 to 5 substitutions per site) (3). Instead, we found that creation of a class of invariant sites is an emergent property of epistasis, and the sequence divergence level dictates whether an equalor a variable-rate model will best fit the data. In de la Paz et al. (2), neutral molecular evolution was simulated for protein families with a fitness function defined by a model inferred via direct coupling analysis (e.g., Fig. 1A). We counted the number of sites that experienced zero, one, two, three, and more substitutions directly from the generation by generation simulated sequences (e.g., Fig. 1B), which avoids confounding effects of estimation errors. This site frequency spectrum (SFS) will follow a Poisson distribution when the data are compatible with a single (S) evolutionary rate across sites (4). SFS will follow a negative binomial distribution if the evolutionary rates among sites are Γ distributed (1). We tested model fits to SFS at different levels of sequence divergence, allowing for a category of invariant sites (I+) with S and Γmodels (5). For each dataset, the best-fit model from four candidates (S, I + S, Γ, and I + Γ models) was selected using corrected Akaike information criterion (AICc) because these models are not all nested (6, 7). The S model fits best for a vast majority of datasets at sequence divergences less than 0.5 substitutions per site (Fig. 1C). That is, epistasis did not create significant evolutionary rate variability among sites at low sequence divergences. For divergences larger than 0.5 substitutions per site, models allowing for an invariant class of sites (I + S and I + Γ) provide the best fits (Fig. 1C). Similar patterns are observed for the other nine protein domain families analyzed by de la Paz et al. (2) (Fig. 2). The Γ model alone was not the dominant model except for one case (Fig. 2). Therefore, at sufficient evolutionary distances, epistasis generates many invariant sites and substitution rate variability, likely because incompatible substitutions at coupled sites are under greater negative selection than substitutions at loosely coupled and uncoupled sites. At relatively low divergences, fewer coupled sites are mutated, likely resulting in less differential negative selection among sites. So, even a single-rate model may fit the data well. Overall, our results suggest that epistasis provides a mechanistic explanation for the abundance of invariant sites beyond what is explained by an equal (S)or a variable (Γ)-rate model. Consequently, models incorporating a separate class of invariant sites best describe evolutionary rates of proteins.