Meta-analysis of dispensable essential genes and their interactions with bypass suppressors

Dispensable essential gene datasets

We retrieved DEGs in S. cerevisiae from two systematic experimental datasets (Liu et al, 2015; van Leeuwen et al, 2020) and from a study that compiled data from the literature (van Leeuwen et al, 2020). Cases of bypass suppression were identified by looking for rare survivors in populations of 100–150 million cells deleted for an essential gene (van Leeuwen et al, 2020; experimental dataset), by following germination of single-deletion mutant spores (Liu et al, 2015), or by a mixture of methods, including directly testing the effect of a mutation on strain viability (van Leeuwen et al, 2020; literature dataset). Because mutation rates are generally low and specific point mutations that can bypass an essential gene are thus unlikely to arise within a single spore, the Liu et al study mainly identified cases of suppression that involved changes in chromosome number, which occur more frequently than point mutations. In contrast, the van Leeuwen et al (2020) study identified essential genes that were bypassed by ploidy changes and by single nucleotide changes in the genome, and as a result identified a higher number of DEGs. The set of tested genes are explicitly mentioned in the systematic studies, whereas for the literature set they are unknown and, therefore, we used all essential genes in S. cerevisiae. The combined dataset contained the DEGs identified in any of the three individual datasets. As tested genes, we considered all tested genes in the systematic studies and the dispensable genes identified in the literature set. We randomly generated 1,000 sets of genes of the same sizes as the individual datasets, sampling from the corresponding set of tested genes.

We calculated the overlap between the different datasets by counting the number of dispensable genes found across two and three datasets (Fig 1A). We repeated the same process in the randomly generated datasets to derive empirical P-values.

Essentiality data

In our analyses, we used essentiality data from S. cerevisiae (van Leeuwen et al, 2020), S. uvarum (Sanchez et al, 2019), C. albicans (Segal et al, 2018), S. pombe (downloaded in November 2021 from PomBase [Harris et al, 2022]), and human cell lines (Hart et al, 2015). We considered human essential genes those that were required for viability in at least three of the five cell lines tested. In C. albicans, genes with essentiality confidence scores above 0.5 were classified as essential and the remaining genes as non-essential.

Orthology mappings

We used PantherDB 16.1 (Mi et al, 2021) to identify orthology relationships (Figs 2, 5, S2, and S5). When indicated, we also used OrthoMCL (Li et al, 2003), SonicParanoid (Cosentino & Iwasaki 2019), based on the popular InParanoid (Sonnhammer & Östlund, 2015), and PomBase (Wood et al, 2012) orthology mappings.

Functional enrichment of dispensable essential genes

For each DEG set and each of the 14 broad functional classes previously defined (Costanzo et al, 2016), we calculated the fold enrichment as the fraction of DEGs annotated to that functional class with respect to the corresponding fraction of core essential genes (Fig S1A). We calculated the statistical significance with two-sided Fisher’s exact tests and corrected for the multiple tests across the 14 functional classes using the FDR. We considered the cases with a P < 0.05 and FDR < 10% as significant enrichments and depletions.

Enrichment for non-essential genes in the Sigma1278b strain

For each dispensable gene set, we calculated the fold enrichment as the ratio of DEGs identified as non-essential in the Sigma1278b strain divided by the analogous ratio of core essential genes (Fig S1D). P-values were calculated using two-sided Fisher’s exact tests.

Complex dispensability bias

For each DEG set, we counted the number of complexes (Meldal et al, 2021) in which all essential subunits were identified either as dispensable or core essential genes. We repeated the same process using the randomly generated datasets to derive empirical P-values (Fig S1B).

Properties of dispensable essential genes

