Performance of a phylogenetic independent contrast method and an improved pairwise comparison under different scenarios of trait evolution after speciation and duplication

Despite the importance of gene function to evolutionary biology, the applicability of comparative methods to gene function is poorly known. A specific case which has crystalized methodological questions is the “ortholog conjecture”, the hypothesis that function evolves faster after duplication (i.e., in paralogs), and conversely conserved between orthologs. Since the mode of functional evolution after duplication is not well known, we investigate under what reasonable evolutionary scenarios phylogenetic independent contrasts or pairwise comparisons can recover a putative signal of different functional evolution between orthologs and paralogs. We investigate three different simulation models, which represent reasonable but simplified hypotheses about gene function (our “trait”) evolution. These are time dependent trait acceleration, correlated changes in rates of both sequence and trait evolution, and asymmetric trait jump. For each model we tested phylogenetic independent contrasts and an improved pairwise comparison method which accounts for interactions between events and node age. Both approaches loose power to detect the trend of functional evolution when the functional trait accelerates for a long time following duplication, with better power of phylogenetic contrasts under intermediate scenarios. Concomitant increase in evolutionary rates of sequence and of trait after duplication can lead to both an incorrect rejection of the null under null simulations of trait evolution, and a false rejection of the ortholog conjecture under ortholog conjecture simulations by phylogenetic independent contrasts. Improved pairwise comparisons are robust to this bias. Both approaches perform equally well to trace rapid shift in traits. Considering our ignorance of gene function evolution, and the potential for bias under simple models, we recommend methodological pluralism in studying gene family evolution. Functional phylogenomics is complex and results supported by only one method should be treated with caution.

simplified hypotheses about gene function (our "trait") evolution. These are time 23 dependent trait acceleration, correlated changes in rates of both sequence and trait 24 evolution, and asymmetric trait jump. For each model we tested phylogenetic 25 independent contrasts and an improved pairwise comparison method which accounts 26 for interactions between events and node age. improved pairwise method using three simple simulation models with different 93 parameters of divergence following speciations and duplications. 94 We simulated ultrametric trees (n = 10000), each with 100 tips generated under a pure 98 birth-death model using the TreeSim R package (Stadler 2011) with a fixed speciation 99 rate of 0.4, and an extinction rate of 0.1. Since our study aims to mimic simple 100 evolutionary scenarios, where standard trait evolutionary models are frequently applied 101 on comparatively smaller phylogenies (Chira & Thomas 2016), we limited our tree size 102 to100 tips. We call these simulated time trees as "Original time trees" (Fig 1A), with 103 the branch lengths corresponding to units of time (i.e. Million Years -My). For 104 consistency, we used the same set of 10000 original time trees for all further analyses. 105

