Phylogenetic placement of whole genome duplications in yeasts through quantitative analysis of hierarchical orthologous groups

Introduction

The potential roles of gene duplication in driving evolutionary innovation have now long been recognized.¹ Whole genome duplication (WGD) is a major event leading to duplicated genes, found several times over the course of evolution in many eukaryotic organisms,² and leads to the formation of an organism (or a cell) with additional copies of its entire genome. This event is often followed by massive gene losses due to non-functionalization, however some duplicated genes might undergo subfunctionalization (gene partitioning) or neofunctionalization (acquire a novel gene function) which will favor their preservation.³

To understand the circumstances leading to WGD in eukaryotes, Marcet-Houben and Gabaldón proposed a hypothesis describing an ancient allopolyploidization event (hybridization event between two different species of the yeast prior to the genome doubling).⁴ In case of an autopolyploid scenario (genome doubling in a single species), benefits from neofunctionalization would only take effect in the long run. In contrast, an allopolyploid scenario would lead to an initial selective advantage, because of the fact that it brings together different physiological properties and isolates the newly formed lineage sexually.⁴

Elucidating the causes of WGD remains a non-trivial task, as such events are often ancient and leave few traces, since over the course of evolution polyploidy is usually not maintained.^5–7 For the same reasons, dating these events becomes challenging. One type of remaining traces in available genomes are conserved paralogs (i.e. ohnologs) which are defined as descendants of an ancestral gene derived from a duplication event. Paralogs are also referred to as ohnologs when resulting from WGD.

Finding homologous relationships (genes derived from the same ancestor) is a key concept in evolutionary biology used by researchers to describe genetic differences and similarities between different species. One way to categorize genes is by inferring so-called hierarchical orthologous groups (HOGs). HOGs are defined as “sets of genes that have descended from a common ancestral gene in a given ancestral species”.⁸ These nested groups of orthologs and paralogs, defined at each taxonomic level, can be used to reconstruct ancestral genomes and infer gene families. This also provides information on evolutionary events such as gene duplications, gains or losses between known or reconstructed genomes. These inferred relationships play a critical role for this study.

Over the past few decades, multiple WGD events on the tree of life have been detected.⁹ In early vertebrate evolution, Susumu Ohno proposed that two rounds of complete genome duplication occurred (often referred as 2R). There is now strong evidence confirming this hypothesis.¹⁰ Additionally, a consensus has been reached in favor of a third round (3R) WGD occurring in the fish lineage and thus also referred to as the fish-specific genome duplication.¹¹ As mentioned earlier, WGD events have also been established in other eukaryotic lineages. For instance, an ancestral genome duplication has been placed before the common ancestor of extant seed plants, followed by another in the common ancestor of extant angiosperms.¹² Moreover, multiple WGD events have been estimated among eudicotyledons and monocotyledons.¹² Over the course of evolution, another well-known case occurred just before the separation of the yeast species Vanderwaltozyma polyspora from the S. cerevisiae lineage.¹³

Although found in vertebrates, plants, fungi, and invertebrates, WGDs remain relatively rare evolutionary events that have far-reaching repercussions on speciation and biodiversity dynamics.^2,14 Being able to place them across the tree of life is crucial to our understanding of species’ evolutionary history. Additionally, polyploidy has important effects on human health. Belonging to the hallmarks of cancers, polyploidization events are found in roughly one third of human cancers. Besides cancer, genome duplications also occur in normal cells such as hepatocytes and in premalignant lesions, impacting other disease conditions.¹⁵ Thus, expanding knowledge of WGD is of great importance for cancer and other disease research.

Currently, the main existing approaches for detecting WGD using genomic data are based on the detection of large conserved syntenic regions, gene counts along phylogeny, and synonymous mutation rates.^16,17 These methods all have their own drawbacks. For instance, investigations through synteny require high-quality genome assemblies providing gene order, and might be biased by small scale duplications or remaining layers of older WGD.¹⁸ In methods based on the number of synonymous substitutions per synonymous site (Ks), signal is lost in older WGDs due to Ks stochasticity. Moreover, these methods suffer from Ks saturation which might lead to the formation of artificial peaks.¹⁹ Such limitations could be mitigated by leveraging a combination of multiple techniques when trying to infer WGD events. Finally, to cope with the rapid increase in available genomic data, scalable methods need to be developed.

