A broadscale analysis of host-symbiont cophylogeny reveals the drivers of phylogenetic congruence

Symbioses exert substantial biological influence, with great evolutionary and ecological relevance for disease, major evolutionary transitions, and the structure and function of ecological communities. Yet, much remains unknown about the patterns and processes that characterise symbioses. A major unanswered question is the extent to which symbiont phylogenies mirror those of their hosts and if patterns differ for parasites and mutualists. Addressing this question offers fundamental insights into evolutionary processes, such as whether symbionts typically codiverge with their hosts or if diversity is generated via host switches. Here, we perform a meta-analysis of host-symbiont phylogenetic congruence, encompassing 212 host-symbiont cophylogenetic studies that include ~10,000 species. Our analysis super-sedes previous qualitative assessments by utilising a quantitative framework. We show that symbiont phylogeny broadly reflects host phylogeny across biodiversity and life-history, demonstrating a general pattern of phylogenetic congruence in host-symbiont interactions. We reveal two key aspects of symbiont life-history that promote closer ties between hosts and symbionts: vertical transmission and mutualism. Mode of symbiosis and mode of transmission are intimately interlinked, but vertical transmission is the dominant factor. Given the pervasiveness of symbioses, these findings provide important insights into the processes responsible for generating and maintaining the Earth's rich biodiversity.


INTRODUCTION
Symbioses are prolonged and intimate associations between organisms of different species (Wilkinson, 2001), typically defined in terms of an interaction between a larger host organism and a smaller symbiont (Estrela et al., 2016).Symbiotic relationships may be mutualistic, where both partners benefit from the interaction, parasitic, where one partner benefits and the other suffers a cost, or commensalistic, where one partner benefits while the other neither benefits nor suffers a loss (Leung & Poulin, 2008).Symbioses exist across all major phylogenetic lineages and are a common feature of life (Figure 1).
A key unanswered question in the study of symbiosis is the extent to which symbiont phylogenies are congruent with those of their hosts.Given an ancestral host-symbiont association, multiple rounds of strict codivergence will result in a symbiont phylogeny that

S Y N T H E S I S A broadscale analysis of host-symbiont cophylogeny reveals the drivers of phylogenetic congruence
Alexander Hayward 1  | Robert Poulin 2 | Shinichi Nakagawa 3  perfectly mirrors host phylogeny (Figure 2A).In practice, perfect congruence is considered to be rare, due to the effects of disruptive events such as host switching (Figure 2B), extinction (Figure 2C), independent divergence (Figure 2D), or failure to diverge (Figure 2E), which all lead to incongruence in host-symbiont phylogenies (Paterson et al., 2003) (Figure 2F).Meanwhile, imperfect congruence can arise through host-switch divergence, also referred to as host-shift speciation or pseudocospeciation (De Vienne et al., 2007;de Vienne et al., 2013), whereby repeated host switches followed by divergence events mirror host phylogeny (congruence is imperfect since symbiont branch lengths will be shorter than host branch lengths).
The idea that symbionts share phylogenetic congruence with their hosts is encapsulated by 'Fahrenholz's Rule', which states that parasite phylogeny reflects host phylogeny (Eichler, 1942;Fahrenholz, 1913).Fahrenholz's Rule was formulated over one hundred years ago and remains untested formally.Yet a persistent assumption of host-specialisation and accompanying host-parasite phylogenetic congruence remains within the field of hostsymbiont interactions (Nylin et al., 2018), despite limited support from the literature (de Vienne et al., 2013).Given the availability of statistical methods to quantify host-symbiont phylogenetic congruence (Extended Data Table 1), and an accumulation of studies that apply these approaches, opportunities now exist to examine the generality of host-symbiont phylogenetic congruence within a formal quantitative framework.
Here, we apply meta-analytical methods to a large body of published studies that quantify phylogenetic congruence for particular symbioses, to test the extent of congruence and factors of influence across host-symbiont diversity.We consider reported test statistics from the two most widely applied approaches for estimating host-symbiont phylogenetic congruence: one based on estimations of shared discrete macroevolutionary events (i.e., 'TreeMap', Page, 1994), and one based on phylogenetic distance comparisons to estimate overall similarity (i.e., 'Parafit', Legendre et al., 2002).Specifically, considering 212 published hostsymbiont cophylogeny studies that include ~10,000 species, spanning a wide spectrum of host and symbiont taxa, we examine the following fundamental outstanding questions on host-symbiont cophylogeny: 1. To what extent are host and symbiont phylogenies congruent?
2. Does host-symbiont phylogenetic congruence vary according to mode of symbiosis?(i.e., parasitism vs. mutualism) 3. Are cophylogenetic patterns influenced by host and symbiont taxonomy?4. Is host specificity an important determinant of cophylogeny? 5. Does host environment affect cophylogeny?(i.e., endosymbiosis vs. ectosymbiosis) 6.How does mode of transmission influence cophylogeny?(i.e., vertical transmission vs. horizontal transmission) We find that both mutualists and parasites show phylogenetic congruence with their hosts, but the pattern is stronger for mutualists.Furthermore, considering a selection of traits predicted to be of importance, vertical transmission (i.e., 'parent-to-offspring transmission') emerges as the main trait that promotes host-symbiont cophylogeny.

