Revisiting the influence of learning in predator functional response, how it can lead to shapes different from type III.

Predator/parasitoid functional response is one of the main tools used to study predation behavior, and in assessing the potential of biological control candidates. It is generally accepted that predator learning in prey searching and manipulation can produce the appearance of a type III functional response. Holling proposed that in the presence of alternative prey, at some point the predator would shift the preferred prey, leading to the appearance of a sigmoid function that characterized that functional response. This is supported by the analogy between enzyme kinetics and functional response that Holling used as the basis for developing this theory. However, after several decades, sigmoidal functional responses appear in the absence of alternative prey in most of the biological taxa studied. Here, we propose modeling the effect of learning on the functional response by using the explicit incorporation of learning curves in the parameters of the Holling functional response, the attack rate (a), and the manipulation time (h). We then study how the variation in the parameters of the learning curves causes variations in the shape of the functional response curve. We found that the functional response product of learning can be either type I, II, or III, depending on what parameters act on the organism, and how much it can learn throughout the length of the study. Therefore, the presence of other types of curves should not be automatically associated with the absence of learning. These results are important from an ecological point of view because when type III functional response is associated with learning, it is generally accepted that it can operate as a stabilizing factor in population dynamics. Our results, to the contrary, suggest that depending on how it acts, it may even be destabilizing by generating the appearance of functional responses close to type I.

INTRODUCTION

Most biological systems involve an array of intricate relationships among organisms which are of paramount importance to understand the patterns of stability which according to May (1973), Gillman and Hails (1997), is the tendency of the system to stay or move away from the equilibrium, and biodiversity of communities (McCann, 2000). Therefore, it is crucial to count on reliable methods to have the best predictions and understanding of population and community dynamics, and eventually, to support wildlife management decisions (Pettorelli et al., 2015).

Functional response (Holling, 1959; Solomon, 1949) is a mathematical framework used to describe the ability of organisms to consume resources based on their availability. In this contribution, we will refer to carnivores, herbivores, parasites, parasitoids, hyperparasitoids, and some herbivores that consume the whole plant (such as phytophagous plankton that eat algae) as “predators.” The term “prey” here includes all different types of living organisms or food resources being consumed by the predator. As the survival of predators depends on their ability to exploit variable densities of prey, these organisms must be able to detect, process, and assimilate the prey as a function of its abundance; this ability is influenced by several factors. According to Holling (1966), the three basic components of the response of predators are (i) the attack rate (linked to the ability to find prey: a); (ii) the time prey is exposed: t; and (iii) the handling time (how fast a prey is consumed: h). In a classical paper, Holling (1959) characterized three types of functional response: type I response, in which the predator consumes its prey at a constant rate regardless of the prey density, and therefore it results in a linear relationship between prey density and consumption rate; the handling time is zero or near zero. Type II response (Holling’s disk equation), in which saturation occurs mostly because the handling time imposes a limit to the rate at which the prey is consumed, therefore it results in a rectangular hyperbola in which the rate of prey consumed asymptotically approaches 1/handling time. Finally, a type III response is a sigmoidal curve. The mathematical reason for the change of shape is that Holling’s disc equation now is a quadratic function of the prey density; the result is an “acceleration” of the attack rate, but it keeps the limitation caused by handling time. The component that produces this effect is learning (Holling, 1966). Real (1977) incorporated the possibility to shift between types II and III functional responses by using the enzyme kinetic models of Barcroft and Hill (1910). In Real’s approach, the attack rate depends on a Power‐Law of resource density as a = bN, where b is the attack coefficient, N is the number of preys in the environment, and q is an exponent that influences the shape of the functional response from a hyperbolic type II functional response (q = 1) to a strict type III functional response (q = 2) and beyond these bounds.

The population consequences of each type of response are different, for instance, the stability of predator and prey populations strongly depends on whether predator consumption rates increase linearly (type I functional response) or follow a saturating function (type II and III functional responses) with prey densities (Hastings, 2013). Type III functional response is assumed to be able to stabilize predator‐prey systems (Hassell, 1978; Hassell et al., 1977; Murdoch & Oaten, 1975; Rall et al., 2008) since its lower efficiency at low prey densities would allow the prey population to recover from population bottlenecks and, in consequence, avoid local extinctions. While at high densities, it would increase the speed of consumption, helping to avert outbreaking‐type dynamics. However, the relationship between functional response types and stability is not simple. On the other hand, stability predictions differ depending on whether functional response parameters are derived. Several examples in predator‐prey systems were recorded for a type III functional response; however, the influence of learning on attack rate and handling time, and the consequence of these changes on the functional response, are poorly known.

