Multi-Level Determinants of Parasitic Fly Infection in Forest Passerines

The study of myiasis is important because they may cause problems to the livestock industry, public health, or wildlife conservation. The ecology of parasitic dipterans that cause myiasis is singular, as they actively seek their hosts over relatively long distances. However, studies that address the determinants of myiasis dynamics are very scarce. The genus Philornis include species that may be excellent models to study myiasis ecology, as they exclusively parasitize bird nestlings, which stay in their nests until they are fully fledged, and larvae remain at the point of entry until the parasitic stage is over, thus allowing the collection of sequential individual-level infection data from virtually all the hosts present at a particular area. Here we offer a stratified multi-level analysis of longitudinal data of Philornis torquans parasitism in replicated forest bird communities of central Argentina. Using Generalized Linear Models and Generalized Linear Mixed Models and an information theory approach for model selection, we conducted four groups of analyses, each with a different study unit, the individual, the brood, the community at a given week, and the community at a given year. The response variable was larval abundance per nestling or mean abundance per nestling. At each level, models included the variables of interest of that particular level, and also potential confounders and effect modifiers of higher levels. We found associations of large magnitude at all levels, but only few variables truly governed the dynamics of this parasite. At the individual level, the infection was determined by the species and the age of the host. The main driver of parasite abundance at the microhabitat level was the average height of the forest, and at the community level, the density of hosts and prior rainfall. This multi-level approach contributed to a better understanding of the ecology of myiasis.


Introduction
A vast variety of dipterans affects the health of animals, including humans, either by acting as disease vectors [1,2] or through parasitism [3]. Some flies with parasitic immature stages cause myiasis, a tegument or cavity infection with larvae, which may represent great economic losses for livestock industry (e.g. species of Lucilia [4]), public health concern (e.g. Dermatobia hominis [3]), or contribute to wildlife species declines (e.g. Philornis downsi [5]).
Parasites that cause myiasis are among those that actively seek their hosts over relatively long distances, which results in a particular life-cycle and natural history. Despite this, the ecology of myiasis has seldom been studied in detail. Some research has been conducted, but each study only sheds light on few selected factors that influence myiasis occurrence. Current data available include associations of parasitic dipterans with climatic factors [6,7], host density [7,8], and nest and host characteristics [8][9][10][11]. Myiasis ecology studies that jointly encompass the analysis of several factors at different levels of biological hierarchy would provide us with a better understanding of the web of causation behind myiasis dynamics, describing the relative contribution of each variable while enabling the assessment of interactions between factors at different levels.
Philornis Meinert (Diptera: Muscidae) is a genus of flies that includes several parasitic species, whose larvae parasitize New World bird nestlings [12]. Most parasitic Philornis spp. cause subcutaneous myiasis, with burrowing larvae that feed on nestling blood, tissue and fluids [13]. These parasites harm nestlings, causing mortality and reduced fitness and growth [7,10,14]. Subcutaneous Philornis larvae constitute a very interesting model to study the ecology of myiasis, as they principally parasitize bird nestlings, which stay in their nests until they are fully fledged, and larvae are easily spotted and remain at the point of entry until the parasitic stage is over, thus allowing the collection of sequential individual-level infection data from virtually all the hosts present at a particular area.
Here we offer a stratified multi-level analysis of longitudinal data of Philornis torquans (Nielsen, 1913) parasitism in two forest bird communities of central Argentina. The effect of host-, microhabitat-and community-level factors and some of their relevant interactions were assessed. Although many host and environmental variables might be associated with the occurrence of this parasite, it is expected that the integrative approach employed will enable us to pinpoint the major drivers of infection risk. While studies that place the focus on a limited group of risk factors are appropriate to demonstrate causal associations between a putative cause and a parasite, only integrative approaches allow us to assess what variables explain most of the variance when all of them are interacting together. We predict that, although many host and environmental factors will be found to be associated with Philornis infection, only a few truly drive the dynamics of parasitism.

