Unbiased generation of phylogenetic profiles for 31406 human orthogroups

Unlike the examples highlighted in Figure 1, the majority of human genes arose through duplication. Following a duplication event, genes can take different evolutionary trajectories, frequently sharing the functions of the ancestor (subfunctionalization) or acquiring new functions (neofunctionalization) (Figure 2A) (Conant and Wolfe, 2008; Kensche et al., 2008). It is therefore necessary to correctly connect orthologs in each species to either daughter or to a group of homologous human genes (co-orthologs) that reflects the common ancestor (Figure 2A). Assigning human genes to such groups (orthogroups) was a central goal of our analysis. The importance of this step becomes clear in a comparison to other widely used approaches to generating phylogenetic profiles (Figure 2B). The reciprocal best BLAST match criterion (Best Bidirectional Hit or BBH)(Altenhoff and Dessimoz, 2009) applied directly to the daughter genes makes incorrect connections in species that split off before the duplication event (Dalquen and Dessimoz, 2013), while a homology search generates identical profiles for both daughters, linking only to the ancestral function (Li et al., 2014; Tabach et al., 2013a). Separate orthogroup profiles tracking each duplication event provide the most unbiased functional predictions (Figure 2B).

Open in a separate window

Figure 2

Binary eukaryotic phylogenetic profiles for 31406 human orthogroups

(A) Schematic to demonstrate the impact of human gene history on the distribution of ortholog functions across species. Genes are represented by filled circles with colors denoting independent functional trajectories. The duplication event is highlighted using a dotted line. (B) A comparison of three different approaches to generating phylogenetic profiles (blue box=present, white box=absent) for the gene family and species tree from (A). The specific functional predictions that each profile can generate are listed to the right. BBH=Best Bidirectional Hit. Dark blue is used to denote a separate phylogenetic profile (Gene X) identified through a similarity search. (C) The top BLAST hit against human TTC21B in each species colored according to a BBH match for the corresponding protein against TTC21B (black), TTC21A (green), or any other human gene (red). Gray lines highlight inferred orthogroup thresholds. Below, resulting phylogenetic profiles (blue/white) and inferred orthogroup tree. See also Figure S2. (D) Extension of algorithm to entire genome. (E) The full hOP-map, containing binary phylogenetic profiles for 31406 orthogroups, listed in decreasing order of estimated first appearance (see Experimental Procedures). Gray lines denote key branch points. See also Figure S2.

For this purpose we needed to obtain unconstrained hierarchical orthogroups defined exclusively by relationships between human genes, on a genome-wide scale. Despite the extensive literature on orthology inference (Kristensen et al., 2011), we could not find a suitable resource meeting these requirements (see Supplemental Experimental Procedures). This fact and the observation that a simple pairwise BBH approach can outperform sophisticated tree-based algorithms, especially when large genomes and many species are involved (Kristensen et al., 2011), led us to develop the modified BBH strategy outlined below. Briefly, we used BLASTp (protein BLAST)(Altschul et al., 1990) bit scores to iteratively join genes into orthogroups and simultaneously identify a score threshold for orthogroup presence (using the longest protein encoded by each gene, see Experimental Procedures). For a given ‘query’ gene or orthogroup (TTC21B in Figure 2C), we used BLAST to identify its closest candidate homolog and the corresponding bit score in each species. We also made use of a reverse-BLAST strategy to define for each of these candidate homologs, its best hit in the human genome (BBH, Figure 2C). We then defined a threshold for identifying a true ortholog of the query as the highest BLAST bit score between the query and any candidate homolog whose BBH in the human genome (‘target’, TTC21A in Figure 2C) was not the query. The target and query were then joined into a new orthogroup and their BLAST scores pooled (TTC21A/B, tree in Figure 2C). This procedure was applied sequentially until the new threshold fell below background levels (see Experimental Procedures), leading to a unique phylogenetic profile for each orthogroup (Figure 2C, Figure S2A). Since other species branching off from the human lineage all evolve their copy of the target (duplicated) gene at approximately the same rate (Figure 2C, after correcting for the higher divergence of a subset of species- such as parasites, see Supplemental Experimental Procedures), this thresholding strategy enables a robust tracking of duplications without the need for explicit models of evolution.