According to Holling (1966), type III curves (S‐shaped) are indicative of organisms that show some form of learning behavior. These organisms have developed general responsiveness to many stimuli and can filter out irrelevant stimuli. Likewise, they can learn and separately channel information from different stimuli. These channels are not permanently established since the learned association will disappear unless it is reinforced or undergoes different experiences. The three key features of this behavior are associative learning, information channeling, and forgetting (Holling, 1966). Such features give organisms great flexibility which allows them to focus on a few stimuli and still retain the ability to take advantage of changes in the environment. In dynamic populations, when the prey density is very low, the predator might not associate this stimulus with a reward because the prey is so rare. Conversely, if prey density increases, the predator could become more responsive to the specific stimuli of the prey through learning. Tinbergen (1960) called this behavior the development of a specific searching image.

Holling’s model also reproduces prey switching, where the predator will consume preferentially (or more than proportionally), the most abundant prey. Thus, the predator will “switch” to another prey once the relative abundance of the different prey species reaches a critical threshold, which usually is near the inflection point in the sigmoidal functional response curve. Based on the enzyme kinetics equation, the shape of the curve is mediated by an N parameter which is the number of encounters that a predator must have with its prey before the predator is maximally efficient in consuming that prey. As the N term multiplies, the curve becomes increasingly more like a switch function (Real, 1979).

However, many predators do not have access to alternative prey due to their specificity, like some biological control agents (Byeon et al., 2011), or under laboratory conditions, they are exposed to only one type of prey; however, they exhibit a type III functional response quite frequently (Dunn & Hovel, 2020; Van Lenteren et al., 2016; Yazdani & Keller, 2016). Consequently, a different type of learning should take place, not mediated by the presence of alternative prey, but by the accumulated experience of the organism when searching, manipulating, and consuming prey. As in any learning process, the organism should then exhibit a learning curve (Shaw & Alley, 1985), in which the accumulated experience would translate into a modification of the functional response as a result of the experience.

Learning has been found extensively in almost all animal taxa (Manning & Dawkins, 1998; Shettleworth, 2001). These phenomena have long been described directly in parasitoid or predatory insects (Haverkamp & Smid, 2020; Little et al., 2019; Turlings et al., 1993; Vet et al., 1995). However, very few authors have studied how learning alters the parameters of the functional response curves, for example, Mendes et al. (2018) have found that, in egg predatory mites, the experienced females have significantly smaller manipulation times when compared to naive ones, but their attack rate is the same. Other authors investigated how pesticides affect predator efficiency, either because the predator attacks less prey or because of a decreasing searching time (He et al., 2012; Martinou & Stavrinides, 2015). These results show that learning can occur separately in attack rate or handling time, which makes the enzyme kinetic approach not fully compatible with the results of laboratory experiments or monophagous insects. Therefore, to explore the relationship between learning and functional response, an alternative model is necessary.

To test the hypothesis that the learning in predators produces the appearance of a type‐III functional response, in this study, we propose to explicitly incorporate learning curves in the parameters of the Holling’s disk equation (functional response type II) and to analyze what changes are produced in the functional response shape by applied learning on its fundamental parameters, the attack rate, and the handling time.

METHODS

Model

As a starting point, we used Holling’s disk equation of functional response type II: where D is the prey density, a the attack rate, h the handling time, T the length of the experiment, and N the consumed prey. The total consumed preys after a certain amount of time was called N.

An exponential learning curve, based on Thurstone (1919) model and Estes (1950) statistical theory of learning in a and h, was added to the functional response type II model. Exponential is a very simple model in which the learning is proportional to the time (T) that take to do a given task b, in a differential equation form is dT/dN = −lT, as a consequence, according to Newell and Rosenbloom (1981), the learning occurs at a constant rate l as a function of the experience (here the number of consumed preys, N). By including this model, both parameters are allowed to improve as a function of preys attacked, resulting in a monotonic increase and decrease of a and h, respectively as shown in Figure 1.