Study area
The data were collected from two native forest patches located in the center of Santa Fe Province (Argentina). One is a reserve belonging to Universidad Nacional del Litoral (centre at 60u 559 W, 31u 239 S), the other one, a private field known as ''Mihura'' (centre at 60u 479 W, 31u 309 S) (Mr. Kling and Mr. Gimenez kindly allowed us to carry out our fieldwork in their property). The distance between the two sites is 20km. The area sampled within each site was 40 ha. Both sites represent relicts of the biogeographic province 'El Espinal' and are located alongside the Salado River. The climate in the region is Pampean Temperate, with an average annual temperature of 18uC (mean minimum = 12uC; mean maximum = 23uC) (extracted from www.climayagua.inta.gov.ar). These forests are breeding grounds for around 100 bird species, mostly Passeriformes [15]. More than 20 species have been found to be parasitized by Philornis torquans larvae in this region, but only three account for over 90% of the records [7].

Data collection
The fieldwork to obtain the data spanned the breeding seasons of 2008-2009 (September -April) and 2009-2010 (September -May). As previous observations showed that the parasitic larvae do not occur until around 10 weeks after the first eggs of the season are laid, the observations from the first weeks (September and October) were not used in the analysis. Every week, each area was exhaustively examined for active nests. All broods detected were thoroughly examined on a weekly basis. Nestlings from species known to be the predominant hosts in the region (i.e. Pitangus sulphuratus L., Phacellodomus ruber Vieillot, Phacellodomus sibilatrix Sclater) [7], were sampled three times a week (every 2-3 days), for a more detailed collection of information on individual level variables. All nestlings were thoroughly examined to determine the presence of Philornis larvae, which were recorded and classified into developmental stages from instar one to three. In the two breeding seasons (totaling 70 weeks) and from both sites, we collected 6629 observations from 3255 nestlings (1400 broods) of 57 bird species. Of them, 2291 observations belonged to predominant hosts (43 broods of Ph. ruber, 79 broods of Ph. sibilatrix and 70 broods of Pi. sulphuratus). Other species that contributed with considerable amount of data (1423 observations) were Furnarius rufus Gmelin (85 broods), Paroaria coronata Miller (89 broods) and Sicalis flaveola L. (105 broods). The data will be made freely available upon request.
The factors evaluated were chosen for their biological sense on the basis of current knowledge on the bio-ecology of birds, dipterans and ectoparasites. A preliminary exploratory analysis allowed us to identify independent variables that were highly correlated between them, of which we chose the most strongly associated with the response to avoid co-linearity problems.
The data obtained at the individual level were species, age (days from hatching), body mass (g), tarsus length (mm) and body mass index (body mass/tarsus length). Additionally, once a week, a small blood sample (,10 ml) was taken from a nail cut into a heparinized capillary tube to carry out haematological analyses as described by Lucas and Jamroz [16].
Microhabitat was referred to as the environment surrounding the brood, including the more immediate (e.g. nest material) and less proximal (e.g. vegetation structure) factors. Nest variables were the species which built the nest, nest external material (grass, mud, sticks, vegetable fiber or wood), chamber material (feathers, grass, sticks, vegetable fiber or wood), degree of environmental exposure of the chamber (open/closed), brood size, nest height (meters above the ground), and nest support (tree species, ground, etc.). The immediate surroundings of the nest were classified according to proximity to water bodies and the vegetation structure, which consisted of determining presence/absence and average height of three different vegetation strata (trees, bushes and grass), and their type (dominant species for each vegetation strata). To this aim, the study sites were characterized in situ, measuring the proximity to water bodies and assessing the vegetation structure for reference nests, so that all the area of each site was covered. Nests that were in close proximity to a reference nest (,30 m) received the same environmental values.
The variables recorded at the community level were weekly precipitation (mm; from the closest meteorological station), weekly mean temperature and relative humidity (using HOBO Pro v2 data loggers placed in the middle of each study area), and total and specific bird nestling density (number of nestlings per Ha in total and by group of bird species).
Finally, from one of the sites (the reserve), we have systematically recorded data on brood density and Philornis parasitism from the breeding season of 2006-2007 through 2011-2012 (6 years), which enables a preliminary analysis to explore inter-annual fluctuation trends in parasitism and their associations with bird community structure and climatic variables, giving a different time-scale context to the results.
For a detailed description of all variables see Table S1.

