Geographic variation in body size and plumage colour according to diet composition in a nocturnal raptor

selection, we showed that predator–prey interactions can affect spatial phenotypic variation by promoting the evolution of local adaptations.


Introduction
Interactions between prey and predators are amongst the major selective forces that drive a diverse suite of functional phenotypes in both counterparts (Kerfoot andSih 1987, Tollrian andHarvell 1999). Feeding relationships can indeed shape elaborate behavioural, physiological and morphological adaptations to prevail in the evolutionary arms race (Edmunds 1974, Surlykke and Miller 1985, Reimchen 1980, West et al. 1991, Brodie III and Brodie Jr et al. 1999, as well as have a prominent and pervasive impact on the stability of food webs, on the structure of animal communities 2 and on the dynamics of natural populations (Emmerson and Raffaelli 2004, Brose et al. 2006b, Petchey and Dunne 2012.
One of the most recognized patterns of predator-prey interactions, which is consistent across aquatic and terrestrial biomes, is that they are intrinsically size-structured (Emerson et al. 1994, Brose et al. 2006a, Hatton et al. 2015. Indeed, it is generally well acknowledged that a predator's size correlates with its prey size (Vézina 1985, Radloff and Du Toit 2004, Owen-Smith and Mills 2008, Costa 2009). According to the optimal foraging theory (Stephens et al. 2007), large animals should select food items which are generally larger than those usually eaten by small species. However, it is also recognized that predators of limited size are usually forced to hunt small prey items (Brown andMaurer 1989, Barclay andBrigham 1991), partially because their gape is smaller comparing to larger animals (Shine and Sun 2003), while large ones are able to detect, capture and consume prey of more variable sizes (Brown andMaurer 1989, Scharf et al. 2000). In addition, larger species usually live in larger home and geographic ranges than smaller ones, thus potentially encountering a wider variety of prey. According to these premises, a positive association between a predator's body size and its diet niche breadth is expected (Gittleman 1985, Scharf et al. 2000. Such an allometric scaling in the relationships between predator size and prey size and diversity is a crucial aspect of food web and metabolic theories (Emmerson and Raffaelli 2004, Brose et al. 2006b, Otto et al. 2007, Petchey and Dunne 2012, Kalinkat et al. 2013. It is important also to note that the strength of predatorprey interactions, however, is not constant across time and space. As a consequence, the magnitude of selection on phenotypic traits of predators and prey is expected to fluctuate spatially and temporally, thus promoting the local emergence of functional phenotypes (Van Buskirk 2002, Kishida et al. 2007). Associations between predators and prey phenotypes, especially size, have been widely studied, but most of the current knowledge is based on comparative studies among species (Vézina 1985, Radloff and Du Toit 2004, Owen-Smith and Mills 2008, Costa 2009, but see Török 1993 or withinpopulation analyses (Simpfendorfer et al. 2001, Ingram et al. 2011. However, intraspecific spatial covariation in phenotypic traits between predators and prey has been rarely investigated (Schwaner 1985, Erlinge 1987, Forsman 1991. This is especially the case for species living across wide geographical ranges. This is unfortunate because its examination is fundamental to understanding how local adaptations emerge and potentially drive the evolution of spatial polymorphism (Huey et al. 2000).
To partly fill this gap, in the present study we examined the covariation between prey composition and some phenotypic traits of a widely-distributed nocturnal predator, the western barn owl Tyto alba. This species typically hunts small mammals across its entire range of distribution (Romano et al. 2020b), which is comprised between the British Isles and south Scandinavia to South Africa, including Middle East (easternmost point in Iran), Arabian Peninsula, Madagascar, most of the Mediterranean islands and most of the African archipelagos in the Atlantic and Indian Oceans (Romano et al. 2019(Romano et al. , 2020a for details). In particular, we investigated whether wing length, a good proxy of body size, bill length and plumage colouration vary geographically according to some diet features, such as percentage of mammals, prey size and diversity. In this taxon, body size (i.e. wing length) does not vary along geographic or climatic gradients according to Bergmann's rule (Romano et al. 2021a), thus suggesting that other selective pressures may drive its spatial variation. Conversely, bill length has been shown to vary consistently with Allen's rule, with populations living at higher latitudes and altitudes (and smaller temperatures) showing a smaller beak than those inhabiting closer to the equator (Romano et al. 2020a). However, considering the obvious importance of this trait in consuming prey, it might be also locally affected by diet. Finally, the degree of brownish colouration has been shown to be not only associated to climatic conditions (Romano et al. 2019), but also to foraging behaviour within population (Roulin 2004, Charter et al. 2012, San-Jose et al. 2019, therefore leaving open the possibility that it is potentially affected by local diet.
According to the previous literature, we predicted that wing length and bill length should be positively associated with prey size. Indeed, larger individuals with a larger bill should be better able to hunt and process larger prey. Because small predators are often restricted in hunting small prey, we also expected that wing length and bill length should be positively associated with prey diversity. Finally, given the lack of scientific literature on these topics, we refrain from proposing specific predictions for the other relationships between phenotypic traits and diet features.

Diet data collection
Diet data were collected using published information extracted from scientific papers, grey literature and PhD/ master theses reporting description of pellets content (details in Romano et al. 2020b). We then selected the papers/reports collecting at least 90 vertebrate prey (Marti et al. 2007), and in which at least of 80% of mammalian prey were identified at the genus level. Following Romano et al. (2020b), we also included diet data of a small island in Cape Verde archipelago (Fogo island, Siverio et al. 2008) with a smaller sample size (50 identified prey). However, the results are qualitatively very similar excluding this datum. We retained papers including information collected in a single site and in geographical small regions (e.g. island, county, district), as well as studies including single-year and multiple-year data. Following Romano et al. (2020b), in order to reduce the temporal variability in diet composition, analyses were limited to data collected after 1940.
For every location, we reported the geographical coordinates, and the proportion of each terrestrial mammalian prey genus i (number of individuals of genus i divided by the total 3 number of terrestrial mammals contained in the pellets) was recorded. When reported, the proportion of the terrestrial mammalian prey over the total amount of vertebrate prey was also computed. Because many studies focused on non-volant mammals only and considering that they represent a minor component of the barn owl's diet (Roulin and Christe 2013), bats were included among the other vertebrates, with birds, reptiles and amphibians. Information on invertebrates was discarded because a very small part of the papers reported reliable information about this food source and because it contributes to a minimal part to the diet (Obuch and Benda 2009, Muñoz-Pedreros et al. 2016, Romano et al. 2020b for details).
Because many studies did not include information at the species level and in order to account for variability in the species diversity among genera, Shannon diversity index of prey was calculated at the genus level (i.e. using the proportion of each mammalian prey genus). The Shannon diversity index was calculated only on mammals because they are the major component of barn owl diet and because in many cases the other vertebrates were generically reported as birds, reptiles or amphibians (Romano et al. 2020b for details). We used Shannon diversity index because it is a good proxy for summarizing diet diversity, irrespectively of prey taxonomy, as reported by recent studies in raptors, including the species under investigation (Milana et al. 2018, Janžekovič and Klenovšek 2020, Romano et al. 2020b). Finally, average body mass of each mammalian prey genus was also estimated. Then, prey were coded into four categories depending on their size (small = maximum 25 g, medium = between 26 and 50 g, large = between 51 and 100 g, very large = more than 100 g). Further detailed information are reported in Romano et al. (2020b).

Collection of phenotypic data
Information on wing length, bill length and degree of reddish colouration of the breast was collected on thousands of T. alba specimen across Europe, Middle East and Africa, which were deposited in museums or collected by private citizens (Romano et al. 2019, 2020a for details). All measurements were made by the same experimenter (A. Roulin). Wing and bill length were measured to the nearest 1 mm and 0.1 mm, respectively, and their measurements are highly repeatable (Romano et al. 2020a). Importantly, we used wing length as a predictor of body size because information on other traits associated with body size (e.g. body mass) was not available on museum specimens and because the vast majority of the studies investigating body size variation in birds used this trait as a surrogate for size (Ashton 2002). Also, this trait was used for the same purpose in previous analyses of the species under investigation (Romano et al. 2021a, Romano et al. in press). The plumage colour of the underside region of the body, which varies from white to reddish-brown both among and within populations (Supporting information; Roulin 2003, Romano et al. 2019, was categorized on an 8-level scale, ranging from 1 for dark brownish to 8 for white. Hereafter, we refer to dark/darker and pale/paler individuals according to this white-to-brown categorization of plumage colour. This categorization is a reliable information on the feather colour because it strictly mirrors the reflectance in the range of visible light as measured by a spectrometer (Dreiss and Roulin 2010). For all the specimens, we also collected geographical coordinates of the location where they were collected.

Association between diet and phenotypic data
We then combined the above-described datasets with the following procedure. Firstly, we divided the distribution range of T. alba into cells of variable side length, from 100 to 400 km. Within each cell we pooled all the diet information provided by different studies. For example, when in a 200 by 200 km cell there were three sites where diet information was available, they were pooled in order to combine all the information about diet data in the same region. We simply averaged the values of different prey categories per site within each cell (e.g. mean proportion of small prey, mean Shannon diversity index, mean proportion of mammals). However, considering that different studies collected information on different sample size (range: 50-233 540 prey items per site, median: 904 prey items per site), with a second procedure, we summed up all the data from different sites (i.e. sum of small mammals, sum of vertebrates, sum of mammal prey, etc…) and then computed a global value for every diet parameter within each cell. This procedure allowed us to weigh more diet data collected in sites where a larger sampling effort was done. Analyses were performed using data extracted with both procedures, but in the Results section we only show the results obtained using the mean values of each diet parameters, while the others are reported in the Supporting information.
To each cell we also associated the barn owl specimens whose coordinates of the recovery site were contained within the cell. However, when diet information was available for islands smaller than a 200 by 200 km cell, we associated to each datum only phenotypic data of the barn owl inhabiting that given island. Again, we averaged values of phenotypic traits among specimens within each cell. Considering that the barn owl shows moderate phenotypic variability within populations, we only considered cells that had data on at least four specimens to increase the accuracy on phenotypic information.
We repeated this procedure for cells of different size (100, 200 and 400 km side length). After this step, the sample size of locations was adequate to perform analyses for cells of 100 and 200 km side length. For the analyses, we considered cells of 200 km side length as the best compromise between a large geographical coverage of the diet data, a large sample size and an accurate estimate of barn owl's phenotype. This is why analyses using such a spatial scale were reported in the main text. However, analyses were repeated on cells of 100 km to check whether results were consistent using different scales. Unfortunately, we could not do the same for cells of 400 km side length because of a considerable reduction of the sample size (a reduction of ca 30% of cells) and because the geographic coverage of data was extremely biased towards Europe.
On the whole, our analyses included 116 cells containing on average 17.27 barn owls (range: 4-229, Supporting information). However, sample sizes of different phenotypic traits slightly vary because of missing data (e.g. bill was broken, wingless individuals, etc…). The total number of barn owl specimens in the analyses is 3100 and the total mammal prey items identified in the locations included in the analyses is 723 556.

Statistical analyses
To examine variation in phenotypic traits according to diet parameters, we used linear models using the 'glmmTMB' package (Brooks et al. 2017) in R (ver. 3.5.1). Variation in wing length, bill length and plumage colour was analysed separately (Romano et al. 2019(Romano et al. , 2020a, including Shannon diversity index, and proportion of small prey (up to 25 g) as predictors. We focused on proportion of small prey rather than on average size of the prey for two important reasons. First, prey size was estimated at the genus level and this information is therefore not precise enough to properly obtain an accurate estimate of average prey size. This is not the case for proportion of small prey because each genus considered in the analyses is included in a single size category. Second, average prey size can be considerably affected by the presence of a small number of prey items of extreme size, while this is not the case for proportion of small prey.
We note that we chose a different cut-off of prey size (25 g rather than 50 g) compared to the one used in a previous study that focused on different barn owl lineages living in other continents (Romano et al. 2020b) for a specific reason: prey size is generally much smaller in the regions inhabited by T. alba than those inhabited by the other barn owl lineages (i.e. T. furcata and T. javanica), and prey smaller than 50 g comprised more than 90% of the prey almost everywhere in its distribution range (Romano et al. 2020b). In practice, T. alba almost entirely relies on prey smaller than 50 g, and using such a cut-off the distribution of this predictor is highly biased toward higher values. This is not the case for a 25 g cut-off, which results in a normal distribution of this predictor. In the present study, we could not include analyses on the other barn owl lineages because of small sample size.
In addition, we repeated the analyses also including proportion of mammals in the diet. However, because proportion of mammals was never a significant predictor (Results) and considering that this information was not available for all the cells (i.e. some diet studies were focused on mammals only), thus resulting in a reduction of sample size, the final analyses were repeated excluding this term. Although within single population, darker owls seem to consume more cricetids and paler more murids (Roulin 2004, Charter et al. 2012, an analysis at such a low taxonomic level (i.e. the family level) were prevented because of the large diversity of small mammal assemblages across the distribution range of the barn owl. For example, no cricetids have been reported in the diet of African barn owls (Romano et al. 2020b), thus preventing us to perform proper analyses on this.
Since both phenotype and diet vary between islands and mainland (Romano et al. 2019, 2020a, b, Romano et al. in press, Janžekovič and Klenovšek 2020, in the models we also included a dichotomic factor indicating if the cell under investigation was on an island (coded as 1) or a mainland (coded as 0). However, analyses performed on continental populations only provided very similar results to those shown in the Result section (details not shown).
Diet and phenotypic traits are not randomly distributed across space (Romano et al. 2019(Romano et al. , 2020a. In order to account for non-random distribution of locations, in all the models we also accounted for spatial autocorrelation, by adding an exponential correlation structure considering the distances between all the pairs of latitude-longitude coordinates of the centre of each cell. The phenotypic traits under investigation have been shown to be predicted by climatic factors, such as temperature and precipitation (Romano et al. 2019(Romano et al. , 2020a, potentially affecting our results. The models described above were therefore repeated with the inclusion of mean annual temperature and annual rainfall of the centroid of each cell, extracted from the WordClim data repository (<http://www.worldclim.org/>, Fick and Hijmans 2017) for the period 1970-2000 (for details, Romano et al. 2019Romano et al. , 2020a, as additional predictors. In order to check whether significant diet predictors were included in all the best supported models (in terms of Akaike information criterion; AIC hereafter), we also performed a model selection analysis for each dependent variable, using function dredge of the R package 'MuMIn' (Barton 2009).
For all the models, we checked for collinearity (i.e. variance inflation factor), residual diagnostics, normal distribution and homoscedasticity, by using the 'performance' package (Lüdecke et al. 2020). No collinearity (variance inflation factor was always smaller and 2.24) and heteroscedasticity was detected, and the distribution of the model residuals was normal. We finally checked for the presence of outliers (Grubb's test) which were removed from the final models (this was the case for two datapoints in the analyses of wing length).

Results
Wing length is negatively predicted by the proportion of prey smaller than 26 g (Table 1, Fig. 1a). This result implies that body size is larger when the diet is composed mainly of large prey. In addition, a similar pattern was observed when in the same model the proportion of small prey was replaced by mean prey size (slope ± SE = 0.146 ± 0.045; t = 3.24; p = 0.012). No effect of diet diversity (Table 1) and proportion of mammals (slope ± SE = −2.924 ± 3.176; t = −0.92; p = 0.36) was observed on wing length.
Bill length was not affected by any of the predictors (Table  1). This was also the case when wing length was added to the models as an additional predictor in order to test variation in relative bill size (details not shown).
Variation in plumage colour is also significantly affected by the proportion of small prey (Table 1). Indeed, plumage colour becomes darker at increasing proportion of small prey (Fig. 1b). However, mean prey size did not predict variation in this plumage trait (slope ± SE = 0.013 ± 0.009; t = 1.50; p = 0.13). No effect of diet diversity (Table 1) and proportion of mammals (slope ± SE = −0.059 ± 0.549; t = −0.11; p = 0.91) was shown on plumage colour.
Very similar results were obtained when diet parameters were collected in cells of 100 km side, but the signs of the models of bill length were opposite (Supporting information), and when they were summed up from different sites within cell rather than averaged (Supporting information). In the latter analysis, however, wing length also increases with Shannon index. The same statistically significant relationships were obtained when climatic variables were included in the models as additional predictors (Supporting information), thus indicating that the above results are not determined by a covariation between diet, phenotype and climate. This result was confirmed by model selection analyses because all the best supported models in terms of AIC always included proportion of small prey among predictors (Supporting information). The other predictors that were always present in the best supported models of wing length and plumage colour were annual precipitation and island, respectively (Supporting information).
Finally, significant trends in both wing length and plumage colour according to percentage of small prey were found also when the analyses were restricted to prey which are rodents (wing length: slope ± SE = −5.149 ± 2.042; t = −2.52; p = 0.012; plumage colour: slope ± SE = −1.189 ± 0.387; t = −3.07; p = 0.002). Therefore, these patterns are not due to the presence of insectivores in the diet (e.g. Soricidae and Afrosoricidae), which represent a non-negligible fraction of the diet at temperate and boreal regions (Romano et al. 2020b), and which are generally smaller than rodents (Cotgreave and Stockley, 1994).

Discussion
One of the main results found in the present study on a very large geographic scale was that populations of the western barn owl that had longer wings hunted larger prey compared to those showing on averagely smaller wings. This finding is in line with most of interspecific comparative studies (Vézina 1985, Radloff and Du Toit 2004, Owen-Smith and Mills 2008, Costa 2009, but see Török 1993) and within-population analyses (Simpfendorfer et al. 2001, Ingram et al. 2011 showing a predator-prey size allometric scaling, and suggests that large individuals are able to exploit a larger fraction of the encountered prey, including large ones. This pattern can be therefore interpreted as an adaptation of barn owl body size for an efficient use of the available sources of food. Indeed, the net energy gain per prey item consumed is expected to increase with increasing prey size, although larger prey are probably more costly to consume and/or to hunt.
Given the correlative nature of the data, however, we are not in the position of arguing causality in the observed covariation between prey and predator size. Indeed, on the one hand, populations composed of generally large individuals might be specialized in foraging on relatively large prey, and, vice versa, those consisting mainly of small individuals might be constrained to preferentially hunt smaller mammals, possibly because their gape size is small (Shine and Sun 2003). On the other hand, the presence of different prey assemblages might have driven the evolution of large-bodied barn owls in regions where prey are generally larger by exerting a negative selection on smaller individuals and/or because larger prey confer more nutrients thus potentially promoting a larger body size growth. In addition, we cannot exclude the possibility that the observed association would be on wing length per se, and not via an effect on body size. Indeed, wing shape and length are known to considerably affect flight ability and foraging strategies (Gamauf et al. 1998, Corvidae et al. 2006. For example, it is possible that long-winged birds are more specialized in hunting large prey, which should generally cover longer distances in a unit of time. Body size, wing morphology and hunting strategy may thus have coevolved in different populations in order to maximize success in hunting the most available prey. However, irrespectively of the direction of causality and the phenotypic target of selection, our study is the first, to the best of our knowledge, to show an intraspecific spatial variation in body size of a predatory bird according to variation in diet size composition. It adds to the limited known cases of adaptations of predator body size according to the size distribution of prey, which was previously observed in other vertebrate taxa, such as mammals (Erlinge 1987) and snakes (Schwaner 1985, Forsman 1991.   Table 1, and other predictors were fixed as their mean values. This result also helps to understand and interpret better our recent finding about the lack of any clear variation in wing length along geographic and climatic gradients in T. alba, a contrasting result when compared to other species of the same genus (Romano et al. 2021a). Wing length seems therefore more strongly affected by diet composition, and thus on food availability, than climatic factors, at least in this taxon.
Contrary to our expectations, we cannot find any significant trend in (relative and absolute) bill length according to diet parameters, even if the direction of the relationship between bill length and prey size is negative. Although such a result seems surprising, this is not if we consider the hunting behaviour of the barn owl. Indeed, this predator attacks prey using its claws, rather than the bill, and then consumes them after reducing the prey into small pieces (Csermely et al. 2002). While bill shape certainly evolved in response to the typical carnivore diet of the barn owl, intraspecific variation in bill size has therefore been mainly driven by thermoregulation reasons (Romano et al. 2020a), rather than on local prey composition.
Another interesting result disclosed by our analyses is that larger prey are present at a larger proportion in the diet of populations showing paler plumage. How to explain this result is matter of speculation. Considering that foraging behaviour but also hunting habitats have been shown to vary according to plumage colour morphs within population (Charter et al. 2012, San-Jose et al. 2019, it is possible that paler-plumaged barn owls might be better specialized in capturing averagely larger prey and/or mainly hunt in habitats where larger prey are more abundant. An alternative scenario rests on the possibility that paler plumage colour might have evolved as a byproduct of selection towards a large body size, which in turn have emerged in response to prey size composition. Indeed, it has been shown that within-population paler individuals are generally heavier than more reddish ones (Roulin 2006). The same relationship between plumage colour and body size (i.e. wing length) also emerges among-population in the sample of barn owls studied here (slope ± SE = 0.048 ± 0.019; t = 2.59; p = 0.010), thus making such a possibility valid. Finally, we note that a possible causal link between diet and plumage melanisation exists. Indeed, quality and quantity of food ingested seem to have direct effects on melanin plumage traits (McGraw et al. 2002, Poston et al. 2005, McGraw 2007, Galván et al. 2019, with individuals consuming a larger content of dietary amino acids and minerals showing more melanic feathers. Following the above reasoning and considering that larger prey contain more nutrients, brownish, rather than paler, feathers should be typical of individuals eating large prey. Our result is therefore in contrast with this previously suggested link between diet and plumage colour. However, we note that in the barn owl plumage colouration is mainly determined by genetic, than environmental, factors (Roulin and Dijkstra 2003), thus making such a link less likely at least in this species.
Finally, contrarily to our expectations, we did not find any effect of diet diversity or taxonomic composition on any of the phenotypic traits under investigation. A possibility is that our Shannon diversity index computed on mammal genera only did not properly mirror the entire range of the prey consumed by different populations. Taken together, these results thus indicate the main driver of local phenotypic variation is the size of the prey, rather than its taxonomic composition.
In conclusion, we found a covariation between prey size and predator size and melanisation across the distribution range of in a widely-distributed bird predator. Irrespectively of the direction of causality and the phenotypic target of selection, the present study shows that predator-prey interactions are important drivers of spatial phenotypic variation by promoting the evolution of local adaptations. Future studies combining information on prey size and composition at different spatial and temporal scales are needed to test for the generality of the patterns documented here for the barn owl.

Transparent Peer Review
The peer review history for this article is available at https:// publons.com/publon/10.1111/jav.02716