FIGURE 1

Theoretical learning curves proposed for the handling time (h) and the attack rate (a), both curves are of exponential learning type, where they tend exponentially to a final asymptotic value from an initial point. h 0 and a 0 are respectively the initial values of the handling times and the attack rate, while h and a are the final values of both variables. The learning amplitudes (Δ, Δ) are defined by the distances between the initial and final values of h or a

For a, the model is a monotonically increasing function is: where l is the learning rate of a per attacked prey, and a is the maximum possible attack rate for this species, with 0 < a ≤ 1. Here the attack rate a increases asymptotically from the initial a at the beginning of the experiment (a 0) to a at a rate of l. If l = 0, there is no learning, but if l = 1, the learning is maximum, and a is approached after a single prey is consumed. How much it is possible to learn is Δ = a–a 0.

A similar model was proposed for the handling time h, with the difference that in the case of h, it decreases with experience and asymptotically tends to zero instead of one as in the case of a, so where l is the learning rate of h per attacked prey, and h is the minimum handling time for this species. The handling time h tends asymptotically from h 0 (the handling time of the inexperienced predator) to h at a rate of l. How much it is possible to learn is Δ = h 0–h. Both Δ (Δ, Δ) are called learning amplitude.

Finally, if l, l, Δ, or Δ are equal to zero there is no learning and the functional response function becomes Holling’s type II disk equation. The resulting dynamic of this model is shown in Figure 1, both a and h depend on the initial condition (a 0, h 0), the asymptotic values (a, h), the learning rate (l, l), and how much the predator can learn (Δ, Δ).

Prey depletion

The effect of prey depletion was also tested using the Rogers (1972) approach, in which the preys are a fixed pool in the experimental arena and are removed without reposition, so the differential equation is modified as: where N 0 is the initial number of prey available, and N is the consumed prey, the available prey (N) at time T is N = N 0–N. As a consequence, the available prey is constantly removed and its number decreases asymptotically toward zero.

Analysis

Under a learning context, the function would be expected to be convex at low prey densities, because the improvement in the predator’s ability to consume prey as a result of its experience is greater than its limitation in the ability to consume prey at a handling time greater than zero. At the inflection point, the improvement in the ability of the predator to consume prey as a result of learning is exactly compensated at that point by the limitation in the ability to consume prey, so the function becomes purely limited by the manipulation time.

Therefore, a characteristic that allows identifying the type III functional response is its sigmoid convex‐to‐concave shape. There, the function of consumed preys (N) as a function of density (D) such as in any sigmoid function has an inflection point, so it is the second derivative zero. Before the inflection point, the function is convex (positive second derivative), and after this point, it is concave (negative second derivative). In type II functional response, the second derivative is always negative and asymptotically approaches zero, and on type I response, it is always zero, as it is a linear function (Figure 2). Therefore, the shape of the functional response was defined as the function of the second derivative of the functional response as a function of the offered preys. As explained above, the cases in which there was an inflection point in which the second derivative was zero while decreasing from positive to negative values were classified as type III, if it was always negative, they were classified as type II, and if it was zero or near‐zero, as type I. Additionally, if the slope of the second derivative was too low (near‐zero), and it was barely noticeable, they were considered near‐type II, and if the slope was low and always near zero, they were considered near‐type I. The three functional response functions with their corresponding first and second derivatives are shown in Figure 2.

FIGURE 2

Theoretical functional response curves according to Holling’s (1959) paper, and its first and second derivatives, with the diagnostic characteristics of each type. Preys attacked as a function of preys offered for Holling’s three functional responses (a), (b) is the first derivative of the function, and below (c), the second derivative. The type III functional response can be characterized by the presence of an inflexion point in the slope of the curve, while the other two do not contain critical points of any kind

To test whether learning in terms of improvement of attack rate and manipulation times generates type III functional responses, the analysis was aimed at identifying the different shapes of the functional response curve in different conditions of learning, here identified as parameters l and l. Other parameters affecting the behavior of the functional response curve were a 0, a, h 0, and h. To simplify the analysis, we assumed that a = 1 and h 0 = 1, so there are only four parameters to analyze, l, a 0, l, and h.

