Unveiling the functional fate of duplicated genes through expression profiling and structural analysis

Least diverged ortholog conjecture analysis across the Tree of Life

This study investigated the least diverged ortholog conjecture, which implies that the gene copy with the shortest branch in the gene tree tends to retain the ancestral function after a duplication event. Analysis was performed on a data set comprising 143 organisms representing different lineages of the Tree of Life, across 15,693 gene trees. This conjecture relies solely on observed branch length and classifies gene copies with the shorter branch as the least diverged ortholog and those with the longer branch as the most diverged ortholog.

However, when gene families have different evolutionary rates, it is not clear whether small differences in branch length estimates can accurately classify gene copies into these two categories (LDO/MDO). Here, in order to detect distinct evolutionary rates on each of the branches, a statistical method was developed that allows for these family-specific evolutionary rates. This method is able to identify branches in a gene tree with accelerated rates of protein sequence evolution using only branch length information. This offers two key advantages: (1) it allows us to leverage precompiled libraries of gene trees such as PANTHER, eliminating the need to construct a custom data set; and (2) it enables the analysis of ancient gene duplication events, where nucleotide substitutions have reached saturation and the usual approach of using K_a/K_s cannot be reliably applied (Smith and Smith 1996; Gharib and Robinson-Rechavi 2013). Additionally, this method assigns an expected branch length to each gene tree branch, providing a model-based score that reflects deviations from the typical evolutionary rate of the gene family. As a result, this method attempts to classify a branch as typical, short, or long relative to the expected rate for the gene family.

To capture differences in evolutionary rates across branches in each gene tree, the method estimates gene family-specific rates of evolution (family-specific factors) for each gene tree (Fig. 1A). A factor greater than one indicates that a gene family is evolving more rapidly than average, potentially reflecting relaxed functional constraints, whereas a factor below one suggests slower evolution, which may reflect stronger functional constraints. The distribution of these factors has a mean of 1.06 and a standard deviation of 0.76 (Supplemental Fig. S1; Supplemental Table S1). A slight skew toward gene families with factors below one suggests that a considerable number may be evolving under stronger evolutionary constraints.

Figure 1.

Analysis and calculation of expected gene tree branch lengths. (A) Graphical illustration of the calculation of the expected branch lengths (X). On the left, a species tree depicts the phylogenetic relationships of four species with branches highlighted by different colors and branch numbers (b₁ to b₆). On the right, two gene trees show a duplication event (node marked with a star). The first gene tree is complete, and the second tree has a gene loss represented by a dashed line. f_R denotes the family-specific rate. When no losses occur (Gene family 1), there is an exact correspondence between the speciation branch in the gene tree and the species tree (X₁ to X₁₁). When there is a loss (Gene family 2), the speciation branch in the gene tree represents the sum of the corresponding branches in the species tree (X₁₃). Note that in both gene families, when a duplication occurs, all paralogous branches are calculated, and the expected values remain the same. (B) Distributions of the differences from expected branch lengths in branches with duplications (pink) and without duplications (light blue). The sample size (n) is indicated in the figure. The P-value for the Mann–Whitney U test is shown between parentheses (P = 0.0).

Open in new tabDownload PowerPointLink to figure

Using the family-specific factors, the expected branch length was computed for all branches in every gene tree (Fig. 1A) before computing the Z-score of the differences between the observed and expected branch lengths. A significant increase in evolutionary rate was observed between branches with no duplication (average −0.02) and those with duplications (0.07; P < 0.001 Mann–Whitney U test) (Fig. 1B), suggesting that paralogs tend to evolve faster than orthologs. Similar results were observed when each major lineage was analyzed independently (Supplemental Fig. S2).