This study aims to propose a novel approach, based on hierarchical orthologous groups of duplicated genes (henceforth referred to as duplication HOGs), to locate WGD events on the given phylogenetic tree of a set of organisms. The method relies on the organisms’ proteomes, i.e. protein sequence annotations, for the detection of WGD resulting from allopolyploidy or autopolyploidy. This study is focused on the yeast lineage as it contains a single well-established WGD.¹³ Additionally, yeast species have relatively small and less complex genomes, facilitating computation and analysis, thus favoring this choice. Applying our proposed method on two datasets comprising yeast species allowed proper identification of previously reported WGD events (details in the Results section). It is worth mentioning that the method could be used for newly sequenced genomes as well.

Methods

Our method uses the HOGs in addition to the species tree as inputs. In the following section, we described how one can estimate HOGs for the genome of interest. Then, in the section “Detection of whole genome duplications (WGD)”, our method to detect WGD is presented in detail. Of note, we evaluated our method using 32 yeast genomes from the Saccharomycetales clade (details in the Results section).

Orthology inference

We used the OMA Standalone software²⁰ to infer HOGs and to reconstruct ancestral genomes. To obtain HOGs, OMA Standalone performs multiple steps. First, the data undergoes quality and consistency checks. Second, all-against-all sequence alignments are computed in order to infer pairwise orthologs. We exported precomputed all-against-all protein sequence alignments from the OMA browser.^21,22 Using precomputed all-against-all alignments as the input of the OMA standalone software significantly lowers the needed computational resources. In contrast to most other orthology inference resources, during the all-against-all alignment phase, OMA compares multiple isoforms in order to identify and utilize the isoform with the best cross-species matches.²¹ Next, pairwise orthologs are used to build an orthology graph. Finally, the GETHOGs algorithm (“Graph-based Efficient Technique for Hierarchical Orthologous Groups”) utilizes this graph and the given species tree for HOG inference. If no species tree is given, OMA estimates it using a least-squares distance approach based on OMA groups²³ (Figure 1). OMA standalone was run with the selected sets of yeasts to obtain an estimated species tree in Newick format and HOGs in an OrthoXML file.²⁰

Figure 1. Quantitative duplication hierarchical orthologous group (HOG) analysis pipeline flowchart.

The main steps of the developed method are illustrated on this flowchart. (For the complete description of the method, see the Methods section).

Results

A diverse range of yeast organisms were used in this project in order to assess the performance of our WGD placement method. We chose yeast for a proof of concept because this clade is well documented and the species have simple genomes relative to other eukaryotes. These genomes were extracted from the Orthologous MAtrix (OMA) database.²⁰ This database currently contains 2,496 genomes (Nov. 2022) from across the entire Tree of Life and is continuously being updated. As part of method evaluation, we first consider a small dataset, “small yeast”, of five species including Schizosaccharomyces pombe, Kluyveromyces lactis, Zygosaccharomyces Rouxii, Vanderwaltozyma polyspora and Saccharomyces cerevisiae. In this dataset, the well-known model organism S. pombe acts as a rooting species. It is known that K. lactis and Z. rouxii are potential candidates of species that might have led to the WGD in the hybridization theory⁴ and thus branching off before the known WGD placed before V. polyspora’s branch. In other words, Vanderwaltozyma polyspora and Saccharomyces cerevisiae are known to have undergone polyploidization.⁴ Another dataset, “large yeast”, was used as proof of concept and is composed of 32 species from the Saccharomycetales order found on the OMA database in addition to the species Penicillium chrysogenum for rooting.