Statistical analysis
The data analyses were conducted at three hierarchical levels, the individual, the microhabitat (nest and surroundings), and the community. Instead of approaching the analysis using a single hierarchical model, which is the usual method used for multi-level analyses [17], we utilized a stratified analytical approach, examining different study units separately, as we judged it more informative and robust. The analyses at lower levels placed the focus on variables acting at that level (e.g. host age), but conveniently included variables from higher levels to avoid potential confounding phenomena (e.g. if more rainfall increases both Philornis intensity and nestling body growth, ignoring rainfall in the data analysis at the individual level will produce a distorted association between tarsus length and parasitism). Relevant interaction terms were also assessed.
The analyses were conducted using generalized linear mixed models (GLMM) and generalized linear models (GLM) with negative binomial responses. The software used was R (The R Project for Statistical Computing; http://www.r-project.org) and the specific package was glmmADMB (function 'glmmadmb').
Model selection was carried out using an information theory approach, following Burnham & Anderson [18]. Starting from a global (maximum) model, model selection and comparison was carried out in a stepwise manner using Akaike information criteria (AIC) or, when overdispersion was present, Quasi-Akaike information criteria (QAIC) [19]. Our approach to model selection was as follows: firstly, all models with AIC values no greater than 5 units compared to the best model were considered. Then the DAICi between models was calculated as: where AIC min denotes the AIC value of the best model [20]. DAIC i was used to estimate the relative likelihood (l i ) The latter allowed the calculation of the weight (w i ) of each model Finally, the models were ordered according to their w i , and the top models that summed 0.9 (cumulative w i ) were selected [18]. The multimodal inference was done using the weighed mean of the b coefficient and its standard error. The terms that were considered significant were those with a coeficient's 95% confidence interval that did not include 0. The analysis at the individual level was done with GLMM and included only the data collected from predominant hosts (Pi. sulphuratus, Ph. ruber, Ph. sibilatrix). The study unit was the individual chick at a given time (t i ). The response variable was larval abundance per nestling [21], i.e. the number of larvae that were parasitizing the chick at t i (the value for a non-parasitized host was '0'). The independent variables of interest were host species, age with linear and quadratic terms, previous measures of condition (taken at previous visits: 2-3 days previously, t 21 ; or 4-5 days previously, t 22 ), which included body mass, body mass index, and total red blood cells and white blood cells. The potential confounders and interactions included are described in Table S2. Briefly, they were nest variables such as brood size and nest material, week and its quadratic term to account for differences along the breeding season, and the community variables site, year, minimum and maximum weekly mean temperatures, and weekly cumulative precipitation, each at different time lags (t 0 through t 26 ). We also controlled for the density of all or each of the three predominant host species at different time lags (t 0 through t 24 ). In addition, densities of nestlings of unsuitable host species (nonpasserines) and an index of biodiversity (Shannon's) [22] were assessed at different time lags, to explore a possible 'dilution effect' [23]. For variables that included time lags, the mean or sum estimates for a week were calculated using the values from seven consecutive days. For lags at t 0 , the values included were those of the day of the observation (day 0 ) plus the ones of the preceding 6 days (day 21 , day 22 , day 23 , day 24 , day 25 , and day 26 ); for lags at t 21 we used day 27 through day 213 , and so on. Relevant two-way interactions (between two variables of interest or a variable of interest and a potential confounder) were assessed (see Table S2). To account for the lack of independence between observations from the same nest and nestling we included the nested random effect 'nest id/individual id'.
At the microhabitat level, the analysis was also done with GLMM and included data collected from the six most abundant passerine species (the predominant ones plus F. rufus, P. coronata and S. flaveola). The study unit was a single brood at t i . The response variable was the mean larval abundance per bird, i.e. the number of larvae/number of nestlings in the brood. The independent variables of interest were the brood species, the species that built the nest, brood size, presence of nestlings of other species (i.e. parasitic: cowbirds, Molothrus bonaerensis Gmelin and striped cuckoos, Tapera naevia L.), environmental exposure of the nest chamber (open/closed), nest height, external and internal nest material (see above), and presence/absence, type and average height of different vegetation strata (trees, bush and grass). Quadratic terms were added to account for non-linearity of continuous independent variables (e.g. nest height). The potential confounders included were the same community level variables assessed for the individual level models (see Table S2). Also, relevant two-way interactions were assessed. The random effect included was 'nest id'.
A GLM was used to analyze the data at the community level. In this case, the study unit was the whole bird community (either from the reserve or Mihura) at a given week (t i ). The response variable was the mean larval abundance per bird (total larvae t i / number of nestlings t i ) 610 (i.e. the mean number of larvae per ten nestlings at a given week). The independent variables of interest were week, minimum and maximum weekly mean temperatures and humidity at different time lags (t 0 through t 26 ), weekly cumulative precipitation (sum of mm fallen at t i ) at the same time lags, and density of hosts at different time lags (t 0 through t 26 ), as follows: potential hosts (22 spp.; all host species that were found to be parasitized by P. torquans), predominant hosts (the 3 spp. named above), unsuitable hosts (non-passerines), and single species: Pi. sulphuratus, Ph. ruber and Ph. sibilatrix. Past abundance per bird of third instar larvae (t 21 through t 26 ) was also assessed. The corresponding quadratic terms were also included. The potential confounders assessed were site and year. The relevant two-way interactions included were those between year, site, temperature, precipitation and host densities at different time lags (see Table S2).
For inter-annual comparisons we used data from seasons 2006-2007 through 2011-2012. Each season, the sampling started the week of September 21 st , and ended when no more breeding was detected, except in 2006-2007, when the fieldwork was interrupted at the beginning of March because the study area was flooded, and in 2010-2011, when it had to be interrupted by the end of January due to health problems of the field workers. To account for this we used data from October to February. The analysis consisted in Spearman's Rho correlation tests between annual measures of parasitism (annual prevalence and mean parasite abundance per nestling) and climatic and host abundance data.