Following this analysis, duplication events were identified across all gene trees and classified into six categories (Fig. 2A; Supplemental Table S2): (i) normal-normal: both branches have a typical rate of evolution and are consistent with their expected lengths: (ii) short-short: both branches are significantly shorter than expected; (iii) long-long: both branches are significantly longer than expected; (iv) normal-long: only one branch is significantly longer than expected; (v) short-normal: only one branch is significantly shorter than expected; and (vi) short-long: one branch is significantly shorter than expected and the other is significantly longer than expected.

Figure 2.

Branch length analysis following duplication events. (A) Depiction of the six categories indicating differences in branch length following a duplication event (node marked with a star). The two rows divide the symmetric and asymmetric evolution of paralogs following duplication. The light silhouette (gray mouse) indicates the least diverged ortholog (LDO), and the dark silhouette (black mouse) indicates the most diverged ortholog (MDO). (B) Species tree showing the major lineages analyzed. The numbers in parentheses indicate the number of species included from the lineage. Gray bars represent the percentage of total duplications contributed by each lineage. Colored bars indicate the proportion of duplications within each category (color scheme as in A), relative to the total number of duplications in that lineage. The absolute numbers of duplications are provided in Supplemental Table S3.

Open in new tabDownload PowerPointLink to figure

The analysis gave a total of 926,814 duplication events across the Tree of Life, with an average of 68 duplications per family. When examining the major lineages (Fig. 2B, gray bars), the Viridiplantae lineage exhibits the highest percentage of duplications from the total number of duplications (68%) (Fig. 2B; Supplemental Table S3), followed by Deuterostomia (16%) and Protostomia (4%). Internal nodes within this tree contributed only 5% of the total number of duplications (Supplemental Table S3). These findings highlight the significant role of gene duplication in the evolutionary history of plants, in agreement with previous studies (Jiao et al. 2011; Soltis et al. 2015).

Further analysis of the distribution of duplications across the six categories within each major lineage (Fig. 2B, colored bars) revealed that the vast majority of duplications fell into the normal-normal group (∼82%). The second most abundant category was the normal-long group (∼13%), which includes paralogs with significantly different branch lengths. The remaining four categories (short-short, long-long, short-normal, and short-long) accounted for only ∼6% of duplications. When grouping these categories based on their evolutionary pattern, 86% displayed symmetric evolution (both branches following similar rates: normal-normal, short-short, and long-long), with 14% exhibiting asymmetric evolution (one significantly longer: normal-long, short-normal, short-long). These results indicate that following duplication events, the majority of paralogous pairs have similar evolutionary rates despite observed variations in absolute branch lengths.

Gene duplications in internal branches represent more ancient duplications compared to species-specific duplications (terminal branches). To investigate any differences, the distribution of duplications between these two types of duplications was examined within each major lineage subtree, excluding three major lineages which consist of only one species (Trichoplax adhaerens, Nematostella vectensis, Monosiga brevicollis). Across all cases, a large proportion of duplication events were present in terminal branches (∼70%) (Supplemental Fig. S3A), which was similar in all major lineages (Supplemental Fig. S3B). When analyzing the distribution of the duplications between the two evolutionary models, a larger proportion of duplications followed symmetric evolution (Supplemental Fig. S3C) in all major lineages (Supplemental Fig. S3D). Furthermore, the proportion of duplications following asymmetric evolution was higher in internal than in terminal branches, suggesting that ancestral duplications may have more time to evolve divergent gene copies.

LDOs show greater conservation of protein structure

Because protein structure is associated with protein function (Orengo et al. 1999), structural similarity can be used to test for functional divergence. However, structural similarity can be influenced by sequence similarity, because structural predictions often rely on sequence information. To investigate this relationship, structural similarity (measured by Foldseek LDDT) and sequence similarity (percentage of identity) were calculated between paralogs classified into symmetric and asymmetric categories. A stronger correlation between structural and sequence similarity was observed in symmetric duplications (Pearson correlation coefficient [PCC] = 0.75), whereas asymmetric duplications showed a weaker correlation (PCC = 0.55) (Supplemental Fig. S4). This substantial difference between the two indicates that different evolutionary constraints or mechanisms act following symmetric versus asymmetric duplication.