Starting with the 19973 human protein coding genes (NCBI), the algorithm yielded 31406 orthogroups (Figure 2D; Supplemental Experimental Procedures for orthogroup naming conventions) and 31406 matching hOP-profiles (Figure 2E), also available on our website. This analysis generated a comprehensive one-to-one orthology dataset (useful for experiments across multiple species, Figure S2B) and, importantly, confirmed that the majority of human genes can be assigned to gene families (Figure S2C), highlighting both the utility and necessity of an orthogroup-based approach to phylogenetic profiling. We detected a significant fraction of human genes tracing back to early eukaryotes (approximately 25%), with major innovation in both vertebrates and later mammals (>50%, Figure S2D- top)(Blomme et al., 2006). Importantly, our strategy generated phylogenetic profiles for ancestors of many of these vertebrate and mammal-specific genes (Figure S2D- middle), gaining predictive power from the higher inferred loss frequencies in early-branching lineages like the fungi (Figure S2D- bottom). An additional benefit of our approach was the ability to reconstruct gene trees across all gene families, which, though restricted to the human lineage, provide a close match to the literature for well-studied examples (Berg et al., 2001; Boureux et al., 2007) (Figure S2E, S2F, Supplemental Experimental Procedures).

Creation of a genome-wide pair-wise phylogenetic co-occurrence matrix

The next major challenge was to define a scoring metric that effectively extracts the predictive value of similarity between pairs of phylogenetic profiles, and then to generate a pairwise phylogenetic co-occurrence matrix between all human orthogroups (hOP-matrix, Figure 3A). Existing methods for comparing phylogenetic profiles fall into two broad categories: linear metrics based on correlation or mutual information and tree-based methods (Kensche et al., 2008). Linear metrics do not account for the interdependence of related species, leading to distorted similarity scores (Figure 3A), while tree-based methods do not scale well and also require the specific fitting of gain/loss models that must also account for noise in homology measurements. A third option, combining the strengths of both approaches, is to use joint binary transitions between presence (1s) and absence (0s) in a pair of linear profiles ordered in accordance with a 1D projection of the species tree (keeping closely related species adjacent) as a proxy for shared loss and gain events (Cokus et al., 2007). Based on the premise that more shared loss events would support a stronger case for co-evolution (Figure 3A), we defined a phylogenetic co-occurrence score (PCS) as the sum of shared transitions weighted to distinguish matches of higher and lower confidence (black bars in Figure 3A) and included a penalty for mismatches (see Experimental Procedures).

Open in a separate window

Figure 3

Generating a genome-wide phylogenetic co-occurrence matrix

(A) Schematic for generating a genome-wide pairwise co-occurrence matrix for human orthogroups. The inset highlights 3 pairs of hypothetical orthogroups indistinguishable if independence of species is assumed. However the profiles in case 1 and 2 can be reconciled with likely models of a single gain event (green dot) and a single loss event (red dot) respectively while case 3 involves 3 independent loss events. We used weighted transitions (black bars of varying thickness) to construct the phylogenetic co-occurrence score (PCS, see Experimental Procedures). (B) The PCS was benchmarked against STRING (Experimental Procedures). A range of mismatch penalty (mp) values were compared to each other and the original ‘runs’ algorithm (Cokus et al. 2007) (C) The upper panel shows the cumulative fraction of co-evolving pairs confirmed by STRING for each PCS threshold (black) compared to a random bootstrap (gray), with the lower panel showing the total number of pairs at each corresponding threshold value, at mp=0.6. Red line indicates the threshold used for further analysis. (D) Table highlighting a selected list of enriched and depleted functional terms at PCS >=10 (Experimental Procedures). (E) All genes contained in orthogroups involved in interactions with PCS>=10 were grouped into 3 categories (both poorly studied, citation count <3; one poorly studied, citation count <3 for gene 1 and >10 for gene 2; all other genes; see Supplemental Experimental Procedures for further details). See also Figure S3.