Measures of host-symbiont phylogenetic congruence
We considered studies that examine host and symbiont phylogenetic congruence using knowledge of the associations that exist between symbiont and host tips.We included reported test statistics (p values and the number of randomisations; see below) from the two most widely applied analytical approaches for estimating host-symbiont phylogenetic congruence: TreeMap (Charleston, 1998;Charleston & Perkins, 2003;Page, 1994), based on estimates of shared discrete macroevolutionary events, and ParaFit (Legendre et al., 2002), based on phylogenetic distance comparisons.Both approaches adopt a null hypothesis of independent evolution between host and symbiont phylogenetic trees, but they interrogate the hypothesis using different forms of phylogenetic data.In TreeMap, congruence is assessed by mapping symbiont phylogeny onto host phylogeny and quantifying apparent host-symbiont codivergence events (out of a maximum of n -1 for a phylogeny of n symbiont taxa), vs. other events (i.e., independent speciations, host switches, and sorting events).Host and symbiont phylogenies are considered congruent if the total is significantly greater than the distribution of scores observed when symbiont phylogeny is randomised a large number of times and mapped back onto host phylogeny, with the resultant p value indicating the level of significance (Charleston, 1998;Charleston & Perkins, 2003;Page, 1990Page, , 1994)).In contrast, instead of proposing evolutionary scenarios, ParaFit adopts a distance matrix approach to quantifying congruence between host and symbiont phylogenies (Legendre et al., 2002).Specifically, the global fit between host and symbiont phylogenies is calculated as the sum of squares of the cross-products of principal coordinates describing each tree, weighted by a matrix describing host-symbiont associations (i.e., 'H-P links'), This is compared to the distribution of values obtained by permutating the dataset many times to determine significance (i.e., the global test of host-symbiont cophylogeny 'ParaFitGlobal') (Legendre et al., 2002).
A potential advantage of ParaFit over TreeMap is that it offers a test of the relative contribution of each individual host-symbiont link to the overall level of congruence observed, which TreeMap does not.Additionally, ParaFit can utilise any distance data, and does not require an estimate of phylogeny, since the distances used are explicitly considered to be functions of host and symbiont phylogeny (Legendre et al., 2002).ParaFit does not correct for phylogenetic non-independence, and host-symbiont pairs with deeper most recent common ancestors (MRCAs) are effectively weighted over those with shallower MRCAs, since deeper MRCAs represent multiple descendent pairs of taxa (thus terminal MRCAs are represented once, but deeper MRCAs are represented multiple times in pairwise distance matrices) (Schardl et al. 2008;de Vienne et al., 2013).Thus, under certain circumstances distance-based methods such as ParaFit may be more sensitive to dependence between trees.Event-based approaches such as TreeMap are not affected by this issue, but they can be prohibitively computationally intensive for large datasets containing many host and symbiont taxa.
Several more recent cophylogenetic programs have been published since TreeMap and ParaFit (Supporting Information: Appendix S1, Table S1.1).Our approach was to include the most widely used events-based program and the most widely used global fit program in this study (i.e., TreeMap and ParaFit respectively).However, we briefly summarise other major programs for assessing phylogenetic congruence in date order below.Similarly to TreeMap, Tarzan (Merkle & Middendorf, 2005), is an event-based approach for assessing host-symbiont congruence, but Tarzan allows the use of symbiont divergence time estimates, permitting temporal teasing apart of codivergence events, and facilitating consideration of the mechanisms responsible for generating congruence (Merkle & Middendorf, 2005).However, the accuracy of Tarzan has been questioned (Conow et al., 2010), and very few published studies have utilised the programme.Jane, another event-based approach, builds on Tarzan's shortcomings (Conow et al., 2010).Additionally, Jane incorporates time zone estimates for both host and symbiont trees, the use of upper bounds on host-switch distances, and flexibility in costs associated with hostswitches over time and phylogenetic distance.However, Jane has been used in much fewer published studies of cophylogeny than either ParaFit or TreeMap, and as with other events-based approaches, it is relatively slow and codivergence is considered to be the most parsimonious hypothesis, which is likely to be an unrealistic assumption (de Vienne et al., 2013).PACo (Balbuena et al., 2013) is a global fit method, but unlike ParaFit, PACo assumes dependence between phylogenies.PACo includes a more nuanced implementation than ParaFit, facilitating more detailed analyses into individual hostsymbiont links, and is considered more statistically robust (Balbuena et al., 2013).However, PACo has a much lower uptake in published studies than ParaFit to date.COALA (Baudet et al., 2014) is an event-based approach, which estimates event frequencies using an approximate Bayesian computation approach, with the aim of increasing confidence in a set of reconciliation costs.However, COALA allows each symbiont leaf to be associated with a maximum of just one host leaf, and very few published studies have utilised the approach.Random TaPas (Balbuena et al., 2020) uses a global-fit approach, and provides a test statistic that is inversely proportional to the cophylogenetic signal present, as well as a summary of codivergence events involved in the construction of the hypothesised tanglegram between host and symbiont phylogenies.Random TaPas permits identification of the individual host-symbiont interactions, phylogenetic tips and nodes that contribute most strongly to overall congruence.It can also incorporate the effect of phylogenetic uncertainty, use time-calibrated trees to distinguish codivergence, and it allows the rapid analysis of extremely large datasets.While Random Tapas does not fully account for phylogenetic nonindependence, the methodology it employs aims to ameliorate its influence, and it represents a very promising program upon which to base future meta-analyses.
Both ParaFit and TreeMap presume fully resolved and accurate host and symbiont phylogenies, and neither method takes phylogenetic uncertainly into account.The vast majority of existing host and symbiont cophylogeny studies are based on phylogenies estimated using single or just a few markers.The field is changing rapidly, and the use of multi-gene and genomic-scale datasets will improve phylogenetic accuracy.However, until a large body of such studies is available, analyses based on single or small numbers of markers represent the main source of data for inclusion in meta-analyses.In the meantime, there is no evidence to suggest that the use of small numbers of markers for phylogeny estimation affects assessments of host-symbiont phylogenetic congruence in a systematic manner, other than to potentially reduce observed signals of congruence somewhat (since incongruence is much more likely than congruence among randomly generated phylogenies).Meanwhile, due to the influence of variation in time and space, we acknowledge that host-symbiont interactions can be more nuanced than the snapshot typically captured in published studies, resulting in some oversimplification of the patterns.
Phylogenetic congruence arising as a result of vertical transmission and codivergence, and phylogenetic congruence arising via other means (such as successive host-switching across host phylogeny), can produce a similar phylogenetic signature.The mechanisms responsible for generating congruence are not differentiated by Parafit or TreeMap, as this requires accurately dated host and symbiont phylogenies, of which relatively few pairs are available currently.However, our analyses are unaffected by this issue, since they were not directed at unravelling alternative scenarios underlying host and symbiont speciation, but rather towards considering broadscale patterns in phylogenetic congruence.