We queried a panel of 17 numeric and four binary features to characterize the set of DEGs (Fig 1C). Protein-level numeric properties included abundance (Ho et al, 2018), half-life (Belle et al, 2006), degradation rate (Christiano et al, 2014), number of structural domains identified by Pfam 34.0 (Mistry et al, 2021), fraction of structurally disordered residues calculated by VLS2b (Peng et al, 2006) downloaded from d2p2.pro (Oates et al, 2013), sequence length, and number of protein–protein interactions degree (Koch et al, 2012). For members of protein complexes (Meldal et al, 2021), we counted the number of different complexes in which they were found and calculated their cocomplex degree as the number of protein partners present in those complexes. Gene-level numeric properties included transcript count (Lipson et al, 2009), expression variance under different environmental conditions (Gasch et al, 2000), coexpression degree calculated as the number of genes with similar expression profiles (i.e., MEFIT scores > 1) (Huttenhower et al, 2006), multifunctionality calculated as the number of GO SLIM (Ashburner et al, 2000) biological process annotations downloaded from SGD (Cherry et al, 2012), number of paralogs (i.e., copy number) (Koch et al, 2012), ratio of non-synonymous to synonymous substitutions (dN/dS) to quantify sequence evolution (Koch et al, 2012), and number of yeast (i.e., yeast conservation) and distant (i.e., broad conservation) species in which the gene is conserved (Koch et al, 2012). We also defined four binary features to describe if proteins had a structural domain (Mistry et al, 2021), localized to a membrane (Babu et al, 2012), or belonged to a protein complex (Meldal et al, 2021), or if genes had any duplicate according to YeastMine (Cherry et al, 2012).

For each numerical feature, values that corresponded to DEGs were z-score normalized using the median and SD of the values that corresponded to core essential genes. A resulting positive median z-score identifies a feature in which DEGs tend to have higher values than core essential genes. Conversely, a negative median z-score identifies a feature in which DEGs tend to have lower values than core essential genes. Dot size in plots is proportional to the absolute median z-score value and the dot color identifies the set of genes with higher feature values. We calculated the statistical significance by means of Mann–Whitney U tests. For each binary feature, we calculated the fold enrichment as the ratio of DEGs with that particular feature divided by the equivalent ratio for core essential genes. Dot size is proportional to the absolute value of the log₂ of the fold enrichment and the dot color identifies the set of genes with a higher ratio for a particular feature. We calculated the P-values with two-sided Fisher’s exact tests. We corrected for the multiple tests across the 21 features by calculating the FDR. All shown dots correspond to features with significant differences (P < 0.05, Mann–Whitney U test or Fisher’s exact test, and FDR < 10%) between the gene sets. For visualization purposes, all significant enrichments with median z-scores or fold enrichment values below 0.1 were shown with dots of the same size. We followed the same approach to characterize (i) dispensable essential versus non-essential genes (Fig 1C); (ii) essential versus non-essential genes (Fig 1C); (iii) DEGs with multiple suppressors versus DEGs with a single suppressor (Fig S3F).

Analyses on S. cerevisiae strains

We downloaded gene presence/absence data for a large panel of S. cerevisiae strains (Li et al, 2019a) and defined several core pangenome gene sets at different stringency levels (see x-axis in Fig S1E). For instance, a threshold of 10 identifies the core pangome composed of all genes, absent only in 10 strains or less. For each DEG dataset and pangenome, we calculated the fraction of DEGs missing from the pangenome and the corresponding fraction for core essential genes, from which we calculated the fold enrichment. We also calculated fold enrichments for core essential genes versus the complete set of DEGs and for the non-essential genes versus essential genes. P-values were calculated with two-sided Fisher’s exact tests.

We retrieved precomputed LOF data for S. cerevisiae strains (Peter et al, 2018) from http://1002genomes.u-strasbg.fr/files/, including frameshift mutations and missense mutations predicted to be deleterious by SIFT (Ng & Henikoff, 2001). We calculated the number of strains in which these mutations affected each gene and aggregated the results per gene set (i.e., dispensable essentials, core essentials, and non-essentials). P-values were calculated using two-sided Fisher’s exact tests (Fig S2A).

