Novel ACE2 protein interactions relevant to COVID-19 predicted by evolutionary rate correlations

Taxon selection and data collection

Our evolutionary rate correlation (ERC) approach requires orthologous protein sequence data across a large number of taxa with well-defined phylogenetic relationships. Calculation of evolutionary rates requires a resolved phylogeny of the taxa analyzed that is scaled to evolutionary time. Our ERC calculations utilize the TimeTree (Kumar et al., 2017) to generate a time-scaled phylogenetic tree using the mammalian taxa that are represented in OrthoDB sequence data (Fig. S1). The tree generated is in units of millions of years and is based on a compilation of many phylogenetic-dating related studies. The tree is utilized as a base topology in phylogenetic analysis and its branch lengths are used to measure time for calculating evolutionary rates from the resultant individual protein trees (Fig. 1). Additional details on the data set are provided in the Supplemental Text. Large files are deposited in FigShare (https://doi.org/10.6084/m9.figshare.14637450) and are listed in Table S1.

Well-defined orthologous sequence data is sourced from OrthoDB (Kriventseva et al., 2019) at the “mammalia” (taxonomic id: 40674) taxonomic level. Since OrthoDB sequence data is gathered from a variety of sources and clustered algorithmically (unsupervised), primarily based on sequence similarity (Kriventseva et al., 2015), related paralogous proteins are often clustered with each other even if canonically annotated as functionally distinct proteins (Table S2). Additionally, since the data sources for sequences can have varying levels of completeness, most ortholog groups on OrthoDB are missing sequence data for one or more taxa represented in the database. So, a majority of the data we initially selected was from single-copy ortholog groups with at least 90 of the 108 possible mammalian taxa present. In addition, some proteins with a possibly relevant function to COVID-19 pathologies (such as XCR1, and IFNAR2) or other relevant pathways in ortholog groups that did not meet the initial criteria, but that had minimal paralogy issues, were included. Paralogous sequences were manually disambiguated based on published protein annotations and phylogenetic analysis. If a taxon in a given sequence had duplicate sequences that could not be disambiguated, the taxon was excluded in phylogenetic and ERC calculations for the specific proteins involved. In total, 1,953 orthologous protein groups are used in analyses.

Phylogenetic calculations and protein alignments

To prepare orthologous sequence data for ERC calculation, each set of protein sequences are first aligned using the MAFFT software package (Katoh & Standley, 2013) using the following arguments: “––maxiterate 1000 ––localpair ––anysymbol”. Since the sequences come from data sources with varying levels of quality and multiple alignment programs can be imperfect, the aligned sequences must then be trimmed. The alignments are trimmed using the trimAl software package (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) using the “-automated1” argument to remove poorly aligned regions. These final prepared alignments are then used to generate maximum-likelihood phylogenies. The IQ-TREE software package (Minh et al., 2020) is used to estimate protein branch lengths (equivalent to average substitution counts per site). Specifically, the “LG+F+G+I” model (which utilizes an empirically derived amino acid substitution matrix) is used with the following additional parameters: “-B 1000 -st AA -seed 1234567890” and the TimeTree phylogeny is provided to constrain output tree topology to reduce possible branch length estimation errors with the “-g” option. These trees are the basis of ERC calculations. Protein branch lengths are based on the average number of changes in amino acids at each residue in the alignment. The resultant branch lengths are paired with corresponding branches in the TimeTree to quantify branch-specific rates to be used for ERC calculations (described below). ERCs calculated with the more complete phylogeny (108 species) had short branch problems in oversampled taxonomic groups (described below and in the Supplemental Text). We therefore used a reduced phylogeny consisting of 60 taxa for subsequent ERC analyses.

Calculation of ERCs

Our evolutionary rate correlation (ERC) method is designed to predict protein-protein interactions using evolutionary data (Yan, Ye & Werren, 2019), and is based on protein evolutionary rates on terminal branches of the mammalian phylogeny (Fig. 1). We found that the more complete phylogeny (108 species) had short branch problems that inflate ERC spearman rank correlations (discussed in Supplemental Text). Most notably, there was an association between branch time and protein rate for many proteins, with oversampling in some taxonomic units (e.g. in Primates and Rodentia) leading to many ERCs being driven by relatively short branches (Supplemental Text). We attempted to control for these effects initially by using a partial correlation method, but found that it was not sufficient due to correlations between residuals and branch time (Supplemental Text). We then removed taxa that contributed short branches in our phylogeny based on either a 20MY or 30MY divergence time threshold (Supplemental Text) and recalculated branch rates for all proteins. We found that the 30MY threshold short branch removal eliminated significant branch time to protein rate correlations for the majority of proteins (87.5%). The resultant rate data no longer has branch time to branch rate as a confounding cofactor, and the ERCs themselves are no longer biased by extremely short branches and taxonomic oversampling (Supplemental Text). The resulting data set is composed of 60 taxa and is used in our subsequent analyses of ERCs.

Using the adjusted data set, ERCs are calculated for every possible combination of protein pairs for which a tree has been generated. Every protein pair for which an ERC is calculated has each respective tree and the TimeTree topology is pruned to only include the shared taxa between the two, using the “ETE3” Python package (Huerta-Cepas, Serra & Bork, 2016). ETE3 is also used to extract the terminal branch lengths of each pruned tree. Evolutionary rates are calculated by dividing the terminal protein-tree branch lengths (average substitutions per site) by the corresponding branch in the TimeTree (measured in millions of years). Terminal branches are used for calculations as they do not have shared evolutionary histories and are therefore independent. The resulting rates have the unit of average substitution per site per millions of years. Given the resultant rates, evolutionary rate correlations are then calculated by performing a Spearman’s rank correlation test (Yan, Ye & Werren, 2019) using the Python package “SciPy” (Virtanen et al., 2020).

Multiple test corrections

P-values are corrected using the Benjamini–Hochberg FDR multiple-test correction procedure implemented in the Python package “statsmodels” (Seabold & Perktold, 2010). The FDR correction is applied to each respective protein’s set of ERCs. Correlation test results are non-directional, but FDR corrections are dependent on the rank of each correlation’s p-values. Since the rank of each correlation test value on respective protein lists vary, the FDR-corrected p-values of a given protein pair can differ depending on the specific ERC set. An ERC is considered significant if the FDR-corrected p-value is less than 0.05.

ERC set enrichment analysis

To summarize the common biological function of proteins that tend to have strong ERCs, gene set enrichment analysis is performed on the top 2% of ERCs (by ρ) of each protein(s) of interest (including the protein itself), including only proteins with ERCs that are significant following an FDR correction at a significance level of 0.05. At most, a protein of interest will have 41 proteins included for enrichment analysis (2% of the total 1,953 proteins plus itself). Protein set enrichment analyses are performed using the Enrichr service (Xie et al., 2021) via the Python bindings provided by the “GSEApy” Python package (Fang et al., 2021) given the background of the full set of 1,953 proteins. We calculate enrichment results for ACE2 and all of its top 20 ERC partners. Additional enrichment analyses were also performed on a case-by-case basis based on relevance, including the reciprocal rank networks. Enrichments are performed using selected relevant term databases: KEGG_2019_Human, GO_Biological_Process_2018, GO_Cellular_Component_2018, GO_Molecular_Function_2018, Reactome_2016,WikiPathways_2019_Human, Tissue_Protein_Expr_from_Human_Proteome_Map, Tissue_Protein_Expr_from_ProteomicsDB, and Jensen_TISSUES.

Enrichment results for terms that are significant at FDR-adjusted p < 0.05 for all analyses are placed into a single table, organized by the enrichment term database (File 3). The outputs from different databases can contain redundant terms to each other, so only the most significant of the redundant terms are reported for any enrichment analysis in the main text.

Reciprocal rank network (RRN) generation

To evaluate and visualize the strongest ERCs centered around proteins of interest, “reciprocal rank networks” (RRNs) are produced. Reciprocal ranks refer to the fact that a significant ERC between two proteins can have different ranks in the two respective protein ERC lists because some proteins have more and higher ERCs than others. To focus on networks of proteins with strong reciprocal rank correlations, we have constructed networks based on proteins with reciprocal ranks of 20 or less (RR20), which is the top 1% in each protein’s highest ERCs based on ρ values. Specifically, we have developed an ACE2 centric reciprocal rank network by the following steps (1) for ACE2, select its top 1% (20) proteins, (2) for each of those proteins, select additional proteins in their ERC list with reciprocal rank 20 or less, and then (3) Given the core set of proteins generated in the previous two steps, connect proteins which have a unidirectional rank of 20 or less.

The resultant network represents the strongest ERCs centered around a protein of interest (in this case ACE2), along with the immediate neighborhood of the strongest ERCs surrounding the protein of interest. The ACE2 Core Reciprocal Rank (CRR) was initiated with the four proteins to which ACE2 has RR20 ERCs (CLU, TMEM63C, FAM3D, and L1CAM), with GPR141 added due to its RR1 strong connection to CLU and unidirectional connection to ACE2. ACE2 also has highly significant ERCs to proteins that do not rank it in their top one percent. Therefore, a similar approach has been used to generate an ACE2 reciprocal network initiated with the top 10 proteins to which ACE2 has highly significant ERCs, but are not reciprocally RR20 ranked, with a subsequent one cycle RR20 built upon these. This ACE2 Unidirectional Reciprocal Rank Network (URR) contains strong network connections to ACE2 through its high unidirectional ERCs. Steps 2–3 were omitted as the network becomes extremely large following just the first step, and our focus is on examining close connections to ACE2 based on ERC analysis.

ERCs within and between protein complexes

To compare whether calculated ERCs are stronger between known interactions versus non-interactions, the protein complex database, CORUM (Giurgiu et al., 2019), was used to retrieve known complexes. The “Core Complex” dataset was downloaded and filtered for human complexes to eliminate redundancy, resulting in 233 protein complexes from this CORUM data set which have two or more components present in our 1,953 protein ERC set, representing 258 pairwise ERC comparisons. As these protein complexes have redundancy (i.e. some complexes contain overlapping protein pairs), the set was further restricted to complexes containing unique protein components—resulting in 139 effective unique complexes considered. To test whether ERCs within complexes are higher than between complexes, all pairwise ERCs within complexes were compared to the median ρ value for each pair to proteins present in non-redundant CORUM set that are not in complex with either of these proteins. A Wilcoxon matched signed-rank test was performed using the “wilcox.test” function in base R (version 3.6.1; with parameters “paired” and “exact” set to “TRUE”) on the in-complex ρ values and the median out-of-complex ρ values, to test if the in-complex ρ values were significantly greater than the median out-of-complex ρ values. In addition, as there were many complexes with a majority of subcomponents not present in our 1,953 datasets, the likelihood of individual pairs directly interacting within the complex decreases with the increasing number of proteins in a complex. Therefore, an additional Wilcoxon matched signed-rank test was performed on members of protein complexes composed of five or fewer proteins.

Testing for taxonomic effects

We use three methods to test for taxonomic effects on the calculated ERCs, (1) multiple linear regression, (2) analysis of covariance (ANCOVA), and (3) non-parametric independent contrasts (Garland, Harvey & Ives, 1992). For the regression and ANCOVA approaches, rate data was grouped by mammalian taxonomic orders accessed via ETE3 (Huerta-Cepas, Serra & Bork, 2016) and treated as an independent variable. The independent contrasts test uses the mammalian topology previously created with TimeTree (Kumar et al., 2017) to generate independent contrasts within the phylogeny. Statistical tests for each method are performed using base R (version 3.6.1). See the Supplemental Text for details.

Testing whether branch rates increase when extending branch time within clades

To test whether increasing branch length results in increasing protein evolutionary rate, we selected separate phylogenetic groups (clades) from the full phylogeny (Fig. S1) that contain short branch lengths. Protein evolutionary rate was calculated for each protein on the short branch, and then sequentially recalculated after removing adjacent taxa to extend the branch internally (Fig. S5). In this way, the protein evolutionary rate was examined as branches are extended internally in independent clades within the tree. Comparing original branches to the 20MY correction resulted in 12 clades for which time scales change between 20MY and 30MY corrections, and 16 clades for which time scales change between 0MY and 30MY. Tests on each branch’s rate against the respective adjusted rate were performed using two-tailed Wilcoxon Matched Signed Rank Tests (Base R v3.6.1), first for proteins of interest (e.g. ACE2) and then for the full protein set. Results are described in the Supplemental Text.

Basic approach

The basic methods are outlined here to provide context for the results which follow. To identify candidate protein interactions using evolutionary rate correlation, we utilized the consensus TimeTree phylogenetic reconstruction for mammalian species (Kumar et al., 2017). A total of 1,953 proteins (including ACE2) were aligned and evolutionary rates for each protein were then calculated for terminal branches of the tree (Fig. 1). This was determined by dividing the protein-specific branch length on each terminal branch by terminal branch time from the consensus tree (Yan, Ye & Werren, 2019). Maximum likelihood branch lengths were estimated in IQ-TREE (Minh et al., 2020) using an empirical amino acid substitution matrix (see methods for details). To investigate evolutionary rate correlations (ERCs) among proteins, Spearman rank correlations were calculated for every protein pair using terminal branch rates (Fig. 1). Due to the large number of comparisons, a Benjamini–Hochberg false discovery rate (FDR) correction was calculated for each protein’s ERC set (significance threshold α = 0.05). We subsequently found that many proteins show a positive correlation between terminal branch time and evolutionary rate, and observed that short branches in relatively oversampled taxa significantly contributed to this correlation (Supplemental Text). We, therefore, removed species that accounted for short branches, which eliminated the protein evolutionary rate to branch time correlation (see Methods and Supplemental Text for details). ERCs were then recalculated, and our ERC analyses are based on this set of 60 taxa.

In addition, we tested whether the observed lower rates of protein evolution for short terminal branches in the phylogeny are due to rates actually increasing over evolutionary time, versus a taxonomic effect. This was accomplished by examining changes in protein rates in independent clades as terminal branches were effectively extended by selective removal of flanking taxa. The analysis shows that evolutionary rates for many proteins increase as branch length is increased (described in more detail in Supplemental Text, Fig. S4). A possible explanation for the pattern is that protein coevolution is mostly episodic, and short branches in a phylogeny are less likely to capture such events. In additional analyses, we tested for but did not find significant confounding effects of taxonomy on the ERC results (Supplemental Text).

Our analyses are focused on candidate protein interactions involving ACE2 using evidence of highly significant ERCs. For this purpose, we first examine proteins in ACE2’s highest 2% of ERCs (top 40 proteins), all of which are highly significant after FDR correction (Table 1). Some of these ACE2 ERC proteins have been previously implicated in severe COVID-19 or SARS-CoV-2 gene expression effects on infected cells. However, while they have not been previously identified as having protein interactions with ACE2, this is predicted by our ERC analysis.

X-C Motif Chemokine Receptor 1 (XCR1) provides an illustrative example. XCR1 is a cytokine signaling receptor and ACE2’s 2nd highest ranked ERC, with a highly significant evolutionary rate correlation. XCR1 is in a small genomic region that is implicated in severe COVID-19 by genome-wide association studies (Severe Covid-19 GWAS Group, 2020; Fricke-Galindo & Falfán-Valencia, 2021). Another example is Interferon alpha/beta receptor 2 (IFNAR2) which, in a genome-wide association study (GWAS) and multi-omic analysis by Pairo-Castineira et al. (2021), was implicated in severe COVID-19. We therefore added it to our analysis, and surprisingly found it to be highly ranked (5th) among ACE2 ERCs. Clusterin (CLU) is the 3rd strongest ERC of ACE2 and the ACE2-CLU pair show high reciprocal ranks to each other (3rd in ACE2’s set, 8th in CLU’s set). CLU prevents the aggregation of misfolded proteins in the blood and delivers them to cells for degradation in lysosomes (Sánchez-Martín & Komatsu, 2020). CLU connects to key proteins in the coagulation pathway based on its reciprocal rank network (“ERC Reciprocal Rank Networks Implicate Coagulation Pathways and Immunity”, Fig. 2). CLU has been implicated in coronavirus infections, as one of only two proteins showing significant expression changes in cells infected by three different coronaviruses tested, including SARS-CoV-2 (Singh et al., 2021). The examples above lend credence to the proposition that the ERC approach is detecting ACE2 protein interactions that have implications to COVID-19.

Differences in ERC rank between protein pairs for the same correlation can occur because some proteins have higher and more extensive ERC connections than others. As a result, while two proteins can have a significant ERC with each other, each one’s rank may differ in their respective ERC lists, as illustrated for ACE2 and GEN1 (Table 1). GEN1 (Flap endonuclease GEN homolog 1) is ACE2’s top-ranked ERC, and is a DNA nuclease whose primary functions are resolution of DNA Holliday junctions and DNA damage checkpoint signaling (Chan & West, 2015). This protein shows high ERCs and is ranked highly in the ERC sets for many other proteins, suggesting central connectivity. As described further in “The ACE2 Unidirectional Reciprocal Rank (URR) Network”, GEN1 shows unexpected enrichments for immune functions, perhaps related to its role in DNA damage checkpoint signaling.

Because our focus is on identifying strong candidate interactions involving ACE2 and its predicted partners, we utilize the rank information to identify proteins with high reciprocal ranks. Specifically, we focus on the strongest reciprocal ranks (RR) defined by ranks of less than or equal to 20 (RR20), which is the highest one percent of each protein’s ERCs, and use these to develop reciprocal rank networks (“ERC Reciprocal Rank Networks Implicate Coagulation Pathways and Immunity”). Although speculative, we posit that protein pairs with high reciprocal ranks are likely to be strongly coevolving (i.e. both partners evolving reciprocally due to selective pressures acting on interacting domains between them). In contrast, protein pairs with a significant evolutionary rate correlation only one ranks highly (e.g. within the top two percent) in the ERC set of the other, are more likely to be due to “unidirectional” evolution. The rationale is that proteins with many significant ERC partners are under selective pressures primarily from their top evolutionary partners, whereas other interactors evolve primarily in response to the forces shaped by their stronger partner(s). We emphasize that this interpretation is speculative, and requires further exploration to determine what factors shape reciprocal ERC ranks between proteins.

The view that ERCs are detecting protein interactions relevant to COVID-19 is further supported by the analysis of ACE2 reciprocal rank ERC networks (“ERC Reciprocal Rank Networks Implicate Coagulation Pathways and Immunity”). Noteworthy in this regard are additional proteins in the coagulation pathway, such as Coagulation Factor V (F5), Fibrinogen Alpha Chain (FGA), Fibrinogen Beta Chain (FGB), and Fibrinogen Gamma Chain (FGG). Thrombosis (blood clotting) is a major pathology of COVID-19 (Gupta et al., 2020). Connections of ACE2 with the proteins above could relate to severe blood clotting problems in COVID-19 infections. ACE2 networks also show strong enrichments of cytokine signaling, viral (and pathogen) infections, and inflammatory response terms (File 3), which are clearly relevant to COVID-19 pathologies such as cytokine storms and systemic inflammation.

In yet other cases, we have found proteins with significant ACE2 ERCs or ACE2 network connections, but for which there is little functional information, such as GPR141. We can use their ERCs to suggest possible functions for future investigation. Finally, ERCs for proteins of known function (such as F5 and GEN1) indicate likely additional roles, suggesting these proteins have unrecognized “moonlighting” functions (Jeffery, 1999).

Below, we first describe proteins of interest to which ACE2 has significant ERCs, summarize aspects of their known biological functions, and examine significantly enriched functional categories for these ERCs. We then build and evaluate two different networks for ACE2 interacting proteins (“ERC Reciprocal Rank Networks Implicate Coagulation Pathways and Immunity”), one of which reveals connections to coagulation pathways and the other to cytokine-mediated signaling, viral response, and immunity. Finally, we discuss the potential implications of these predicted ACE2 interactions to COVID-19 pathologies and propose some specific hypotheses that emerge from this analysis.

Top ERC interactions link ACE2 to COVID pathologies

To investigate protein associations of ACE2, we first determined the protein enrichment categories for its top 2% ERC proteins (based on Spearman rank correlation coefficients, ρ) using the gene set enrichment package Enrichr (Xie et al., 2021) (Table 2). The top two KEGG_2019_Human enrichments are for complement and coagulation cascade related (FDR = 2.0E−03) and cytokine-cytokine receptor interaction related (FDR = 2.0E−03) terms. This finding is consistent with two hallmarks of COVID-19 pathology, abnormal systemic blood-clotting (thrombosis) and cytokine storms (Coperchini et al., 2020; Fei et al., 2020). Additionally, several terms related to viral/bacterial-specific infection are significantly enriched, such as Tuberculosis (FDR = 1.4E−02), HPV infection (FDR = 1.4E−02), measles (FDR = 2.4E−02) and Hepatitis C (FDR = 3.1E−02). Gene Ontology Biological Process also shows enrichment for tumor necrosis factor (TNF) pathways, including the signaling pathway (FDR = 3.9E−03) and cellular responses (FDR = 1.6E−02). Additional terms are shown in Table 2.

The ACE2 ERC analysis indicates that ACE2 is “coevolving” with proteins involved in the complement and coagulation pathways, cytokine signaling, TNF, and pathogen response pathways. Here, we summarize results and background information on some of the key proteins among ACE2’s ERCs (more extended summaries of each protein are in the Supplemental Text).

Among ACE2’s strongest ERCs are proteins involved in immunity. For example, XCR1 (X-C Motif Chemokine Receptor 1) is ACE2’s 2nd top-ranked ERC (ρ = 0.67, FDR = 6.2E−05). It is a chemokine XCL1 receptor involved in immune response to infection and inflammation (Lei & Takahama, 2012). Strikingly, the Severe Covid-19 GWAS Group (2020) detected a small genomic region containing six genes that significantly associate with severe COVID-19, one of which is XCR1. Our finding that XCR1 is ACE2’s 2nd highest ERC interactor lends independent support for a relationship between COVID-19 and XCR1. Furthermore, it suggests that an interaction between ACE2 and XCR1 could be involved in COVID-19 pathologies. To our knowledge, there are no other reports of interactions between these two proteins.

Another striking connection of ACE2 ERC to immunity is through IFNAR2 (Interferon alpha/beta receptor 2), which has a highly significant ACE2 ERC correlation (ρ = 0.62, FDR = 6.1E−04). IFNAR2 forms part of an important receptor complex with IFNAR1 (Thomas et al., 2011) involved in interferon signaling through the JAK/STAT pathway to modulate immune responses. IFNAR2 has been implicated in severe COVID-19, based on mendelian randomization, genome-wide associations, and gene expression changes (Liu et al., 2021; Pairo-Castineira et al., 2021). Our data provide independent support for a role, possibly mediated through ACE2 interactions. Interferon pathways are important in antiviral defense, but also can contribute to cytokine storms and COVID-19 pathologies (McKechnie & Blish, 2020). Other immune-related proteins with high ERC connections to ACE2 include TLR8 (Toll-like Receptor 8), FAM3D (FAM3 metabolism regulating signaling molecule D), and PLA2G7 (phospholipase A2 group VII).

Coagulation pathway proteins figure prominently in ACE2 ERC-predicted protein interactions (Table 3, Fig. 2). This is reflected both in significant enrichment for coagulation cascade proteins in the top 2% strongest ACE2 ERCs (Table 2) and the strong reciprocal rank network for ACE2 (“ERC Reciprocal Rank Networks Implicate Coagulation Pathways and Immunity”, Fig. 3). The finding has obvious potential implications to a hallmark pathology of COVID-19, systemic coagulopathy (Wright et al., 2020; Medcalf, Keragala & Myles, 2020). A list of coagulation and blood-related proteins associated with ACE2 is presented in Table 3. Among ACE2’s top 2% ERCs associated with coagulation pathway are Coagulation Factor V (F5), Protein C inhibitor (SERPINA5 aka PCI), and Thrombin Receptor 2 (F2RL2) (Table 1).

Also relevant to coagulopathy are Clusterin (CLU) and the orphan G protein-coupled receptor 141 (GPR141). The chaperone protein CLU has a soluble form that circulates in the blood and is part of the “cleaning squad” that clears misfolded extracellular proteins for delivery to lysosomes and degradation (Itakura et al., 2020; Sánchez-Martín & Komatsu, 2020). It is the 3rd highest ACE2 ERC (ρ = 0.63, FDR = 1.5E−04), and these two proteins show strong reciprocal ranks (3, 8), likely supporting biological interactions. Relevant to this point is that both ACE2 and CLU have soluble forms that circulate in the blood (Itakura et al., 2020). Of direct relevance to COVID-19 and possible ACE2-CLU protein interactions, Singh et al. (2021) found in cells infected with different coronaviruses (SARS-CoV-2, SARS-CoV, and MERS-CoV), only two genes were found to be differentially expressed in all three, with CLU being one.

CLU’s top 2% strongest ERCs show highly significant enrichment for terms relating to coagulation cascades and clot formation (File 3, e.g. “Complement and coagulation cascades”, FDR = 6.3E−12), as well as significant terms that are relevant to immunity, such as “Immune system” (FDR = 4.8E−03) and “activated immune cell type” (FDR = 3.4E−05). Among its top ERC proteins relevant to coagulation process are Coagulation Factor V (F5, ρ = 0.67, FDR = 9.1E−06, rank 3), Fibrinogen Gamma chain (FGG, ρ = 0.59, FDR = 1.7E−04, rank 18), Coagulation Factor XIII B chain (F13B, ρ = 0.63, FDR = 2.8E−05, rank 19), and Fibrinogen Alpha chain (FGA, ρ = 0.57, FDR = 2.9E−04, rank 27) (Fig. 3, File 1). Notably, fibrinogen is a major binding “client” of Clusterin in stressed plasma (Wyatt & Wilson, 2010). Little is known about GPR141; however, the ERC analysis suggests an important role in blood coagulation. Among GPR141’s top ERC proteins relevant to coagulation process are Kininogen 1 (KNG1, ρ = 0.60, FDR = 9.3E−04, rank 5), Plasminogen Activator (PLAT, ρ = 0.58, FDR = 6.5E−04, rank 6), Thrombin (Coagulation Factor II or F2, ρ = 0.58, FDR = 6.5E−04, rank 7), Fibrinogen Beta chain (FGB, ρ = 0.57, FDR = 6.5E−04, RR 11, 11), Complement C1s (C1S, ρ = 0.54, FDR = 1.6E−03, rank 22), F2R-like thrombin (also called trypsin receptor 3; F2RL3, ρ = 0.52, FDR = 2.6E−03, rank 37), and Coagulation Factor V (F5, ρ = 0.52, FDR = 1.7E−03, rank 39) (Fig. 2, File 1).

GPR141 has a highly significant ERC to CLU, with these two proteins being each other’s first ranking ERCs (ρ = 0.68, FDR = 9.1E−06, RR 1,1). The pattern suggests a strong biological interaction, although none is described in the literature. The result supports investigating functional interactions between CLU and GPR141, based upon their high ERC and reciprocal ranks. Our network analysis below (“ERC Reciprocal Rank Networks Implicate Coagulation Pathways and Immunity”) further supports extensive interconnections among ACE2, Clusterin, GPR141, and coagulation pathway proteins, implicating the protein interaction pathway as a possibly significant contributor to disruption of coagulation in COVID-19 disease. Coagulation cascade proteins found in the ACE2’s top 2% ERCs, ACE2 reciprocal rank network, and Clusterin-GPR141 associated proteins are highlighted in Fig. 2.

Androgen Receptor (AR, ρ = 0.57, FDR = 8.8E−04, rank 13) is the receptor for the male hormone androgen. It plays a major role in reproductive system development, somatic differentiation, and behavior (Matsumoto et al., 2008). Androgen-AR signaling induces ACE2 (Wu et al., 2020), while knockdowns of AR result in downregulation of ACE2 (Samuel et al., 2020). AR agonists also reduce SARS-CoV-2 spike protein-mediated cellular entry (Deng et al., 2021). Additionally, AR is associated with COVID-19 comorbidities (Dolan et al., 2020), and recently implicated in the severity of COVID-19 in women with polycystic ovarian syndrome, a disorder associated with high androgen levels and androgen sensitivity (Gotluru et al., 2021). Our ERC finding indicating ACE2 and AR coevolution suggests regulatory feedback between these two proteins, which could be relevant to COVID-19 severity and other sex differential pathologies, such as cardiovascular disease (Viveiros et al., 2021).

Other notable significant ACE2 ERCs (Table 1) include Metabolism regulating signaling molecule D (FAM3D), Transmembrane-protein 63C (TMEM63C); Collagen Type IV Alpha 4 (COL4A4), L1 cell adhesion molecule (L1CAM), and ITPRIP-like 2 (ITPRIPL2). More detailed information on these and other proteins mentioned in this section is provided in “ERC Reciprocal Rank Networks Implicate Coagulation Pathways and Immunity” and the Supplemental Text.

ERC reciprocal rank networks implicate coagulation pathways and immunity

As mentioned previously, two proteins with a significant evolutionary rate correlation (ERC) may often “rank” each other differently in their respective top ERC connections. This occurs because some proteins have more extensive ERC connections than others. High reciprocal ERC ranks between protein pairs may be more indicative that they are under strong coevolutionary pressure in their sequence and function. We have thus found it useful to evaluate these reciprocal rank connections as a network. The rationale is that such proteins are likely to be reciprocally evolving (“coevolving”). To build reciprocal rank networks, we use protein pairs that reciprocally share ranks less than or equal to 20 (RR20), which are the top one percent for each protein’s respective ERC set.

A core ACE2 reciprocal rank network was generated by building reciprocal rank connections (RR20) outward of ACE2, to provide a backbone set of RR20 protein connections. The backbone was expanded on by adding the RR20 connections of the non-ACE2 backbone proteins. Unidirectional ERCs (≤ rank 20) were then added between proteins within the RR set to produce an ACE2 Core Reciprocal Rank (CRR) Network (Fig. 3). The network is designed to capture features of ACE2’s protein interactions as revealed by the strong reciprocal evolutionary correlations among proteins.

ACE2 also has highly significant ERCs to proteins that do not rank ACE2 within their top 1% of ERCs, due to those proteins having more protein interactions with higher ERCs. A second network was therefore generated using ACE2’s top ten unidirectional ERCs, followed by calculating the RR20 associations for those proteins. This second network is referred to as the ACE2 Unidirectional Reciprocal Rank (URR) Network (Fig. 4).

These are presented below. In general, the reciprocal ranks analysis lends credence to our proposition that ERCs reveal real biological interactions, as well as providing predictions for novel protein interactions possibly of importance to COVID-19 pathologies and protein-interaction networks.

The ACE2 unidirectional reciprocal rank (URR) network

ACE2 also has highly significant ERCs with interacting proteins that are unidirectional, meaning that ACE2 ranks these proteins in its top 2%, but the partner protein does not rank ACE2 within its top 2% due to higher ERC correlations with other partners (Table 1). Some of ACE2’s highest-ranking proteins fall into this category, including GEN1 (rank 1), XCR1 (2), IFNAR2 (5) KIF3B (6), and ITPRIPL2 (7), FAM227A (8), TLR8 (9), COL4A4 (10), F5 (12), and AR (13). To focus on strong protein connections in this set, we took the top ten proteins with unidirectional ERCs for ACE2 and then added their reciprocal rank 20 (RR20) partners. The resulting ACE2 Unidirectional Reciprocal Rank (URR) Network contains 69 proteins (Fig. 4).

Notable in the network are many proteins involved in immunity and cytokine signaling, such as IFNAR2 (Interferon alpha/beta receptor 2), XCR1 (X-C Motif Chemokine Receptor 1), and ICOS (Inducible T Cell Costimulator). There are also Toll-Like Receptors TLR8 and TLR9, which stimulate innate immune activity (Forsbach et al., 2011), and Tumor Necrosis Factor related proteins such as TNSFS18, TNTSF15, TNFRSF9, and TNRRSF1A.

Enrichment analysis of the URR network generates 72 significant terms (File 3). The network is highly enriched for cytokine-cytokine receptor interaction (FDR = 6.5E−06), I-kappaB kinase/NF-kappaB signaling (FDR = 1.6E−06), necroptosis (FDR = 3.3E−03), viral infections, such as Human Papillomavirus (FDR = 5.7E−04) and Herpes virus (FDR = 3.5E−03), JAK-STAT and PI3K-AKT signaling pathways, Toll-like receptor signaling, and immune system Homo sapiens (FDR = 3.7E−03).

XCR1 is the 2nd highest ACE2 ERC. It is the receptor for chemokine XCL1, which is produced in response to infection and inflammation, and during development of regulatory T cells (Lei & Takahama, 2012). Furthermore, XCR1 maps to a region implicated in severe COVID-19 by a genome-wide association study (Severe Covid-19 GWAS Group, 2020). As seen in Fig. 4, XCR1 forms a RR subnetwork with six other proteins (ICOS, CCR5, WDR66, TNSFS15, PRSS38 and FAM227A), three of which are known to be involved in immunity. ICOS (Inducible T Cell Costimulator) is reciprocally evolving with XCR1 based on their ERC interaction. It is an inducible T Cell stimulator that is essential for T helper cell responses (Hutloff et al., 1999; Tafuri et al., 2001). In addition, ICOS signaling is impaired in COVID-19 patients requiring hospitalization (Hanson et al., 2020). The high ERC between ACE2 and XCR1 and high reciprocal ranks of XCR1 to ICOS suggests that the disruption of an ACE2-XCR1 interaction could have a contributory role in COVID-19. C-C Motif Chemokine Receptor 5 (CCR5) forms a significant RR ERC with XCR1 as well. Several studies have implicated CCR5 variation and expression to be associated with COVID-19 severity (Gómez et al., 2020; Hubacek et al., 2021; Kasela et al., 2021), while others have not (Bernas et al., 2021). TNFSF15 is a third immune response protein in the XCR1 RR subnetwork that shows elevated expression in patients with severe COVID-19 (Jain et al., 2021). We recognize that the involvement of these immune-related proteins in COVID-19 does not require an effect mediated through ACE2. Instead, their protein evolutionary correlations suggest that ACE2 may play a contributory role to COVID-19, possibly through XCR1-related pathways, as suggested by the network analysis.

IFNAR2 is another protein that is highly correlated with ACE2 (ρ = 0.62, FDR = 6.1E−04) and is also implicated in severe COVID-19 by GWAS and expression data (Liu et al., 2021; Pairo-Castineira et al., 2021). It has RR20 ERCs with ten other proteins and is embedded in a complex web of interactions with members of the ACE2 network. Here we draw attention to a few key features. Notably, IFNAR2 and IFNAR1 are RR partners, as expected given that they combine to form the IFN-alpha/beta receptor, which is the receptor for both alpha and beta interferons. IFNAR2 forms a high RR relationship with TNFRSF1A (ρ = 0.84, FDR = 4.8E−12, 1,1 reciprocal ranks). This protein is the receptor for TNFα and the pathway affects apoptosis and inflammation regulation. Jin et al. (2015) found that ACE2 deletion increases inflammation through TNFRSF1A signaling, lending further support to a functional association between ACE2 and this protein.

GEN1 is the highest-ranking ACE2 ERC protein (ρ = 0.67, FDR = 4.2E−05), and it functions as a resolvase of Holliday junctions and a DNA damage checkpoint signaling (Chan & West, 2015). Frankly, we are perplexed by the functional significance of ACE2-GEN1 correlated evolution. As observed in the ACE2 network, GEN1 is a highly interconnected protein, with 14 RR20 connections in the network. This result suggests that GEN1 may have additional functions beyond DNA replication. Indeed, although its second-highest RR is to CC2D1B (2,1), a protein involved in mitosis, its highest RR is to Interferon Lambda Receptor 1 (IFNLR1), with an impressive Spearman correlation of ρ = 0.89 (FDR = 6.2E−17). As IFNLR1 binds cytokine ligands and stimulates antiviral response, this suggests some feedback mechanism between GEN1 and the immune system, possibly related to its functional role in DNA damage checkpoint signaling. Indeed, its top 2% ERCs show enrichment for multiple viral infection terms (File 3). Therefore, it appears that GEN1 has a “hidden life” that ERC analysis suggests warrants exploration.

The Collagen Type IV A4 subnetwork (Figs. 4, 5) lends further credence to the view that ERCs can detect proteins with likely binding partners. COL4A4 is a component of the Collagen Type IV protein complexes in basement membranes in the extracellular matrix of various tissues, including the kidney glomerulus and vascular endothelial cells, and lung alveoli (Myllyharju & Kivirikko, 2001). COL4A4, COL4A3, and COL4A5 complex with each other in the basement membranes of kidney glomeruli – mutations in these COL4A proteins are known to cause different kidney disorders (Torra et al., 2004; Wiradjaja, DiTommaso & Smyth, 2010). Consistent with their expected binding, COL4A4 and COL4A3 are each other’s reciprocal best partners (ranks 1,1) and highly correlated with each other (ρ = 0.88, FDR = 4.4E−16). Both show highly significant ERCs to COL4A6 (rank 6,5 for COL4A4 ρ = 0.83, FDR = 2.1E−12; rank 22,30 for COL4A3 ρ = 0.78, FDR = 1.6E−10). Thus, evolutionary rate correlations show highly significant ERCs among Collagen Type IV proteins known to physically interact. A future direction is to use ERCs to more precisely define predicted coevolving protein segments, which could be used to inform docking simulations and experimental studies.

COL4A5 also has significant ERCs to COL4A3 (ρ = 0.71, FDR = 2.2E−08) and COL4A4 (ρ = 0.71, FDR = 1.7E−08), but these do not qualify as RR20 due to the large number of high ERCs for COL4A5. Interestingly, COL4A5-MUC15 are top-ranking partners (ranks 1,1) with a very high ERC (ρ = 0.89, FDR = 3.2E−16). MUC15 is a cell surface protein that is believed to promote cell-extracellular matrix adhesion and it is implicated in affecting influenza infection (Chen et al., 2019), which may increase its relevance in the context of COVID-19 infection. ERCs may help to inform candidate domains within each protein that are involved in their expected binding affinity.

Coagulation Factor V (F5) is known for its role in the coagulation cascade. However, F5 is a highly ERC-connected protein, with 43 proteins ranking it in their respective top 5 highest ERCs. This connectedness is also reflected in the RR20 network shown below (Fig. 6). F5 has 16 RR20 connections out of the maximum 20 possible. Although F5 is a vital protein in the coagulation cascade, its top 16 RR connections indicate immune functions, including Interferon λ receptor 1 (IFNLR1; RR 4,10) and Oncostatin M Receptor (OSMR; RR 1,4). This is reflected in the enrichments among its 16 RR proteins for the JAK-STAT signaling pathway (FDR = 8.7E−03) and response to cytokine (FDR = 2.5E−02). Similarly, the F5 top 2% ERC show enrichments for 54 terms (File 3); notably many related to inflammatory response (FDR = 1.1E−03) and the complement system (FDR = 8.4E−03). The functions of several of F5’s RR20 partners are not well known, such as C14orf140 and C5orf34. Their top 2% enrichment suggests cytokine receptor activity (FDR = 2.7E−02) for C14orf140, and Human Complement System (FDR = 1.9E−03) and cytokine receptor activity (FDR = 2.1E−02) for C5orf34. In conclusion, F5 appears to have a “secret life” of strong protein interactions reflecting moonlighting functions with extensive signaling or modulation roles beyond coagulation regulation.

ERCs and protein interactions

We postulate that ERCs detect proteins that are coevolving due to functional interactions. Furthermore, we propose that physical binding is an important mechanism contributing to significant ERCs between proteins. This is consistent with anecdotal observations from this study of high reciprocal rank ERCs among the fibrinogen components FGA, FGB, & FGG, the Collagen Type IVA proteins COL4A4, COL4A3, and COL4A6 proteins, and Interferon alpha/beta proteins IFNAR2 and IFNAR1.

To further investigate the role of binding affinity, we examined the mammalian protein complex database CORUM (Giurgiu et al., 2019) to determine whether significantly higher Spearman rank correlations (ρ values) are found among proteins within known protein complexes. A set of 139 protein complexes (excluding those with overlapping proteins) were identified which contain at least two members from our ERC data set, for a total of 258 pairwise comparisons. We compared the ρ values of within complex proteins to the median values for proteins outside the complex and found that Spearman rank correlations of within complex proteins were significantly higher than its between complex values according to Wilcoxon matched signs rank tests (WMRST) under a significance level of α = 0.05 (p = 5.2E−04), with a median increase of 6.3% (File 11). Many of the complexes contain large numbers of proteins, reducing the probability of direct physical contact between individual members. We therefore also analyzed only proteins from complexes with 5 or fewer members (96 pairs). In this case, the median ρ value increase is 15.8% (WMSRT p = 6.2E−03). The results support the view that proteins within known complexes show higher ERCs than between complexes, and further implicate physical contact as a contributor to ERCs. However, other studies have found ERCs between proteins that do not bind to each other, but are involved in shared function, such as metabolic pathways (Clark, Alani & Aquadro, 2012). Thus, future research is needed to better understand the different biological drivers of ERCs between proteins.