To test the least diverged ortholog conjecture, controlling for sequence divergence, structural similarity was computed between each paralogous gene copy and a corresponding co-orthologous outgroup gene. The outgroups og1 and og2 were selected to have similar levels of sequence divergence to the MDO and the LDO, respectively (see Methods). To balance the need for comparable sequence divergence with the availability of suitable cases, the absolute difference in sequence identity between MDO-og1 and LDO-og2 was required to be at most 20%.

In asymmetric duplications, the median structural similarity for MDO-og1 pairs (0.723, dark orange) remained consistently lower than that for LDO-og2 pairs (0.7662, light orange) (Fig. 3A). In contrast, differences were minor for symmetric duplications (MDO-og1: 0.8506; LDO-og2: 0.8557). These findings suggest that, although sequence similarity influences structural similarity, structural differences may persist beyond what sequence alone can explain, potentially reflecting functional divergence.

Figure 3.

Structural similarity between paralogs and orthologs. (A) Structural similarity of each paralog to two different orthologs (og1 and og2) (n = 466,476). Stars on the tree nodes indicate duplication events. Box colors correspond to the arrows in the figure, indicating the comparisons made. The P-value is indicated on top of the figure for the Mann–Whitney U test (**) P < 0.01. All comparisons and P-values are indicated in Supplemental Table S4. (B) Comparison of the excess of structural divergence and the asymmetry strength (n = 466,476). The shaded area represents the 95% confidence interval around the mean. The x-axis represents increasing bins of asymmetry strength.

Open in new tabDownload PowerPointLink to figure

To further explore this, the asymmetry strength (measured as the observed branch length differences between paralogs) was compared to the excess structural divergence (calculated as LDDT[MDO, og1] – LDDT[LDO, og2]) (Fig. 3B). The analysis revealed that greater asymmetry was associated with a stronger structural divergence of the MDO copy. A positive correlation was observed up to an asymmetry strength of approximately 1.8 (Fig. 3B). The decline in the excess of structural divergence may reflect the retention of conserved structural motifs despite the large sequence divergence.

To ensure that these results were not influenced by the choice of sequence identity difference threshold, additional analyses were performed using different distance similarity thresholds between 5% and 50%. In all cases, similar trends were observed (Supplemental Figs. S5, S6). Taken together, these results indicate that increased asymmetry in sequence evolution between paralogs is associated with greater structural divergence from the outgroup, even after controlling for sequence similarity. This supports the use of protein structure as a proxy for functional divergence. In line with the LDO conjecture, asymmetric duplications exhibit greater functional divergence, with the LDO copy more likely to retain the ancestral function.

Differences in branch lengths support differences in gene coexpression

To test functional differences across genes, gene coexpression was also used as a proxy for functional similarity. For animals, expression atlases of 16 species across 32 anatomical entities were obtained from the Bgee database (Bastian et al. 2021). For plants, publicly available expression atlases from 20 species across 72 anatomical samples were used (Supplemental Table S5; Julca et al. 2021; Koh et al. 2023). In order to assess if the observed differences in branch length following gene duplication have functional implications, gene coexpression between paralogs resulting from symmetric and asymmetric evolution was compared using the Pearson correlation coefficient. Subsequently, the direction of change was investigated by comparing the expression profiles of each paralog to their co-orthologous outgroup gene.

When analyzing pairs of paralogs in each category, within both animals and plants, paralogs from the asymmetric category showed lower correlation values compared to those from the symmetric category (Fig. 4A). Furthermore, because most duplications were located in terminal branches and represent more recent duplication events, their gene expression profiles were evaluated separately (Supplemental Fig. S7). In both cases, terminal and internal branches present consistent results, with paralogs in the symmetric category showing slightly higher coexpression values. Notably, symmetrically evolving paralogs arising from duplications on the terminal branches within plants exhibited the highest coexpression (Supplemental Fig. S7). This observation may suggest that recent duplications have undergone less evolutionary divergence due to a shorter time frame. However, the same pattern is not observed in animals.