For each strain, we counted the genes affected by CNL events in a panel of S. cerevisiae strains (Peter et al, 2018) and aggregated the result per gene set (Fig S2C). P-values were calculated using two-sided Fisher’s exact tests. Finally, we retrieved dN/dS data for the same panel of S. cerevisiae strains and grouped them by gene set (Fig S2B). P-values were calculated using Mann–Whitney U tests.

Orthology relationships of dispensable essential genes

For each gene, we calculated its orthology relationships in C. albicans, S. pombe, and human (Figs 2A and B and S2D and E). Specifically, we considered gene absence, gene duplication (including 1:N and N:M orthology relationship), N:1 relationships, and 1:1 orthologs. For 1:1 orthologs, we evaluated the essentiality in the target species. For each species and property, the fold enrichment was calculated as the fraction of DEGs with respect to the fraction of core essential genes with that property. P-values were calculated by two-sided Fisher’s exact tests. We used the same approach to compare dispensable essential to non-essential genes (Fig S2H) and non-essential to essential genes (Fig S2G and I).

Gene age

For each gene, we calculated its age by identifying the farthest species from S. cerevisiae with a present ortholog. We used orthology relationships for 98 species from PantherDB (Mi et al, 2021). The phylogenetic tree to calculate species relationships was downloaded from UniProt (UniProt Consortium, 2021), and for each species, we calculated the distance to S. cerevisiae as the number of main branches separating them. Thus, genes with age 0 are specific to S. cerevisiae and not present in any other of the 98 species, whereas age five corresponds to genes present in the most distantly related species. We grouped gene ages for each gene set (core, dispensable, and non-essentials) and calculated P-values with Mann–Whitney U tests (Fig 2C).

Gene loss

For each gene of age X, we calculated the fraction of species closer to S. cerevisiae (distance < X) in the phylogenetic tree with that gene absent from their genome. For instance, for a given gene of age 3, we calculated the fraction of species at distance 1 or 2 to S. cerevisiae with the gene of interest absent. We aggregated data for each gene set (core essentials, dispensable essentials, and non-essentials) and calculated P-values by means of Mann–Whitney U tests (Fig 2E).

Cancer cell lines

We used fitness data from genome-scale CRISPR–Cas9 knockout screens in 1,070 cancer cell lines from DepMap (Meyers et al, 2017). For each gene, we calculated the median effect of gene knockout on cell proliferation and the associated SD across all cell lines. P-values were calculated using Mann–Whitney U tests (Fig 2F).

Sequence analysis

For all 1:1 ortholog pairs between S. cerevisiae and S. pombe, we calculated their protein sequence identity (Fig 2G). Sequence length similarity was calculated as the length ratio between the shortest of the sequences with respect to the longest (Fig S2J). Thus, values closer to one describe sequence pairs of similar length, whereas values closer to 0 correspond to sequences of very different lengths. P-values were calculated using Mann–Whitney U tests. We followed the same approach to compare S. cerevisiae and C. albicans sequences (Fig S2K and L).

Interaction data

We combined suppression interactions from our recent study (van Leeuwen et al, 2020) with interactions found in the literature (van Leeuwen et al, 2020) including only deletions of essential genes suppressed under standard conditions. We generated 1,000 randomized networks respecting the topology (i.e., maintaining the total number of connections of each gene) using the BiRewire R package (Iorio et al, 2016). We calculated the number of bypass suppression pairs present in both datasets and compared that value to the number of overlapping pairs in randomized networks to derive an empirical P-value (Fig 3A).

Functional overlaps

We calculated the fraction of bypass suppression gene pairs that coded for proteins localized to the same subcellular compartment (Huh et al, 2003), had MEFIT (Huttenhower et al, 2006) coexpression scores above 1.0, were annotated to the same biological process GO term (Myers et al, 2006; Costanzo et al, 2016), and coded for members of the same complex (Meldal et al, 2021) and molecular pathway (Kanehisa et al, 2016). We repeated this calculation with the non-interacting gene pairs in the bypass suppression network and derived fold enrichments and P-values using two-sided Fisher’s exact tests. We applied this approach to the individual and the combined datasets (Fig S3A).