Literature search
A literature search was performed using Google Scholar on 19th March 2019, which identified 368 citations of the original paper for TreeMap (Page, 1994), and 332 citations of the original paper for Parafit (Legendre et al., 2002), resulting in a total of 700 articles that were screened to extract metrics for inclusion in our meta-analysis.Articles that did not contain cophylogenetic analyses were immediately excluded.Studies focussing at the population level were also excluded, as these do not represent true cophylogenetic analyses at the macroevolutionary level.Additionally, studies that included less than four taxa were excluded from consideration, as these do not provide sufficient power for inclusion in the meta-analysis.Studies that did not report the test statistic for congruence were also necessarily excluded.After exclusions were applied, a sample of 212 separate host-symbiont cophylogenetic analyses remained, with a mean of ~23 host-symbiont links per study, upon which our meta-analysis is based (Appendix S1, Table S1.2).

Data extraction
Metrics collected for this study are described below, with reference to Table S1.1 in Appendix S1.A short citation of each study was recorded under 'authors', and the year of publication was recorded in 'year'.Hosts and symbionts were classified broadly according to Linnean taxonomy for 'host_tax_broad' and 'symbiont_tax_broad' as either: invertebrate, vertebrate, plant or microbe (i.e., microscopic symbionts such as fungi, protozoa, bacteria, viruses).
We adopted the mode of symbiosis and mode of transmission between host species specified by the authors in each individual study for 'symbiosis' and 'mode_ of_transmission_broad'.In cases where either mode of symbiosis or mode of transmission were not directly specified by authors, we consulted the literature for clarification.In a small number of studies restricted to bacterial intracellular symbionts, the mutualism-parasitism distinction was not defined by the authors and either no further information was available, or a symbiont was cited in the literature as being both a mutualist or a parasite, depending on which study was considered.The nature of the relationship between bacterial intracellular symbionts and their hosts is complex, and in some cases they may display both beneficial and detrimental effects simultaneously (Zug & Hammerstein, 2015).In a few cases of conflict or where authors did not explicitly state mode of transmission for bacterial intracellular symbionts, we assumed a mode of transmission in line with the majority of available references.We only encountered one study where authors categorised the mode of symbiosis as commensalism (McFrederick & Taylor, 2013).On the continuum of symbioses from pure parasitism (fitness losses for the host) to mutualism (fitness gains for the host), commensalism represents a single point where losses and gains for the host precisely equal zero.Consequently, commensalism is an unlikely and unstable state, easily tipped to one side or the other with any small change in external conditions.Thus, the lack of widely recognised groups of commensals is the likeliest explanation for the scarcity of studies on commensalism in our data (note that we did not include this category, commensalism, in our analyses).
The total number of host tips that were linked to a symbiont taxon were summed to provide 'host_tips_ linked', which in a very few cases was corrected to remove multiple sampling of the same host species, to provide 'host_tips_linked_corrected'.The total number of symbiont tips with a link to a host taxon were summed to provide 'symbiont_tips_linked', while the total number of individual links between hosts and symbionts was recorded as 'total_host_symbiont_ links'.If all symbionts in a phylogeny were strict specialists, such that each one had a single link to a single host, 'total_host_symbiont_links' would simply equal 'symbiont_tips_linked'.However, because symbionts are often associated with more than one host, the value of 'total_host_symbiont_links' was often higher than the total number of symbionts included in a study.Thus, a measure of symbiont generalism was captured using 'host_range_link_ratio', defined as 'total_host_ symbiont_links' divided by 'symbiont_tips_linked', providing the mean number of host-symbiont links observed per symbiont taxon, with the measure increasing with increasing generalism.An alternative estimate of symbiont host specificity was captured using 'host_ range_taxonomic_breadth', which considers Linnean taxonomic rank, and was calculated by assigning an incremental score to successive host taxonomic ranks per symbiont in turn (i.e., single host species = 1, multiple host species in the same genus = 2, multiple host genera = 3, multiple host families = 4, multiple host orders = 5), summing the total score across all symbionts, and dividing by 'symbiont_tips_linked' (i.e., the total number of symbionts).Consequently, 'host_range_tax-onomic_breadth' increases with symbiont generalism, such that symbiont phylogenies containing symbionts capable of infecting hosts from a wide range of taxonomic ranks are assigned a greater score.
The number of phylogenetic permutations performed by authors during cophylogenetic analyses was recorded as 'no_randomizations', which poses a unique problem in our meta-analysis (discussed in the section 'Publication bias and sensitivity analysis').The resultant p value from each study was recorded as 'p_value', whereby observed p values decrease with a decreasing likelihood of observing host-symbiont cophylogeny by chance alone (i.e., as calculated during permutation tests performed by authors during TreeMap or ParaFit analyses).