Figure 4.

Functional analysis of paralogs using expression data. (A) Gene coexpression values of paralogs following symmetric (n = 275,988) and asymmetric (n = 31,671) duplications, calculated using the Pearson correlation coefficient (PCC). The P-values are indicated on top of the figures for the Mann–Whitney U test. (*) P < 0.05, (**) P < 0.01. (B) Gene coexpression (PCC) between the LDO (dashed lines) and MDO (solid lines) genes and their closest outgroup co-ortholog following symmetric (blue) and asymmetric (orange) duplications. The sample sizes for this plot are as follows: plants symmetric n = 114,702, plants asymmetric n = 13,533, animal symmetric n = 17,493, and animal asymmetric n = 7490. (C) Violin plot showing the tissue specificity differences (delta Tau) between the paralogs following symmetric (animals n = 77,889, plants n = 275,988) and asymmetric (animals n = 17,205, plants n = 31,671) duplications. (D) Tissue specificity differences (delta Tau) between the LDO (dashed lines) and MDO (solid lines) genes and their closest orthologs following symmetric (blue) and asymmetric (orange) duplications. The sample sizes are the same as in B. In plots B and D, the x-axis indicates the divergence time between the species compared in time intervals. Plant and animal data sets were compared independently. The comparisons used are shown above each panel, with duplication events marked by stars.

Open in new tabDownload PowerPointLink to figure

The results may be influenced by the different species’ evolutionary history (e.g., polyploid plants). Therefore, each species was analyzed independently and showed consistent results (Supplemental Fig. S8). Among the animals, 13 species showed similar patterns, indicating lower coexpression values between paralogs within the asymmetric category. However, three species, Felis catus, Drosophila melanogaster, and Gorilla gorilla gorilla, did not show significant differences between the two categories. For plants, 19 species supported the previous pattern, and only one (Amborella trichopoda) did not. Moreover, the polyploid plants (Brassica napus, tetraploid; Triticum aestivum, hexaploid; Glycine max, tetraploid) exhibited major differences between symmetric and asymmetric paralogs. In summary, these results reveal a general trend where paralogs following asymmetric evolution exhibit less coexpression, suggesting that sequence divergence between paralogs correlates with differences in gene expression. Nevertheless, discrepancies observed in a few species highlight the importance of including many different species in this kind of study.

To investigate the direction of functional changes following duplication—whether the least or most diverged ortholog conserves the ancestral gene function—coexpression patterns were analyzed between each paralogous gene copy and their co-orthologous outgroup gene. For plants, 112,015 duplications were analyzed, with 99,995 symmetric and 12,020 asymmetric duplications. For animals, 22,885 duplications were examined, of which 16,273 were symmetric and 6612 were asymmetric duplications. Then, within each category, gene coexpression was computed (PCC) between each copy and the outgroup gene (Fig. 4B). Differences were observed in both symmetric and asymmetric categories, with greater differences observed between the asymmetric category (Fig. 4B). Moreover, in both data sets—animals and plants—the gene following the longest branch (MDO in the asymmetric category) consistently displays low coexpression values (Fig. 4B), suggestive of a lower functional similarity with the ancestral gene. Similar results were observed by analyzing each species independently (Supplemental Fig. S9).

To study the effect of time since duplication, coexpression was compared across species' divergence times (Fig. 4B). In both animals and plants, more ancient paralogs exhibited a slight decrease in coexpression.

In summary, these results suggest that paralogous genes descending from duplications with significantly different branch lengths have more distinct coexpression values. In particular, the most diverged orthologs display the greatest change in gene expression, thus suggesting that they have a greater functional divergence from the ancestral gene.

MDOs show more functional specificity