Complex monochromaticity by suppression mode

We selected a non-redundant set of 17 protein complexes with at least two dispensable essential subunits in the bypass suppression network. We only kept one representative complex when several complexes had the same set of DEGs. For each complex, we calculated if all dispensable essential subunits could be suppressed by the same suppressor mode (LOF or GOF). Note that in one complex, all subunits could not only be suppressed by LOF suppressors but also by GOF suppressors (indicated by “LOF & GOF” in the panel). We counted all complexes with this monochromaticity in suppression mode and compared that value to the number of monochromatic complexes in a set of 1,000 randomized bypass suppression networks to derive an empirical P-value (Fig 3C). We applied the same approach to the two individual suppression networks to discard a bias in the literature dataset (Fig S3D).

Network modularity based on cocomplex relationships

We counted the number of DEGs within the same protein complex (Meldal et al, 2021) that shared at least one suppressor. We repeated the same calculation using pairs of DEGs belonging to different complexes to derive fold enrichment and a P-value calculated with a two-sided Fisher’s exact test. We followed the same approach querying for interactors of bypass suppressors instead of the interactors of DEGs (Fig 3D).

Functional preferences

We annotated the dispensable essential and suppressor genes in the network using 14 broad functional classes previously defined (Costanzo et al, 2016). We then calculated the number of bypass suppression gene pairs within each pair of classes and repeated the process in 1,000 randomized bypass suppression networks to derive empirical P-values. We used the median frequencies in the randomized set to calculate the fold enrichments (Fig S3B, top). Only fold enrichments of significant associations are shown (P < 0.05 and FDR < 10%).

We calculated fold enrichments for suppressors as the fraction of those genes in each class with respect to the corresponding fraction of background genes (Fig S3B, bottom). We calculated P-values using two-sided Fisher’s exact tests and corrected for multiple testing using the FDR.

Overlap with genetic interactions

To evaluate the overlap between bypass suppression interactions and genetic interactions, we used a global genetic interaction map (Costanzo et al, 2016), which includes data for most gene pairs in S. cerevisiae. LOF alleles of the gene pairs were screened using hypomorphic and deletion alleles to query essential and non-essential genes, respectively. We used the standard cutoffs to identify negative and positive genetic interactions. For gene pairs screened more than once (for instance, using different hypomorphic alleles of the same gene), we implemented a consensus approach in which we considered the gene pair to genetically interact if 50% or more of the corresponding allele pairs showed an interaction. We split the bypass suppression interaction pairs by their suppression mode and counted for each set the number of overlapping pairs with negative and positive genetic interactions. Next, to calculate the expected overlap by chance with genetic interactions, we generated a background set including all possible DEG and suppressor gene pairs, and removed the pairs present in the bypass suppression interaction network. Then, we calculated the overlap of this background gene pairs with negative and positive interactions. We compared the corresponding ratio to the overlaps obtained using the bypass suppression interaction network, and calculated the statistical significance by two-sided Fisher’s exact tests (Fig 3B).

Clustering of genetic interaction profiles

We downloaded the genetic interaction profile similarities computed using the complete genetic interaction profiles from https://thecellmap.org/ (Usaj et al, 2017). For genes with multiple alleles, we averaged the similarity values across alleles. We performed hierarchical clustering of the genetic interaction profiles using the R function hclust, and used a cutoff of 0.2 to define a network of genes with high-profile similarities (Fig 3F).

Agreement in copy number changes and suppression mode across S. cerevisiae strains