We employed a database of known interactions to optimize and validate the PCS (STRING; see Experimental Procedures). We found that varying the mismatch penalty (mp) value had a strong effect on PCS behavior (Figure 3B), namely that higher stringency improved the known interaction recall (true positive) rate but came at the cost of fewer total interactions recovered (increased false negatives). We found higher stringency (mp=0.6) useful for a genome-wide exploration of functional connections (Figures 4-​-5),5), though predictions using a lower penalty value could provide better coverage for specific functional clusters (available on website). We also confirmed that orthogroup-based profiling increases predictive power over a pair-wise BBH approach restricted to genes not assigned to gene families (Figure S3A). Figure 3C plots the cumulative fraction (upper panel) and number of predicted interactions (lower panel) as a function of PCS at mp=0.6, after filters were applied to exclude orthogroups that either appeared too recently or contained too many genes to produce useful functional predictions (Figure S3B and Supplemental Experimental Procedures).

Open in a separate window

Figure 4

Unbiased clustering of human genes into co-evolving modules

(A) Schematic illustrating the module expansion strategy (see Experimental Procedures). (B) Histogram representing the final sizes of 1074 modules. (C) Modules linked to the Fanconi Anemia (FA) pathway. Red stars indicate predicted interactions. (D) Schematic illustrating aspects of the FA pathway with genes in hOP-modules color-coded according to (C). (E) Module linked to splicing, with each orthogroup annotated to specific splicing complexes. See also Figure S4.

Open in a separate window

Figure 5

Global analysis of hOP-modules

(A), (C), (D) The first two principal components obtained from the PCA of consensus profiles corresponding to modules with 3 or more components (334 hOP-modules, Experimental Procedures) plotted against each other superimposing module size (A), functional enrichment (hypergeometric, FDR<0.1) in green (C), or certain highlighted functional categories (D). (B) Examples of modules with consensus profiles, also marked in (A). (E) Examples of uncharacterized hOP-modules, also marked in (C). (F). Examples of functionally enriched hOP-modules, also marked on the map in (D). Red stars indicate refractory genes. See also Figure S4; Table S2 for details of all hOP-modules.

We found that that the strongest co-evolving pairs (2101 unique genes, PCS>= 10) were strongly enriched for large protein complexes, metabolic pathways and some organelles but devoid of genes involved in canonical signaling, immune responses and transcriptional control (Reactome Pathways, see Figure 3D and Table S1), closely mirroring trends in bacteria (Campillos et al., 2006). This observation argues for a generalizable tendency for cellular networks with strong internal coupling (protein complexes and metabolic pathways) to form evolutionary modules over interlinked signaling and transcriptional pathways. Despite these constraints, the finding that the set of strongly co-evolving pairs includes more than 700 genes for which little or nothing is known about its function (Figure 3E and Figure S3C) highlights the power of our approach as a discovery tool. We also noted a strong enrichment for genes involved in monogenic diseases (hypergeometric p-value <1×10⁻⁴, Supplemental Experimental Procedures) at the same PCS threshold, highlighting the potential for identifying new functional connections relevant to genetic disorders.

Unbiased genome-wide assignment of co-evolving gene pairs to functional modules

We identified co-evolving modules (hOP-modules) in an unbiased genome-wide approach starting with pairs of proximal orthogroups to which, in sequential steps, a single additional orthogroup was added as the top hit from the weighted PCS average against existing module members (Figure 4A) until a PCS threshold was reached (see Experimental Procedures). The first striking observation was that the size distribution was skewed heavily towards smaller modules with a large number of single pairs failing to pick up additional orthogroups (Figure 4B). These modules often map to one or more sections of highly conserved cellular pathways and multi-protein assemblies (Figure 4C-E), highlighting an unanticipated degree of modularity within these pathways (Figure 4D) and suggesting novel connections (red stars in Figure 4C; for example, SPIDR is a recently identified scaffold protein that has been linked to DNA repair (Wan et al., 2013) but not directly to the Fanconi Anemia pathway). It is worth noting that, in cases like the splicing-related module #41 (Figure 4E) where genes are largely annotated purely based on homology or high-throughput proteomics (Cvitkovic and Jurica, 2013), co-evolution can serve as an independent confirmation of shared function.