Effect size
We used p values obtained from the randomisation tests implemented in TreeMap and ParaFit as measures of phylogenetic congruence.These p values were converted into r and its transformation Zr.We can calculate r equivalent via t values with df = N − 2 from p values (one-tailed) (Rosenthal & Rubin, 2003), following the formula: where N is sample size, with a sampling variance estimate for r, Var(r) (Borenstein et al., 2011).Note that Equation 1 is exactly the same as the formula used to convert t values from regression (or correlation) analysis to r (Pearson's correlation coefficient), commonly used in a meta-analysis.In Rosenthal and Rubin (2003), r equivalent was proposed to deal with situations where one only has a p value and a corresponding df, assuming this p value comes from a t test, when comparing two independent groups.Therefore, r equivalent is equivalent to a pointbiserial correlation coefficient between measurements from two groups, and a dummy variable (0 and 1, indicating two groups), and is analogous to Pearson's correlation coefficient.
To perform a meta-analysis and to enable an evaluation of cophylogenetic patterns across individual studies, we made two assumptions.First, t-distributions and permuted (randomised) distributions can be likened; a clear and strong association certainly exists between these two distributions (Proschan et al., 2014).Under this assumption, although TreeMap and ParaFit use different randomisation procedures, if their p values are obtained from the same cophylogenies, they should be identical or very similar.We tested this by correlating effect sizes of both TreeMap and ParaFit based on the same cophylogenies from the 11 studies that employed both methods; they were closely, although not perfectly, correlated (r = 0.781, t 9 = 3.757, and p = 0.0045; for more, see Appendix S2).Second, we supposed that congruence between two phylogenies could be conceptualised as correlation.Under Equation 1, N is the number of pairs of two variables used to calculate a correlation.In our case, we cannot get the exact number of pairs, but we could obtain the number of 'effective' pairs (N ep ): where n h and n s are the sample size for hosts and symbionts respectively, and when n h = n s = n, N ep = n; a version of Equation 3 is used to calculate 'effective sample size' for an unbalanced design (ie., n h ≠ n s ).It may not be dif- ficult to see N ep as an analogue to the number of pairs.Nonetheless, it is possible that our two assumptions would not completely hold especially when sample sizes are small.Importantly, r equivalent is currently the only suitable method by which to conduct a meta-analysis using p values (Rosenthal & Rubin, 2003), and this is the case for our meta-analysis.Furthermore, there is no a priori reason to believe that the conversion from p values to t-values would bias results either way (but see the section 'Publication bias and sensitivity analysis').
We did not use raw values of r equivalent for our metaanalysis, but we employed Fisher z-transformed r, Zr (Zr equivalent ) and its sampling variance, Var(Zr), which are calculated here as follows: ( For meta-analysis, Zr has two major advantages over r: (1) Zr is approximately normally distributed, and thus, it does not have boundaries, unlike r (i.e., −1 and 1), which could violate assumptions of normality; and, (2) the sampling variance for r, Var(r), includes the term r in addition to N (Equation 2), while Var(Zr) only requires N (Equations 4).The observed r required in Equation 2 itself has error so that Var(r) is more prone to error than Var(Zr) (Borenstein et al., 2011).Where necessary and for ease of interpretation, we back-transformed Zr to r (i.e., r equivalent ) using this formula: An additional statistical consideration is that in many studies the number of host-symbiont phylogenetic permutations (randomisations) performed by authors limits the lowest observable p value statistic (e.g., often the number of phylogenetic permutations was set to 999 or 1000, which limits p ≤ 0.001).Thus, when p = 0, we used p = 0.001 for 999 and 1,000, p = 0.0001 for 9999 and 10,000, and so on (Appendix S1, Table 1.2).This created ceilings for ~32% of the effect sizes included in our analysis (i.e., these effect sizes could have been larger if more permutations were performed by the authors of a study).Since the true likelihood of achieving the observed number of codivergences by chance alone can be considerably lower than, for example, p = 0.001 (with randomisations = 1,000 and p = 0), particularly in cases of perfect or near-perfect host-symbiont congruence involving many taxa, our analyses consequently represent a conservative measure of host-symbiont phylogenetic congruence, and should be considered to be the lower bound of congruence that exists for the considered studies (sensitivity analyses relevant to this point are discussed below).
For meta-analytic models, we quantified 'heterogeneity' measures (I 2 ), which indicate the amount of variance unexplained after controlling for sampling variance (i.e.Equation 5) (Borenstein et al., 2017;Higgins et al., 2003;Nakagawa & Santos, 2012), while for meta-regression, we calculated the 'variance explained' by moderators (fixed effects), with marginal R 2 (sensu Nakagawa & Schielzeth, 2013).Notably, we calculated not only 95% confidence intervals for all of the estimates associated with our meta-analysis, but also 95% prediction intervals (presented in figures).The prediction interval indicates the probability of obtaining an effect size value from a new study (with no sampling error) within this interval; 95% of the time, a new effect size would fall between the 95% prediction interval (IntHout et al., 2016;Nakagawa et al., 2020).Recently, reporting prediction intervals has been strongly encouraged in medical sciences (IntHout et al., 2016).
All model specifications, model selection procedures and associated coding are provided in our Supporting Information (Appendix S1-S5), where we also provide statistical justification for merging datasets for both 'TreeMap' and 'ParaFit'.

Publication bias and sensitivity analysis
We acknowledge the possibility of bias in the literature.For example, there may be a bias towards studies showing stronger congruence, if there is: (a) a preference towards testing associations where authors suspect cophylogeny a priori, or (b) a lower likelihood of authors publishing results if no significant pattern of cophylogeny is identified.However, cases of non-congruent cophylogenies are just as informative and worthy of publication as cases of congruence.In our meta-analyses, we have assumed that there is no such publication bias where non-significant results are less likely to be published, and related biases such as time-lag bias, or the decline effect where the magnitude of effect sizes declines over time (Koricheva & Kulinskaya, 2019;Rothstein et al., 2005).
To test for publication bias, we examined funnel plots, which are scatterplots of effect size against study precision (e.g., here the inverse of standard error, which is proportional to sample size) (Egger et al., 1997).In the absence of publication bias, studies with high precision should lie close to the mean, while studies with low precision should lie evenly distributed on both sides of the mean, resulting in an approximately funnel-shaped distribution.Funnel plot asymmetry indicates a relationship between effect size and study precision, which can suggest publication bias (e.g., asymmetry due to missing effect sizes with high precision; see below), or true heterogeneity in a dataset.Since residuals of effect sizes are less likely to be influenced by heterogeneity (Nakagawa & Santos, 2012), we applied funnel plots that utilise residuals from meta-regression (after removing predicted fixed effects, i.e. marginal residuals).Furthermore, we utilised contour-enhanced funnel plots, an improvement over the original funnel plot which include contour lines that indicate different significance levels (Peters et al., 2008).Since a funnel plot is a visual method that does not offer a test of statistical significance, we also applied 'Egger's regression' test, which statically tests funnel asymmetry, and is often used in conjunction with funnel plots.The original Egger's regression test assumes data independence; because in our dataset certain effect sizes came from the same study, we applied a version that takes into account non-independence (Nakagawa & Santos, 2012).This version of Egger's regression can be implemented by adding sampling standard error (the square root of sampling variance, i.e., Equation 5) as a moderator, to a multilevel meta-regression model (Moreno et al., 2009;Rodgers & Pustejovsky, 2020).Finally, to test for time-lag bias (a decline effect), we fitted the year of publication ('year') as a moderator in a multilevel metaregression (Appendix S4).
Aside from potential issues of publication bias, it is notable that our meta-analysis had a unique problem.This is because the number of randomisations ('no_randomizations') used by the original authors sets an upper limit on the observed effect size (e.g., with randomisations of 100 and 1000, the smallest p values obtainable are 0.01 and 0.001 respectively; see above).We predicted two potential outcomes due to this ceiling effect.First, the ceiling would result in funnel asymmetry.However, asymmetry due to the ceiling effect would be distinctive and distinguishable visually from the asymmetry due to publication bias (for the details see the section 'Caveats and recommendations' in Results and Discussion).Second, the ceiling would lead to underestimates of the overall effect size, category-wise effects (effects of different groups), and contrasts between categories (e.g., parasites versus mutualists).Therefore, we conducted two sensitivity analyses to show that our underestimation of effect sizes did not influence our main conclusion.For the first test, we ran linear mixed models to examine if there was any statistical difference in the number of randomisations (fitted as log('no_randomations')) among categories of moderators.If no statistically significant differences were detected, this would confirm that this issue is spread across the three categories identified as important by the model selection procedure (i.e., 'symbiosis', 'host_tax_broad', 'mode_of_transmis-sion_broad'; Appendix S4, Table S4.12-S4.13).For the second test, we ran binomial family generalised linear mixed models with the logit link of whether a ceiling was reached or not (e.g., p = 0.01 for 100 randomisations).We intended to show that categories (groups) with high effect sizes were more affected by the ceiling effect.The second analysis also indicated that the higher the effect size estimate for a particular category (e.g., the mutualism group or the vertically transmitted symbionts group), the more underestimated the estimate.Both analyses were fitted in the lme4 package, version 1.1-23 (Bates et al., 2015) (Appendices S4 and S5).

Phylogenetic congruence
Firstly, we examined the extent to which symbiont phylogeny is congruent with host phylogeny, by performing a broadscale analysis of cophylogeny across biodiversity.We found support that host and symbiont phylogeny show a strong signal of congruence (r [all] = 0.490, 95% confidence interval CI = [0.446:0.532], I 2 = 61.2%; Figure 3A, and Appendix S2, Table S2.1).Since the interpretation of r equivalent (Equation 1) is not straightforward, we provide a benchmark based on our mean sample size of ~23 host-symbiont pairs.Specifically, when N = 23, and r = 0.3, 0.5, or 0.7, the resulting p values correspond to p = 0.08, 0.008, or 0.0001, respectively.Thus, our observed value of r equivalent = 0.49 corresponds to a p value of just above 0.008, indicating strong support for hostsymbiont phylogenetic congruence across biodiversity.We note that the strength of the association identified may be an underestimate of the true pattern, given the limits of current approaches to quantify congruence between host and symbiont phylogenies (for further discussion, see 'Caveats and recommendations' below).

Parasitism and mutualism
We then tested the more specific prediction that parasite phylogeny mirrors host phylogeny, across multiple host and parasite taxa.Parasitism is an antagonistic symbiotic interaction, defined in terms of benefits for the parasite and losses for the host (Poulin, 2007).Consequently, parasitism is only in the interests of the symbiont.Nonetheless, if parasites are specialised and tightly coevolved to single host species (Summers et al., 2003), this may favour evolutionary histories that result in strong congruence between host and parasite phylogenies.A signature of congruence is also expected if parasites make preferential host switches to closely related host species.Alternatively, if parasites commonly switch to novel and distantly related host species, either with parasite speciation following host-switching events (Giraud  , 2010;Rundle & Nosil, 2005), or without speciation (Clayton et al., 2015), we expect a weak or absent signal of phylogenetic congruence.Our results show that host and parasite phylogenies display greater similarity than expected by chance, with parasite phylogeny tending to mirror host phylogeny (r [parasite] = 0.451, CI = [0.396:0.503]; Figure 3B and Appendix S3, Table S3.1).While the pattern for host-parasite studies is weaker than that observed for all host-symbiont studies, it is consistent with a broad interpretation of Fahrenholz's Rule, to include significant congruence greater than expected by chance, rather than just the extremely unlikely case where host and symbiont trees are perfect mirror images of each other.
Next, we compared our findings for parasites to those for mutually beneficial symbiotic interactions.Mutualisms are defined in the context of 'reciprocal benefits' (sensu Bronstein, 2015).However, mutualism may also be viewed within a framework of 'reciprocal exploitation' (Herre et al., 1999).Within either framework, the incentive to cheat (Leigh, 2010;Trivers, 1971) may destabilise the longevity of mutualisms over protracted evolutionary timescales, potentially resulting in no greater overall phylogenetic congruence between hosts and mutualists than that observed for parasites.Alternatively, mutualisms may be robust to such challenges and thus be evolutionarily persistent (Douglas, 2008;Ferriere et al., 2007), in which case a pattern of tight host-mutualist phylogenetic congruence may be more prevalent.To explore patterns, we tested hostmutualist phylogenetic congruence across diversity, while simultaneously providing a general comparison of host-symbiont phylogenetic congruence between mutualism and parasitism.We found that host and mutualist phylogenies show statistically significant greater congruence than that observed between host and parasite phylogenies (r [mutalist] = 0.559, CI = [0.488:0.621], and, the difference, i.e., contrast, between mutualist and parasite symbionts, b [mutualist-parasite] = 0.144, CI = [0.028:0.260], R 2 = 0.054; Figure 3B and Appendix S3, Table S3.1).An immediately obvious potential explanation for the observed pattern is that reciprocal benefits provide a more cohesive force to unite host and symbiont evolutionary trajectories, than the antagonistic interactions involved in parasitism, but see continued discussion below.

Host and symbiont taxonomy
To further dissect cophylogenetic patterns between symbiont and host phylogenies, we examined whether taxonomy exerted an effect by partitioning studies according to broad taxonomic groups for hosts and symbionts.We found that effect size was positive across all host taxonomic groupings considered, with hostsymbiont associations involving microbial hosts showing the greatest congruence, followed by invertebrate hosts, vertebrate hosts, and lastly plant hosts (r [microbe] = 0.708, CI = [0.564:0.810]; r [invertebrate] = 0.557, CI = [0.479:0.627]; r [vertebrate] = 0.474, CI = [0.410:0.534]; r [plant] = 0.360, CI = [0.250:0.460]; and R 2 = 0.162; Figure 3C and Appendix S3, Table S3.2).For symbiont taxa, effect size was also positive across all taxonomic groups, but the explained variation (R 2 ) in phylogenetic congruence was lower, meaning that host taxonomy has more explanatory power for variation in congruence in our dataset (r [microbe] = 0.497, CI = [0.436:0.554]; r [invertebrate] = 0.467, CI = [0.399:0.531]; r [vertebrate] = 0.414, CI = [−0.189:0.790]; r [plant] = 0.814, CI = [0.585:0.923]; and R 2 = 0.0785; Figure 3D and Appendix S3, Table S3.3).It is unclear why associations involving microbial hosts should show the greatest congruence.However, by definition, symbionts of microbial hosts are intracellular symbionts, and it is possible that this results in greater host specificity and stricter vertical transmission (see below), than for extracellular symbionts of multicellular hosts.Meanwhile, the finding that host-symbiont cophylogenies involving symbionts that are plants show markedly higher congruence is partially due to the combined influence of strong observed congruence between mutualistic algal symbionts and their fungal hosts, and the small sample size in this category (reflecting a scarcity of studies on symbionts that are plants in relevant cophylogenetic literature).

Host specificity
Host specificity is a central aspect of symbiont life history with widespread ecological and evolutionary implications (Poulin, 2007;Thrall et al., 2007), which may exert an independent effect on host-symbiont phylogenetic congruence beyond the mutualist-parasite divide.If specialist symbionts are evolutionarily or ecologically 'tied-in' to their host associations, undergoing few or no host switches (Clayton & Johnson, 2003;Hafner et al., 2003;Hafner et al., 1994), they are expected to display high phylogenetic congruence with their hosts.However, if specialists possess an ability to make host-switches of varying frequency and magnitude, while remaining associated with one or few hosts (Hall et al., 2016;Krumbholz et al., 2009), congruence may be low.Meanwhile, truly generalist symbionts are expected to show low phylogenetic congruence with their hosts, unless dominant interactions occur within a subset of total host range (Charleston & Robertson, 2002).Although observed point estimates were in the expected direction (i.e., lower congruence with decreasing host specificity), no statistically significant association was identified for either measure of host specificity (the slope, b [ln(host range)] = −0.037,CI = [−0.187:0.114], R 2 = 0.002, and b [ln(taxonomic breath)] = −0.036,CI = [−0.216:0.143], R 2 = 0.001; Appendix S3, Figure S3.1-S3.2;Table S3.4-S3.5),suggesting that host specificity may not be an important general determinant of host-symbiont cophylogeny (at least, given the variation in these measures of host generalism observed in our data set).

Mode of transmission
Lastly, we examined the relevance of mode of transmission, to test the general expectation that vertical transmission (i.e.'intergenerational transmission'), as opposed to horizontal transmission (i.e.'infectious transmission'), should promote greater phylogenetic congruence (Ebert, 2013;Moran et al., 2008;Nieberding & Olivieri, 2007).If a symbiont is vertically transmitted (whether it has positive, neutral, negative effects on the host), this should automatically decrease the likelihood of host-switching, and increase the likelihood of codivergence.Therefore, phylogenetic congruence should emerge as a consequence.In contrast, horizontal transmission should facilitate exposure to novel hosts, potentially leading to host-switching, and decreased host-symbiont phylogenetic congruence.In line with this prediction, we found that mode of transmission was a significant predictor of host-symbiont phylogenetic congruence, with the strongest effect observed for vertical transmission, followed by the effect observed for mixed modes of transmission (symbionts are transferred via either route), and with horizontal transmission showing the lowest effect size (r [vertical] = 0.636, CI = [0.561:0.701]; r [both] = 0.521, CI = [0.432:0.600]; r [horizontal] = 0.419, CI = [0.361:0.474], and R 2 = 0.187; Figure 4B and Appendix S3, Table S3.7).An issue affecting our data is that all of the strictly vertically transmitted symbionts included in our analysis are classified as mutualists.However, many of the parasitic symbionts included show a mixed mode of transmission, and splitting by mode of symbiosis reveals that parasites with mixed modes of transmission show greater phylogenetic congruence with their hosts than horizontally transmitted parasites (r [both(parasite)] = 0.511, CI = [0.415:0.595], and r [horizontal(parasite)] = 0.418, CI = [0.354:0.595]; Figure 4C and Appendix S3, Table S3.8).Similarly, mutualistic symbionts show a stepwise increase in phylogenetic congruence with hosts, from horizontal to mixed, and from mixed to vertical transmission, with vertically transmitted mutualists displaying considerably greater phylogenetic congruence than horizontally transmitted mutualists (r [horizontal(mutualist)] = 0.42 7, CI = [0.289:0.548], r [both(mutualist)] = 0.607, CI = [0.332:0.787], and r [vertical(mutualist)] = 0.633 CI = [0.556:0.699]; Figure 4C and Appendix S3, Table S3.8).The pattern of changes seen in modes of transmission were similar between parasites and mutualists, highlighting that vertical transmission may be the key variable that explains the effect on phylogenetic congruence (further associated results including model selection procedures and modelaveraged results are provided in Appendix S3, Figure S3.4-S3.5, and Table S3.9-S3.13).

Caveats and recommendations
The validity of our results relies on the assumption that considerable publication bias is not present among the studies included in our meta-analysis, meaning that our dataset is a representative sample of available evidence on phylogenetic congruence between hosts and symbionts.If publication bias were present, it would result in funnel asymmetry with missing small effect sizes of low precision (Egger et al., 1997).Although significant asymmetry was observed for a few studies, which can be seen in the residual funnel plot (Figure 5A S4), this asymmetry was caused by missing small effect sizes of high precision.As mentioned above, this is the signature of the ceiling effect set by the limited numbers of randomisations performed by authors.Notably, once standard error is modelled as a moderator (i.e.Egger's regression, see above), no sign of asymmetry was apparent (Figure 5B), indicating potentially no publication bias.Also, we found little sign of a decline in effect size over time (Appendix S4, Figure S4.2-S4.3 and Table S4).
In addition to publication bias tests, we investigated whether the ceiling effect results in systematic bias in our results, which would invalidate our main findings.We found that the number of randomisations among categories of moderators were distributed randomly (Appendix S5, Table S5.1-S5.3).Therefore, we were in general more likely to underestimate categories with a larger mean effect (i.e., greater phylogenetic congruence), than a smaller overall effect although not always (Appendix S5, Table S5.4-S5.6).This finding indicated that our results were likely to be conservative in terms of overall congruence, as well as differences in congruence among categories.Related to this point, we used r equivalent as our effect size by converting p values from the randomisation tests in TreeMap and ParaFit primarily because this is the only effect size statistic that can be obtained given the information available.We reiterate that our use of r equivalent relies on a number of assumptions including the two mentioned above, and consequently our results need to be interpreted carefully.This is especially so for the overall effects, yet the resulting contrasts (differences) should not be biased (at least not overestimated).We quote Rosenthal and Rubin (2003)-'We think of r equivalent as a first-aid kit to be used for the time being until we can get to a highly sophisticated medical center.The medical center would be better, but it may be a long way away', with which we concur.Further, we provide the following recommendations for future (and past) studies of codivergence using TreeMap, ParaFit and related software packages.Researchers should set the number of randomisations much higher (e.g., >100,000).Given considerable recent increases in computational speed, even with large sets of hosts and symbionts, this is realistic, especially when applying newer methods (e.g., Balbuena et al., 2020).We should conduct reanalyses of published studies using more robust molecular phylogenies based on larger numbers of markers, and where necessary, incorporate phylogenetic uncertainty into analyses.Wherever possible, authors should seek to use dated phylogenies, so that mechanistic hypotheses relating to codivergence or the longevity of patterns can be teased apart (e.g., does host-symbiont phylogenetic congruence break down over time?).Additionally, future studies should make all data readily available, including phylogenetic tree files in Nexus or Newick format, in the spirit of open science (Nakagawa et al., 2020).Therefore, future meta-analyses can directly incorporate such data.This will also enable data to be readily checked and reanalysed in a consistent manner, and reused as methods continue to develop.Collectively, these recommendations will greatly help in determining a more accurate estimate of the overall magnitude of congruence, and will allow more detailed analyses on cophylogeny and the mechanisms that drive it.

CONC LUSION
Employing a large quantitative and systematic review of the cophylogeny literature, we tested major hypotheses regarding the drivers of host-symbiont phylogenetic congruence.We found support for a general pattern of congruence between host and symbiont phylogenies.Significant congruence was observed both for hosts and their parasites, and for hosts and their mutualists, but it was significantly stronger in the latter case.We also identified an effect of mode of transmission on host-symbiont phylogenetic congruence, finding that vertical transmission is correlated with greater congruence.Thus, we suggest that: (a) mutualism promotes greater host-symbiont phylogenetic congruence, and, (b) this is driven primarily by a predominance of vertical transmission among mutualists.
Our study emphasises the apparent universal importance of codivergence in generating diversity in interspecific interactions.We acknowledge that we adopt a classic cophylogenetic framework in this study, considering cases where symbionts are typically associated with very few or even single host species.Different processes may operate in diffuse symbioses, where symbionts are associated with a large number of different host species (Stanton, 2003).However, evidence of cophylogenetic patterns among gut microbiota and their hosts is emerging (Groussin et al., 2020;Groussin et al., 2017;Moeller et al.,  have also been reported across taxa and scales in more diffuse mutualisms, such as plant-pollinator associations.For example, a recent broadscale analysis of flowering plants and their pollinators found evidence of commonplace cophylogenetic signal across taxa and ecological scales (Hutchinson et al., 2017).However, more research is required to evaluate how widespread phylogenetic congruence is in diffuse mutualisms.
The observation that strict vertical transmission is uncommon in parasites fits with predictions that vertical transmission should favour decreased harm to hosts, since reducing host fitness ultimately reduces a symbiont's own opportunities for transmission (Ewald, 1987;Yamamura, 1993).Accordingly, parasite lineages forced to express strict vertical transmission have been demonstrated to shift along the parasitism-mutualism continuum towards greater host benevolence, both empirically (Bull et al., 1991;Stewart et al., 2005), and comparatively (Clayton & Tompkins, 1994;Herre, 1993).A key question remains whether mutualism facilitates the evolution of vertical transmission, or if the process typically proceeds as argued above (from a starting point of parasitism).A study of bacterial symbionts concluded that horizontal transmission was the most primitive mode of transmission, and that vertical transmission is an inescapable evolutionary end point, given the mutational processes affecting symbionts that adopt it (Sachs et al., 2011).However, these conclusions may be taxon-specific, and require further examination across host-symbiont taxonomic diversity.
Elucidating the processes that drive symbiont diversification is fundamental to understanding the factors responsible for generating the Earth's biodiversity, given that symbionts represent a considerable proportion of total species diversity (Poulin, 2014;Windsor, 1998).Meanwhile, unravelling the factors that influence hostsymbiont phylogenetic congruence is crucial for improving knowledge of host-switching, a major research focus in host-parasite interactions, with direct implications for the prediction and mitigation of zoonoses, emerging infections, and the control of agricultural pests (Gortazar et al., 2014;Morens et al., 2004;Woolhouse et al., 2005).
Our results show that non-random cophylogenetic patterns are widespread, but the relative contribution to diversity of codivergence vs. host switches remains unclear.The importance of host switches in generating diversity does not appear to be universal, and may be context or taxon dependent, or vary over time (Hay et al., 2020;Navaud et al., 2018).Testing the relative contribution of codivergence versus host switches requires correlating the number and magnitude of host switches in given host-symbiont cophylogenies with the extent of overall cophylogeny detected.Currently, estimation of the magnitude of host switches in individual studies is problematic as the necessary phylogenetic data are typically missing.However, as access to authors' data improves, this represents a promising avenue for future studies to disentangle the relative contribution of these important processes.
To the best of our knowledge, this study represents the first quantitative appraisal of host-symbiont cophylogeny, which is a central aspect of host-symbiont evolution.With this work, we hope to initiate a new direction in the study of symbiosis, towards formal quantitative and systematic analyses that seek to address diverse questions regarding the nature of host-symbiotic evolutionary relationships.Importantly, considerable variation in congruence exists among host phylogenies and both parasite and mutualist phylogenies.Therefore, a core challenge is to identify which factors are of importance in fostering close cophylogeny.As the number and resolution of individual cophylogenetic studies increases, the power of meta-analyses will grow, thereby permitting increasingly detailed interrogations of patterns and underlying mechanisms.

AC K NOW L E DGE M E N T S
The authors thank Dr Ben Longdon, Professor Chris Bass, and Professor Angus Buckling for providing comments on an early version of the manuscript, Dr Hanne Løvlie for discussion, and three anonymous reviewers for useful comments.AH was awarded a travel grant from the Nilsson-Ehle-donation fund to support research in theoretical and applied genetics, by the Royal Physiographic Society in Lund, to establish the collaboration leading to this work.AH is currently supported by a Biotechnology and Biological Sciences Research Council (BBSRC) David Phillips Fellowship (grant number: BB/N020146/1).SN was supported by an Australian Research Council (ARC) Discovery Grant (DP200100367).

AU T HOR S H I P
RP conceived the study, all authors designed the study.SN performed the analyses and described methodologies (wrote the Supporting Information, commented by AH and RP).AH collected the data and drafted the manuscript.RP and SN provided comments and contributed to the final manuscript.

COM PET I NG I N T E R E ST S
The authors declare no competing interests.

PE E R R EV I EW
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ele.13757.

DATA AVA I L A BI L I T Y STAT E M E N T
A copy of the source data and R code used in this study have been deposited at figshare.com:https://doi.org/10.6084/m9.figshare.14393309

OPE N R E SE A RC H BA DGE S
This article has earned Open Data and Open Materials badges.Data and materials are available at: https://figsh are.com/artic les/datas et/Datas et_for_ELE_EV_ELE13 757/

F
I G U R E 1 (A) Examples of mutualism and parasitism from across biological diversity: (a) Ectomycorrhizal fungi (white) gain carbohydrates, host spruce tree gains increased root (brown) water and nutrient absorption; (b) Figwasps gain habitat and food for young, host fig tree gains pollination; (c) Rumen bacteria gain habitat and food, host cow gains nutrition from digesting excess bacteria; (d) Clownfish gain protection from fish predators, host anemone gains protection from anemone predators.Examples of parasitism: (e) HIV uses host cellular machinery to replicate; (f) Tongue-eating isopods feed on host tissue; (g) Human head lice feed on host blood; (h) Wheat leaf rust fungi steal host resources to reproduce.Credits: (a) André-Ph.D. Picard, (CC BY-SA 3.0).(b) Mike Gordon, (CC BY-NC-ND 2.0).(c) Daniel Schwen, (CC BY-SA 3.0).(d) Nick Hobgood, (CC BY-SA 3.0).(e) CDC/C.Goldsmith, P. Feorino, E. L. Palmer, W. R. McManus, Public Domain.(f) Elkin Fricke, (CC BY-SA 4.0).(g) Gilles San Martin, (CC BY-NC-ND 2.0).(h) USDA/James Kolmer, Public Domain.(B) Hypothetical host and symbiont evolutionary trees, illustrating examples of perfect phylogenetic congruence between host and symbiont trees (top), and incongruence between host and symbiont trees (bottom).Narrow grey lines indicate which hosts harbour a particular symbiont.Modified fromPoulin (2007)

F
Congruence and incongruence in host-symbiont phylogenies: (A) Hypothetical host and symbiont evolutionary trees that perfectly mirror each other, resulting in perfect phylogenetic congruence.Grey lines indicate which hosts harbour a particular symbiont.(B) Host switching of a symbiont lineage to a new host lineage.(C) Extinction of a symbiont lineage.(D) Independent divergence of a symbiont lineage without host.(E) Failure of a symbiont lineage to diverge with its host lineage.(F) Hypothetical host and symbiont evolutionary trees that do not mirror each other, as a consequence of disruptive events (i.e.processes b-e listed above), resulting in phylogenetic incongruence.(G) Repeated host switches followed by divergence events that mirror host phylogeny, resulting in imperfect congruence

F
I G U R E 3 Phylogenetic congruence between hosts and symbionts (having grey and black backgrounds, respectively, in cartoons).Orchard plots show group-wise mean(s) with their 95% confidences intervals (thick lines) and 95% prediction intervals (thin lines) and observed effect sizes based on various sample sizes split by: (A) overall; (B) mode of symbiosis; (C) host taxa; and (D) symbiont taxa.The confidence interval indicates uncertainty around our overall estimate, while the prediction interval shows a possible range for a new effect size (without sampling errors;Nakagawa et al., 2021) ; Egger's regression: b [sampling se] = 1.389,CI = 0.944: 1.833]; see also Appendix S4, Figure S4.2-S4.3 and Table

F
I G U R E 4 Phylogenetic congruence between hosts and symbionts (having grey and black backgrounds, respectively, in cartoons).Orchard plots show group-wise mean(s) with their 95% confidences intervals (thick lines) and 95% prediction intervals (thin lines) and observed effect sizes based on various sample sizes split by: (A) site of symbiosis; (B) transmission mode; (C) type of symbiosis × transmission mode.The confidence interval indicates uncertainty around our overall estimate, while the prediction interval shows a possible range for a new effect size (without sampling errors;Nakagawa et al., 2021) Residual funnel plots: (A) the residuals from the meta-regression with the four moderators (symbiosis, host_tax_broad, mode_ of_tansmission_broad, and log(host_trange_link_ratio); (B) the residuals from the meta-regression model which is the same as before but with sampling standard error (SE) as an extra moderator (i.e.,Egger's regression)