We defined CNL and copy number gain (CNG) events as having a copy number below 1 or above 1, respectively, per haploid genome as defined by the 1,011 genomes project (Peter et al, 2018). The study used a combination of flow cytometry, sequencing coverage analysis across 1-kb windows, and allele frequency distributions to determine copy number changes, including both aneuploidies and segmental duplications (Peter et al, 2018). For each bypass suppression gene pair, we calculated the number of strains in which both genes had a CNL (i.e., co-loss events, CNL–CNL) and a CNL event for the DEG and CNG for the suppressor gene (i.e., loss-gain events, CNL–CNG). We disregarded 18 hypermutated strains with copy number changes in >33% of the genes in the bypass suppression network and aggregated co-loss and loss-gain events for all bypass suppression gene pairs after splitting pairs by their suppression mode (LOF or GOF). We repeated the same calculation with a background set of gene pairs composed of all possible DEG-suppressor pairs, after removing the gene pairs in the bypass suppression network. We compared the proportion of co-loss versus loss-gain events for the LOF and GOF bypass suppression pairs, for LOF bypass suppression pairs and background pairs, and for GOF bypass suppression pairs and background pairs (Fig 4A). We calculated the statistical significance by two-sided Fisher’s exact tests. Finally, we followed the same approach to count the number of strains in which each set (LOF and GOF bypass suppression and background pairs) overlapped more often with co-loss than loss-gain events (and vice versa). We compared the resulting proportions as explained above and calculated the statistical significance with two-sided Fisher’s exact tests (Fig S4A).

Co-mutation in S. cerevisaie strains

We defined as deleterious mutations the missense mutations predicted as damaging by SIFT, indel mutations, and frameshift mutations (see the section “Analyses on S. cerevisiae strains” above). For each DEG, we retrieved the strains in which it had a deleterious mutation and checked if any of its bypass suppressor genes was also mutated in any of those strains. We counted the number of DEGs co-mutated in any strains with any of their suppressor genes, the number of DEGs mutated alone, and the number of DEGs not mutated in any strain. We repeated the same process using 1,000 randomized bypass suppression networks. We performed this calculation using (1) LOF bypass suppression pairs and haploid strains (Fig 4B); (2) LOF bypass suppression pairs and diploid strains (Fig S4B); and (3) GOF bypass suppression pairs and haploid strains (Fig S4C).

Phenotypic changes across species and presence of bypass suppressor mutations

We hypothesized that the relationships between DEGs and bypass suppressor mutations identified in S. cerevisiae should be reflected in the evolutionary landscape of other species. To test this hypothesis, we identified DEGs that were non-essential or absent in a given target species and evaluated if the bypass suppressor mutations were fixed in the given target genome. To determine if a bypass suppressor mutation was fixed in another species, we took into account the effect of the suppressor mutation on gene function. Briefly, for LOF suppressors, we evaluated if mutations in the target species would reduce the gene activity with respect to S. cerevisiae. Conversely, for GOF suppressors, we evaluated if the mutations would increase the gene activity.

First, we annotated the orthology relationship of each DEG in S. pombe. We only considered DEGs absent in S. pombe or with a 1:1 ortholog. For genes with 1:1 orthologs, we annotated the essentiality of the ortholog in that species. We also annotated the orthology relationships of bypass suppressor genes in S. pombe. For suppressors with 1:1 orthologs, we performed a sequence alignment between the protein sequences of both species.

We next describe the set of rules that we evaluated to identify cases with equivalent bypass mutations in S. pombe. Briefly, in LOF suppressors, we looked for orthologs with decreased activity with respect to the S. cerevisiae gene, whereas in GOF suppressors, for orthologs with increased activity. The first set of rules was based on orthology relationships. We considered S. pombe to have a LOF bypass mutation if the suppressor gene was absent. Also, if it had an N:1 ortholog, which could be similar to a copy number decrease and, thus, a decrease in activity. Suppressors with more than one ortholog in S. pombe or with a 1:1 ortholog were considered non-equivalent LOF bypass mutations because their copy number did not decrease. Conversely, we considered as GOF bypass mutations cases in which the suppressor gene had more than one ortholog, similar to increasing their copy number and their activity, and non-equivalent GOF bypass mutations cases in which there was a N:1, 1:1, or absent ortholog in S. pombe.