Annotation of internal node events 106
We annotated 99 (number of tips -1) internal node events ("speciation" or 107 "duplication") so that each tree had at least one speciation and one duplication node 108 events. We considered 3 fixed proportions of duplication events (number of duplication 109 events / total number of internal nodes): 0.2, 0.5, and 0.8 in this study (S1 Fig). At first, 110 we randomly annotated "duplication" events to the internal nodes based on the 111 proportion of duplication events we considered for the tree set. We then assigned 112 "speciation" events to the rest of the nodes of the tree. We, thereby, had 3 sets of 10000 113 original time trees with 3 different proportions of speciation and duplication events (S1 114  simulated trait values at tips. We used a simulated tree with 20 tips for illustration. 120 Branches following gene duplication are asymmetrically painted using the phytools R 121 package (Revell 2012) to introduce different sequence evolutionary rates to those 122 branches. This demonstrates how the scale (i.e. the branch lengths) of the calibrated 123 tree can differ from the original time tree. 124 all our main text analyses to compare the performances of different approaches in 130 testing the ortholog conjecture on different simulation models. To verify that our results 131 did not depend on the number of duplication events in a tree, we reproduced results 132 using proportions of duplications of 0.5 and 0.8 (supporting materials). 133

Simulation of trait 134
In this study, we considered a trait with a range between 0 to 1. We call it, ! like the 135 tissue-specificity score used in several previous studies of the ortholog conjecture 136 For modeling of "Time dependent trait acceleration", we wrote our own functions to 148 paint tree with different colors to simulate trait on each original time tree using a BM. 149 In this model, we used a 5 times higher trait evolutionary rate over a specific period of 150 time following gene duplication. The trait evolutionary rates do not vary between events 151 under the null hypothesis. 152 "simulate.rates" function to "simulate.rates_heterogeneous" function to introduce 155 different sequence evolutionary rates to the different node events. To avoid any bias in 156 our inference, we used the sequence evolutionary rates, K of 0.02, as used for 157 simulation of the trait ! (( ! = 0.02) in a null scenario. These trees are "Substitution 158 trees" (Fig 1B), where the branch lengths of each tree representing substitution rates. 159 To generate "Rates of sequence and trait evolution model", we considered a multiplier 160 of rate parameter (K = 0.02) to a randomly chosen daughter branch of a chronogram 161 following gene duplication. 162 An advantage of simulations is that we have the original time trees. With real data, we 163 never have these original time trees. Rather we use "Substitution trees" (Fig 1B) to 164 generate pseudo time trees, i.e. "Calibrated time trees" (Fig 1C), for hypothesis testing.  However, due to simulated tree structure, very few trees passed time calibration step 170 when we fix focal speciation nodes, and their ages across all trees. Hence, we used the 171 ages of all the speciation nodes of a tree to time calibrate them. The scale of the original 172 tree drastically changes after calibration if we do not fix the age of the root node. It may 173 impact the inference of phylogenetic method due to its branch length dependence. 174 Hence, we also used the age of the root of each simulated original time tree to maintain successfully time calibrated for 0.2, 0.5, and 0.8 proportion of duplication events, 177 respectively. 178

Phylogenetic analyses 179
For each internal node of a gene tree, we calculated the PIC of simulated trait ! from 180 the APE R package using the following formula: 181 where PICn is the node of which PIC is to be computed, ! # , and ! ! are the trait values

Pairwise analyses 188
We considered all possible pairwise comparisons for each tree to obtain Pearson 189 correlation coefficients R between ! values of events, as was done by Dunn et al. distinguish the effects of events on correlation coefficients. In contrast, our improved 198 pairwise comparison considers an ANCOVA analysis (i.e. aov(R ~ My + Event + My*Event)) to find out the significance level. In all tests of our improved pairwise 202 approach, we included the interaction terms between events and node age, and 203 compared with model without interaction term to determine the best model among 204 them. Since none of our models supported linear fit as best fit model, we used 205 polynomial (quadratic or cubic) regression in our improved pairwise analyses. 70% of 206 our data is used to train the model, and then we used them to predict the value of rest 207 of the 30% test set data. Based on the obtained root mean square error, we chose the 208 best fit (quadratic or cubic) polynomial model for each simulation study to avoid curve 209 over fitting problem. In addition, we considered repeated 10 fold cross validation 210 approach to choose best polynomial model. 211

Statistical significance 212
For any method, we considered a result is significant if P < : = 0.05. Since we used a 213 large number of tests (n = 84) in our study, we additionally applied a Bonferroni 214 correction where a result is considered to be significant when P < :′, where :′ = 0.05 / 215 84 = 5.95e -4 . 216

matters-in-gene-family-evolution. 219
Immediately after duplication, both duplicates might experience accelerated evolution 222 for a period of time due to relaxed selection (Rogozin 2014). This can be represented 223 by a simple model of acceleration of trait evolutionary rate for a given time after 224 duplication (Fig 2). After that period of acceleration, the evolutionary rates return to 225 the pre-duplication level. In this model, speciation does not have any impact on 226 evolutionary rates: a speciation during the time period of the duplication-caused 227 acceleration does not revert the gene to a "speciation" rate. To study our capacity to 228 recover evolutionary signal under this model, we simulated with different durations of 229 accelerated evolution, and different proportions of duplications to speciations in the 230 gene trees (Fig 2). 231 In this comparative study, we used an improved pairwise comparison due to the low 232 power of a simple pairwise comparison approach as used by Kryuchkova-Mostacci and 233 Robinson-Rechavi (2016) (S1 Table). First, both methods recover the proper 234 evolutionary signal under the null, with no difference of PIC between duplication and 235 speciation branches (Fig 3A), and no difference between ortholog and paralog pairs 236 ( Fig 3B). Second, when there are few duplications in proportion to speciations, both 237 approaches recover the evolutionary signal, namely faster evolution of paralogs or of 238 duplication branches (Fig 3C-F, S1 Table). This is despite the acceleration lasting 50% 239 to 90% of the time since the oldest duplication. It is possible because the rarity of 240 duplications provides many speciation branches without acceleration (Fig 2). While the 241 performance of both methods at recovering the evolutionary signal decreases when the 242 proportion of duplication nodes increases, the PIC is more robust than the pairwise an acceleration lasting 70% to 90% of old duplication age, neither approaches can 246 recover the signal of the ortholog conjecture (S1 Table). Under such scenarios, most 247 speciation branches or ortholog pairs have evolved mostly or entirely under the 248 influence of a duplication-caused acceleration, and there is no signal to recover. on the other hand is that the null is not always properly recovered. When the trait 276 evolutionary rate is not impacted by duplication, but the sequence evolutionary rate is 277 (scenario 2 of Table 1), the PIC approach wrongly rejects the null for trait evolution. 278 This is due to time tree calibration using sequence evolution as a reference. Indeed there 279 is no difference in PIC when the original time tree is used (A in S2, C in S3 and S4 280 Figs). However, the higher sequence evolutionary rates after duplication lead to higher 281 expected variance, and thus the calibrated time trees produce lower contrasts for 282 duplication than for speciation, with a significant difference in PIC (A in S2, C in S3 283 and S4 Figs). Moreover, if both the trait and the sequence evolutionary rates accelerate 284 after duplication events, but the sequence acceleration is higher (scenario 4c of Table  285 1), the PIC can again estimate lower contrasts for duplication branches with calibrated 286 time trees (Fig 4E, K in S3 and S4 Figs). This can lead to wrongly accepting the null if 287 a one-sided test is used, or rejecting the null in the wrong direction. This bias does not 288 impact the pairwise approach, under neither of these scenarios (Fig 4F, B Table 1. C. 300 speciation in all 3 tree models, and rejects the null hypothesis. D. Higher correlation 302 coefficients for orthologs than for paralogs also support the conjecture in this case. E-303 F. The ortholog conjecture scenario 4c of Table 1, where both the rates ratios varies but 304 with an excess of sequence evolutionary rates ratio than trait ratio following duplication 305 than speciation. E. 2 folds higher trait evolutionary rates following duplication supports 306 the conjecture for the original time tree. 8 folds higher sequence evolutionary rates after 307 duplication generate 8 times higher expected variances for duplications, and produce 308 lower contrasts for duplications than speciation in the substitution rates trees. This led 309 us to obtain 2/8 times lower contrasts for duplications than for speciations in the this model (S1 Table). Thirdly and most worryingly, when both sequence and trait 371 produced the expected signal under the null (S1 Table) Table 1), the PIC method fails to detect the signal of ortholog conjecture 403 ( Fig 4E, K in S3 and S4 Figs). In this case, improved pairwise comparisons seem to be 404 a better choice due to its branch lengths independence. 405 The only model for which the ratio of duplication to speciation events impacts the 406 results is with time-dependent trait acceleration. When there are many duplications, and 407 the duplication-induced acceleration lasts for 70% to 90% of old duplication age, most 408 speciation branches also evolve under the influence of trait acceleration like duplicates. 409 Hence, both the approaches fail to recover the signal of the ortholog conjecture in a real 410 ortholog conjecture scenario (S1 Table). Under intermediate scenarios, PIC allows to 411 better recover the signal of the ortholog conjecture than pairwise comparisons. 412 Overall the main message of this study is that our ignorance of gene function evolution 413 makes it difficult to choose one method as the gold standard. Each method has its own 414 limitations, and we call for methodological pluralism, at least in the present state of our 415 knowledge. Relying on a single approach to interpret result can be problematic in many 416 cases of comparative functional genomics study involving complex gene family 417 evolution. These results apply in principle not only to the ortholog conjecture, but to 418 any other cases of gene trait evolution (e.g., horizontal gene transfer), where evolution 419 might be gradual or by jumps, affect a more or less large proportion of branches of the 420 gene tree, and be confounded with changes in sequence evolutionary rates. In  Table 1: Hypothetical scenarios of the rates of sequence and trait evolution model. 438 439 rT: Ratio of trait evolutionary rates, and rK: ratio of sequence evolutionary rates 440 between speciation and duplication. Ratio as 1 refers to no change, while ratio > 1 refers 441 to higher rates following gene duplication than speciation in the corresponding rate 442 parameter(s). 1: Both the rate parameters remain same under a null, 2. The trait 443 evolutionary rate is not affected by duplication, but the sequence evolutionary rate is in     The branch lengths of the trees are in million years, My. For A-C, the proportion of 609 duplications are 0.2, 0.5, and 0.8, respectively. This means that out of 99 internal nodes 610 20 nodes, 50 nodes, and 80 nodes are duplication events, respectively in A-C. Rest of 611 the nodes are annotated as speciation nodes. We used an asymmetric model, where after duplication one among the two duplicates 615 experience acceleration. A quadratic polynomial is used in our improved pairwise 616 approach. A-B. Scenario 2 of Table 1, A. Due to no variation in trait evolutionary rates 617 following speciation and duplication events, the original time trees recover an expected 618 signal of no difference in contrasts. An excess of sequence evolutionary rates following 619 duplication than speciation produces lower PICs for duplications compared to 620 speciations for the substitution rate trees, and for the calibrated time trees. This leads to 621 a false rejection of the null in a real null condition for the calibrated trees. B. The 622 improved pairwise approach shows the expected trend by showing no difference in 623 correlation coefficients between events in this case. C-F. The ortholog conjecture 624 scenario 4a and 4b of Table 1. When the trait evolutionary rates ratio is higher than (C-625 D), or equals to (E-F) the sequence evolutionary rates ratio of events, higher changes 626 in trait rates following duplications make both the approaches to recover the expected 627 signals in support of the conjecture.