Ethics Statement
The animal care and treatments used in this study were approved by the Bioethical Committee of the School of Veterinary Medicine of Universidad Nacional del Litoral (PROTOCOL No. 02/08).

Results
Larvae of P. torquans were found in 22 out of 57 species examined. At the individual level, the highest levels of parasitism were observed in Pi. sulphuratus, with 41.2% parasitized and a mean infection intensity of 11.1 larvae/nestling, followed by Ph. ruber with a prevalence of 12.6% and a mean intensity of 13.4 and Ph. sibilatrix, with 10.6% prevalence and a mean intensity of 3.6. At the brood level, the prevalence and mean intensity were, respectively, 44.9% and 45.0 larvae/brood for Pi. sulphuratus; 35.3% and 8.9 for Ph. sibilatrix; 24.2% and 53.8 for Ph. ruber; 24.2% and 15.2 for F. rufus; 16.1% and 8.7 for Pa. coronata and 10.7% and 3.2 for S. flaveola. For a summary of the mean parasite abundance and prevalence by year, see Table 1.

Individual-level variables associated with infection
The individual variables that consistently remained as important predictors of P. torquans infection were host species, age (linear and quadratic terms) and the interaction between them ( Table 2).
The pattern of infection varied with the age of the nestling, but it was different in each of the three species. The larval abundance per nestling increased immediately after hatching to reach a peak and decrease afterwards, but the peak occurred at different times for each host species: at 2 days of age for Ph. sibilatrix, 6 for Ph. ruber and 12 for Pi. sulphuratus (Figure 1).
There was a strong interaction between species and Pi. sulphuratus density two weeks previously. While the larval abundance per nestling in Ph. sibililatrix was consistently lower than in the other two species, P. torquans parasitized Ph. ruber mainly two weeks after the densities of Pi. sulphuratus had been low (Figure 1), indicating there is plasticity in host preference dependent on host availability.
All other individual variables assessed were not significantly associated with Philornis larva abundance.

Microhabitat-level variables associated with infection
The analysis at this level was conducted using 643 observations from 312 broods of the six most abundant species. The models with best goodness of fit are shown in Table S3. Broods of Pi.  Figure 2). Regarding the vegetation surrounding the nests, there was a significant interaction between dominant tree species and forest height (with both linear and quadratic terms), indicating that the occurrence of P.torquans declined with the height of the forest, except in those dominated by Gledittsia triacanthos L., the only exotic tree species (Table 3; Figure 2). In low to medium height forests, the mean larval abundance was greater when the dominating species were Acacia caven Molina, Geoffrea decorticans (Gill. ex Hook. & Arn.) Burkart and Celtis tala Gillies ex Planch than in forests of Prosopis aff. alba, Aspidosperma quebracho-blanco Schltr. and G. triacanthos. Nests surrounded by bush had significantly more larvae than those built in the absence of this vegetation cover, while the association with the herbaceous stratum height was negative (Table 3; Figure 2). The relationship between the microhabitat variables and Philornis did not depend on the brood species (the relevant interaction terms were not significant).