When available, orthogonal sources of data helped link the observed evolutionary modularity to its underlying causes. Analyzing a curated compendium of splicing genes (Hegele et al., 2012; Wahl et al., 2009) revealed two clear evolutionary trajectories for ancient splicing genes (Figure S4A), one that is persistent with a few losses in protists, and a second that is characterized by a gradual pruning of the splicing machinery in fungal lineages (Figure S4B-C). These observations parallel previous work on intron losses in fungi (Nielsen et al., 2004) and explain the wide range of complexes observed in hOP-modules enriched for splicing (Figure 4E).

We took advantage of our unbiased approach to ask if a global analysis of hOP-modules would reveal general principles underlying the modular co-evolution of human genes in addition to generating a set of high-confidence functional predictions. We applied principle component analysis (PCA) to ‘consensus profiles’ (Figure 5A-B, Experimental Procedures) for 334 hOP-modules with 3 or more components. The first two principal components (Figure S4D) were dominated by the evolutionary age of hOP-modules and overall loss frequency (examples from Figure 5B define the corners of the PC space in Figure 5A, 5C and 5D). This analysis revealed clear constraints upon the potential range of phylogenetic patterns and also separates out large, poorly-resolved hOP-modules corresponding to mammal-specific (box 1 in Figure 5A) and paneukaryote modules (box 2 in Figure 5A) with a shortage of informative loss events. Analysis of functional enrichment (Reactome/COMPLEAT, Figure 5C) helped identify completely uncharacterized modules (Figure 5E) that may reflect undiscovered cell functions, as well targeted predictions for refractory genes clustering within functionally enriched hOP-modules (Figure 5F). Modules belonging to related functional categories often occupied neighboring areas of the map (Figure 5D). Nevertheless, a large number of consensus profiles were isolated, belonging to modules with unique loss patterns often linked to metabolic functions (Figure 5D, Table S2). Modules occupying sparser areas can be analyzed with a lower threshold to discover more refractory genes without blurring functional boundaries (Figure S4E, S4F). Finally, modules enriched for more than one function may imply novel connections between them, such as Module #18 that is linked to two apparently distinct metabolic pathways, tyrosine/phenylalanine catabolism and the peroxisomal degradation of branched fatty acids (Figure 5F). Interestingly, older studies suggest that these two pathways might be linked (Vamecq and Van Hoof, 1984).

Identifying novel interactors of the WASH complex

As a first proof of principle, we focused on the regulation of Arp2/3, an ancient multi-protein machine that nucleates actin filaments and organizes them into branched networks (Campellone and Welch, 2010). Arp2/3 can be activated by WASP-WAVE-WASH-family nucleation-promoting factors (Class I NPFs, Figure 6A-B). The Arp2/3 complex as well as the WAVE and WASH nucleation-promoting complexes emerged from our analysis with distinctive loss patterns (Figure 6A) (Kollmar et al., 2012). The WASH complex was of particular interest to us because it is incompletely characterized and we identified promising candidates in its corresponding hOP-module and neighborhood (Figure 6C, Supplemental Experimental Procedures): CCDC22 has been identified in pull-downs with FAM21A (Harbour et al., 2012) and COMM Domain (COMMD) proteins (Starokadomskyy et al., 2013); DSCR3 is an ortholog of the WASH-interacting retromer component VPS26A (Gomez and Billadeau, 2009); the poorly-studied RAB21 localizes to early endosomes (Simpson et al., 2004). The list also provides a close match with very recent predictions from an orthogonal approach (Li et al., 2014).

Open in a separate window

Figure 6

Investigation of genes with a predicted link to WASH complex function