As described in the method section, once HOGs were inferred with OMA standalone, the estimated number of gained, duplicated, retained or lost genes on each branch was obtained and visualized with a python package pyHam (Figure 2A,B). We notice a striking gene loss on the branch leading to the extant Saccharomyces cerevisiae South African wine strain (strain AWRI796), which may be due to its domestication. However, many of the gene losses might be false positives, which can happen when the extant genome is of lower quality or incomplete. Branches leading to Kluyveromyces lactis and Zygosaccharomyces rouxii have slightly less duplication than the other branches, but this data on its own is not sufficient to discriminate against possible WGD candidates (Figure 2A,B).

Figure 2. Hierarchical orthologous groups (HOGs) reveal evolutionary events in the yeast lineage.

(A) Inferred phylogenetic tree of the “small yeast” dataset using OMA standalone pipeline. Numbers in black correspond to branch length and numbers in purple to the branch id found in the bar plots. Colored stacked bars on branches illustrate the non-adjusted number of gained, duplicated, retained or lost genes from one node to another (these species have roughly 5000 genes in average). The same information can be found as a barplot on panel B. The red-colored branch corresponds to the placement of the whole genome duplications (WGD) detected by our method. (B) Bar plots represent the number of gained genes, duplication events, duplication HOGs, retained genes, and lost genes for the “small yeast” dataset. Branch ids correspond to the ids found on the phylogenetic tree of panel A. (C) Same bar plots as panel B adjusted by the distance (branch length plus pseudocount of one) and the average genome size for each branch (id).

However, when counting the number of duplication HOGs adjusted by genome size and branch length, a distinctive larger relative quantity (seen as an outlier on the distribution) appears on the branch of the known WGD (branch id 5; Figure 2C). The significance threshold is determined by adding the third quartile to the multiplication of the interquartile range by 1.5 (Figure 3A). This allows clear placement of the WGD event on the species tree. Note that this significant value remains with a pseudocount of five (see Figure 3B). Indeed, the next highest quantity (branch id 8) lies at only 52% of this significant value.

Figure 3. Distinctive peak in accordance with the known whole genome duplications (WGD) branch.

(A) Bar graph of the relative (by min-max normalization) adjusted (by branch length plus one and average genome size) number of duplication hierarchical orthologous groups (HOGs) for each branch of the “small yeast” dataset (see Figure 2 for branch id correspondence). On the right side, the box plot represents the distribution of these values. The distribution outlier is considered as a significant value and thus as a WGD candidate. Red colored bars correspond to branches that are significant over a range of different pseudocounts. (B) Replication of panel A barplot with a pseudocount of five. (C) Adjusted number of duplication HOGs for each branch of the “large yeast” dataset. Pseudocount added to the distance changed to five. Red colored bars correspond to branches that remain significant across a range of different pseudocounts (e.g. 1,5 and 10). See Figure 4 for branch id correspondence.

Applying the same approach on the bigger data set containing 32 species from the Saccharomycetales order resulted in two possible WGD candidates. The first corresponds to the aforementioned WGD branch with the branch id 25 (Figure 4). The second, with an adjusted number above 0.8, corresponds to the branch leading to the extant species Pichia sorbitophila with the branch id 42. Interestingly, it was shown that Pichia sorbitophila results from allopolyploidization and has a genome of about twice the size of its ancestor.²⁷ In this data set, we can see smaller differences among adjusted HOG counts which makes it more difficult to set the threshold to distinguish WGD branches, i.e. those with the highest duplication HOGs (Figure 4).

Figure 4. Quantitative duplication hierarchical orthologous group (HOG) analysis results on the “large yeast” dataset.

(A) Inferred phylogenetic tree of the “large yeast” dataset. Established post-whole genome duplications (WGD) species under the light-green box. (B) Graph represents the number of duplication HOGs adjusted by the distance (branch length plus one) and the average genome size for each branch (id). On the right side, the box plot represents the distribution of these values. Keeping significant values that remain with higher pseudocounts (Figure 3C) correspond to found WGD candidates (branch id 25 and 42). Filled blue colored circles on the phylogenetic tree correspond to remaining significant values with a pseudocount of one.

Discussion