Community-level variables associated with infection
The best models showing the community-level variables associated with mean parasite abundance are shown in Table S3. Past precipitation and current and past host density appeared to be the most important variables (Table 4; Figure 3).
There was a strong positive association between mean larval abundance and the precipitation 4 and 5 weeks previously. The mean minimum temperature 5 weeks in the past was also positively associated with mean Philornis abundance. Past maximum temperature (6 weeks previously) showed a weak, though significant, positive association with mean parasite abundance during the season 2009-2010, but such an association was not observed during the first breeding season.
The density of suitable hosts was also a very important driver of Philornis abundance (Table 4; Figure 3). The most strongly positively associated variables were current Pi. sulphuratus density, past predominant host density (6 weeks previously) and past potential host species (3 weeks previously). Increases in such densities resulted in strong increases in mean parasite abundance at the community level. On the contrary, very significant negative associations were observed with past densities of non-passerine birds in one of the study sites. In Mihura, increases in the density of non-passerine nestlings in that site were followed by decreases in the mean larval abundance in the community, but this association was not present in the reserve. Less strong but consistently significant in both sites was a negative association with Ph. ruber two weeks previously.

Interannual comparisons
The data of the interannual comparisons are shown in Table 1. The only statistically significant associations observed in this preliminary analysis were one between cumulative precipitation during the first two thirds of the season and mean parasite abundance (Rho = 0.943; p = 0.017) and between Pi. sulphuratus abundance and annual prevalence (Rho = 0.943; p = 0.017) (Figure 4).

Discussion
Host-parasite interactions may be influenced by mechanisms operating at a range of biological levels. Studies at different scales of hierarchy contribute to a better understanding of the ecology of wildlife diseases [24]. Such an approach allows a more holistic perspective, thus unveiling the relative size of the effects and facilitating the elucidation of the mechanisms underlying the dynamics of parasites. However, because of its difficulty, such an integrative approach has seldom been used. Although there are prior attempts [25], it was not until the late 1990's that multi-level (or hierarchical) disease studies were seriously undertaken (reviewed by Diez Roux [26]). These hierarchical approaches mostly consisted of including variables from different hierarchical scales in the same statistical models. While this is an appropriate procedure to control for the potential confounding effect of higherlevel variables, the inference about the effect of higher level variables may be flawed [27,28]. We used a stratified analytical approach using different units of analysis, as it can yield more informative and robust results. We evaluated the factors determining the variation in P. torquans abundance using four different study units at increasing hierarchical scales, from the intimacy of an individual host, to the dynamic complexity of the bird community. Changing the level of the study unit allows us to truly establish how interactions at lower levels of the biological hierarchy might affect the parasitism dynamics at higher levels. Moreover, examining the effects across the hierarchical gradient provides the opportunity to assess the true biological significance of an observed association. For example, as P. torquans frequently parasitizes Ph. ruber one may expect that the dynamics of this parasite are positively influenced by the abundance of broods of this bird species. However, the hierarchical model with the individual as the study unit shows a very strong negative influence of Ph. ruber density on the number of larvae per nestling (Table 2). When following the community in time, however, this negative association between Ph. ruber density and Philornis abundance persists, but at a much lower magnitude (Table 4; Fig. 3; e.g. a 10fold increase in previous Ph. ruber density results in a 200-fold decrease in current larval abundance per bird in the individual level model, but only a 50% decrease in the community level model). Yet at a longer time scale, the influence of Ph. ruber density was not strong enough to be detected. This clearly shows that using a single hierarchical model with the individual as the study unit would overestimate the negative influence of Ph. ruber density on Philornis dynamics.
We found several significant associations at all levels, but only few variables truly governed the dynamics of this parasite. At the individual level, the infection was determined by the species and the age of the host. The main driver of parasite abundance at the microhabitat level was the average height of the forest, and at the community level, the density of hosts and prior rainfall.