To avoid the curse of dimensionality, the parameters were tested in pairs, with the ones influencing the attack rate being analyzed separately from those influencing manipulation time. So one analysis was performed manipulating a 0 and l, keeping h constant (l = 0, h = h 0) and another with h and l keeping a constant (l = 0, a = a 0 ).

Time also influences the results, as the experiment became longer, the number of prey consumed and therefore the experience increased, so to keep the dimensions as low as possible, the length of the experiment was fixed to one time unit. As a consequence, T is removed from the equations as it is equal to one in a multiplication. Substituting the a and l parameters from the functional response Equation 1 by the Equations 2 and 3 will result in an equation where the N variable is present on both sides, as it happens with the Rogers random predator model from Equation 4, so we decided instead to solve numerically the models from Equations (1), (2), (3), (4) in the form of differential equations, as in Rosenbaum and Rall (2018) until the time equal one to find the results of this study. Equation 1, then became

Equation 2 expressed as a differential equation is:

Equation 3 also expressed as a differential equation is:

The Rogers random equation is simply the Holling’s disk equation with an exponentially decaying N substituted into it, so in its differential equation form it became: and where as in Equation 4, N is the available prey, so N = N 0 – N.

The final form of the differential equations used for this study was for the parameters a 0 and l, the composition of Equations 1 and 2 into a system of ordinary differential equations: so, a is updated as a function of time using the number of prey consumed. In the prey‐limitation case, Equation 4 replaces Equation 1, and by replacing dN with −dN, the sign of the equations is changed, as a result, the system changes to:

Similarly, the effect of h and l parameters were tested by composing Equations 1 and 3 into an ODE system:

And the effect of prey‐depletion became:

Numerical methods

The tests were simulations at 1,002,001 combinations of parameters (a 1001 × 1001 matrix). For the attack rate a, the values used in the simulations were a range of a 0 between 0.01 and 1.0 at intervals of 0.01, and with a range of l between 0 (no learning) and 1 (maximum learning) at intervals also of 0.01. For the manipulation time h, the range of h values was between 0.1 and 0 at intervals of 0.01, always starting from 1, and a range of learning rates l as in l, with a range between 0 and 1, also at intervals of 0.01. All the analyses were performed without and with prey limitation (Holling and Rogers models respectively). Each simulation was run to a fixed length of one time unit, and the prey densities were between 1 and 200 at intervals of 0.01. When the simulations were performed iterating over l and a 0, l was considered zero (no learning in h), and h was fixed to 0.1; on the other hand, when the iterations were over l and h, l was considered zero (no learning in a), and a was fixed to 0.1.

For each combination of parameters, the differential equations were solved by integrating numerically the differential equations using the fourth‐order Runge‐Kutta (RK4) integration method (Süli & Mayers, 2003), at steps of 1/1024 time units. All the simulations were performed with a series of routines written in the Python programming language which are available in Appendix S1.

RESULTS

Learning attack rate

The simulations showed that under learning that improves attack rates, the functional response was predominantly of type II, especially at low rates of learning (Figure 3), and at high levels of a 0 (initial attack rate). On the other hand, at low levels of a 0, the response began to transition from type II to type III responses, as the learning rate increased.

FIGURE 3

Functional response as a function of the combination between the logarithm of the learning rate (l) on the Y‐axis, and the initial value of the attack rate (a 0) on the X‐axis, without limitation by prey (a, c) and with limitation (b, d). As explained in the methods, the remaining parameters were l = 0.0 (no learning), h = 0.1, length of the experiment = 1 time unit, and the range of prey densities used to calculate the derivatives were between 1 and 200 at steps of 0.01. In graphs a and b, the grayscales and contour lines show the prey density (N) at which the first inflexion occurs. Graphs c and d, on the other hand, show the maximum value of the second derivative of the functional response curve. Maximum values of second derivative greater than zero indicate type III functional responses, values close to zero, but positive, are near‐type II functional responses, finally, negative values indicate type II functional responses

In an intermediate zone, the functional response was characterized as near‐type II, this response has characteristics of a type III response, such as a positive second derivative or close to zero, but very attenuated, so it is visually indistinguishable from type II (Figure 4).

FIGURE 4