(A) Consensus profiles for the Arp2/3 complex, SCAR/WAVE complex and WASH complex (module ID on the left). (B) Schematic highlighting the role of branched actin in various cellular locations and specificity provided by the different Class I NPFs, with inset showing recruitment of Arp2/3 followed by nucleation of a new branch. (C) hOP-profiles from module #6 and extended neighborhood (Supplemental Experimental Procedures) containing WASH complex components (gray) and candidates (red stars). (D), (E) mTurquoise-tagged DSCR3 or RAB21 (blue in merge) were co-expressed with mCitrine-tagged KIAA1033 (green in merge) in HeLa cells also immunostained with an antibody to EEA1 (D, red in merge) or p62 (E, red in merge) and imaged at 63X on a confocal microscope. Images are maximum intensity projections. In the DSCR3 images, arrows highlight characteristic bright co-localizing EEA1-negative foci. In both cases, white boxes represent magnified regions on right and bottom. Scale bar=10 μm. (F) HeLa cells treated with siRNA for 48 hours were subjected to acute stress (10 μM Bortezomib for 30 minutes), and immunostained for pericentrin (green) and p62 (red) and labeled with Hoechst (blue) to mark nuclei. Sample images from 3 treatment conditions obtained at 20X are shown. Scale bar=20 μm. (G) Mean p62 puncta intensity per cell (Supplemental Experimental Procedures) for each condition. Bars represent the mean and standard deviation of two replicates with >1000 cells per replicate from a single experiment (experiments could not be pooled due to the variable fraction of control cells expressing p62, but trends were reproducible across >3 independent experiments). Red stars mark a significant deviation from control (p-value<0.05, estimated using multiple pairwise comparison after one-way ANOVA). (H) The average number of puncta per cell, analyzed as in (G). See also Figure S5.

The WASH complex has been implicated in regulating endosome morphology (Derivery et al., 2009), endosome-to-Golgi transport (Gomez and Billadeau, 2009) and autophagic flux (Xia et al., 2014; Zavodszky et al., 2014) in mammalian cells with data from Dictyostelium (Park et al., 2013) suggesting that the autophagy role might be conserved. To test for such a potential link, we co-expressed fluorescently tagged candidates DSCR3 and RAB21 together with the WASH complex component SWIP/KIAA1033, followed by immunostaining for early endosomes (EEA1) or autophagosomes/aggresomes (p62/SQSTM1). Both DSCR3 and RAB21 colocalized strongly with the WASH complex marker in puncta that partially overlapped with both compartment markers (Figure 6D, 6E). Inspection of larger EEA1-positive rings revealed the WASH marker also in small surface domains, consistent with previous reports (Figure 6D, (Derivery et al., 2009)) while RAB21 was excluded from these domains (Figure 6D, 6E). While multiple COMMD proteins also colocalized strongly with the WASH marker (Figure S5A-B), we found that, unlike the other candidates, they also aggregated when expressed at high levels (data not shown) (Vonk et al., 2014) and did not analyze them further. We confirmed that knockdown of WASH1, which destabilizes the entire complex, led to a visible collapse of the endosomal membrane network (Figure S5C) and also a strong increase in the number and intensity of p62 autophagosome-like puncta visible after acute stress (Figure 6F-6H). Interestingly, knockdown of DSCR3 or RAB21 phenocopied the latter effect (Figure 6F-6H), indicating a likely parallel disruption of autophagic flux. While correlative, these experiments imply a shared function with the WASH complex and provide a starting point for future investigation.

Generation of hOP-profiles

The species used for phylogenetic profiling were restricted to high-quality verified RefSeq genomes to prevent propagation of annotation errors. We first identified the top-scoring homolog (BLASTp bit scores) for the longest protein encoded by each human gene (19973 unique genes, NCBI, March 2013) in each fully annotated genome (177 eukaryotic species, NCBI, March 2013), as well as the reciprocal best match in the human genome for each contributing protein in another species (BBH). In the first step of our iterative algorithm, the third highest BLASTp bit score (Supplemental Experimental Procedures) against a human gene (query) belonging to a gene with a different human BBH match (target) was used to combine the target and query into a single orthogroup and simultaneously set the threshold for orthogroup presence. BLASTp scores for each new orthogroup were pooled and the algorithm iterated until a bit score threshold (50) was reached or the orthogroup size exceeded 100, resulting in 31406 orthogroups for 19973 genes. The orthogroup thresholds were then applied to a matrix of the top BLASTp bit scores (corrected for species evolving significantly faster than their neighbors) against each orthogroup to generate 31406 binary phylogenetic profiles (available in File S1, see Supplemental Experimental Procedures for details).