Individual level
The abundance of P. torquans larvae was largely influenced by host species. Even though a previous study in the same region showed that three species were parasitized above the rest in prevalence and intensity [7], here we find a very marked difference between these predominant hosts, with Pi. sulphuratus the most frequently and intensively parasitized. However, this host preference was highly dependent on the density of Pi. sulphuratus two weeks previously. Adult flies seek the preferred host, Pi. sulphuratus, but when its broods are scarce they can adapt to other hosts, principally Phacellodomus in the region studied. This host-preference plasticity determines a temporary host-straggling process that allows the parasite to persist when the main host is not present. However, the 2-week delay is puzzling. It is noteworthy that this 2week lag was consistently found at different scales of analysis in both negative and positive associations with densities of different host species. It might be related to the behavior of adult flies, which is largely unknown. One possible explanation is that after emerging from the puparium, the flies seek an active nest to develop their adult stage, and it may take some days until gravid females are ready to lay their larvae on the nestlings.
Nestlings are susceptible to parasitism by Philornis for the whole duration of the nestling period, but previous reports showed that the number of larvae increases with age, reaching a peak by the end of the second third of the period, after which it declines [10,29,30]. The explanation offered for this pattern is that body surface available, the tolerance to parasitism and chances of being found by an adult fly increases as the chick ages, but older nestlings will fledge before the larva can complete development and therefore the parasite would risk not finding a suitable substrate for pupating [10,30]. In our study, the relationship between nestling age and larval abundance was consistent with this pattern, but the peak took place at very different ages depending on the host species. These different temporal patterns of infection are particularly striking considering that it is the same parasite species, and suggests there might be a plasticity in the parasite synchronization [31], or perhaps a differential immunocompe- Table 3. Variables of interest of the best models describing microhabitat (nest and surrounding) factors associated with mean P. torquans abundance in a brood. tence, immune response, or a dissimilar maturation of the immune system between species [32]. Further studies should be aimed at elucidating the mechanisms behind this pattern.

Microhabitat level
The influence of vegetation structure surrounding the nest on Philornis abundance was remarkable. We found significant associations with every vegetation stratum. While nest height had no effect on brood infection, the consistent and strong general pattern was a reduction in mean larval abundance as the average forest height increased (except in forests dominated by the only exotic species). This suggests that what is important for the parasite is not the height at which its hosts are, but the microenvironment associated with differential forest height. The presence of medium stratum and the lower stratum height were also associated with P. torquans abundance. Segura and Reboreda [8], studying Pa. coronata broods, did not find an association between P. seguyi and vegetation structure. On the other hand, Le Gros et al. [33] did find an association with vegetation structure, but contrary to our findings, as they observed a greater proportion of nests parasitized by Philornis porteri in grasslands than in bushlands or forests, indicating that the influence of the vegetation structure on the dynamics of Philornis parasitism depends on the parasite species involved. This strong association with vegetation structure alerts us to the impact that land use change, deforestation and other disruption of the forests caused by humans may have on the dynamics of this hostparasite interaction.
Regarding the more intimate microhabitat, the nest factors studied were not found to be important (only associations explained by host species were found). Species that build the nest with similar materials and structure (Ph. sibilatrix and Ph. ruber) showed significant differences in their levels of Philornis parasitism. Likewise, S. flaveola, which uses abandoned F. rufus nests to breed, had much less larval abundance than the latter, indicating that the species is more important than the nest characteristics.