Consumed preys as a function of offered preys under different learning (l) rates for attack rate (a, d), and first (b, e) and second derivatives (c, f). Plots a, b, and c are the models without prey depletion, so the X‐axis is the prey density, while d, e, and f are with depletion according to the Rogers model, therefore the X‐axis is the initial number of prey. With high learning rates, the functional response approaches type III functional response (positive second derivative at low prey density as a consequence of learning, and then negative as a consequence of handling time limitation). As explained in the methods, the remaining parameters were l = 0.0 (no learning), h = 0.1, length of the experiment = 1 time unit

At low initial attack rates (<0.3) and learning rates >0.1, the functional response is a Holling type III. Interestingly, when the learning rate becomes very high (>0.5), the near type II functional response again requires lower values of a 0, since its inflection point becomes very close to zero in very low prey densities, so the curve becomes closer to type II at similar values of a 0, but with lower learning values.

Under the Rogers model with prey depletion, the results are similar to those with constant density, the major difference is that the area in the parameter space in which the functional response is type III and/or close to II, is slightly larger, reflecting the effect of the reduction in the number of preys available in the shape of the functional response. On the other hand, the maximum second derivative under the Rogers model is lower than with the absence of prey limitation.

Learning handling time

Simulations carried out with manipulation time learning showed results where most of the combinations of parameters produce type II functional responses. Only in combinations of parameters with very high learning rates and very low minimum manipulation times the functional response begins to differ from type II to resemble type I (Figures 5 and 6), within the upper left corner of the said graph (l = 1, and h = 0), the functional response is type I.

FIGURE 5

Functional response as a function of the combination between the logarithm of the learning rate (l) on the Y‐axis, and the minimum handling time (h) on the X‐axis without prey depletion (a, c) and with prey depletion (b, d). As explained in the methods, the remaining parameters were l = 0.0 (no learning), a = 0.1, length of the experiment = 1 time unit, and the range of prey densities used to calculate the derivatives were between 1 and 200 at steps of 0.01. In graphs a and b, the grayscales and contour lines show the prey density (N) at which the first inflexion point occurs. Graphs c and d, on the other hand, show the maximum value of the second derivative of the functional response curve. Maximum second derivative values greater than zero indicate type I functional responses, values close to zero, but negative, are near‐type II functional responses, finally, negative values indicate type II functional responses

FIGURE 6

Consumed preys as a function of offered preys under different learning (l) rates for handling time (a, d), and first (b, e) and second derivatives (c, f). Plots a, b, and c are the models without prey depletion, so the X‐axis is the prey density, while d, e, and f are with depletion according to the Rogers model, therefore the X‐axis is the initial number of preys. With high learning rates, the functional response approaches type I functional response (positive second derivative at low prey density as a consequence of learning, and then negative as a consequence of handling time limitation). As explained in the methods, the remaining parameters were l = 0.0 (no learning), a = 0.1, length of the experiment = 1 time unit

As shown in Figure 6, the learning of the manipulation time never generated visually similar functional responses to type III; on the contrary, as the learning rate l improved, the functional response looked more and more like type I. When the minimum handling time (h) was very short, the curve equaled type I, otherwise, the curve resembled an intermediate between a truncated type I and type II. At low prey densities, the second derivative is always negative, only with very high learning rates (l > 0.5), the second derivative becomes slightly positive at higher prey densities, to become slightly negative again. This does not conform to a type III functional response, but rather a near‐type I, because the derivative is very small, and the two critical points are unnoticeable. The effect of prey depletion was again small, with an area of type I functional response smaller than with Holling’s model.

DISCUSSION

Learning can produce all types of functional responses, depending on what parameters it affects. Under conditions of high learning amplitude, the functional response differed from type II. At low learning rates in both a and h, the result is a type II functional response, at high levels of learning in a prey‐predator is type III, while at high learning rates in h, the resulting functional response approaches type I, so if h = 0, as the experiences accumulate, the handling time approaches asymptotically to zero, and the functional response approaches dn/dt = aN, which is the functional response type I. Learning can only produce type III functional responses if it affects the attack rate.