Previous studies have shown that, after gene duplication, paralogs tend to increase the levels of tissue specificity, which occurs promptly after the duplication event (Huminiecki and Wolfe 2004; Huerta-Cepas et al. 2011). To test whether the least diverged ortholog (paralog with a shorter branch) or the most diverged ortholog (paralog with a longer branch) has different tissue specificity, the Tau score (defined in Methods) was computed for each gene in every species.

To compare each category, the difference in tissue specificity was calculated (delta Tau) between the most diverged ortholog and the least diverged ortholog in each group (symmetric and asymmetric duplications, Fig. 4C). A value of zero indicates no difference in sample-specificity between the compared genes, whereas positive values indicate that the MDO exhibits more specificity and negative values signify that the LDO is more specific. The results of this analysis show that paralogs within the asymmetric category, in both animals and plants, have slightly higher positive values, indicating that the MDOs tend to have a more tissue-specific pattern of expression (Fig. 4C). Similar results were observed when analyzing duplication events separately for internal and terminal branches, with terminal branches exhibiting more pronounced differences (Supplemental Fig. S10) and independently within each species (Supplemental Fig. S11). These results suggest that the gene copy which exhibits greater sequence divergence tends to have greater tissue specificity.

To determine whether the most diverged ortholog also displays increased tissue specificity compared to the ancestral gene, the difference was computed between the Tau values for each paralogous gene copy and their co-orthologous outgroup gene (Fig. 4D). A zero value indicates no difference in tissue specificity, whereas positive values suggest increased specificity of the paralogous gene. In both categories, there are differences between the paralogs, with the MDO genes exhibiting higher delta Tau values. Specifically, in both plants and animals, the MDO gene of the asymmetric category consistently showed the greatest increase in tissue specificity over time (Fig. 4D). Moreover, similar results were observed when analyzing individual species (Supplemental Fig. S12). In addition, the delta Tau in animals showed a slight decrease over time, whereas in plants this decrease was not apparent (Fig. 4D). This result indicates that tissue specificity in animals may decrease over time, with recent duplications exhibiting higher specificity. To summarize, these findings suggest that when one gene copy is under selective pressure to retain the ancestral function (LDO), the other with the longer branch (MDO) tends to acquire a more specialized role, which may differ over time.

Well-characterized gene families corroborate the LDO conjecture

To explore the findings of this study, well-characterized gene families from the literature with available experimental data were examined. The analysis revealed that symmetric duplications (normal-normal category) generally retain similar functions between paralogs, whereas asymmetric duplications (normal-long category) are more likely to diverge functionally. In these asymmetric pairs, the least diverged ortholog typically preserves the ancestral function.

A representative example is the RNF113 gene family (Supplemental Fig. S13A), which duplicated independently in human and mouse. The human duplication is classified as normal-long, whereas the mouse duplication is normal-normal. In Drosophila, the ortholog mdlc is broadly expressed and essential for both neuronal development and spermatogenesis (Carney et al. 2013; Brattig-Correia et al. 2024). In humans, the LDO gene RNF113A retains broad expression and the ancestral neuronal role (Tsampoula et al. 2022), whereas RNF113B (MDO gene) is testis-specific and essential for meiosis (Brattig-Correia et al. 2024). In mouse, both Rnf113a1 and RNf113a2 retain broad expression, yet Rnf113a2 (MDO gene) is required for spermatogenesis and Rnf113a1 (LDO gene) for neuronal development (Brattig-Correia et al. 2024), indicating some functional partitioning despite the classification as normal-normal.

Several clear examples of functional redundancy in normal-normal pairs were also identified. SEC10a/b in Arabidopsis (Supplemental Fig. S13B) share exocytosis-related functions (Vukašinović et al. 2014). The paralogs Sar1a/b in mouse (Supplemental Fig. S13C) are fully interchangeable for viability (Tang et al. 2024). For the paralogs Cdx1/Cdx2 (Supplemental Fig. S13D), a knock-in of Cdx2 rescues the Cdx1-null phenotype without defects (Savory et al. 2009).