The second set of rules we used to evaluate equivalent mutations was based on protein sequences. We only considered frameshift, nonsense, and missense mutations of suppressor genes with 1:1 orthologs. For the rest of cases, only the orthology rules (see above) were applied. The position of the nonsense and frameshift suppressor mutations identifies the part of the protein that should remain functional. Functionality encoded beyond that residue is compromised. Thus, we considered 1:1 orthologs in S. pombe with a shorter sequence than the position of the nonsense or frameshift suppressor mutation as LOF bypass mutations. Conversely, we considered cases in which the ortholog sequence was equal or longer than the position of the nonsense or frameshift suppressor mutation as non-equivalent LOF bypass mutations. In cases with missense mutations, we performed a sequence alignment between the S. cerevisiae suppressor gene and its 1:1 ortholog in S. pombe. We considered the ortholog to have an equivalent LOF bypass mutation if the same mutated residue or a gap was found in the aligned mutated position of the ortholog sequence. If the aligned mutated residue was the same as in the WT S. cerevisiae sequence (i.e., unmutated), we considered the ortholog to have a non-equivalent LOF bypass mutation. Cases in which the aligned position had different residues in S. pombe (not the WT and not the suppressor mutation) could not be classified as either equivalent or non-equivalent LOF bypass mutations. For GOF suppressors with a missense mutation and a 1:1 ortholog in S. pombe, we also performed a sequence alignment between the suppressor gene and its 1:1 ortholog. We considered the ortholog to have an equivalent GOF bypass mutation if the same mutated residue was found in the aligned mutated position of the ortholog sequence. If the aligned mutated residue was the same as in the WT S. cerevisiae sequence (i.e., unmutated), we considered the ortholog to have a non-equivalent GOF bypass mutation. The rest of cases could not be classified as either equivalent or non-equivalent GOF bypass mutations. We also evaluated missense mutations of suppressors with unknown suppression mode that had 1:1 orthologs in S. pombe. Cases with the exact same mutation in the ortholog were classified as equivalent bypass mutations, whereas cases in which the residue did not change in the 1:1 ortholog were classified as non-equivalent bypass mutations. The remaining suppressors with unknown suppression mode were not evaluated. Importantly, in suppressor genes with a frameshift, nonsense, or missense mutation, and with a 1:1 ortholog in S. pombe, the sequence-based assessment took precedence over the orthology based evaluation.

Finally, we considered a DEG to have an equivalent bypass suppressor in S. pombe if any of its suppressors satisfied that criteria. We grouped DEGs by their essentiality in S. pombe, expecting DEGs with equivalent phenotypes in S. pombe (i.e., absent or 1:1 non-essential orthologs) to have equivalent bypass suppressors more often than DEGs with a 1:1 essential ortholog. We calculated the fraction of genes with equivalent bypass suppressors for both gene sets to derive fold enrichment and the P-value with a one-sided Fisher’s exact test. We compared the fold enrichment of the bypass suppression network to a set of randomized bypass suppression networks, which we used to derive an empirical P-value.

We repeated the exact same process (i) using C. albicans sequences, orthology relationships, and essentiality annotations; (ii) using orthoMCL (Li et al, 2003), SonicParanoid (Cosentino & Iwasaki, 2019), and PomBase (Wood et al, 2012) as alternative orthology mappings; (iii) considering only DEGs with a single bypass suppressor to control for the bias introduced by gene degree; (iv) removing bypass suppression pairs from the literature which may have been potentially identified by phylogenetic approaches; (v) removing cocomplex and copathway bypass suppression pairs which may be more prone to present similar phylogenetic patterns; (vi) switching LOF and GOF annotations to demonstrate the specificity of our sets of rules; (vii) removing every node in the network to discard dependence on a single gene.