Theoretically, a linear functional response is possible when a predator can search and handle different prey simultaneously, or when the handling time is negligibly small (Hassell, 2000; Jeschke et al., 2002, 2004). A consequence of this work is that the type I functional response, which is usually interpreted as typical of filter feeders (Jeschke et al., 2004), or in general is not associated with learning processes, can be also a result of complex behaviors. Examples of type I functional responses have also been found in some parasitoid species (Kaçar et al., 2017; Mills & Lacan, 2004), and in filter‐feeding birds (Arzel et al., 2007), a taxon of animals capable of learning complex behaviors. Arzel et al. (2007) found a switch point between two types I functional responses with different slopes, showing that the complexity of the foraging behavior might imply several tasks that once optimized cannot be improved further and result in a curve with a series of discontinuities. Here we propose that this type of functional response is also a product of learning like type III, only that it is a different type since it has the same shape, but is caused by a different mechanism. So, it can be called type Il if it is generated by learning, and type If it is generated by filter feeding. On the other hand, learning in terms of improvement in the attack rate would produce responses increasingly similar to III, while if the learning occurs in the optimization of handling times, a type Il response would be generated.

Prey depletion did not change the overall pattern of the results. The main difference is that the maximum second derivative under the Rogers model is lower than with the absence of prey limitation, reflecting that as the available preys are reduced with time, so do the opportunities to learn, as a consequence the second derivative is smaller. Therefore, in an arena with a limited number of preys, the animals learn less, because they run out of preys. Here, the appearance of a turning point is earlier due to depletion of the prey and not due to learning. So, the functional response is more often type III, but not because of learning. The effect of handling time was similar with an area in which the functional response is type I or near‐type I, smaller than without prey depletion. The main reason is that as the number of available prey decreases with time, the functional response begins to be limited earlier, and keeps its functional response II shape even though its handling time decreases to a near‐zero value.

As proposed here, only the predator exhibits learning abilities, but prey may also learn, decreasing the effectiveness of predators (Brown, 2003; Zhao et al., 2006). Prey learning by personal experience is also costly because it involves surviving an encounter with the predator, so prey‐learning might resort also to social learning (Galef & Laland, 2005). Aggressive defensive behavior of preys might lead to complex patterns, including type‐IV functional response, in which the number of consumed preys decreases after reaching a maximum, product of prey collective defense actions (Líznarová & Pekár, 2013). To the best of our knowledge, no theoretical model is describing the effect of prey learning on the functional response, and the interactions between it and the predator abilities, so it remains an open topic in predator‐prey theory.

In terms of population dynamics, it is generally accepted that the type III functional response may have stabilizing effects on prey dynamics as proposed by Oaten and Murdoch (1975), the reason is the increase in the probability for a prey to be killed as their density increases, which means the presence of a positive second derivative on the functional response curve, as in the type III curve. Here we observed that this phenomenon (the transition from a type II response to a type III) occurs only in the case of strong learning in attack rate, but not in learning in handling time. However, the handling time learning produces a curve with an asymptote that increases with time and prey density, which can give a different type of stabilization that requires further more specific studies.

Another issue is that, since learning is an accumulated process, it will interact with population dynamics in the form of a delayed effectiveness response. For example, given that the prey population will decrease after a peak or outbreak as some herbivorous insects as described by Berryman et al. (1987). Under the approach used here, with cumulative learning, some predators might remain more effective for a while after a decrease in prey population. In the case of parasitoids, as the lifespan is the same as its hosts it is not an issue, while in the case of vertebrate insectivores it will cause further instability or local extinctions, by consuming more than proportionally when the population is small after an outbreak.

CONCLUSIONS

The results obtained in this study show that learning can change the functional response of predators in different ways since it is generally accepted in the literature by generating either type I, II, or III and intermediate forms in the absence of alternative prey. Therefore, learning can both be a stabilizing or destabilizing factor in the population dynamics, depending on which type of prey consuming behavior it affects.

CONFLICT OF INTEREST

The authors declare no competing interests.

AUTHOR CONTRIBUTIONS

Octavio Augusto Bruzzone: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). María Belén Aguirre: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Jorge Guillermo Hill: Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Eduardo Gabriel Virla: Conceptualization (equal); Supervision (lead); Writing – review & editing (lead). Guillermo Logarzo: Conceptualization (equal); Funding acquisition (lead); Resources (equal); Supervision (lead); Writing – review & editing (lead).

Appendix S1

Click here for additional data file.