Community level and inter-annual comparisons
Our finding that climatic variables, namely precipitation and temperature, were positively associated with mean larval abundance coincide with previous myiasis studies [6,7,34]. The time lags involved might provide clues on the duration of parts of the life cycle of the fly, which is so far ignored. The association with precipitation was much stronger than the one with temperature, and was clearly evident in the preliminary inter-annual comparison.
Other strong driver of Philornis abundance at the community level was host density. High densities of Pi. sulphuratus nestlings resulted in concomitant high mean larval abundance in the community, which merely indicates that the main host is particularly abundant during that week and so the chances of finding larvae are high, as when we re-ran the analysis excluding Pi. sulphuratus hosts in the response (i.e. mean larval abundance in hosts other than Pi. sulphuratus) such an association disappeared, showing that more Pi. sulphuratus did not result in higher parasitism in other hosts (data not shown). The positive association with past densities of predominant hosts was expected, as many suitable hosts in the past meant more resources for larva to proliferate resulting in greater production of adult flies. This delayed-density dependence of parasite abundance has been documented for several host-parasite systems [35]. Significant negative influences of past host abundances were inconsistent. In one of the study sites, at identical densities of predominant hosts, we observed that high past densities of non-hosts were associated with a reduction in the mean larval abundance in the community, suggesting a dilution Table 4. Parameters of the best models describing the community level factors associated with mean P. torquans abundance in a given week.  effect [23]. The dilution effect has been demonstrated for a tickborne disease [36] and a rodent virus transmitted by direct contact [37], but never for a myiasis. However, our results should be interpreted with caution, as the association was found in only one site, and when we re-ran the analysis at the community level focusing on the main host (i.e. the response variable was larval abundance on Pi. sulphuratus), the association with past nonpasserine nestling density disappeared. Densities of Ph. ruber also showed a negative association with Philornis. High densities were followed by a reduction in Philornis abundance two weeks later. A hypothesis that would explain this unexpected result is that the parents and older Ph. ruber nestlings may feed on adult Philornis flies that approach to lay their propagules or as they hatch from the pupae that are lodged in their nests (Ph. ruber is an insectivorous species).
Time lags appeared to very important in this system, especially for the community-level analysis. It should be taken into account that using many time lags may present over-fitting problems [38]. However, the strength of the associations found in our models argues against the probability of incurring in a type 1 error due to over-fitting.
The preliminary interannual comparison showed that Philornis abundance was associated with the precipitation during the first two thirds of the breeding season, and that the annual prevalence depended on the number of Pi. sulphuratus nestlings present in the area. It is remarkable to observe that with such limited statistical power (n = 6) we found two statistically significant associations, which indicates that the relationships are of great magnitude, consistent with what was observed at lower time and hierarchy scales. Moreover, an exploratory multivariable linear model found three variables significantly positively associated with mean parasite abundance: precipitation, Pi. sulphuratus abundance, and mean minimum temperature (adjusted R 2 = 0.9996; see Table S4).
Dudaniec et al. [39] observed P. downsi parasitizing nestlings at higher intensities in years of high precipitation. They interpreted that increases of P. downsi intensity in rainy years were determined by a positive response of host numbers to rainfall, and not by a direct effect of rainfall or temperature. Here we show that the combination of multi-level and multi-variable analyses allows us to dissect the effects of host density from those of climatic variables. Besides the significant influence of host density, we found an additive direct effect of temperature and rainfall.
It is noteworthy that both climatic and host density variables were consistently found to be important predictors in the lower hierarchy models.

Conclusions
The stratified multi-level approach employed permitted an integrative perspective, unveiling the relative size of the effects and the assessment of interactions between variables at different levels, discriminating between critical and less important factors, and enabling us to observe that the direction and strength of the main relationships were consistent across hierarchical levels and in some cases even through different time scales. For instance, rainy weeks were followed by weeks of increased parasite abundance, and rainy years yielded greater overall parasite production. This helps us realize that precipitation is a driver of Philornis abundance both in the short and the long term. Other major driver that consistently appeared strongly associated with Philornis infection across all hierarchy scales was the density of suitable host nestlings. These two findings entail that there is potential for global environmental change to affect the impact of Philornis on their hosts. Forest fragmentation is resulting in high brood densities in the remaining forest patches [40], and in some areas of South America the weather is becoming wetter and warmer [7,41], a scenario that might enhance conditions for Philornis (although flies are likely to do poorly at higher temperatures not measured in our study).
A broader perspective proved to be more informative than a focalized approach. For instance, the vast majority of Philornis studies placed the focus on a single host species. In this study we demonstrated that the infection pattern of the same parasite species along the duration of the nestling period was largely different in each of the three host species examined. At the same time, studying specific relationships at lower levels of the biological hierarchy produces complementary information that helps unveil underlying mechanisms. For example, an investigation of the immune response against Philornis in Phacellodomus spp. and Pi. sulphuratus would establish if there are differences in tolerance and resistance between these host species, which would in turn contribute to understanding why there is such a marked host preference, and why months and years with high Pi. sulphuratus abundance result in months and years with greater P. torquans abundance, respectively.
A stratified multi-level analysis proved useful to improve our mechanistic understanding of ecological processes associated with parasites and to assess at higher levels the net effect of processes interacting at lower levels. Biogeographical studies and long-term data sets will be crucial for characterizing the context under which these net effects may be different.
Further studies using experimental approaches should explore in more detail the associations identified here to determine its causal nature and establish their underlying mechanisms. In addition, this kind of integrative empirical data may prove invaluable to support theoretical research aimed at explaining and predicting how parasitism at the level of individuals scales up into complex ecological processes [24].