In the normal-long category, PRP18A/B in Arabidopsis (Supplemental Fig. S13E) show functional divergence. PRP18A (LDO) retains a role in alternative splicing, similar to its yeast ortholog, whereas PRP18B (MDO) shows no detectable expression in 2-week-old seedlings, and T-DNA insertion mutants exhibit no observable phenotype compared to wild-type plants (Kanno et al. 2018). The paralogs LAZY1/LAZY6 in Arabidopsis (Supplemental Fig. S13F) also illustrate this pattern: LAZY1 (LDO) controls shoot orientation as in rice, whereas LAZY6 (MDO) has a different expression profile and no significant phenotype in mutants (Li et al. 2007; Yoshihara and Spalding 2017; Waite and Dardick 2024). Similarly, studies in knockout mice suggest that Becn1 and Becn2 (Supplemental Fig. S13G) have distinct functions in regulating autophagosome formation and mitophagy (Galluzzi and Kroemer 2013; Quiles et al. 2023).

Overall, these examples are consistent with the broader patterns described in this study. Nevertheless, certain cases, such as the Rnf113 family in mice, highlight the complexity of functional evolution postduplication and suggest that exceptions or intermediate scenarios exist. Whereas general trends across the Tree of Life were the focus of this study, gene-specific dynamics should be considered in future analyses.

Gene tree processing

To study the least diverged ortholog conjecture, the species tree for 143 organisms and 15,693 gene trees were downloaded from PANTHER v18.0 (Thomas et al. 2022). PANTHER gene trees are annotated with evolutionary information in their internal nodes, including speciation, gene duplication, horizontal gene transfer (HGT), or unknown events (nodes for which the type cannot be reliably inferred). Specifically in bacteria, HGT events appear to be underestimated in this data set due to the methodology and the type of data used (single reference genome per bacterial species rather than a pangenome). Some specific branches in the gene trees were excluded from the analysis. These include branches with poorly fitting distances (arbitrarily set to 2.0 substitutions/site in PANTHER) and some duplication events that involve: duplications at the root of the gene tree (no speciation event to compare); those with horizontal gene transfer or “unknown” events annotated on either side of the duplication; as well as those with gene losses immediately following the duplication. Finally, a total of 15,577 gene trees were used for further analysis (Supplemental Table S2).

Calculation of the expected gene tree branch lengths

Expected gene tree branch lengths were calculated for 2,648,252 branches across 15,577 gene trees using a simple evolutionary model: the expected length between two speciation events in a gene tree is equal to the corresponding branch in the species tree scaled by a gene family-specific rate. In the case of gene loss, the model assumes that the branch length in the gene tree represents the sum of all the branches corresponding to the branch in the species tree scaled by the family-specific rate (Fig. 1A).

By expressing this as a system of approximate linear equations, it is possible to optimize the family rates and species tree branch lengths, constrained by the observed branch lengths within the many gene families. That is, $|X|\approx \sum _{b_i\in X}\left|b_i\right|f_R$, where X represents a specific branch in a gene tree (|X| being the length of that branch), which corresponds to a set of branches in the species tree, each denoted b_i (with length |b_i|), and f_R representing the family-specific factor.

The species tree used in PANTHER only provides topological relationships. As such, both the species tree branch lengths and family-specific factors were estimated from the observations of gene tree branch lengths. The method consists of two main steps:

Using these species tree branch lengths and family-specific factors, expected branch lengths were calculated in the gene trees, before identifying statistically significant rate shifts (Supplemental Table S2). The Python scripts generated for this method are available on GitHub (https://github.com/DessimozLab/ldo_study) and in the Supplemental Scripts.

Identification of branch length differences after duplication

Because all paralogs resulting from a duplication event can be classified as either least diverged orthologs (shorter branch) or most diverged orthologs (longer branch), a statistical method was developed to determine whether the differences observed in the branch length postduplication were significant.

For 15,577 gene trees, the difference between observed and expected branch lengths (calculated under the model described in the previous section, from speciation event to speciation event) was computed (Fig. 1A). Then, to assess the statistical significance of the observed excess or depletion, these differences were converted into standardized Z-scores using the z-score function from the SciPy Python package v1.13.0 (Virtanen et al. 2020). An excess substitution can be interpreted as a potential signal of branch-specific adaptive selection and/or relaxed constraints. The resulting distributions of the normalized differences between branches with (876,261) and without (1,771,991) duplication were compared using the Mann–Whitney U test (also using SciPy). All further analysis filtered out gene trees without duplication events, which resulted in 13,597 gene trees and 364,829 duplication events distributed across different lineages of the Tree of Life.

For each duplication event, the two branches (Fig. 1A) were compared using their Z-score of the difference between observed and expected (from speciation to speciation event) by computing P-values and using a two-tailed test (with significance level α = 0.05) to identify significantly shorter and longer than expected branches. When more than two gene copies were present, indicating either successive duplications, multifurcation (i.e., duplication nodes with more than two descendant branches), or both, the analysis was performed at every node starting from the leaves (Supplemental Fig. S14). Then, the least diverged ortholog was identified and compared to the other paralogs. This approach ensures that a long branch is included only once in the analysis.

Based on these comparisons, duplication events were placed into six categories (Fig. 2A): (i) normal-normal: both branches are not significantly different; (ii) short-short: both branches are significantly shorter than expected; (iii) long-long: both branches are significantly longer than expected; (iv) normal-long: only one branch is significantly longer than expected; (v) short-normal: only one branch is significantly shorter than expected; and (vi) short-long: one branch is significantly shorter than expected and the other is significantly longer than expected. Using this information, the data set was divided into two duplication models: symmetric (normal-normal, short-short, long-long) and asymmetric evolution (normal-long, short-normal, short-long). All subsequent analyses were based on this classification into these two categories: symmetric and asymmetric duplication events.

Pairwise comparisons of structural and sequence identities

Structural data for 1,829,120 proteins out of 1,968,858 of the PANTHER data set were downloaded from the AlphaFold Protein Structure Database v4 (Varadi et al. 2022). The 139,738 remaining proteins did not have structural predictions available. Corresponding amino acid sequences for all proteins with available structural data were obtained from the UniProt database using the REST API (downloaded on 19 March 2025).

Pairwise structural similarity was computed for 926,814 paralog pairs using Foldseek version 8.ef4e960 (van Kempen et al. 2024). The Foldseek LDDT (Local Distance Difference Test) score reports the average LDDT of the alignment. LDDT scores express the percentage of interatomic distances and range from 0 (not conserved distances) to 1 (perfect model) (Mariani et al. 2013).

To estimate sequence similarity, multiple sequence alignments were generated for each gene family using MAFFT v7.526 (Katoh and Standley 2013). From these alignments, pairwise percentage identity between paralogs was calculated. The Pearson correlation coefficient was then computed to assess the relationship between sequence and structural similarity, separately for symmetric and asymmetric duplication categories.

In addition, pairwise sequence and structural similarity were calculated for each paralog and their co-orthologous outgroup gene (857,663 comparisons), using the approach described above. Two outgroup genes were selected to account for differences in the absolute branch length of each paralog. For the first outgroup (og1), co-orthologous genes were selected from the closest outgroup species before the duplication event in the gene tree. Then, the closest gene in that species (least diverged copy) was selected.

For the selection of the second outgroup (og2), the aim was to minimize the differences between LDO-og2 and MDO-og1, such that the sequence differences between MDO and og1 would be comparable to the differences between LDO and og2. The branch distance between MDO and og1 was calculated, and then an external ortholog with a similar distance to the LDO copy was selected as og2. To also account for sequence divergence, pairwise percentage identities between MDO and og1 and LDO and og2 were calculated using the multiple sequence alignment generated before (Supplemental Table S6).

For the analysis using the second outgroup (og2), the ideal scenario is that the difference in sequence similarity between MDO-og1 and LDO-og2 would be zero. In such cases, the comparison between MDO-og1 and LDO-og2 would be perfectly balanced in terms of sequence or branch distance. However, applying this strict filter drastically reduced the data set size (from 548,818 to 1319 (0.2%)). To address this, six more relaxed thresholds were used to filter for differences in percentage identity of: <5%, <10%, <20%, <30%, <40%, and <50%.

Distributions of the different data sets were compared using the Mann–Whitney U test.

Asymmetry strength and excess of structural divergence analysis

To test whether the observed branch length differences between paralogs (referred to as asymmetry strength) correlate with structural divergence, the branch length differences for both symmetric and asymmetric duplications were calculated. Then, the differences in structural similarity (LDDT) between the comparisons LDO-og2 and MDO-og1 were computed (Supplemental Table S7). This analysis was performed for all data sets generated using the different sequence similarity and branch length thresholds described in the previous section.

Compilation of gene expression atlases

As another functional characteristic of paralogous genes, expression atlases of 16 animals and 20 plants were used. For animals, RNA-seq gene expression data and anatomical annotations (UBERON terms) were retrieved from Bgee 15.1 (Bastian et al. 2021). This consisted of 2639 RNA-seq experiments across 32 different anatomical samples. For plants, expression atlases were taken from previous studies (The Plant Expression Omnibus) (Julca et al. 2021; Koh et al. 2023), and anatomical entities were obtained from the sample names of RNA-seq experiments in the NCBI BioSample database. Manual curation was conducted to define anatomical entity names (e.g., leaf, flowers), resulting in 72 different anatomical samples for plants. All species included, for both animals and plants, had at least four different anatomical samples (Supplemental Table S5). Also, only genes with TPM ≥ 2 in at least one sample were retained (plants = 446,852 and animals = 167,106 genes). Then, the TPM data of each gene was log-transformed using the hyperbolic arcsine function (Johnson and Krishnan 2022): $\arcsin h\left(x\right)=\ln \left(x+\sqrt{x^2+1}\right)$. Finally, the average value was computed for samples coming from the same anatomical entity.

Gene coexpression analysis

To test for functional differences between pairs of paralogs (animals 95,094, and plants 307,659 pairs), as well as each paralog and the closest orthologous outgroup (animals 24,983, and plants 128,235 pairs), gene-gene coexpression patterns were computed using the Pearson correlation coefficient. To facilitate comparison across different species, only equivalent anatomical entities were used (e.g., leaf for Arabidopsis and tomato) (Supplemental Table S5). When selecting an orthologous outgroup to compare the paralogous genes to, the closest species was used for which data was available (Supplemental Table S8). Additionally, species were removed from the outgroup analysis when there were too few orthologous genes (<50 genes). Finally, distributions of the different data sets were compared using the Mann–Whitney U test.

Analysis of gene expression specificity

Sample-specificity of genes, based on expression data, was calculated using the Tau score (Yanai et al. 2005; Kryuchkova-Mostacci and Robinson-Rechavi 2017). Tau values range from 0 to 1, where 0 indicates a gene is broadly expressed and 1 denotes tissue-specific expression. Gene pairwise comparisons were performed by calculating the difference in tissue specificity (delta Tau) between the two genes. Comparisons were performed for the pairs of extant paralogous genes (animals 95,094, and plants 307,659 pairs), as well as independently for each gene copy with the chosen orthologous outgroup gene (animals 24,983, and plants 128,235 pairs). Differences in their distribution were assessed using the Mann–Whitney U test.

Estimation of divergence times

Divergence times (adjusted times) between species were obtained from the TimeTree database (Kumar et al. 2022) and are expressed in millions of years ago (MYA). These time estimates were used to plot the differences between genes in the expression profile and sample specificity analyses.