Coevolution leaves a weak signal on ecological networks

. One of the major challenges in evolutionary ecology is to understand how coevolution shapes species interaction networks. Important topological properties of networks such as nestedness and modularity are thought to be affected by coevolution. However, there has been no test whether coevolution does, in fact, lead to predictable network structure. Here, we investigate the structure of simulated bipartite networks generated under different modes of coevolution. We ask whether evolutionary processes in ﬂ uence network structure and, furthermore, whether any emergent trends are in ﬂ uenced by the strength or “ intimacy ” of the species interactions. We ﬁ nd that coevolution leaves a weak and variable signal on network topology, particularly nestedness and modularity, which was not strongly affected by the intimacy of interactions. Our ﬁ ndings indicate that network metrics, on their own, should not be used to make inferences about processes underlying the evolutionary history of communities. Instead, a more holistic approach that combines network approaches with traditional phylogenetic and biogeographic reconstructions is needed.


INTRODUCTION
Dynamics of ecological communities are fundamentally shaped by their networks of interacting species (Thompson 2005). The study of these networks, which involves the classification of interspecific relationships and the strength of their reliance upon one another, has important implications for ecology and conservation-informing, for example, the ability of communities to maintain ecosystem function in the face of disturbance. Specifically, interaction patterns are thought to affect the resilience of a network to disturbances such as fluctuating species abundances, the introduction of new species, or the extinction of existing species.
Ecological networks are frequently characterized using two network topology metrics-modularity and nestedness (Bascompte 2010, Th ebault andFontaine 2010). Modular community interactions are more insular, occurring within separate groups or "modules" more often than between modules. Conversely, in nested networks, a few species interact with many species while the other species interact with progressively smaller subsets of those species. Both mutualistic and antagonistic networks can exhibit a high degree of nestedness (e.g., plants and pollinators , hosts and parasites [V azquez et al. 2005] or modularity (e.g., plants and pollinators [Olesen et al. 2007], hosts and parasites [Krasnov et al. 2012], plants and seed dispersers [Donatti ❖ www.esajournals.org 1 April 2017 ❖ Volume 8(4) ❖ Article e01798 et al. 2011]. Both nested and modular communities may also exhibit phylogenetically structured interactions, where closely related species have more similar interaction patterns than distantly related species (Cattin et al. 2004, Rezende et al. 2007, Donatti et al. 2011. A strong phylogenetic signal indicates that network patterns are constrained by past evolutionary history (e.g., as a result of trait conservatism). Coevolutionary dynamics likely play a critical role in determining how communities are structured, because coevolution shapes and maintains the traits involved in species interactions. Coevolution can lead to interaction modules either because phylogenetically closely related species will often have similar trait values (Lewinsohn et al. 2006) or because convergent evolution among phylogenetically distantly related species will lead to similar trait values (Donatti et al. 2011, Danieli-Silva et al. 2012. Coevolution can alternatively lead to nestedness if phenotypes constrain interactions (i.e., some links are "forbidden", Jordano et al. 2003, Santamar ıa and Rodr ıguez-Giron es 2007). Coevolution can also affect the degree of dependence between interacting lineages (the "intimacy" of species associations) which can also affect network structure indirectly (Ollerton 2006, Fontaine et al. 2011. Selection that leads to the avoidance of sharing interaction partners (socalled partner overlap) and biological limits that influence how many partners a species can have (interaction niche breadth) are expected to lead to more intimate interactions. For example, antagonistic coevolution favors parasites that are able to overcome host defenses, which can subsequently lead to parasites infecting only a few mutually exclusive host species (e.g., overcoming one host's defenses might come at the expense of the ability to overcome other host's defenses; Thompson 2006, Fontaine et al. 2011). In contrast, more generalized interactions with considerable partner overlap may evolve if selection favors phenotypes that are compatible with traits common to many species in the community, as might be the case for free-living mutualists (Thompson 2005, Fontaine et al. 2011. Actual networks exhibit levels of interaction intimacy that lie along a continuous spectrum, with varying degrees of interaction niche breadth and partner overlap (Fontaine et al. 2011).
A variety of other processes can also affect species interaction patterns including ecological dynamics (Krause et al. 2003, Bastolla et al. 2009, Th ebault and Fontaine 2010, spatiotemporal species distributions (Pimm et al. 1991, V azquez et al. 2009, Pillai et al. 2011, Encinas-Viso et al. 2012, stochastic events (Krishna et al. 2008, Canard et al. 2012, Jonhson et al. 2013), or combinations thereof (V azquez et al. 2009, Suweis et al. 2013. Understanding which ecological, evolutionary, and stochastic processes have consistent effects on interaction patterns (i.e., in determining the degree to which communities are nested or modular) is critical for elucidating the mechanisms that structure communities, as network topology correlates with stability. Compared to randomly assembled communities, nested mutualistic communities may be more resilient to species extinction and fluctuations in abundance (Memmott et al. 2004, Fortuna and Bascompte 2006, Burgos et al. 2007, Th ebault and Fontaine 2010, Valdovinos et al. 2013, but see Allesina andTang 2012, James et al. 2012). Modular interaction patterns, particularly in food webs, are also thought to be more stable than random interactions (May 1972, Krause et al. 2003, Stouffer and Bascompte 2011 because fluctuations in species abundance are largely contained within modules and are less likely to spread throughout an entire community (Krause et al. 2003, but see Lawton 1980, Th ebault andFontaine 2010). Additionally, modeling work suggests that when interactions are phylogenetically structured, species extinction may trigger extinction cascades of related species and a more pronounced loss of taxonomic diversity (Rezende et al. 2007).
Recently, several authors have used theoretical models to explore the effect of evolution on networks (Loeuille 2010, Nuismer et al. 2013. However, there have been no explicit tests of how different evolutionary assembly mechanisms may lead to predictably different network topologies. Here, we simulate coevolution and community assembly of two interacting lineages of species under four different modes of coevolution, each of which corresponds to a limiting case of the possible coevolutionary mechanisms (coevolution with cospeciation, coevolution without cospeciation, cospeciation without coevolution, and neither coevolution nor cospeciation) to test whether coevolution leads to predictable ❖ www.esajournals.org 2 April 2017 ❖ Volume 8(4) ❖ Article e01798 network structure, specifically nestedness, modularity, and phylogenetic interaction signal. We also examine whether interaction intimacy, measured as the range of interaction niche breadth and partner sharing between species, accentuates the coevolutionary signal of these different assembly mechanisms on the structure of networks.

Community generation approach
Our network generation process proceeds as follows: 1. Simulate two phylogenetic trees for interacting lineages 2. Simulate evolution of traits along those trees 3. Build interaction networks, using species' trait values to determine which pairs of species can interact 4. Evaluate the nestedness, modularity, and interaction phylogenetic signal of networks Depending on the evolutionary scenario under consideration, the mechanisms for tree construction and trait evolution vary in order to simulate the corresponding coevolutionary and cospeciation processes (Steps 1-2). All simulations were conducted in R version 3.3 (R Core Team 2008). All simulation code is available at https:// github.com/lponisio/network_assembly.git.
We do not explicitly make a distinction between predatory, mutualistic, parasitic, or commensal interactions. These scenarios differ in how species interactions affect fitness but, because species' fitnesses are not a component of our data simulation process, our conclusions should apply to networks containing interactions of all types. To facilitate biological intuition, however, we will develop our methods using language most appropriate for mutualistic interactions.

Phylogeny and trait generation
We begin by generating phylogenetic trees for a fixed number (15) of resource-providing species (e.g., plants) and resource-seeking species (e.g., pollinators). We do this using a birth-death process which approximates diversification of clades by assuming homogeneous rates of speciation (birth, k) and extinction (death, l) across taxa and time (Stadler 2012). We next simulate trait evolution along these trees as a random walk through trait space (Brownian motion model; Paradis 2012), and in one scenario (see Coevolution without cospeciation), we modify the Brownian motion model to include selective optima that assert an attracting force on random walk trait evolution (Ornstein-Uhlenbeck model [Gillespie 1996, Nuismer andHarmon 2015]). According to the Brownian motion model, a trait X evolves at a rate r: where W(t) is drawn at random from a normal distribution with mean 0 and variance r 2 . Similarly, the Ornstein-Uhlenbeck model is a random walk in which trait values are pulled toward some "optimal" value, h, with an attraction strength proportional to the parameter a. The model has the following form: The four scenarios under which we simulate trait evolution are as follows: 1. Coevolution and cospeciation ( Fig. 1a): High levels of species-level interaction specificity leading to tight coevolution and cospeciation and, consequently, evolution of congruent phylogenies (Farenholz 1913). The plausibility of this scenario has been extensively criticized for its lack of realism (e.g., Klassen 1992 [Althoff et al. 2012]). Different parasite groups, including macro-parasites (internal and external) and their vertebrate hosts also do not exhibit a strong signal of cospeciation (for a review, see Hoberg and Brooks 2008). However, this scenario provides a useful benchmark, because communities assembled under this scenario represent an extreme outcome of coevolution. We simulate this mode of coevolution by generating communities where the tree topology and trait values (simulated under Brownian motion, Eq. 1) of interacting taxa are identical (Fig. 1a).
2. Coevolution without cospeciation ( Fig. 1b): Trait complementarity is the result of natural selection honing the traits involved in interactions across many lineages (e.g., corolla tubes of Interacting communities are composed of phylogenetically related and unrelated species that have converged on similar traits (Olesen et al. 2007, Donatti et al. 2011, Danieli-Silva et al. 2012, Krasnov et al. 2012. This pattern of trait convergence is often observed in interacting species, such as patterns of fruit design in unrelated plant species (Jordano 1995), oil production in orchids (Ram ırez et al. 2011), and ectoparasites of mammals (Krasnov et al. 2012). To investigate this scenario, we generate communities where speciation of the interacting partners occurs independently (and thus phylogenetic topologies differ between the two lineages) but trait distributions of interacting species are probabilistically matched (Fig. 1b). To link the trait distributions, the traits of one interacting lineage are first simulated under Brownian motion evolution (Eq. 1). These traits are then used as lineage-specific optima (h, Eq. 2) in an Ornstein-Uhlenbeck process to generate the traits in the other interacting lineage. We kept the strength of attraction to the optima (a, Eq. 2) constant at 1. To approximate reciprocal selection, we then use the traits of the second lineage as optima in another Ornstein-Uhlenbeck process to regenerate traits for the first lineage. Each lineage's traits are thus used as optima for the other lineages to evolve toward. 3. Cospeciation without coevolution (Fig. 1c): Concurrent allopatric speciation of lineages across a shared biogeographic and climatic landscape leads to congruence between phylogenies of interacting species, without reciprocal selection on species' traits (Smith et al. 2008, Althoff et al. 2012. For example, the phylogenetic congruence between yuccas and yucca moths is likely due to shared biogeographic history rather than coevolution (Althoff et al. 2012). When species co-occur geographically, it is possible that they also coevolve. Here, however, we present an extreme case where no coevolution takes place so that we can quantify the effect of cospeciation alone. To model this, we generate communities where interacting taxa share the same tree topology, but the traits of the partners' traits evolve independently via Brownian motion evolution (Eq. 1, Fig. 1c).
4. No coevolution, no cospeciation ( Fig. 1d): Speciation and trait evolution proceed independently of one another. Such independent evolution could occur if environmental or developmental constraints on the traits involved in interactions (e.g., the body size of pollinators and the flower size of plants) enable interactions to occur without coevolution taking place. To construct communities under this scenario, we generate trees and trait values of interacting taxa independently, the latter via Brownian motion evolution (Eq. 1, Fig. 1d).
The shape of phylogenetic trees has the potential to affect the structure of networks (Chamberlain et al. 2014), so we generated birth-death trees for a range of extinction and speciation rates (Fig. 2) for all of the evolutionary scenarios. Specifically, we simulated trees with primarily deep or primarily shallow divergences (Fig. 2). We also tested whether the trait distribution affected network topology by altering the variance of the Brownian motion process (r 2 , Eq. 1) while holding the tree age constant (because trait variance and tree age are linearly related under a Brownian motion model).  Pires et al. 2011). Depending on the variance in trait values and the structure of the phylogenetic tree, the distribution of traits was normal to right-skewed, which is consistent with empirical findings (Stang et al. 2006(Stang et al. , 2009). In addition, our formulation of trait ranges reflects the observation that species with more extreme trait values tend to ❖ www.esajournals.org 5 April 2017 ❖ Volume 8(4) ❖ Article e01798 have relatively narrow or wide ranges of partners with which they can interact (Stang et al. 2009, Warren et al. 2010.

Linkage rules and interaction intimacy
The degree of interaction intimacy in a community is determined by the amount of partner overlap in the community and the interaction niche breadth of the constituent species. We varied both components of interaction intimacy by changing the specificity of the interactions which we did by changing the maximum acceptable difference in trait values that would still enable two potential partners to interact (Appendix S1: Fig. S1). The range of acceptable trait values for where h i is the trait value of species i , and d indicates how similar a partner's trait value must be (relative to h i ) in order to enable an interaction. We varied d from 0.1, which indicates that individuals only accept a narrow range of trait values in their partners, to 1, which indicates that they accept a wide range. d above 1 or below 0.1% yielded interaction matrices where nearly all or nearly none of the species interacted and, thus, they were not considered. We also generated weighted networks, in which partners interact in proportion to the amount of overlap in their trait ranges, and unweighted networks, in which partners interact equally frequently with all others whose trait ranges overlap (x i À x j ).
We calculated niche breadth as d, the mean proportion of potential interaction partners a species interacts with, from the total possible ❖ www.esajournals.org 6 April 2017 ❖ Volume 8(4) ❖ Article e01798 (i.e., the mean species degree). We calculated the partner overlap of species belonging to the same interaction lineage as s ij , the mean pairwise similarity of interaction partners using the Gower metric (Gower 1971, Oksanen et al. 2017). Because we assume that interactions require complementarity of traits between partner species (Santamar ıa and Rodr ıguez-Giron es 2007, Stang et al. 2007, Nuismer et al. 2013, traits of interacting lineages are positively correlated in each of the coevolutionary scenarios. For hosts and parasites, this positive correlation might correspond to a match between traits that govern susceptibility in hosts and mode of attack in parasites. Traits such as the depth of nectar tube in a flower and the length of a pollinator's tongue would also be expected to positively correlate when interactions are mutualistic. When simulating phylogenies, we began with a 1:1 ratio of resource-providing species to resource-seeking species. We did this in order to facilitate comparisons between scenarios where there is necessarily a one-to-one congruence (e.g., the coevolution and cospeciation scenario) with those where there is not. However, species that evolve to possess trait values that prohibit them from interacting with any members of the community were not included in the network. Thus, the communities vary in the number of species in each interacting lineage. We would expect this feedback between species' traits and community composition when species with ill-suited trait values are unable to persist in that community.

Characterizing network topologies
For each of the resulting interaction networks, we calculated topological descriptors, focusing in particular on nestedness and modularity. We use NODF (weighted or unweighted) to evaluate network nestedness (Almeida-Neto et al. 2008). The NODF metric is large when species with fewer partners interact with subsets of partners with which more connected species interact (Almeida-Neto et al. 2008). Many methods exist for partitioning networks into sub-communities for modularity computation, and all have potential pitfalls (Fortunato 2010). We therefore considered three community partitioning methods: (1) a dynamic algorithm via a random walk (Pons andLatapy 2005, Csardi andNepusz 2006), (2) a greedy modularity optimization (Clauset et al. 2004, Csardi andNepusz 2006), and (3) a hierarchical clustering algorithm Girvan 2004, Csardi andNepusz 2006).
We estimated the strength of the phylogenetic interaction signal as the correlation between the evolutionary time separating species and the dissimilarity of their interaction partners. The evolutionary distance between two species is proportional to their phylogenetic distance-the sum of the branch lengths connecting those species (Paradis et al. 2004). We measured the dissimilarity of interaction partners by calculating the relative overlap of interaction partners between pairs of species using the Gower metric (Gower 1971, Oksanen et al. 2017. Following Rezende et al. (2007) and Cattin et al. (2004), we then calculated phylogenetic interaction signal as the correlation between the phylogenetic distance matrix and the interaction dissimilarity matrix using a Mantel test. A high correlation indicates that species with older divergence times are less likely to share interaction partners.
To calculate the effect of coevolution, we selected and compared random pairs of communities that shared the same cospeciation history with and without coevolution (Fig. 1a, c and b,  d). By comparing the observed metrics to those derived from null communities, we can quantify how much the observed level of structure can be explained by the simple rules that determine the probability that two species interact in the null communities. To assemble 99 null communities, we reshuffled the interactions between species but fixed the total number of interactions and species and, in quantitative networks, the distribution of the interaction frequencies (Galeano et al. 2009). We then calculated the mean and standard deviation of the null community metrics, X null and r null , respectively. By standardizing the metrics using z-scores, ðX obs À X null Þ=r null , where X obs is the observed metric, we can compare the metrics between networks of different numbers of interactions and species. Simulating >99 for all analyses presented here was too computationally intensive. However, we verified that conclusions did not differ qualitatively for a subset of scenarios when using 999 vs. 99 null communities. We then calculated the difference between metric z-scores between pairs of communities that had cospeciated but differed in whether or not they exhibited coevolution (Fig. 1a-c and b-d). The ❖ www.esajournals.org 7 April 2017 ❖ Volume 8(4) ❖ Article e01798 mean and standard deviation of these differences across replicate pairs of communities was used as an estimate of the effect of coevolution.

RESULTS
We found that the coevolutionary history of a community can have consistent and detectable effects on the structure of interactions ( Fig. 3; Appendix S1: Tables S1-S3). Specifically, when interaction probabilities were not weighted by the degree of overlap between partner traits, communities that coevolved and cospeciated were more modular and anti-nested (less nested than expected by randomly assembled communities; Fig. 3a, c). This is likely because cospeciating In left panels, all species interactions were equally frequent, whereas in right panels the interaction probability depended on the degree of trait overlap between potentially interacting pairs of species. Relative nestedness and modularity are corrected for the null expectation, given the number of interactions and species in the interaction network (z-scores, top two panels). Scores greater than $ 2 or less than $ À2 are significantly more or less structured than randomly assembled communities. Vertical bars denote 95% confidence intervals. The modularity metric depicted here was calculated using a hierarchical clustering algorithm. ❖ www.esajournals.org 8 April 2017 ❖ Volume 8(4) ❖ Article e01798 clades form interacting modules (e.g., see Fig. 1a). The metric used for calculating modularity (edge-betweenness, greedy optimization, and random walk) did not qualitatively affect the relative differences between coevolutionary scenarios (Appendix S1: Fig. S2). Across all coevolutionary scenarios, communities were, on average, significantly modular but not significantly nested (Fig. 3a-d).
In addition, communities that cospeciated had the highest phylogenetic interaction signal (Fig. 3e, f; Appendix S1: Table S3). That a strong phylogenetic interaction signal characterizes cospeciating communities is likely a consequence of the fact that shared tree topologies between interacting lineages constrain trait evolution such that closely related species will have similar traits and thus share interaction partners, with or without coevolution.
The effect of coevolution did not depend strongly on phylogenetic tree structure or the variance in evolved trait values (Table 1). Because these community characteristics did not have distinct effects on the differences between coevolutionary communities, we restrict our focus to community with shallow divergences between species within an interacting lineage and a low variance in evolved trait values (Fig. 2b).
There was no detectable effect of coevolution on network structure in communities with the same cospeciation history. In both communities that cospeciated and those that did not, the mean effect size of coevolution (calculated as the difference in network metrics or phylogenetic interaction signal between communities with and without coevolution; Fig. 1a-c and b-d) was smaller than the standard deviation of the difference (Table 1). Thus, the effect of coevolution was never significantly different from zero.
In addition, at any level of interaction niche breadth or partner overlap, the effect of coevolution was highly variable and not significantly different from zero (Figs. 4 and 5). Modularity increases slightly across niche breadth in the communities that coevolved and cospeciated in comparison with those that only cospeciated. The variance around those trends is large relative to the effect sizes and, thus, it seems unlikely that they would be of empirical relevance, at least under the range of community parameters considered here. Table 1. The mean and standard deviation of the difference in modularity (z-scores), nestedness (z-scores), and phylogenetic interaction signal between coevolved and non-coevolved communities, controlling for whether the interacting lineages cospeciated (see Fig. 2 for details on the parameters governing tree divergence and the variance in evolved trait values). Notes: The effect of coevolution with cospeciation is calculated as the difference between communities that coevolved and cospeciated and those that only cospeciated. Similarly, the effect of coevolution without cospeciation is calculated as the difference between coevolved communities that did not cospeciate and those with interaction partners that did not coevolve or cospeciate. SD represents the standard deviation. The mean differences are always less than the standard deviation.

DISCUSSION
We found that coevolution leaves a weak signal on network topology in four coevolutionary community assembly scenarios. These scenarios represent extreme possibilities, ranging from coevolution and cospeciation to completely independent evolution and speciation. Real-world communities likely fall somewhere between these extremes and are, thus, likely to exhibit intermediate patterns to those reported here. Because the topological differences reported here are small, detecting signals of coevolution will likely be difficult. Thus, caution should be taken when using only the most commonly studied network topology descriptors (nestedness, modularity, and phylogenetic interaction signal) alone to make inferences about the processes that underlie the evolutionary history of a community. Fig. 4. The effect of coevolution on the relative modularity (a-b), nestedness (c-d), and phylogenetic interaction signal (e-f) across a range of levels of mean interaction niche breadth. Interaction niche breadth is measured as the mean proportion of potential interaction partners a species interacts with out of the total possible. The solid curves represent smoothed mean differences between randomly selected pairs of communities, and the filled area corresponds to 95% confidence intervals. Note that we do not directly vary interaction niche breadth, but instead do so indirectly by varying the size of the range of trait values that a species will accept in a partner (see Appendix S1: Fig. S1 for further interpretation of these values).
❖ www.esajournals.org 10 April 2017 ❖ Volume 8(4) ❖ Article e01798 By assembling communities both with and without coevolution in the same context, our approach differs from those that consider coevolutionary dynamics in isolation (e.g., Loeuille 2010, Nuismer et al. 2013. Here, coevolutionary feedbacks between phylogenies, traits, and interactions were modeled implicitly, so that the same framework could be used for communities that had and had not coevolved. Our coevolutionary communities exhibit similar structural properties to communities generated by explicit coevolutionary models of bipartite networks (Nuismer et al. 2013). We found that trait-matching-based species interactions led to anti-nested networks (Fig. 3)-this is consistent with other studies (Nuismer et al. 2013). However, many real-world communities, and in particular mutualistic communities, tend to be significantly nested in comparison with randomly assembled communities (Fontaine et al. 2011). In models where the resource-seeking species (e.g., the pollinator) has to overcome a barrier to access the reward offered by the resource species (e.g., the plant), networks tend Fig. 5. The effect of coevolution on the relative modularity (a-b), nestedness (c-d), and phylogenetic interaction signal (e-f) for simulated communities across a range of partner overlap values. Partner overlap is measured as the mean similarity of pairs of species belonging to the same interacting lineage. Curves were calculated as described in the caption to Fig. 4. The solid curves represent smoothed means, and the filled area corresponds to 95% confidence intervals.
❖ www.esajournals.org 11 April 2017 ❖ Volume 8(4) ❖ Article e01798 to be nested (Santamar ıa and Rodr ıguez-Giron es 2007, Nuismer et al. 2013). We did not, however, explore this mechanism for determining interactions because we assumed a positive correlation between the traits of the two interacting lineages which would not be valid when barrier-type mechanisms determine interactions. Our finding that different coevolutionary processes do not lead to drastically different interaction patterns in networks sheds light on why other studies have found it difficult to distinguish mechanisms (e.g., Machado et al. 2005, Cruaud et al. 2011, Ram ırez et al. 2011, Althoff et al. 2012. Until recently, coevolution and cospeciation were thought to play a critical role in structuring many one-to-one symbiotic relationships, such as those between figs and fig wasps or yucca and yucca moths. Evidence was largely based on the congruence between traits of interacting species and their phylogenies. Recent work using more synthetic approaches and diverse data sources, however, has shown that these apparent patterns were primarily caused by the shared biogeographic history of interacting lineages leading to cospeciation (Althoff et al. 2012) and/or incomplete sampling (Machado et al. 2005, Cruaud et al. 2011).
Because at least four different mechanisms give rise to the same interaction patterns, additional tests have to be devised and undertaken to assess the contribution of different assembly mechanisms. Specifically, determining whether assembly processes structure interactions will require a synthetic approach-combining network metrics with more traditional phylogenetic, trait, and biogeographic reconstructions (e.g., Ram ırez et al. 2011).

ACKNOWLEDGMENTS
We are grateful to P. R. Guimarães Jr., C. Kremen, D. Ackerly, B. Morrison, and B. Brunner for their helpful comments and insights. This work was supported by funding from the National Science Foundation (GRFP to LCP) and the Natural Sciences and Engineering Research Council of Canada (PDF to LKM), and the Berkeley Institute for Data Science (postdoctoral fellowship to LCP) which is funded by the Gordon and Betty Moore Foundation (Grant GBMF3834 to UC Berkeley) and the Alfred P. Sloan Foundation (Grant 2013-10-27 to UC Berkeley).