Generation of the hOP-matrix

A tree-aware score to compare pairs of hOP-profiles was defined as the linear sum of shared transitions between presence and absence with a positive weight for transitions involving doublets (00→11 and 11→00) and a (negative) penalty for mismatches. This score was calculated for each pair of profiles to generate a 31406*31406 sparse co-occurrence matrix. The STRING database of protein-protein interactions was used to benchmark the score and optimize the weight and penalty factors (maximizing the number of identified STRING interactions as a fraction of the total number of co-evolving pairs). Orthogroups first appearing in the vertebrate branch and orthogroups containing more than 4 genes (poor information content) were excluded for the global analyses of pairs and modules, leaving a total of 14412 orthogroups (scores in File S1, see Supplemental Experimental Procedures for details).

Benchmarking and enrichment statistics

We used the STRING protein-protein interaction database (http://string-db.org/, version 9.1, (Franceschini et al., 2013)) to optimize and benchmark the co-occurrence metrics between hOP-profiles. P-values for functional enrichment were estimated using the hypergeometric test and reported using the false discovery rate (FDR). Terms and associated annotations were obtained from the Reactome database (http://www.reactome.org/, version 48 downloaded 07/ 2014 and the COMPLEAT database (http://www.flyrnai.org/compleat/, downloaded 06/2014, (Vinayagam et al., 2013)). See Supplemental Experimental Procedures for details.

Generation of hOP-modules

Agglomerative modules were created in a stepwise fashion starting with high-scoring seed pairs of orthogroups (PCS>=5 for full coverage of the map). At each step, the top-scoring orthogroup from the weighted average of co-occurrence scores between existing module components and the rest of the co-occurrence matrix was added to the module, and the iterative process halted at PCS=5. To maximize the growth opportunity for high-confidence co-evolving seed pairs, modules were created sequentially starting with the strongest seed pair and module components removed from the general pool. See Table S2 and Supplemental Experimental Procedures for additional details.

Consensus profiles and PCA

Binary consensus profiles were generated for each hOP-module by assigning a 1 to each species with an ortholog in 50% or more of the member hOP-profiles and a 0 otherwise. The 334 consensus profiles corresponding to hOP-modules with 3 or more components were then subjected to a principal component analysis (PCA) implemented using inbuilt MATLAB functions. The first two principal components were used for all subsequent analysis.

Cell culture, siRNA and expression constructs, antibodies

HeLa, NIH3T3, and Hs68 cells were cultured in 10%FBS/high-glucose DMEM/PSG medium (GIBCO). For siRNA experiments or transient transfections, cells were cultured on 96-well plastic-bottom (Corning 3904) or glass-bottom (In Vitro Scientific, P96-1.5H-N) plates. siRNA reagents were obtained from Qiagen (pools of 3). Expression constructs were derived from the human ORFeome collection (version 5.1, http://horfdb.dfci.harvard.edu/hv5/) or from PCR. Antibodies against p62/SQSTM1 (mouse, used at 1:400, BD Biosciences 610832), acetylated tubulin (mouse, used at 1:3000, Sigma T6793), pericentrin (rabbit, used at 1:800, Abcam ab4448), and EEA1 (mouse, used at 1:400, BD Biosciences 610457) were employed for immunofluorescence. See Supplemental Experimental Procedures for detailed protocols and Table S5 for siRNA specifications.

We thank Arnold Hayer, Amy Winans, Steve Cappell and Chad Liu for critical reading of the manuscript, Arnold Hayer for providing GATEWAY destination vectors, Mingyu Chung for a MATLAB nuclear segmentation routine and members of the Meyer lab for helpful discussions and feedback. This work was supported by NIH grants GM063702 and GM030179, the Stanford Center for Systems Biology, and the Stanford Graduate Fellowship in Science and Engineering (G.D.).