On the use of orthology relationships to identify gene activity changes between species

Genes duplicated in the S. pombe lineage but not in S. cerevisiae should result in 1:N orthology relationships (one gene in S. cerevisiae, N genes in S. pombe). At the moment of the duplication event, there is usually an increase in transcript and protein levels (Pavelka et al, 2010) causing an increase in gene activity. To evaluate if the initial increase in expression of duplicated genes is maintained in S. pombe, we quantified expression changes between S. cerevisiae (Lipson et al, 2009) and S. pombe (Grabherr et al, 2011; Koch et al, 2012), which highly correlate for 1:1 orthologs (r = 0.62, P < 0.0005, Fig S5K), by aggregating expression levels of orthology groups. For instance, in 1:N orthology relationships, we compared the expression level of one S. cerevisiae gene to the result of adding the expression levels of the N genes in S. pombe. In N:1 orthology relationships, we compared the aggregated expression levels of the N S. cerevisiae genes to the expression level of a single gene in S. pombe. In N:M orthology relationships, we compared the combined expression levels of N genes in S. cerevisiae to the aggregated expression level of M genes in S. pombe. In 1:1 orthologs we compared the expression levels of the ortholog in each species. For each orthology relationship, we calculated the expression ratio by dividing the expression level of the gene/s in S. pombe by the expression level of the gene/s in S. cerevisiae. In agreement with the initial copy number event, expression ratios were higher for 1:N orthology relationships than for 1:1 orthologs (median values 1.345 and 0.791, respectively, P < 0.0005, Mann–Whitney U test, Fig S5L). Conversely, expression ratios were lower for N:1 than for 1:1 orthology relationships (median = 0.568, P < 0.0005, Mann–Whitney U test). That is, for each gene transcript in S. cerevisiae, there were more transcripts in S. pombe for the aggregated 1:N orthologs than for 1:1 orthologs, and fewer transcripts for N:1 orthologs than for 1:1 orthologs. Therefore, in general, the expression level change resulting from a copy number gain or loss is kept in S. pombe.

To account for cases in which the initial copy number changes and the expression changes in S. pombe may not agree, we implemented an alternative approach to identify genes with increased or decreased activity between S. cerevisiae and S. pombe, by replacing the use of orthology relationships (see the “Phenotypic changes across species and presence of bypass suppressor mutations” section) with the quantification of expression changes between both species. For all genes with identified orthologs, we calculated the S. pombe to S. cerevisiae expression ratio between both species as explained above. Genes with an extreme decrease or increase in expression levels between orthologs (genes within the bottom 5% and top 5% expression ratios, respectively) were labeled as genes with decreased and increased activity, respectively. We considered a LOF bypass mutation equivalent in S. pombe if the suppressor gene was within the bottom 5% expression ratios and non-equivalent otherwise. Conversely, we considered a GOF bypass mutation equivalent in S. pombe if the suppressor gene was within the top 5% expression ratios and non-equivalent otherwise. By combining expression ratios with the rules based on protein sequences and gene absence (see the “Phenotypic changes across species and presence of bypass suppressor mutations” section), we calculated if DEGs absent or non-essential in S. pombe were more likely to have the corresponding bypass suppressors in that species, as explained in the previous section.

However, functional divergence may take place after the duplication event and the initial increase in gene activity may be modulated, which is difficult to evaluate without gene-specific functional assays. To control for the impact to our results of the potential functional divergence after copy number changes, we repeated the same analysis shown in Fig 5B but considering only suppressors with a 1:1 ortholog in S. pombe or absent in that species. Thus, we did not consider orthology relationships (i.e., duplicated suppressors or with N:1 orthologs) or transcript changes to identify equivalent GOF or LOF bypass suppressors in S. pombe.