Modeling site-and-branch-heterogeneity with GFmix

Modeling site-and-branch-heterogeneity with GFmix | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Modeling site-and-branch-heterogeneity with GFmix View ORCID Profile Charley G. P. McCarthy , View ORCID Profile Edward Susko , View ORCID Profile Andrew J. Roger doi: https://doi.org/10.1101/2025.08.07.669136 Charley G. P. McCarthy 1 Institute for Comparative Genomics, Dalhousie University , Halifax, NS B3H 4R2, Canada 2 Department of Biochemistry and Molecular Biology, Dalhousie University , Halifax, NS B3H 4R2, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Charley G. P. McCarthy For correspondence: charley.mccarthy{at}dal.ca edward.susko{at}gmail.com Edward Susko 1 Institute for Comparative Genomics, Dalhousie University , Halifax, NS B3H 4R2, Canada 3 Department of Mathematics and Statistics, Dalhousie University , Halifax, NS B3H 4R2, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Edward Susko For correspondence: charley.mccarthy{at}dal.ca edward.susko{at}gmail.com Andrew J. Roger 1 Institute for Comparative Genomics, Dalhousie University , Halifax, NS B3H 4R2, Canada 2 Department of Biochemistry and Molecular Biology, Dalhousie University , Halifax, NS B3H 4R2, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Andrew J. Roger Abstract Full Text Info/History Metrics Supplementary material Preview PDF A bstract Phylogenetic trees are often inferred from protein sequences sampled from diverse taxa across the tree of life. The compositions of these amino acid sequences may be heterogeneous across both sites and branches, particularly if deep phylogenetic divergences are the focus. Under some conditions, failure to model this compositional heterogeneity can lead to phylogenetic artefacts. However, the computational cost of phylogenetic inference with models accounting for compositional heterogeneity can be prohibitive. The originally proposed site-and-branch-heterogeneous GFmix model accounts for changing relative frequencies of G, A, R, and P (GARP) vs. F, Y, M, I, N, and K (FYMINK) amino acids resulting from extreme variation in G+C content among taxa. This GFmix model modifies a fitted site-heterogeneous profile mixture model in a branch-specific manner using parameters that reflect branch-specific amino acid compositions. This approach has been shown to improve likelihoods and reduce compositional artifacts. However, the original implementation of the model includes constraints which may sacrifice accuracy for computability and is limited to modeling variation in GARP/FYMINK composition. Here we investigate the properties of the original GFmix model in greater depth and present several improvements to the model. The improved GFmix models permit fewer constraints on branch-specific composition parameters, allow modeling of user-defined compositional heterogeneity, and provide for full maximum-likelihood optimization of parameters. We have also developed new methods for detecting compositional heterogeneity directly from sequence data. Analyses of simulated site-and-branch-heterogeneous data indicates that the improved GFmix models better estimate branch-specific compositions and branch lengths in heterogeneous trees. We applied the various versions of the GFmix model to a real dataset with known compositional heterogeneity artefacts. We find that the most complex GFmix model with full maximum likelihood parameter optimization consistently supports the correct tree over the artefactual tree with improved likelihoods. All versions of the GFmix model are available from https://www.mathstat.dal.ca/~tsusko/software.html . I ntroduction Deep phylogenies are often inferred from aligned sequences of multiple proteins retrieved from many species with very different lifestyles and histories from across the tree of life. Increasing the number of sequences and/or species sampled for deeper phylogenies may not always be beneficial ( Philippe et al. 2011 ). Resolving tricky relationships like locating the roots of the animals and land plants subtrees has remained difficult, even with improved data quality and sampling ( Feuda et al. 2017 ; Whelan et al. 2017 ; Wickett et al. 2014 ; Puttick et al. 2018 ). Many of these deep-time phylogenetic problems have remained difficult to resolve in part because most phylogenetic models do not capture the complexity of sequence evolution at this deep evolutionary scale and are therefore prone to systematic errors in tree inference. Substitution processes across sequences may be heterogeneous due to site-specific functional or structural constraints ( Halpern and Bruno 1998 ; Goldstein 2008 ). Amino acid composition across species may be heterogeneous due to genomic changes arising from loss of DNA repair mechanisms, genome reduction or adaptation to extreme environments ( Acosta et al. 2015 ; Zeldovich et al. 2007 ; Siglioccolo et al. 2011 ). Failing to model or mismodeling such heterogeneity can induce long-branch attraction (LBA) artefacts and lead to incorrect tree estimation ( Felsenstein 1978 ; Foster 2004 ; Brinkmann et al. 2005 ). The simplest protein sequence evolution models typically do not capture the heterogeneity of the substitution process over sites and over branches. Amino acid substitutions over the tree at a site are modeled as occurring according to a Markov process with rate matrix Q ij = S ij π j . The stationary frequencies, π j , are usually estimated from the data, but the exchangeability rates between amino acids i and j, S ij , are fixed. Common exchangability matrices such as LG and WAG were estimated based on large empirical data sets. WAG was estimated assuming homogeneous evolutionary rates and substitution processes across all sites ( Whelan and Goldman 2001 ), whereas LG was optimized allowing for differing rates across sites ( Le and Gascuel 2008 ). Heterogeneity of process over sites is typically accommodated via mixtures of these base models where some parameter in the base model differs across mixture classes. For instance, heterogeneity of evolutionary rates across sites is typically modeled using discrete rate multipliers derived from the Gamma (Γ) distribution or estimated from data ( Yang 1994 ; Soubrier et al. 2012 ). Heterogeneous substitution processes across sites have been accommodated in several ways. In some cases, mixtures of empirically-estimated rate matrices are used ( Le et al. 2012 ) or alternatively different site classes are modeled by using a common set of exchangeabilities and a mixture of frequency profiles ( Si Quang et al. 2008 ). The latter profile mixture models allow the π j to vary over mixture classes and assume that all substitution patterns across sites can be adequately described using a set of equilibrium amino acid frequency profiles. The frequency profiles may be derived from empirical data as in the C-series ( Si Quang et al. 2008 ) and UDM ( Schrempf et al. 2020 ) models, estimated via a Dirichlet process as in the Bayesian CAT model ( Lartillot and Philippe 2004 ) or estimated by maximum likelihood (ML) using the composite likelihood method MAMMaL ( Susko et al. 2018 ). Combining site profile models with site-rate heterogeneity models - e.g. LG+C20+Γor CAT-GTR+Γ- is common in large-scale phylogenetics studies and is readily implemented in software such as IQ-TREE and PhyloBayes ( Lartillot et al. 2009 ; Minh et al. 2020 ; Wong et al. 2025 ). Heterogeneity of amino acid composition across species has been modeled using several approaches. In the ML framework, the branch non-homogeneous and non-stationary Correspondence and Likelihood Analysis (COaLA) model uses correspondence analysis to capture the most important features of compositional variation observed in the data while minimizing the number of parameters requiring optimization ( Groussin et al. 2013 ). In the Bayesian framework, the node-discrete compositional heterogeneity (NDCH) model fits multiple amino acid composition vectors across a tree ( Foster 2004 ; Foster et al. 2009 ) whereas the CAT-BP model varies the CAT model across a tree to reflect branch-specific amino acid composition ( Blanquart and Lartillot 2006 ; Blanquart and Lartillot 2008 ). Of these three models, only CAT-BP accommodates both across-site-and across-branch-heterogeneity ( Blanquart and Lartillot 2008 ). Applying the CAT-BP model requires substantial time and computational overhead, and because it is implemented in a Bayesian framework, final interpretations of results are dependent on tree and parameter convergence ( Williams et al. 2021 ). Similar caveats apply to NDCH and the site-heterogeneous CAT-GTR+Γmodel ( Williams et al. 2021 ). Alternative ways to address branch-heterogeneity in amino acid composition have been through data curation by recoding amino acids into larger biochemical groups or removing very heterogeneous sites in a dataset ( Embley et al. 2003 ; Susko and Roger 2007 ; Whelan et al. 2017 ; Martijn et al. 2018 ). Both approaches may improve model fit and parameter convergence for heterogeneous data ( Feuda et al. 2017 ; Whelan et al. 2017 ). However, the efficacy of recoding approaches has been questioned ( Li et al. 2021 ; Hernandez and Ryan 2021 ; Foster et al. 2022 ) and site-removal approaches may yield unexpected results ( Francis and Canfield 2020 ). GFmix is a novel site-and-branch-heterogeneous model implemented in a maximum-likelihood (ML) framework. GFmix was initially developed to resolve several conflicting relationships within the mitochondria-Alphaproteobacteria tree ( Muñoz-Gómez et al. 2022 ). Mitochondria and multiple alphaproteobacterial phyla have undergone convergent genome reduction, which is correlated with increased genomic AT-content and elevated FYMINK-content in their proteomes ( Ettema and Andersson 2009 ; Hershberg and Petrov 2010 ; Acosta et al. 2015 ). Other alphaproteobacterial phyla have comparatively GC-rich genomes and GARP-rich proteomes. This GARP/FYMINK heterogeneity was thought to drive conflicting placements of mitochondria with respect to Alphaproteobacteria, and the monophyly of FYMINK-rich alphaproteobacterial phyla ( Martijn et al. 2018 ; Muñoz-Gómez et al. 2019 ; Fan et al. 2020 ). Application of the GFmix model to the mitochondria-Alphaproteobacteria tree favored placing mitochondria as sister to Alphaproteobacteria and multiple independent origins of FYMINK-rich proteomes within Alphaproteobacteria ( Muñoz-Gómez et al. 2022 ). The model has subsequently been extended to other cases of compositional heterogeneity affecting the root of eukaryotes and relationships within Archaea ( Baker et al. 2024 ; Baker et al. 2025 ; Williamson et al. 2025 ). The GFmix model estimates a branch-specific parameter for each branch of a tree and uses this parameter to make branch-specific modifications to the mixture class frequencies of a previously-fitted profile mixture model ( Muñoz-Gómez et al. 2022 ). The likelihood of the tree given this modified model is then re-estimated. Although GFmix improves tree inference in the presence of compositional heterogeneity, there are limitations to the original implementation. No re-optimization of model and branch length parameters occurs before likelihood re-estimation under GFmix. The branch-specific parameters used to modify the mixture model are estimated using the tree and the observed frequencies of amino acids for the taxa under study and impose parameter constraints ( Muñoz-Gómez et al. 2022 ). The approach is thus not a full maximum likelihood implementation. Although this makes it computationally less intensive, which can be valuable for large data sets, it makes parameter estimation less accurate. Finally, in its original implementation, the model only accounts for GARP/FYMINK heterogeneity resulting from GC-content variation ( Muñoz-Gómez et al. 2022 ). Here, we present improved implementations of the GFmix model intended to loosen the constraints of the original implementation. All implementations now allow user-defined compositional heterogeneity to be modelled, and the most extensive implementations allow some or all model and tree parameters to be re-optimized in a maximum-likelihood framework. We assessed the performance of each implementation of the GFmix model using 16-taxa and downsampled 4-taxa datasets simulated under site-and-branch-heterogeneous GARP/FYMINK conditions. We find that the improved GFmix models more accurately estimate branch lengths and compositional shifts in 16-taxa simulations, and that this improved performance is also observed in more information-sparse 4-taxa simulations. Modeling GARP/FYMINK heterogeneity is appealing because of the connection that those groups of amino acids have to GC-content but, in some settings, there is reason to expect that the predominant heterogeneity will be for different groups (e.g. in Baker et al. (2024) the heterogeneity in halophilic Archaea is between DE/IK).Accordingly, we have developed several methods to identify the the most variable groups of amino acids. We find that the identification methods that use model-based optimization criteria more accurately identify the true groups of varying amino acid compositions in simulations, albeit with an associated increase in computational runtime. Finally, we applied each GFmix implementation to a real dataset with known compositional heterogeneity arising from GC-content variation. We chose as our basis a 54-taxon dataset comprised of nuclear-encoded protein sequences from Rhodophyta and Viridiplantae and nucleomorph-encoded protein sequences from Cryptomonada ( Novak et al. 2024 ). Nucleomorph genomes are highly reduced, and thus AT- and FYMINK-rich at the genomic and proteomic level, respectively ( Moore and Archibald 2009 ). We downsampled the original dataset, retaining the Cryptomonas paramecium nucleomorph, and added nucleomorph-encoded protein sequences from the Rhizarian alga Bigelowiella natans ( Gilson et al. 2006 ). We compared support for two hypotheses - one in which the B. natans and Cryptomonas paramecium nucleomorphs branch together at the base of Rhodophyta, and one in which the nucleomorphs branch separately within Rhodophyta and Viridiplantae based on the biological consensus of their respective origins ( Ishida et al. 1999 ; Douglas et al. 2001 ). We find that the most extensive GFmix implementation consistently supports the latter hypothesis, with improved likelihoods under newer implementations of the model. Overall, our results indicate that GFmix improves likelihood and reduces compositional heterogeneity artifacts, and that newer implementations of the model show substantial improvements in performance over the initial implementation. M aterials and M ethods Implementations of GFmix The GFmix Model The GFmix model for the evolution of an amino acid at a site along a tree builds on the usual base model where, conditional upon the ancestral amino acid, evolution along a branch occurs independently of neighbouring or ancestral branches and according to a continuous-time, time-reversible Markov chain. Time reversibility implies that the Markov chain rate matrix can be decomposed as Q ij ∝ S ij π j for i ≠ j , where π j is the stationary frequency of amino acid j and S is a symmetric exchangeability matrix with positive entries. Here Q ii = −∑ j ∣ j ≠ i Q ij and the proportionality constant is determined by −∑ i π i Q ii = 1 which leads to branch lengths being interpretable as expected numbers of substitutions. We use the fixed LG exchangeability matrix of Le and Gascuel (2008) in the GFmix model; no parameters in S require estimation. To emphasize dependence of Q on π we use the notation Q ( π ). Let z denote the amino acids at all of the nodes for a site, including internal nodes, and let p ( e ) and c ( e ) denote the parent and child nodes for a branch e . Under the base model, the probability of z is a product over edges of matrix exponentials, multiplied by the frequency of the amino acid, z r , at the root, We have allowed that the frequencies Π = [ π (1) … π (2 m− 2) ] differ over branches or at the root. In the base model these would be the same for each branch and any root location would give the same probability. Decomposing z as [ x, y ], where y are the amino acids at internal nodes, the probability of the observed data, x , is obtained by summing (1.1) over the unobserved internal node data, The pruning algorithm of Felsenstein (1973) is used to feasibly calculate (1.2). Profile mixture models ( Si Quang et al. 2008 ) and the discretized Gamma rate model of Yang (1994) allow variation of stationary frequencies and evolutionary rates over sites. They assume frequency classes and rates at sites are drawn from discrete distributions, independenty of each other and independently across sites. The frequency distribution usually has fixed frequency class profiles, π ( c ) , that are drawn with probability w c , the latter requiring estimation. We use the C-series frequencies ( Si Quang et al. 2008 ) throughout our examples, but each GFmix implementation allows other profiles to be input. The rate distribution is the discretized Gamma distribution of Yang (1994) where rate r k ( α ) arises with probability 1 /K ; α is the shape parameter of the distribution and K (set to 4 in examples) the number of rates. Under these models, the conditional probability of the data at a site with rate class k and frequency class c is p M ( x | c, k ; π, t, α ) = p B ( x ; Π ( c ) , r k ( α ) t ) where Π ( c ) = [ π ( c ) … π ( c ) ] and p B ( x ; Π ( c ) , t ) is determined from (1.1)-(1.2). Since ( c, k ) are unobserved, the probability of the observed data x is calculated as The model described has fixed stationary frequencies of amino acids over branches for any given profile class c . GFmix extends the model by allowing these frequencies to vary over the tree but in a specific manner where 𝒢 and ℱ are groups of amino acids; 𝒢= GARP and ℱ= FYMINK in the original implementation. The additional branch-specific parameters require estimation but the multiplicative constants are implied by the constraint . The additional parameters and assumes that the biases towards amino acids in𝒢, ℱ or𝒪, biases that are determined by and , applies over all sites regardless of their rate and frequency profile. Let Π ( c ) ( γ ) = [ π ( c 1) … π ( c 2 m− 2) ]; note that the root edge also has a γ ( r ) parameter vector. The site likelihood is calculated similarly as in (1.3): The GFmix methods described below differ largely in how they provide parameter estimates for the model. Given a set of parameters, the end result is the same, a log like-lihood for the model described above. That log likelihood assumes evolution over sites is independent and thus is obtained as follows Maximum-likelihood implementations - GFP and GFF The maximum likelihood implementations of GFmix are conceptually the most straightforward. They simply maximize the log likelihood (1.6) using the L-BFGS-B algorithm implementation of Zhu et al. (1997) . GFF optimizes all parameters. Partial optimization over γ is done by GFP, holding all other parameters fixed. The other parameters are estimated by maximum likelihood using IQ-TREE under the corresponding profile mixture model (equivalently with ). Original GFmix implemention - OGF The original GFmix implementation was first described in Muñoz-Gómez et al. (2022) for 𝒢 / ℱ equal to GARP/FYMINK and modified in Baker et al. (2024) to accommodate arbitrary 𝒢 and ℱ. For a given branch e , assume that we know b e , the ratio of the total expected frequency ratio of ℱamino acids to the total frequency of 𝒢 amino acids; OGF estimates b e from the vectors, over taxa, of observed frequencies of amino acids as described below. The and and then satisfy the following equations This gives C + 1 equations in C + 2 unknowns. The original GFmix implementation imposed the additional model condition This constraint implies that the weighted average multiplier of the 𝒪 class is 1. That constraint is not essential and not present the maximum likelihood implementations GFP and GFF. To estimate the b e , OGF aggregates all sequences that descend from branch e and then estimates b e as the ratio of the total observed frequency of amino acids in 𝒢 to total observed frequency of amino acids in ℱ. It then solves the C +2 equations in C +2 unknowns given in (1.7)-(1.8) to obtain estimates of and and . Sum of squares difference between node and model frequencies - SSF For a given frequency class c and rate class k , the conditional probability of amino acid j at the root is . Given the conditional probabilities at the parental node p ( e ) of branch e probabilities at the child node c( e ) can be calculated through Starting with , (1.9) can be used recursively, traversing edges in the tree from root to tip to calculate for all nodes l . The expected frequency of amino acid j is then obtained as , were we now explicitly indicate dependence on parameters. We use the algorithm just described to get node frequencies in simulations where we know the parameter values. We also used it to obtain estimates of γ using least squares. Fixing at maximum likelihood likelihood values obtained under the corresponding profile mixture model (equivalently with ), we minimize where sums are over terminal nodes s and amino acids j . Determining compositionally heterogeneous groups The original GFmix model assumed fixed groups, 𝒢 / ℱ = GARP/FYMINK. The rationale for this was both empirical, in that overall GARP and FYMINK frequencies often show a lot of variation over taxa, and because these two groups of amino acids are GC (GARP) and AT (FYMINK) rich ( Muñoz-Gómez et al. 2022 ). The improved GFmix models now allow 𝒢 / ℱ group determination to be based on the data at hand using methods that we now introduce. Binomial test of two proportions Compositional heterogeneity is first identified from data using a binomial test of two proportions ( Baker et al. 2024 ; Williamson et al. 2025 ). This approach requires two apriori specified groups of taxa for which compositional heterogeneity between groups is suspected. For instance, Baker et al. (2024) separated taxa according to whether they exist in hypersaline environments or not. For a given dataset, taxa are divided into two groups. For each amino acid, the composition bias between the two taxon groups is computed as a Z-score in the following equation: where X 1 , X 2 are the total numbers of that amino acid and n 1 , n 2 are the total number of all 20 amino acids across the two taxon groups respectively. This test assumes that the proportions of an amino acid across taxa is approximately normal, with the null hypothesis that p 1 = p 2 . | Z | > 1.96 indicates rejection of the null hypothesis at significance level p 1.96), ℱ ( Z < − 1.96) or 𝒪 (| Z | ⩽ 1.96) on the basis of this Z -score. General optimization routine The main reason for the GFmix model adjustment (1.4) is to allow the frequencies of amino acids within groups to vary over taxa. So we seek groups where the average within-group 𝒢 (similarly ℱ) amino acid frequency shows evidence of substantial variation over edges. Since the simplifying assumption in (1.4) is that frequencies within a class are modified by the same scalar factor, we also seek choices of ℱ and ℱ classes so that frequencies of amino acids within the 𝒢 (or similarly the ℱ) class differ for an edge from the overall frequencies in a similar fashion. If, for instance, 𝒢 /ℱ =GARP/FYMINK is a good choice, then when the frequency of G is elevated for a taxon relative to the overall frequencies, we also expect A to be elevated. In other words, the amino acids within a class should covary across edges. Our approach to meeting the objectives of 𝒢 /ℱ class selection outlined above is to choose the classes that minimizes a goodness of fit criterion. We describe three different possible criteria below. The results of the binomial approach give starting classes. Given a criterion, optimization of the class memberships follows a non-greedy approach, in which both the amino acids in the current classes and all permissible single-amino acid swaps between classes at a given point in the algorithm are assessed under some optimality criterion. Permissible swaps follow a bipartite relationship between classes: single-amino acid swaps are permitted between 𝒢 and 𝒪, ℱ and 𝒪, but not 𝒢 and ℱ. Once all swaps have been assessed, the best-performing swapped classes are compared to the current classes’ with respect to their optimality criteria. If the swapped classes improve upon the current classes, the swapped classes replace the current classes and the optimization routine repeats until no further improvement can be made. The resulting 𝒢,ℱ and 𝒪 are considered the final classes for this optimization routine. Criterion 1: χ 2 test of homogeneity The χ 2 test is a standard measure of compositional homogeneity in phylogenetics ( Foster 2004 ). The test statistic is the standard χ 2 test of homogeneity of multinomial frequencies over groups. In our application, each taxon plays the role of a group and the multinomial frequencies are the summed 𝒢, ℱ and 𝒪 for each taxon. In using the χ 2 test as an optimization criterion, we attempt to identify 𝒢 and ℱ which maximizes the χ 2 test statistic for a given alignment – in other words, the most compositionally-heterogeneous 𝒢 and ℱ. For a given 𝒢 and ℱ, the counts of all amino acids within each class are summed and the χ 2 test is performed with the summed counts using the chisq.test() function in R ( R Core Team 2025 ). Summation of the counts of 𝒢 and ℱ ensures that the χ 2 test is always performed with two degrees of freedom, which means the test statistic is comparable for any size of 𝒢 or ℱ. Criterion 2: Sum-of-squares difference using SSF implementation SSF estimates branch-composition parameters that minimize the sum, over taxa and amino acids, of squared (SS) differences between the observed frequency of an observed amino acid for a taxa and the corresponding expected frequency under the model. It does this for a given choice of 𝒢 and ℱ but returns as a by-product, the optimal SS for that choice. That SS can also be used as a criterion. We attempt to identify 𝒢 and ℱ which minimize the sum-of-squares difference. Criterion 3: Likelihood using OGF implementation Finally, the original implementation of GFmix (OGF) can itself be used as optimization criteria. We attempt to identify 𝒢 and ℱ which maximizes the likelihood of a given alignment, tree and mixture model. For the original implementation of GFmix, 𝒢 and ℱ are used to estimate branch-specific composition parameters b e as detailed above. Simulated data analysis Generating 16-taxon datasets Simulated site-and-branch-heterogeneous sequence data was generated using AliSim, a component of the IQ-TREE 3 software package ( Ly-Trong et al. 2022 ; Wong et al. 2025 ). Data was generated using a balanced 16-taxon simulating tree, which was divided into four 4-taxon clades and rooted at the midpoint. Clades were labelled A-D with taxa in each clade labelled A1-A4 etc. Branch lengths are indicated in Figure 1 . One hundred replicate sequence data sets were generated, each with a length of 10,000 sites with no gaps. Site-heterogeneity was simulated using an LG+C20+Γmodel with four discrete rate categories and α = 0.5 ( Le and Gascuel 2008 ; Si Quang et al. 2008 ; Yang 1994 ). For branch-heterogeneity, all stem and tip branches in clades B and D were assigned a target b e = 0.1, whereas all other branches were assigned a target b e = 1 under the assumptions of the original GFmix model ( Muñoz-Gómez et al. 2022 ). As simulating branch-heterogeneity is not a native feature of AliSim, a custom approach was used in which site-and-branch heterogeneous replicates were simulated on a class-by-class basis and these class sub-replicates were concatenated to form the final replicate for downstream analysis. Detailed explanation of this approach is provided in Supplementary Information . Download figure Open in new tab Fig. 1. Tree topologies used for simulating and downsampled site-and-branch-heterogeneous sequence data. Branches and taxa coloured in grey have amino acid composition π ( GARP ) < π ( FY MINK ), branches and taxa coloured in black have amino acid composition π ( GARP ) = π ( FY MINK ). Simulating branch lengths indicated underneath branches. Top: 16-taxon tree. Bottom: 4-taxon tree. Generating 4-taxon datasets Phylogenetic models generally perform better when more taxa are available, particularly if these taxa can break long branches within a tree ( Graybeal 1998 ). To assess model performance in more taxa-sparse cases, all 16-taxon replicates were downsampled to 4-taxa replicates by extracting one taxon from each clade A, B, C and D. Site-heterogeneous model fitting Site-heterogeneous model fitting was performed for all simulated data using IQ-TREE 3 ( Wong et al. 2025 ), with a LG+C20+Γmodel with mixture weight optimization ( Le and Gascuel 2008 ; Si Quang et al. 2008 ; Yang 1994 ) and using the simulating tree as a fixed tree. Branch-heterogeneous model fitting Branch-heterogeneous model fitting was performed for all simulated data using each of the four GFmix implementations outlined above ( Muñoz-Gómez et al. 2022 ; Baker et al. 2024 ). For each implementation, the compositional heterogeneity to model was set to GARP/FYMINK. Assessment of GFmix model performance Estimating branch lengths with and without maximum-likelihood optimization We assessed the performance of GFF in re-estimating branch lengths by comparing the observed branch lengths estimated under LG+C20+Γand LG+C20+Γ+GFF against the expected simulating branch lengths. Branch lengths estimated under LG+C20+Γand LG+C20+Γ+GFF were extracted for each branch across all replicates. Root mean squared error (RMSE) for a given estimation method and model parameter is the square root of the average, over the 1000 simulated data sets, of the squared differences between the estimated parameter and the true one. For the 4-taxon data sets these were calculated separately for each external branch. The lengths of the two root branches were summed prior to calculation of a squared error so that the RMSEs were for branch lengths in the corresponding unrooted tree. For 16-taxon data sets, results were summarized for clades by taking averages of branch-specific RMSEs. For example, an RMSE reported for Clade A was calculated as the average of the RMSEs for the tip branches A1-A4 and the RMSE for the stem. The reason for this is that because of the symmetry of the tree and processes, the same behaviour is expected for the A1 and A4 branch, eliminating the need to consider them separately. Plots comparing the difference between simulating branch lengths and observed branch lengths in both 16-taxon and 4-taxon datasets were generated using ggplot2 ( Wickham 2016 ). Estimating node GARP/FYMINK frequency ratio Under the GFmix model, expected amino acid frequencies vary continuously over the tree. We assessed the ability of methods to estimate the ratio of aggregate GARP/FYMINK frequencies by comparing those ratios calculated with estimated parameters against the true ratios calculated using simulated parameters. Because of the symmetry of the processes and branch lengths for certain branches, long-run RMSEs of estimation are known to be the same, for instance, for a tip in the A clade and a tip in the C clade. Average RMSE were thus reported to simplify presentation. For the 16-taxa datasets, RMSEs were calculated as follows. For the four simulated clades, GARP/FYMINK ratio RMSEs were separately averaged over A and C tips (A+C tips), A and C stems (A+C stems), B and D tips (B+D tips) and B and D stems (B+D stems). For the two nodes closest to the root, the root mean squared error was calculated as the square root of the mean of the squared error of both nodes combined. RMSEs for the root node were taken from the square root of the mean of the squared error of every root node across all 100 replicates. For the 4-taxa datasets, RMSEs were averaged over the A and C branches, and, separately, over the B and D branches. Similarly as for the 16-taxon data sets for the two internal nodes, the average RMSE - averaged over the two internal nodes - was calculated as well as the the RMSE for the root node. Plots comparing the difference between simulating GARP/FYMINK frequency ratios and observed GARP/FYMINK frequency ratios in both 16-taxon and 4-taxon datasets were generated using ggplot2 ( Wickham 2016 ). Assessing methods of identifying compositional heterogeneity We assessed the performance of our four methods for identifying compositional heterogeneity — the binomial test of two proportions and the three optimization criteria detailed above — across all 16-taxon and 4-taxon replicates. As the simulating heterogeneity for all replicates was GARP/FYMINK, our assessment was primarily based on how often a method correctly assigned 𝒢 = GARP and ℱ = FYMINK and, secondarily, how often a method misassigned 𝒢,𝒪 and ℱ. Each method of identification was performed across each of the 100 16-taxon and 4-taxon replicates, and the 𝒢,𝒪 and ℱ assignments of each method were tabulated. Unique patterns of 𝒢 /𝒪 /ℱ assignment across all methods were identified, and the number of times each method had made each assignment was tabulated. This 𝒢 / 𝒪 /ℱ assignment per-method information was visualized in an UpSet-like plot using ggplot2 ( Lex et al. 2014 ; Wickham 2016 ). Additional assessment of each identification method is detailed in Supplementary Information . Applying GFmix to real data Algal nuclear and nucleomorph dataset Our real data analysis utilized a dataset designed to determine the position of Cryptomonad nucleomorphs relative to Rhodophyta ( Novak et al. 2024 ), downsampled and modified to include Chlorarachniophyte nucleomorph sequences. Nucleomorphs are highly-reduced nuclear genomes that exist within plastids acquired through separate secondary endosymbosis events in Cryptomonads (Cryptista) and Chlorarachniophyceae (Rhizaria) ( Moore and Archibald 2009 ). These nucleomorph genomes are AT-rich at the genomic level and FYMINK-rich at the proteomic level. Cryptomonad nucleomorphs are generally thought to be derived from a deep-branching Rhodophyte ( Douglas et al. 2001 ) and Chlorarachniophycete nucleomorphs are thought to be derived from a Ulvophyte donor within Viridiplantae ( Ishida et al. 1999 ). Novak et al. (2024) sampled 180 protein sequences across 54 taxa from Rhodophyta, Viridiplantae and four nucleomorph genomes from Cryptomonada. Novak et al. (2024) found strong support for the placement of Cryptomonad nucleomorphs as sister to extremophilic Cyanidiophytina at the base of Rhodophyta. Adding B. natans nucleomorph data to the Novak et al. (2024) dataset Novak et al. (2024) did not sample any Chlorarachniophycete nucleomorphs in their analysis. To create a data set where the long-branches and composition biases of the two kinds of nucleomorphs may lead to artefactual attraction, we added protein sequence data from the nucleomorph of the Chlorarachniophyceae alga Bigelowiella natans ( Gilson et al. 2006 ) to the Novak et al. (2024) dataset as follows. Homologs from the B. natans nucleomorph genome were first identified using sequence similarity searches against Cryptomonas paramecium nucleomorph proteins in the dataset with BLASTp (e-value = 1 e − 4 ) ( Altschul et al. 1990 ), and double-checked by sequence similarity searches against the non-redundant protein sequence database ( Sayers et al. 2024 ). This retained 66 protein sequences from the Novak et al. (2024) dataset with homologous data from the B. natans nucleomorph. These sequences were re-aligned using MAFFT and all sequences were concatenated into a single 17,696-site dataset using R ( Katoh and Standley 2013 ; R Core Team 2025 ). Downsampling Novak et al. (2024) dataset and tree hypotheses We downsampled the modified Novak et al. (2024) dataset based on phylogenetic diversity, attempting to include at least one member of each major clade from Rhodophyta and Viridiplantae sampled in the original dataset. We retained 16 taxa for our analysis - nine Rhodophyta, five Viridiplantae, and the nucleomorphs of Cryptomonas paramecium and B. natans . We chose two tree hypothesis to test for our downsampled dataset following the phylogeny in Novak et al. (2024) . The first tree (“NMs Apart”) places the two nucleomorphs as evolving independently: the C. paramecium nucleomorph as sister to Cyanidioschyzon merolae at the base of Rhodophyta, and the B. natans nucleomorph sister to the closest sampled relative of Ulvophyceae - Volvox carteri - within Viridiplantae. The second tree (“NMs Together”) places the two nucleomorphs as sister taxa adjacent to C. merolae at the base of Rhodophyta. We assume that “NMs Apart” is the correct tree based on biological consensus ( Ishida et al. 1999 ; Douglas et al. 2001 ; Moore and Archibald 2009 ; Novak et al. 2024 ) whereas “NMs Together” is an artefact of long-branch attraction and shared compositional bias. Site-heterogeneous model fitting Site-heterogeneous model fitting was performed for each dataset using IQ-TREE 3 ( Minh et al. 2020 ) with tree search constrained to the topologies detailed above. The LG+C20+Γmodel was fit for each tree using mixture weight optimization (i.e., the -mwopt flag). The estimated frequency class (+F) was not included in model fitting, due to its potential compromising effects on tree inference ( Baños et al. 2023 ). Branch-heterogeneous model fitting Branch-heterogeneous model fitting was performed for each dataset using each of the four GFmix implementations outlined above. Trees were rooted between Viridiplantae and Rhodophyta. Initial model fitting was performed using the default 𝒢 / ℱ = GARP/FYMINK. Additional fitting was performed using 𝒢 / ℱ as determined by the bin optimization approach detailed above with the criterion function set to the OGF likelihood (Criterion 3). R esults New implementations of GFmix GFmix is a site-and-branch-heterogeneous model which modifies a site-heterogeneous profile mixture model to reflect branch-specific amino acid composition ( Muñoz-Gómez et al. 2022 ). The general GFmix model assumes three sets of amino acids 𝒢, 𝒪 and ℱ.𝒢 and are ℱ amino acids which are assumed to be heterogeneous in composition across a tree, whereas 𝒪 is assumed to be homogeneous in composition. In this study, we assessed the performance of various implementations of the GFmix model: Original GFmix - OGF: the original implementation of GFmix first outlined in Muñoz- Gómez et al. (2022) and extended to model user-defined 𝒢 / ℱ compositional heterogeneity in Baker et al. (2024) . SSFreq - SSF: estimating branch-composition parameters that minimize the sum-of-squares differences between observed and model amino acid frequencies over taxa. GFmix Partial - GFP: estimating separate branch-specific composition parameters for 𝒢 and ℱ amino acids, γ 𝒢 and γ ℱ respectively, using ML estimation GFmix Full - GFF: estimating (or re-estimating) all model, branch-specific composition and branch-length parameters using ML estimation A brief comparison of each implementation is provided in Table 1 and further information can be found in Materials and Methods and Supplementary Information . View this table: View inline View popup Download powerpoint Table 1. Summary of GFmix implementations assessed in this study, and abbreviations used throughout this study. 𝒢 and ℱ: classes of compositionally heterogeneous amino acids; b e : the estimated frequency ratio b of 𝒢 and ℱ amino acids at branch e of a tree; : the scaling constant for 𝒢 amino acids at branch e of a tree, : the scaling constant for ℱ amino acids at branch e of a tree, w c : weights of a profile mixture model; α : shape parameter of the Γdistribution of site-rate categories. Assessing GFmix performance To assess the performance of each GFmix implementation, we simulated 100 16-taxon alignments under an LG+C20+Γsite-heterogenous model using AliSim ( Ly-Trong et al. 2022 ) with a custom workflow for simulating GARP/FYMINK branch-heterogeneity (see Materials and Methods and Supplementary Information ). We also downsampled these alignments to 4 taxa to assess model performance under more data-sparse conditions ( Graybeal 1998 ). Information on the simulating and downsampled trees for all 16- and 4-taxon replicates is provided in Figure 1 . Estimation of branch lengths Three of the four GFmix implementations compared in this study use branch lengths as estimated under the site-heterogeneous LG+C20+Γmodel by IQ-TREE, whereas the full GFmix model (GFF) is able to re-estimate these branch lengths through maximum-likelihood optimization. Thus, our comparison of branch length estimation for our simulated 16-taxon trees was between branch lengths estimated under an “uncorrected” maximum-likelihood model (hence referred to as UML) and branch lengths re-estimated under GFF. For 16-taxon datasets ( Figure 2 ), UML consistently underestimates branch lengths for tip and stem branches in the high-FYMINK clades B and D, and for the two internal branches closest to the root. UML also overestimates branch lengths for stem branches in the equal-frequency clades A and C ( Figure 2 ). GFF improves branch length estimation across all tip, stem and internal branches - with the tradeoff that there is a larger variance among branch length estimates for the internal branches ( Figure 2 ). These trends can also be observed in the root mean squared error (RMSE) calculations for the 16-taxon trees, with the largest error for GFF observed at the internal branches ( Table 2 ). View this table: View inline View popup Download powerpoint Table 2. Comparison of branch length estimation root mean squared error for 16-taxon trees simulated under LG+C20+Γwith GARP/FYMINK heterogeniety. Download figure Open in new tab Fig. 2. Performance of GFmix implementations across 100 16-taxon replicates simulated under LG+C20+Γ+ { GARP/FYMINK } . Comparisons are between observed estimations and expected estimation based on simulating values. Left: Estimation of branch lengths. Comparison made between estimates under LG+C20+Γmodel by IQ-TREE and estimates under LG+C20+Γ+ { GARP/FYMINK } model with maximum-likelihood optimization of all model and tree parameters by GFF. Right: Estimation of node GARP/FYMINK composition ratios. Comparison made between estimates by each GFmix implementation under LG+C20+Γ+ { GARP/FYMINK } model. Dots: mean value, errorbars: ± 1 σ , dotted line: observed and expected are identical. For the downsampled 4-taxon datasets, UML overestimates the branch lengths for the equal-frequency taxa A and C ( Figure S1 ). Estimation markedly improves for these taxa across replicates when branch lengths are re-optimized using GFF, alongside small improvements in estimation for the high-FYMINK taxa B and D ( Figure S1 ). As with the 16-taxon datasets, there is a tradeoff in improving tip branch length estimation that results in larger variance in length estimates for internal branches ( Figure S1 ). The larger degree of error observed in the 4-taxa RMSEs is likely a reflection of less information in the 4-taxa datasets relative to the 16-taxa datasets ( Table 3 ). View this table: View inline View popup Download powerpoint Table 3. Comparison of branch length estimation root mean squared error for downsampled 4-taxon trees. Estimation of branch GARP/FYMINK composition ratio We assessed the performance of each of the four GFmix implementations at modelling compositional shifts across heterogeneous 16-taxon trees. This was done by comparing the GARP/FYMINK frequency ratio implied by the simulating LG+C20+Γ+GFmix model (where 𝒢 /ℱ = GARP/FYMINK) with those estimated under the alternative GFmix implementations. Performance varied for each implementation. OGF consistently underestimated GARP/FYMINK composition across the equal-frequency nodes, and overestimated GARP/FYMINK composition for all nodes in the high-FYMINK clades B and D ( Figure 2 ). SSF displays improved performance for tip nodes but with substantially larger variance in estimates for stem, internal and root nodes ( Figure 2 ). GFG and GFF estimate GARP/FYMINK compositions better and more consistently across all node types, with larger variance estimates for deeper nodes ( Figure 2 ). These trends can also be observed in the RMSE calculations ( Table 4 ). OGF performs the worst for all node types, SSF performs best in tip composition estimation, whereas GFG and GFF exhibit consistently lower errors across the tree ( Table 4 ). View this table: View inline View popup Download powerpoint Table 4. Comparison of GARP/FYMINK frequency ratio estimation root mean squared error for 16-taxon trees simulated under LG+C20+Γwith GARP/FYMINK heterogeniety. Similar trends are observed in the downsampled 4-taxon datasets. OGF consistently gave biased branch GARP/FYMINK composition estimates across tip, internal and root nodes ( Figure S1 ). SSF performs well at estimating tip compositions, which suggests better performance with less taxa, but fares poorer with internal and root composition estimation ( Figure S1 ). GFG and GFF estimate GARP/FYMINK compositions better and more consistently across all node types, with GFF performing better for equal-frequency taxa A and C and both internal and root node compositions ( Figure S1 ). These trends can also be observed in the RMSE calculations ( Table 5 ). OGF performs worst for all node types, SSF performs best in tip composition estimation whereas GFP and GFF exhibit consistently lower errors across the tree ( Table 5 ). View this table: View inline View popup Download powerpoint Table 5. Comparison of branch length estimation root mean squared error for downsampled 4-taxon trees. Identifying compositional heterogeneity GFmix has been applied in cases of well-known amino acid compositional heterogeneity arising from GC-content variation ( Muñoz-Gómez et al. 2022 ; Williamson et al. 2025 ; Baker et al. 2025 ) or adaptation to hypersaline environments ( Baker et al. 2024 ). In the latter case, the amino acids most important to compositional heterogeneity were DE/IK rather than GARP/FYMNK. There may be other cases where underlying groups that are most important to compositional heterogeneity are either not known or not immediately apparent. In such cases, it would be desirable to identify the groups of amino acids directly from sequence data. We investigated how to identify these groups from data using our simulated datasets, which allows comparisons with truth, since GARP/FYMINK are known to be the heterogeneity-determining amino acids in the simulating model. We applied four methods of compositional heterogeneity identification - the binomial test of two proportions previously applied in conjuction with GFmix analysis ( Baker et al. 2024 ; Williamson et al. 2025 ; Baker et al. 2025 ) and a new optimization-based procedure using three different statistical criteria. Figure 3 outlines the workflow of each identification method, and further details are given in Materials and Methods . We tabulated the assignments of 𝒢, ℱ and 𝒪 for each method across our 16-taxon and 4-taxon datasets to assess how often each method correctly assigned 𝒢 /ℱ as GARP/FYMINK. Download figure Open in new tab Fig. 3. Workflow of 𝒢 /𝒪 /ℱ assignment from simulated sequence data. Left: Assignment using binomial test of two proportions. Given an input alignment, and two groups to which taxa are assigned, tabulate amino acid counts and calculate a Z -score for each amino acid across the two groups. Amino acids are then assigned into 𝒢, 𝒪 or ℱ on the basis of their Z -score. Right: Assignment using non-greedy optimization procedure. All possible 1-amino acid “swaps” are performed between 𝒢 / 𝒪 / ℱ, except any 𝒪 -to- ℱ swaps. Each new 𝒢 / 𝒪 / ℱ swap is assessed using one of three optimization criteria, and the best-performing swap is compared with the current 𝒢 /𝒪 /ℱassignment. If this swap improves the criterion, this swap is now the 𝒢 /𝒪 /ℱ assignment and the optimization procedure is repeated. If this swap does not improve the criterion, the optimization procedure is concluded and the current 𝒢 /𝒪 /ℱ assignment is considered optimized. Refer to Methods and Materials for further information. Comparing method accuracy The binomial test of two proportions is a simple onetime test for the equality of two proportions, in this case the proportion of each amino acid between two groups of taxa ( Baker et al. 2024 ). For our 16- and 4-taxon datasets, this test was performed between the equal-frequency taxa (A* and C*) and the high-FYMINK taxa (B* and D*). Although the method does generally assign GARP to 𝒢 and FYMINK to ℱ it frequently assigns additional amino acids to the two classes which should be assigned 𝒪 to ( Figures 4 and S2 ). Across almost all 16-taxon replicates, Leucine and Valine are alwaysassigned to 𝒢 and Serine is always assigned to ℱ ( Figure 4 ). For 29 of the 100 4-taxon replicates, the method assigns Tyrosine to 𝒪 instead of ℱ ( Figure S2 ). Download figure Open in new tab Fig. 4. Comparison of 𝒢 /𝒪 /ℱassignment for four compositional heterogeneity identification methods across 100 16-taxon alignments simulated under a LG+C20+Γ+ { GARP/FYMINK } model. Left: Assignments of 20 amino acids to𝒢, 𝒪andℱ, with tiles shaded by observed assignment and expected assignment of amino acids underneath the x-axis. Correct assignment - i.e. 𝒢= GARP, 𝒪= DCQEHLSTWV and ℱ= FYMINK - indicated by green ticks. Assignments ordered by from most to least correct. Right: Number of assignments by identification method across 100 16-taxon alignments. Blank cells indicate no assignment by that method. The χ 2 optimization procedure tries to select the 𝒢 and ℱ that minimize the χ 2 statistic for an alignment, using the assignments of the binomial test as a starting point. This method noticeably improves upon the binomial test in terms of closeness to 𝒢 /ℱ = GARP/FYMINK, but still fails to correctly assign 𝒢, 𝒪 and ℱ for any 16-taxon or 4-taxon replicate ( Figures 4 and S2 ). In the majority of both sets of replicates, Phenylalanine and Tyrosine are placed in 𝒪 instead of ℱ ( Figures 4 and S2 ). In smaller numbers of 16-taxon replicates, Histidine and/or Tryptophan is assigned to 𝒢 ( Figure 4 ). Leucine and Serine are assigned to 𝒢 and ℱrespectively in some of the 4-taxon replicates ( Figure S2 ). The final two optimization procedures use two of the GFmix implementations discussed in this study - SSF and OGF meaning that these optimization procedures incorporate model and tree information in their assessment. For SSF optimization, the procedures tries to select 𝒢 and ℱ which minimizes the sum-of-squares difference (SS) between observed and model amino acid frequencies for an alignment given a tree and model. For OGF optimization, the procedure tries to select 𝒢 and ℱ which maximizes the log-likelihood statistic estimated by the original GFmix model for an alignment, given a tree and underlying siteheterogeneous model ( Muñoz-Gómez et al. 2022 ). Both procedures show substantial improvement in 𝒢, 𝒪 and ℱ assignment over the other two methods ( Figures 4 and S2 ). OGF optimization assigns 𝒢/ ℱ = GARP/FYMINK in all 100 16-taxon replicates ( Figure 4 ). SSF optimization correctly assigns 𝒢/ ℱ = GARP/FYMINK in 51 of 100 replicates, and makes a single misassignment in another 25 16-taxon replicates ( Figure 4 ). For the 4-taxon replicates, OGF makes the correct assignment in 81 out of 100 replicates and SSF makes the correct assignment in 45 out of 100 replicates ( Figure S2 ). In another 36 replicates, SSF misassigns one amino acid ( Figure S2 ). Applying GFmix to real data We applied each GFmix implementation to a modified version of a real dataset of protein sequences derived from algal nuclear and nucleomorph genomes ( Novak et al. 2024 ). Nucleomorphs are vestigial nuclei found in membrane-delimited subcellular compartments that also contain plastids in Cryptomonada and Chlorarachniophyceae algae ( Ishida et al. 1999 ; Douglas et al. 2001 ; Moore and Archibald 2009 ). These compartments are relics of “secondary symbionts”; that is eukaryotic algal symbionts that took up residence within a different eukaryotic host and became integrated organelles. Nucleomorphs have undergone extensive genome reduction, and are thus AT-rich at the genomic level and FYMINKrich at the proteomic level ( Moore and Archibald 2009 ). Novak et al. (2024) constructed a dataset of 180 protein sequences derived from 54 genomes from Rhodophyta, Viridiplantae and four nucleomorphs from Cryptomonada. Using this dataset, Novak et al. (2024) placed Cryptomonad nucleomorphs as sister to extremophilic Cyanidiophytina at the base of Rhodophyta. We modified the Novak et al. (2024) dataset to include homologs from the genome of the Bigelowiella natans nucleomorph ( Gilson et al. 2006 ), and downsampled the dataset to 66 genes and 17,696 sites across 16 representative taxa - 9 Rhodophyta, 5 Viridiplantae and the nucleomorph genomes of Cryptomonas paramecium and B. natans (see Materials and Methods for further details). We tested two hypotheses with this dataset - one in which both nucleomorphs branch sister to their closest sampled relative from the literature ( Ishida et al. 1999 ; Douglas et al. 2001 ), and one in which the two nucleomorphs branch as sister taxa at the base of Rhodophyta as a result of long-branch attraction and shared amino acid composition bias. We designated these hypotheses as “NMs Apart” and “NMs Together”, respectively ( Figure 5 ). Download figure Open in new tab Fig. 5. Tree topologies tested for modified Novak et al. (2024) dataset. The modified dataset contains 16 taxa - 9 Rhodophyta (red), 5 Viridiplantae (blue) and two nucleomorph genomes (black). Left: “NMs Apart” tree, in which the two nucleomorphs branch sister with their closest sampled relatives from the literature. Right: “NMs Together” tree, in which the two nucleomorphs branch as sister taxa at the base of Rhodophyta. “NMs Apart” is considered the “correct” tree on the basis of biological consensus (see references in text), whereas “NMs Together” is artefact of long-branch attraction due to the FYMINK-rich composition of both nucleomorph proteomes. Root position represented by a dot. We first tested support for both trees under each GFmix implementation with the default 𝒢/𝒪 classes GARP/FYMINK ( Table 6 ). We found consistent support for the “NMs Together” hypothesis under both the branch-homogeneous LG+C20+Γmodel and two of the four GFmix implementations tested. The correct “NMs Apart” hypothesis was favoured under the implementation of GFmix — LG+C20+Γ+GFF — that uses full maximum like-lihood estimation and under LG+C20+Γ+SSF where and parameters are deter-mined by minimizing a sum of squares. Log-likelihoods improved by at least 1000 units for both hypotheses under each implementation of GFmix relative to LG+C20+Γ( Table 6 ). Calculating the Akaike and Bayesian information criteria for each model shows that LG+C20+Γ+GFF provides the best fit to data for each topology ( Tables S1 and S2 ). Despite favouring the correct tree, LG+C20+Γ+SSF exhibits poorer fit to the data than the other GFmix models ( Tables 6 , S1 and S2 ). View this table: View inline View popup Download powerpoint Table 6. Phylogenetic analysis of two topologies for modified Novak et al. (2024) dataset. “NMs Apart” : nucleomorphs branch on either side of Rhodophyta/Viridiplantae split, “NMs Together” : nucleomorphs branch together as sister taxa at base of Rhodophyta. See Fig. 5 for more details. LG+C20+Γfitting was performed using IQ-Tree 3 with mixture weight optimization. All other fittings performed using various implementations of GFmix with default 𝒢/ℱ (GARP/FYMINK) or 𝒢/ℱchosen by OGF optimization procedure. We repeated this analysis with 𝒢/ℱ classes selected by the optimization procedure previously outlined. We selected OGF as the optimization criterion based on the results of our simulation analysis, and this selected 𝒢/ℱ= GARPVMTHQ/FYCINK. We observed a similar trend to the first analysis - the correct “NMs Apart” hypothesis is only favoured under LG+C20+G+GFF ( Table 6 ). This suggests that optimizing the tree and model parame-ters - e.g. α , mixture weights and branch lengths - in the context of a branch-heterogeneous model is as important as accounting for branch heterogeneity itself. Log-likelihoods improved for both hypotheses with 𝒢/ℱ classes optimized over the default GARP/FYMINK assignments, with the exception of LG+C20+Γ+SSF which produces worse likelihoods than the base model. AIC and BIC also improved with 𝒢/ℱ classes optimized over the default GARP/FYMINK assignments, with the exception of LG+C20+Γ+SSF which produces worse fit than the base model ( Tables S1 and S2 ). This poor performance is likely due to LG+C20+Γ+SSF assigning extreme b e values to both tip branches and deep branches in both trees (see Discussion below for further discussion of this point). D iscussion How well does GFmix perform on simulated data? We simulated site-heterogeneous alignments with branch-heterogeneous GARP/FYMINK compositions, under the assumptions of the GFmix model, to assess the performance of each GFmix implementation discussed in this study. Given that we know the true simulating values, how well did each implementation perform at estimating branch lengths and node-specific GARP/FYMINK composition? What effect does downsampling from 16 to 4 taxa have on GFmix performance? Only the most complex GFmix implementation assessed here — GFF — allows for re-estimation of branch lengths. Thus we compared the branch lengths estimated under the LG+C20+Γ+GFF model with those originally estimated by IQ-TREE under the branch-homogeneous LG+C20+Γmodel (or UML for “uncorrected ML” in our notation). Notably, the UML model had large biases in branch length estimates across most branches in heterogeneous trees. This includes underestimating the lengths of branches undergoing a shift from equal-frequency to high-FYMINK composition, and overestimating the lengths of adjacent branches which are more compositionally-homogeneous. This is an error which LG+C20+Γ+GFF largely corrects, although this results in larger estimation variance at branches close to the root. Assessment of node-specific GARP/FYMINK composition estimation was performed by comparing the expected GARP/FYMINK frequency ratio at each node under the simulating model with those estimated at each node by each GFmix implementation. OGF shows a substantial bias in estimation. SSF shows little bias in estimation but often has the largest standard deviation. Because it reduces the data to tip frequencies prior to estimation, some increase in variance was to be expected relative to ML methods. GFP and GFF perform well across all node types, with larger estimation variance at branches close to the root. These two methods are comparable in their biases and variance of frequency estimation. They do, however, tend to have biases that are larger than SSF. The reasons for this is not clear. The performance trends we observe for 16-taxon datasets largely hold for 4-taxon datasets. The noteworthy exception is the larger variance in composition estimation at internal nodes for 4-taxon datasets relative to the 16-taxa datasets. This can be explained by the relative lack of information available for parameter estimation for models in 4-taxon datasets relative to 16-taxon datasets. Overall, LG+C20+Γ+GFF outperforms LG+C20+Γin estimating branch lengths for heterogeneous trees and the two ML-based GFmix implementations (GFP and GFF) give better root mean squared errors of node-specific composition estimation than the non-ML GFmix implementations (OGF and SSF). However, because of their computational costs, these approaches, particularly GFF, may not always be feasible for taxon-rich or site-rich datasets. The “simpler” OGF or SSF implementations may have utility in these cases. How well can we identify compositional heterogeneity? Although originally implemented to model GARP/FYMINK heterogeneity, the GFmix implementations assessed in this study are capable of modeling user-defined 𝒢/ℱ heterogeneity. For some datasets, the optimal heterogeneity to model may be known a priori -e.g. GARP/FYMINK as a result of GC-content variation - but a method of identifying the groups of amino acids most important to heterogeneity is desirable. We devised four methods for identifying compositional heterogeneity directly from sequence data and assessed their performance in identifying 𝒢/ℱ = GARP/FYMINK from simulated site-and-branch-heterogeneous data. The least computationally-intensive methods — the one-time binomial test of two proportions and 𝒢/ℱclass optimization using the χ 2 test statistic — are also the least accurate in identifying the true heterogeneity from simulated data. The more computationally-intensive, GFmix-based methods —𝒢/ℱclass optimization using either the sum-of-squares difference between observed and model amino acid frequencies estimated by SSF, or calculating the likelihood of the tree under OGF with𝒢/ℱ— can use additional information from the model or tree for 𝒢/ℱ assignment. Perhaps unsuprisingly, these methods perform far better at identifying true heterogeneity from simulated data. Comparing OGF and SSF, OGF had a far larger frequency of correct 𝒢/ℱ identification. It is worth noting, however, that for SSF the most frequent mistake was not correctly assigning a single amino acid, either C or W, both of which are relatively rare in terms of overall frequencies. The nature of common mistakes in 𝒢/ℱ estimation led us to recognize some unexpected behaviour of the GFmix model. GFmix models sequence evolution along a branch according to a Markov process with stationary amino acid frequencies of mixture classes determined by branch-specific 𝒢, ℱ and rate multipliers associated with these classes. Although the model along a branch is stationary, the marginal frequencies of amino acids at the ancestral node for the branch will not be the same as the stationary frequencies for the branch except in special cases such as when the 𝒢 and ℱ rate multipliers are constant throughout the tree. Consequently, the marginal frequencies change continuously along a branch from what they were at the ancestral node and tend toward the stationary frequencies for the branch as it gets longer. However, since branch lengths are finite, frequencies at the end node of the branch can differ substantially from the stationary frequencies that would be approximated with very large branch lengths. Figure 6 illustrates how amino acid frequencies averaged over mixture classes behave under the default LG+C20+Γ+GFmix model, as the GARP/FYMINK frequency ratio shifts from those expected under the standard LG+C20+Γmodel to 0.1 times that ratio, as a function of branch length. As expected, the frequencies of GARP decrease as the frequencies of FYMINK increase, with all frequencies stabilizing as branch length increases. Download figure Open in new tab Fig. 6. Average amino acid frequencies as a function of branch length under LG+C20+Γ+GARP/FYMINK. As the GARP/FYMINK frequency ratio decreases from those of the LG+C20+Γmodel to 0.1 times that ratio – meaning π ( GARP ) < π ( FY MINK ) – the frequencies of GARP decrease and the frequencies of FYMINK increase. This effect stabilizes as branch lengths increase. At shorter branch lengths ( < 2.5), amino acids in the homogeneous 𝒪class also increase or decrease in frequency before returning to equilibrium at longer branch lengths. The exceptions (in bold ) are Leucine, Serine and Valine. The behaviour of amino acids under the GFmix model at shorter branch lengths may explain misassignments of𝒢, 𝒪and ℱfor some identification methods. The behaviour of 𝒪 amino acids in Figure 6 is more surprising. Many 𝒪 amino acids also increase/decrease in frequency at shorter branch lengths, before returning to their original frequency at branch lengths above 2.5. Leucine, Serine and Valine do not follow this behaviour at longer branch lengths: L and V remain at decreased frequencies and S remains at an increased frequency ( Figure 6 ). This is mixture-model-dependent behaviour — C20 includes classes in which, across classes, the frequency of Leucine and Valine are positively correlated with the the frequency of some FYMINK amino acids (particularly Phenylalanine and Isoleucine). Similarly, Serine usage is postively correlated with the use of GARP amino acids (particularly Alanine). Thus, when the frequencies of FYMINK amino acids in these classes are shifted this causes an inverse shift of L and V to compensate. Similarly, decreasing GARP causes an increase in S. This behaviour likely explains the consistent misassignment of L, S and V under more simplistic class membership optimization criteria like the binomial approach that chooses 𝒢/ℱ based on differences in observed amino acids across groups without adjusting for branch-lengths. The behaviour of the other 𝒪 amino acids, which need not be close to limiting (large branch length) values, particularly at branch lengths similar to those used to simulate heterogeneous data in this study, is also a likely cause of misassignment. A possible explanation is that as GARP amino acids shift towards increasing usage of FYMINK, these transitions occur via intermediates in the 𝒪 class, leading to temporal increases (and compensatory increases) in some frequencies. In terms of required runtime, the tree/model-dependent methods of identification are slower than the χ 2 method as is to be expected. For 4-taxon datasets SSF is faster than OGF. This is due to the nature of the sum-of-squares calculation, which is faster for smaller datasets than the standard likelihood calculation. SSF, however, was slower for 16-taxon data sets. Likelihood calculation is slower than sum-of-squares calculation but, for each bin choice, OGF calculates a single likelihood using b e parameters determined by a hierarchical clustering procedure, SSF chooses b e parameters to minimize a sum of squares which can require many sum-of-squares and derivative calculations. It is possible that many parameter combinations give similar sum-of-squares, thus leading to flatter surfaces and longer optimization times. Overall, optimization using likelihood via OGF is the best performing method of identifying amino acids belonging to 𝒢 and ℱ classes from sequence data and has acceptable runtime and scaling. The other methods assessed here may have some utility for less heterogeneous datasets or smaller datasets, but should be used with some caveats regarding their performance and/or compared with a priori knowledge of compositional heterogeneity for a given dataset if available. How well does GFmix perform on real data? GFmix is intended for use in phylogenetic analysis where compositional heterogeneity is suspected. Artefacts arising from compositional heterogeneity have been suspected in a number of deep-time phylogenetic analyses ( Williams et al. 2021 ). However, artefacts can occur in any dataset with an appreciably diverse sampling of taxa ( Muñoz-Gómez et al. 2019 ) or particularly divergent protein sequences ( Foster and Hickey 1999 ). We applied each GFmix implementation to a modified version of a dataset comprised of protein sequences derived from algal nuclear and nucleomorph genomes, adding data from the B. natans nucleomorph genome and downsampling to 16 representative taxa including the nucleomorphs of C. paramecium and B. natans ( Gilson et al. 2006 ; Novak et al. 2024 ). The Novak et al. (2024) dataset was selected as a basis for our real data analysis as nucleomorph genomes are highly-reduced and thus FYMINK-rich at the proteomic level ( Moore and Archibald 2009 ), whereas most of the nuclear genomes retained in the downsampled dataset were comparatively homogeneous with respect to GARP and FYMINK. We tested two hypotheses with this dataset - one in which the two sampled nucelomorphs in this dataset branch with their closest sampled relatives according to the literature ( Ishida et al. 1999 ; Douglas et al. 2001 ), and one in which the two sampled nucleomorphs branch as a sister taxa at the base of Rhodophyta due to long-branch attraction. We found that only full ML estimation with GFmix (LG+C20+Γ+GFF) consistently favoured the independent placement of the two nucleomorphs within Rhodophyta and Viridiplantae respectively. This result is consistent with the accepted biological consensus for the origin of nucleomorphs in Cryptomonada and Chlorarachniophyceae ( Moore and Archibald 2009 ). The log-likelihoods evaluated for each tree and the fit to data according to AIC and BIC improve under each GFmix implementation relative to the base LG+C20+Γmixture model, with one exception, and substantially so for LG+C20+Γ+GFP and LG+C20+Γ+GFF. This suggests that directly estimating branch-specific amino acid composition parameters improves log-likelihood evaluation over the indirect estimation as used in the original implementation of GFmix tested here as LG+C20+Γ+OGF ( Muñoz-Gómez et al. 2022 ; Baker et al. 2024 ). The further improvement of LG+C20+Γ+GFF over LG+C20+Γ+GFP in favouring the “correct” tree is evidence of the utility of re-estimating additional parameters such as branch lengths and mixture weights in the context of a branch-heterogeneous evolutionary process. Finally, we see improved likelihoods and fit to data when using 𝒢 and ℱ classes determined by our optimization procedure over the default GARP/FYMINK assignments. This further supports the validity of inferring compositional heterogeneity from sequence data as we have attempted with our optimization procedure. The exception to this trend is the unexpected behaviour of the LG+C20+Γ+SSF model. When 𝒢/ℱ = GARP/FYMINK, the model favours the correct tree but when 𝒢/ℱ = GARPVMTHQ/FYCINK the model favours the incorrect tree with a greater likelihood difference than the base model. Furthermore, when = GARPVMTHQ/FYCINK the model also produces worse likelihoods than the base model for both trees. The fit of the model to data is the worst of each GFmix model when 𝒢/ℱ= GARP/FYMINK, and worse than the base model when 𝒢/ℱ= GARPVMTHQ/FYCINK. Examining the b e parameters estimated by the model ( Supplementary Information, Supplementary Material ), it appears that this may be a result of the model estimating extremely large b e parameters at the tips and extremely small b e parameters along internal branches in both trees. The SSF implementation of the GFmix model attempts to minimize the differences between the amino acid frequencies observed across all terminal nodes with the stationary frequencies implied by the base model. As such, it is sensitive to branch lengths estimated under the base model. As many of the internal branches of the tree have small estimated branch lengths, extreme b e values are sometimes assigned by SSF in order to change frequencies over edges to best fit taxon frequencies. It is possible that penalizing large or small b e values could lead to better estimation. At present, we would caution that although use of the SSF will give frequencies that better match observed frequencies, it can also lead to poor tree likelihoods. Other considerations and future directions There are limitations to the GFmix model which apply to all implementations tested here. For example, our software implementation of GFmix does not perform tree searching. A set of candidate trees for likelihood evaluations must be determined by the user using IQ-TREE or other means. Additionally, GFmix assumes rate heterogeneity categories are derived from the Gamma (Γ) distribution ( Yang 1994 ) and thus is not currently compatible with the FreeRate model of rate heterogeneity ( Soubrier et al. 2012 ). Finally, it cannot yet be used in combination with protein evolution models which accommodate other complex features, such as the GHOST model of heterotachy ( Crotty et al. 2019 ) or MAST for non-treelike evolution ( Wong et al. 2024 ). GFmix was originally implemented in Muñoz-Gómez et al. (2022) to model amino acid composition biases arising from genome-wide GC-content variation amongst species which drives changes in protein sequence composition such as GARP vs. FYMINK alterations. However, the model can also be applied to cases of selection biases in preferred amino acids as shown in Baker et al. (2024) . GFmix assumes that whatever amino acid composition bias exists for a given taxa (or group), it applies to the entire protein sequence(s) under consideration. For certain trees of interest, compositional heterogeneity may involve several different sets of amino acids in different parts of the tree due to lineage-specific adaptations. This can be observed across Archaea, where different amino acid groups are associated with lineages that have undergone genome reduction (resulting in GARP/FYMINK changes), adaptation to hypersaline environments (leading to DE/IK changes) or adaptation to extremely hot environments (the proteomic fraction of IVYWREL or ILVWYGERKP amino acids) ( Baker et al. 2024 ; Baker et al. 2025 ; Eme et al. 2023 ; Zeldovich et al. 2007 ). For such trees of interest, identifying the optimal composition bias to model using the methods outlined here may be difficult. Further curation of data, such as subsampling trees or partitioning based on shared amino acid composition or other information may be advisable. Another alternative, as used in Williamson et al. (2025) , is to partition sites on the basis of compositional similarity and to apply the GFmix model separately to each partition. Future elaborations of the GFmix approach could investigate the use of different 𝒢 and ℱ amino acid groups for different parts of a tree. C onclusions Accurate inference of deep phylogenies requires models that can accommodate heterogeneous amino acid composition at both the site and branch levels. In this study, we presented new implementations of the site-and-branch-heterogeneous GFmix model. The simpler implementations estimate branch-specific model parameters independently and then use the GFmix model for likelihood inference with these parameters. More complex implementations of the model can estimate some or all of the model and tree parameters in a maximum-likelihood framework. We assessed the performance of these implementations when applied to data simulated under GARP/FYMINK compositional heterogeneity and to real data subject to similar composition biases. Applied to simulated data, we showed that our newer implementations of GFmix perform substantially better at estimating branch-specific amino acid composition changes, and the most extensive implementation is capable of correcting branch length estimation errors made by branch-homogeneous profile mixture models. We also demonstrate that GFmix can be used to identify the underlying amino acid groups most responsible for compositional heterogeneity from sequence data. Applied to an algal nuclear and nucleomorph dataset derived from Novak et al. (2024) , GFmix improves likelihoods and fit to data over branch-homogeneous models and does so substantially for the most extensive implementations. We find that the most extensive implementation supports the independent origins of the nucleomorphs of Bigelowiella natans and Cryptomonas paramecium , in Viridiplantae and Rhodophyta respecitvely, over the artefactual long-branch attraction tree placing the two nucleomorphs as sister taxa. This demonstrates the efficacy of the GFmix model in handling known compositional heterogeneity artefacts within real phylogenetic data. This study demonstrates the utility of modelling both site- and branch-level heterogeneity in phylogenetics. The GFmix model should be considered for use by researchers in cases where compositional heterogeneity may adversely affect phylogenetic tree inference. A uthor contributions A.J.R., E.S. and C.G.P.M. designed and led the study. E.S. designed and implemented all models and parameter estimation methods. C.G.P.M. implemented the analysis code and performed all simulation, downsampling, phylogenetics and data analysis. C.G.P.M. drafted the initial manuscript and all authors contributed to editing the paper. S upplementary M aterial Additional information for generating simulated site-and-branch-heterogeneous sequence data and additional assessments of compositional heterogeneity identification methods is provided in Supplementary Information . Tree parameter estimations for real dataset analysis is provided in Supplementary Material . Simulated and real datasets, and the code used to perform the analysis and generate figures provided in this study is provided in Supplementary Data available on Dryad: Link TBD. A cknowledgements This work and C.G.P.M. were supported by the Moore-Simons Project Call on the Origin of the Eukaryotic Cell, Simons Foundation Grant 735923LPI ( https://doi.org/10.46714/735923LPI ) and by NSERC Discovery Grants awarded to A.J.R. and E.S.. The authors thank Trong Hahn Ly and Minh Bui (Australian National University) for their guidance on simulating site-and-branch-heterogeneous data with AliSim, and Ryo Harada (Dalhousie University) for identifying and adding Bigelowiella natans nucleomorph protein sequence data to the algal dataset analyzed in this study. C.G.P.M. thanks Hector Baños (Dalhousie University and California State University, San Bernardino), Kelsey Williamson (Dalhousie University) and Brittany Baker (Université Paris-Saclay and Institut Pasteur) for advice and feedback in other projects which informed this work. Funder Information Declared Simons Foundation , 735923LPI Natural Sciences and Engineering Research Council, https://ror.org/01h531d29 , Discovery Grants R eferences ↵ Acosta , S. , M. Carela , A. Garcia-Gonzalez , M. Gines , L. Vicens , R. Cruet , and S. E. Massey ( 2015 ). DNA Repair Is Associated with Information Content in Bacteria, Archaea, and DNA Viruses . Journal of Heredity 106 , pp. 644 – 659 . doi: 10.1093/jhered/esv055 . OpenUrl CrossRef PubMed ↵ Altschul , S. F. , W. Gish , W. Miller , E. W. Myers , and D. J. Lipman ( 1990 ). Basic local alignment search tool . Journal of Molecular Biology 215 , 403 – 410 . doi: 10.1016/s0022-2836(05)80360-2 . OpenUrl CrossRef PubMed Web of Science ↵ Baker , B. A. , A. Gutiérrez-Preciado , Rodríguez del Río , C. G.P. McCarthy , P. López-García , J. Huerta-Cepas , E. Susko , A. J. Roger , L. Eme , and D. Moreira ( 2024 ). Expanded phylogeny of extremely halophilic archaea shows multiple independent adaptations to hypersaline environments . Nature Microbiology 9 , 964 – 975 . doi: 10.1038/s41564-024-01647-4 . OpenUrl CrossRef PubMed ↵ Baker , B. A. , C. G. P. McCarthy , P. López-García , R. B. Leroy , E. Susko , A. J. Roger , L. Eme , and D. Moreira ( 2025 ). Phylogenomic analyses indicate the archaeal superphylum DPANN originated from free-living euryarchaeal-like ancestors . Nature Microbiology 10 , 1593 – 1604 . doi: 10.1038/s41564-025-02024-5 . OpenUrl CrossRef PubMed ↵ Baños , H. , E. Susko , and A. J. Roger ( 2023 ). Is Over-parameterization a Problem for Profile Mixture Models? Systematic Biology 73 , 53 – 75 . doi: 10.1093/sysbio/syad063 . OpenUrl CrossRef ↵ Blanquart , S. and N. Lartillot ( 2006 ). A Bayesian Compound Stochastic Process for Modeling Nonstationary and Nonhomogeneous Sequence Evolution . Molecular Biology and Evolution 23 , pp. 2058 – 2071 . doi: 10.1093/molbev/msl091 . OpenUrl CrossRef PubMed Web of Science ↵ Blanquart , S. and N. Lartillot ( 2008 ). A Site- and Time-Heterogeneous Model of Amino Acid Replacement . Molecular Biology and Evolution 25 , pp. 842 – 858 . doi: 10.1093/molbev/msn018 . OpenUrl CrossRef PubMed Web of Science ↵ Brinkmann , H. , M. van der Giezen , Y. Zhou , G. P. de Raucourt , and H. Philippe ( 2005 ). An Empirical Assessment of Long-Branch Attraction Artefacts in Deep Eukaryotic Phylogenomics . Systematic Biology 54 , pp. 743 – 757 . doi: 10.1080/10635150500234609 . OpenUrl CrossRef PubMed Web of Science ↵ S. Smith Crotty , S. M. , B. Q. Minh , N. G. Bean , B. R. Holland , J. Tuke , L. S. Jermiin , and A. V. Haeseler ( 2019 ). GHOST: Recovering Historical Signal from Heterotachously Evolved Sequence Alignments . Systematic Biology . Ed. by S. Smith . doi: 10.1093/sysbio/syz051 . OpenUrl CrossRef PubMed ↵ Douglas , S. , S. Zauner , M. Fraunholz , M. Beaton , S. Penny , L.-T. Deng , X. Wu , M. Reith , T. Cavalier-Smith , and U.-G. Maier ( 2001 ). The highly reduced genome of an enslaved algal nucleus . Nature 410 , 1091 – 1096 . doi: 10.1038/35074092 . OpenUrl CrossRef PubMed Web of Science ↵ Embley , M. , M. v. der Giezen , D. S. Horner , P. L. Dyal , and P. Foster ( 2003 ). Mitochondria and hydrogeno-somes are two forms of the same fundamental organelle . Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 358 , 191 – 203 . doi: 10.1098/rstb.2002.1190 . OpenUrl CrossRef PubMed Web of Science ↵ Eme , L. , D. Tamarit , E. F. Caceres , C. W. Stairs , V. De Anda , M. E. Schön , K. W. Seitz , N. Dombrowski , W. H. Lewis , F. Homa , J. H. Saw , J. Lombard , T. Nunoura , W.-J. Li , Z.-S. Hua , L.-X. Chen , J. F. Banfield , E. S. John , A.-L. Reysenbach , M. B. Stott , A. Schramm , K. U. Kjeldsen , A. P. Teske , B. J. Baker , and T. J. G. Ettema ( 2023 ). Inference and reconstruction of the heimdallarchaeial ancestry of eukaryotes . Nature 618 , 992 – 999 . doi: 10.1038/s41586-023-06186-2 . OpenUrl CrossRef PubMed ↵ Ettema , T. J. and S. G. Andersson ( 2009 ). The -proteobacteria: the Darwin finches of the bacterial world . Biology Letters 5 , 429 – 432 . doi: 10.1098/rsbl.2008.0793 . OpenUrl CrossRef PubMed ↵ Fan , L. , D. Wu , V. Goremykin , J. Xiao , Y. Xu , S. Garg , C. Zhang , W. F. Martin , and R. Zhu ( 2020 ). Phylogenetic analyses with systematic taxon sampling show that mitochondria branch within Al-phaproteobacteria . Nature Ecology and Evolution 4 , 1213 – 1219 . doi: 10.1038/s41559-020-1239-x . OpenUrl CrossRef PubMed ↵ Felsenstein , J. ( 1973 ). Maximum Likelihood and Minimum-Steps Methods for Estimating Evolutionary Trees from Data on Discrete Characters . Systematic Biology 22 , 240 – 249 . doi: 10.1093/sysbio/22.3.240 . OpenUrl CrossRef ↵ Felsenstein , J. ( 1978 ). Cases in which Parsimony or Compatibility Methods will be Positively Misleading . Systematic Biology 27 , pp. 401 – 410 . doi: 10.1093/sysbio/27.4.401 . OpenUrl CrossRef ↵ Feuda , R. , M. Dohrmann , W. Pett , H. Philippe , O. Rota-Stabelli , N. Lartillot , G. Wörheide , and D. Pisani ( 2017 ). Improved Modeling of Compositional Heterogeneity Supports Sponges as Sister to All Other Animals . Current Biology 27 , 3864 – 3870.e4 . doi: 10.1016/j.cub.2017.11.008 . OpenUrl CrossRef PubMed ↵ Foster , P. G. ( 2004 ). Modeling Compositional Heterogeneity . Systematic Biology 53 , pp. 485 – 495 . doi: 10.1080/10635150490445779 . OpenUrl CrossRef PubMed Web of Science ↵ Foster , P. G. , C. J. Cox , and T. M. Embley ( 2009 ). The primary divisions of life: a phylogenomic approach employing composition-heterogeneous methods . Philosophical Transactions of the Royal Society B: Biological Sciences 364 , pp. 2197 – 2207 . doi: 10.1098/rstb.2009.0034 . OpenUrl CrossRef PubMed ↵ Foster , P. G. and D. A. Hickey ( 1999 ). Compositional Bias May Affect Both DNA-Based and Protein-Based Phylogenetic Reconstructions . Journal of Molecular Evolution 48 , 284 – 290 . doi: 10.1007/pl00006471 . OpenUrl CrossRef PubMed Web of Science ↵ Foster , P. G. , D. Schrempf , G.J. Szöllăsi , T. A. Williams , C. J. Cox , and T. M. Embley ( 2022 ). Recoding Amino Acids to a Reduced Alphabet may Increase or Decrease Phylogenetic Accuracy . Systematic Biology 6 72 , pp. 723 – 737 . doi: 10.1093/sysbio/syac042 . OpenUrl CrossRef ↵ Francis , W. R. and D. E. Canfield ( 2020 ). Very few sites can reshape the inferred phylogenetic tree . PeerJ 8 , e8865 . doi: 10.7717/peerj.8865 . OpenUrl CrossRef PubMed ↵ Gilson , P. R. , V. Su , C. H. Slamovits , M. E. Reith , P. J. Keeling , and G. I. McFadden ( 2006 ). Complete nucleotide sequence of the chlorarachniophyte nucleomorph: Nature’s smallest nucleus . Proceedings of the National Academy of Sciences 103 , pp. 9566 – 9571 . doi: 10.1073/pnas.0600707103 . OpenUrl Abstract / FREE Full Text ↵ Goldstein , R. A. ( 2008 ). The structure of protein evolution and the evolution of protein structure . Current Opinion in Structural Biology 18 . Theory and simulation /Macromolecular assemblages, pp. 170 – 177 . doi: 10.1016/j.sbi.2008.01.006 . OpenUrl CrossRef PubMed Web of Science ↵ Graybeal , A. ( 1998 ). Is It Better to Add Taxa or Characters to a Difficult Phylogenetic Problem? Systematic Biology 47 , pp. 9 – 17 . doi: 10.1080/106351598260996 . OpenUrl CrossRef GeoRef PubMed Web of Science ↵ Groussin , M. , B. Boussau , and M. Gouy ( 2013 ). A Branch-Heterogeneous Model of Protein Evolution for Efficient Inference of Ancestral Sequences . Systematic Biology 62 , pp. 523 – 538 . doi: 10.1093/sysbio/syt016 . OpenUrl CrossRef PubMed ↵ Halpern , A. L. and W. J. Bruno ( 1998 ). Evolutionary distances for protein-coding sequences: modeling site-specific residue frequencies . Molecular Biology and Evolution 15 , pp. 910 – 917 . doi: 10.1093/oxfordjournals.molbev.a025995 . OpenUrl CrossRef PubMed Web of Science ↵ Hernandez , A. M. and J. F. Ryan ( 2021 ). Six-State Amino Acid Recoding is not an Effective Strategy to Offset Compositional Heterogeneity and Saturation in Phylogenetic Analyses . Systematic Biology 70 , pp. 1200 – 1212 . doi: 10.1093/sysbio/syab027 . OpenUrl CrossRef PubMed ↵ Hershberg , R. and D. A. Petrov ( 2010 ). Evidence That Mutation Is Universally Biased towards AT in Bacteria . PLoS Genetics 6 , e1001115 . doi: 10.1371/journal.pgen.1001115 . OpenUrl CrossRef PubMed ↵ Ishida , K , B. Green , and T Cavalier-Smith ( 1999 ). Diversification of a Chimaeric Algal Group, the Chlo-rarachniophytes: Phylogeny of Nuclear and Nucleomorph Small-Subunit rRNA Genes . Molecular Biology and Evolution 16 , pp. 321 – 321 . doi: 10.1093/oxfordjournals.molbev.a026113 . OpenUrl CrossRef Web of Science ↵ Katoh , K. and D. M. Standley ( 2013 ). MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability . Molecular Biology and Evolution 30 , 772 – 780 . doi: 10.1093/molbev/mst010 . OpenUrl CrossRef PubMed Web of Science ↵ Lartillot , N. , T. Lepage , and S. Blanquart ( 2009 ). PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating . Bioinformatics 25 , pp. 2286 – 2288 . doi: 10.1093/bioinformatics/btp368 . OpenUrl CrossRef PubMed Web of Science ↵ Lartillot , N. and H. Philippe ( 2004 ). A Bayesian Mixture Model for Across-Site Heterogeneities in the Amino-Acid Replacement Process . Molecular Biology and Evolution 21 , pp. 1095 – 1109 . doi: 10.1093/molbev/msh112 . OpenUrl CrossRef PubMed Web of Science ↵ Le , S. Q. , C. C. Dang , and O. Gascuel ( 2012 ). Modeling Protein Evolution with Several Amino Acid Replacement Matrices Depending on Site Rates . Molecular Biology and Evolution 29 , pp. 2921 – 2936 . doi: 10.1093/molbev/mss112 . OpenUrl CrossRef PubMed Web of Science ↵ Le , S. Q. and O. Gascuel ( 2008 ). An Improved General Amino Acid Replacement Matrix . Molecular Biology and Evolution 25 , pp. 1307 – 1320 . doi: 10.1093/molbev/msn067 . OpenUrl CrossRef PubMed Web of Science ↵ Lex , A. , N. Gehlenborg , H. Strobelt , R. Vuillemot , and H. Pfister ( 2014 ). UpSet: Visualization of Intersecting Sets . IEEE Transactions on Visualization and Computer Graphics 20 , pp. 1983 – 1992 . doi: 10.1109/TVCG.2014.2346248 . OpenUrl CrossRef PubMed ↵ Li , Y. , X.-X. Shen , B. Evans , C. W. Dunn , and A. Rokas ( 2021 ). Rooting the Animal Tree of Life . Molecular Biology and Evolution 38 , pp. 4322 – 4333 . doi: 10.1093/molbev/msab170 . OpenUrl CrossRef ↵ Ly-Trong , N. , S. Naser-Khdour , R. Lanfear , and B. Q. Minh ( 2022 ). AliSim: A Fast and Versatile Phylogenetic Sequence Simulator for the Genomic Era . Molecular Biology and Evolution 39 , msac092 . doi: 10.1093/molbev/msac092 . OpenUrl CrossRef PubMed ↵ Martijn , J. , J. Vosseberg , L. Guy , P. Offre , and T. J. G. Ettema ( 2018 ). Deep mitochondrial origin outside the sampled alphaproteobacteria . Nature 557 , 101 – 105 . doi: 10.1038/s41586-018-0059-5 . OpenUrl CrossRef PubMed ↵ Minh , B. Q. , H. A. Schmidt , O. Chernomor , D. Schrempf , M. D. Woodhams , A. von Haeseler , and R. Lanfear ( 2020 ). IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era . Molecular Biology and Evolution 37 , pp. 1530 – 1534 . doi: 10.1093/molbev/msaa015 . OpenUrl CrossRef PubMed ↵ Moore , C. E. and J. M. Archibald ( 2009 ). Nucleomorph Genomes . Annual Review of Genetics 43 , 251 – 264 . doi: 10.1146/annurev-genet-102108-134809 . OpenUrl CrossRef PubMed Web of Science ↵ Muñoz-Gómez , S. A. , S. Hess , G. Burger , B. F. Lang , E. Susko , C. H. Slamovits , and A. J. Roger ( 2019 ). An updated phylogeny of the Alphaproteobacteria reveals that the parasitic Rickettsiales and Holosporales have independent origins . eLife 8 , e42535 . doi: 10.7554/eLife.42535 . OpenUrl CrossRef PubMed ↵ Muñoz-Gómez , S. A. , E. Susko , K. Williamson , L. Eme , C. H. Slamovits , D. Moreira , P. López-García , and A. J. Roger ( 2022 ). Site-and-branch-heterogeneous analyses of an expanded dataset favour mitochondria as sister to known Alphaproteobacteria . Nature Ecology and Evolution 6 , 253 – 262 . doi: 10.1038/s41559-021-01638-2 . OpenUrl CrossRef PubMed ↵ Novak , L. V. F. , S.A. Muñoz-Gómez , F. van Beveren , M. Ciobanu , L. Eme , P. López-García , and D. Moreira ( 2024 ). Nucleomorph phylogenomics suggests a deep and ancient origin of cryptophyte plastids within Rhodophyta . bioRxiv . doi: 10.1101/2024.03.10.584144 . OpenUrl Abstract / FREE Full Text ↵ Philippe , H. , H. Brinkmann , D. V. Lavrov , D. T. J. Littlewood , M. Manuel , G. Wörheide , and D. Baurain ( 2011 ). Resolving Difficult Phylogenetic Questions: Why More Sequences Are Not Enough . PLOS Biology 9 , pp. 1 – 10 . doi: 10.1371/journal.pbio.1000602 . OpenUrl CrossRef ↵ Puttick , M. N. , J. L. Morris , T. A. Williams , C. J. Cox , D. Edwards , P. Kenrick , S. Pressel , C. H. Wellman , H. Schneider , D. Pisani , and P. C. Donoghue ( 2018 ). The Interrelationships of Land Plants and the Nature of the Ancestral Embryophyte . Current Biology 28 , 733 – 745.e2 . doi: 10.1016/j.cub.2018.01.063 . OpenUrl CrossRef PubMed ↵ R Core Team ( 2025 ). R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing . Vienna, Austria . ↵ Sayers , E. W. , J. Beck , E. E. Bolton , J. R. Brister , J. Chan , R. Connor , M. Feldgarden , A. M. Fine , K. Funk , J. Hoffman , S. Kannan , C. Kelly , W. Klimke , S. Kim , S. Lathrop , A. Marchler-Bauer , T. D. Murphy , C. O’Sullivan , E. Schmieder , Y. Skripchenko , A. Stine , F. Thibaud-Nissen , J. Wang , J. Ye , E. Zellers , V. A. Schneider , and K. D. Pruitt ( 2024 ). Database resources of the National Center for Biotechnology Information in 2025 . Nucleic Acids Research 53 , D20 – D29 . doi: 10.1093/nar/gkae979 . OpenUrl CrossRef ↵ Schrempf , D. , N. Lartillot , and G. Szöllăsi ( 2020 ). Scalable Empirical Mixture Models That Account for Across-Site Compositional Heterogeneity . Molecular Biology and Evolution 37 , pp. 3616 – 3631 . doi: 10.1093/molbev/msaa145 . OpenUrl CrossRef ↵ Si Quang , L. , O. Gascuel , and N. Lartillot ( 2008 ). Empirical profile mixture models for phylogenetic reconstruction . Bioinformatics 24 , pp. 2317 – 2323 . doi: 10.1093/bioinformatics/btn445 . OpenUrl CrossRef PubMed Web of Science ↵ Siglioccolo , A. , A. Paiardini , M. Piscitelli , and S. Pascarella ( 2011 ). Structural adaptation of extreme halophilic proteins through decrease of conserved hydrophobic contact surface. en . BMC Structural Biology 11 , p. 50 . doi: 10.1186/1472-6807-11-50 . OpenUrl CrossRef PubMed ↵ Soubrier , J. , M. Steel , M. S. Lee , C. Der Sarkissian , S. Guindon , S. Y. Ho , and A. Cooper ( 2012 ). The Influence of Rate Heterogeneity among Sites on the Time Dependence of Molecular Rates . Molecular Biology and Evolution 29 , pp. 3345 – 3358 . doi: 10.1093/molbev/mss140 . OpenUrl CrossRef PubMed ↵ Susko , E. , L. Lincker , and A. J. Roger ( 2018 ). Accelerated Estimation of Frequency Classes in Site-Heterogeneous Profile Mixture Models . Molecular Biology and Evolution 35 , pp. 1266 – 1283 . doi: 10.1093/molbev/msy026 . OpenUrl CrossRef PubMed ↵ Susko , E. and A. J. Roger ( 2007 ). On Reduced Amino Acid Alphabets for Phylogenetic Inference . Molecular Biology and Evolution 24 , pp. 2139 – 2150 . doi: 10.1093/molbev/msm144 . OpenUrl CrossRef PubMed Web of Science ↵ Whelan , N. V. , K. M. Kocot , T. P. Moroz , K. Mukherjee , P. Williams , G. Paulay , L. L. Moroz , and K. M. Halanych ( 2017 ). Ctenophore relationships and their placement as the sister group to all other animals . Nature Ecology and Evolution 1 , 1737 – 1746 . doi: 10.1038/s41559-017-0331-3 . OpenUrl CrossRef PubMed ↵ Whelan , S. and N. Goldman ( 2001 ). A General Empirical Model of Protein Evolution Derived from Multiple Protein Families Using a Maximum-Likelihood Approach . Molecular Biology and Evolution 18 , pp. 691 – 699 . doi: 10.1093/oxfordjournals.molbev.a003851 . OpenUrl CrossRef PubMed Web of Science ↵ Wickett , N. J. , S. Mirarab , N. Nguyen , T. Warnow , E. Carpenter , N. Matasci , S. Ayyampalayam , M. S. Barker , J. G. Burleigh , M. A. Gitzendanner , B. R. Ruhfel , E. Wafula , J. P. Der , S. W. Graham , S. Mathews , M. Melkonian , D. E. Soltis , P. S. Soltis , N. W. Miles , C. J. Rothfels , L. Pokorny , A. J. Shaw , L. DeGironimo , D. W. Stevenson , B. Surek , J. C. Villarreal , B. Roure , H. Philippe , C. W. dePamphilis , T. Chen , M. K. Deyholos , R. S. Baucom , T. M. Kutchan , M. M. Augustin , J. Wang , Y. Zhang , Z. Tian , Z. Yan , X. Wu , X. Sun , G. K.-S. Wong , and J. Leebens-Mack ( 2014 ). Phylotranscriptomic analysis of the origin and early diversification of land plants . Proceedings of the National Academy of Sciences 111 , E4859 – E4868 . doi: 10.1073/pnas.1323926111 . OpenUrl Abstract / FREE Full Text ↵ Wickham , H. ( 2016 ). ggplot2: Elegant Graphics for Data Analysis . Springer-Verlag New York . ↵ Williams , T. A. , D. Schrempf , G.J. Szöllăsi , C. J. Cox , P. G. Foster , and T. M. Embley ( 2021 ). Inferring the Deep Past from Molecular Data . Genome Biology and Evolution 13 , evab067 . doi: 10.1093/gbe/evab067 . OpenUrl CrossRef ↵ Williamson , K. , L. Eme , H. Baños , C. G. P. McCarthy , E. Susko , R. Kamikawa , R. J. S. Orr , S.A. Muñoz-Gómez , B. Q. Minh , A. G. B. Simpson , and A. J. Roger ( 2025 ). A robustly rooted tree of eukaryotes reveals their excavate ancestry . Nature 640 , pp. 974 – 981 . doi: 10.1038/s41586-025-08709-5 . OpenUrl CrossRef ↵ Wong , T. , N. Ly-Trong , H. Ren , H. Baños , A. Roger , E. Susko , C. Bielow , N. De Maio , N. Goldman , M. Hahn , G. Huttley , R. Lanfear , and B. Q. Minh ( 2025 ). IQ-TREE 3: Phylogenomic Inference Software using Complex Evolutionary Models . EcoEvoRxiv . doi: 10.32942/x2p62n . OpenUrl CrossRef ↵ M. Matschiner Wong , T. K. F. , C. Cherryh , A. G. Rodrigo , M. W. Hahn , B. Q. Minh , and R. Lanfear ( 2024 ). MAST: Phylogenetic Inference with Mixtures Across Sites and Trees . Systematic Biology 73 . Ed. by M. Matschiner , 375 – 391 . doi: 10.1093/sysbio/syae008 . OpenUrl CrossRef ↵ Yang , Z ( 1994 ). Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. en . Journal of Molecular Evolution 39 , pp. 306 – 314 . doi: 10.1007/BF00160154 . OpenUrl CrossRef PubMed Web of Science ↵ Zeldovich , K. B. , I. N. Berezovsky , and E. I. Shakhnovich ( 2007 ). Protein and DNA Sequence Determinants of Thermophilic Adaptation . PLOS Computational Biology 3 , pp. 1 – 11 . doi: 10.1371/journal.pcbi.0030005 . OpenUrl CrossRef ↵ Zhu , C. , R. H. Byrd , P. Lu , and J. Nocedal ( 1997 ). Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization . ACM Trans. Math. Softw . 23 , 550 – 560 . doi: 10.1145/279232.279236 . OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted August 09, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Modeling site-and-branch-heterogeneity with GFmix Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Modeling site-and-branch-heterogeneity with GFmix Charley G. P. McCarthy , Edward Susko , Andrew J. Roger bioRxiv 2025.08.07.669136; doi: https://doi.org/10.1101/2025.08.07.669136 Share This Article: Copy Citation Tools Modeling site-and-branch-heterogeneity with GFmix Charley G. P. McCarthy , Edward Susko , Andrew J. Roger bioRxiv 2025.08.07.669136; doi: https://doi.org/10.1101/2025.08.07.669136 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13895) Bioinformatics (41951) Biophysics (21456) Cancer Biology (18594) Cell Biology (25520) Clinical Trials (138) Developmental Biology (13381) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24323) Genetics (15612) Genomics (22510) Immunology (17738) Microbiology (40401) Molecular Biology (17184) Neuroscience (88622) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7644) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)