Geographic Variation in Network Structure of a Nearctic Aquatic Food Web

ABSTRACT


INTRODUCTION
7 and definitions of climate variables available in appendix S1) that are highly correlated with latitude 137 across our study sites. In the Sarracenia food web, total species richness is greater at higher latitudes that 138 generally experience lower and more variable temperatures along with lower and more variable amounts 139 of precipitation (Buckley et al., 2003;2010). We predict that food chain length and complexity will 140 increase with latitude due to the greater probability of omnivore and top predator presence in high latitude 141 species-rich webs. Our final prediction is based on the previous finding that within-site compositional 142 turnover shows an inverse relationship with latitude (Buckley et al., 2010). We predict that within-site 143 variability in web structure will track compositional turnover and increase at lower latitudes due to the 144 lack of trophic redundancy in the low latitude species-poor pools (Baiser &Lockwood 2011). 145

METHODS 146
The Sarracenia food web 147 The food web inhabiting the aquatic microhabitat in the leaves of the northern pitcher plant is comprised Feeding interactions in the Sarracenia food web center on a detritus "processing chain" 155 (Bradshaw and Creelman, 1984;Heard, 1994). Prey items that are captured by the plant are shredded by 156 the midge and the sarcophagid fly into particulate organic matter (POM). Bacteria directly decompose 1). We determined the abundances of invertebrates, rotifers, protozoa, and bacteria in each pitcher. We 163 counted and identified all invertebrates in each pitcher and rotifers and protozoa in 0.1-mL sub-samples 164 using a phase-contrast scope at 100×. Protozoa were identified to genus where possible and unidentified 165 protozoa were not used in food web calculations (there were 16 unknown protozoa, 13 of which occurred 166 in less than 2% of pitchers and three which occurred in 6-18 % of pitchers). Bacterial abundances were 167 estimated using serial dilutions (10 -5 and 10 -7 ) for each leaf and plating out samples on half-strength Luria 168 broth agar. Thus, only plate culturable bacteria were included and identified by colony morphotypes. We 169 calculated the density of the aforementioned taxonomic groups as abundance/mL in each pitcher. In total, 170 75 taxa were included in web calculations (see food web metrics below). We determined latitude for each 171 site using the satellite global positioning system and recorded the total volume of pitcher fluid for each 172 pitcher (see Buckley et al., 2003 and 2010 for details on site selection, leaf selection, sampling protocol, 173 and a complete list of species found in all food webs). 174

Food web metrics 175
Feeding interactions (hereafter links) between the species of the Sarracenia food web were assigned 176 based on previous studies (Addicott, 1974;Forsyth & Robertson, 1975;Heard, 1994;Cochran-Stafira & 177 von Ende, 1998;Miller et al., 2002;Butler et al., 2008) and direct observation. We constructed an n × 2 178 matrix for each of the 780 food webs, where the n rows are the number links; the first column of the 179 matrix contains the predator species identity and the second column contains the prey species identity for 180 each link. We used Network3D (Williams, 2010) to calculate a suite of 15 metrics that characterize 181 complexity, chain length, type of taxa, and variation in trophic strategy for the Sarracenia food web 182 (Table 1). Because metrics for most well studied food webs co-vary to some degree (Vermatt et al.,

Nitrogen data 186
Pitcher plants receive atmospheric N (in the form of NH 4 and NO 3 ) from rain and snowmelt that fill the 187 pitchers. This atmospheric deposition can affect pitcher morphology and habitat structure for the food 188 web (Ellison & Gotelli, 2002), and pitcher plant population dynamics (Gotelli & Ellison, 2002). We 189 estimated deposition levels at each sampling site during the year of the survey to investigate these 190 potential effects on food web structure. We used nitrogen deposition from National Atmospheric 191 Deposition Data (NADP) monitoring stations that were closest to our sample sites in the United States 192 sites. Therefore, for consistency, we estimated N deposition (total N = NH 4 + NO 3 as precipitation-194 weighted mean concentration in mg/L) at all our sites (i.e. United States + Canadian sites) in the summer 195 quarter (July-September) using the AURAMS model (Moran et al,. 2008) and used this estimate as a 196 predictor variable for all sites in our analyses of Sarracenia food webs The estimates for United States 197 sites were well-correlated with empirical NADP data (r = 0.66, p < 0.0001), and we assumed similar 198 accuracy for Canadian sites. Further details on modeling N deposition are given in appendix S2. 199 and PC2) as predictor variables. Climate variables and latitude were transformed to standard deviation 208 units for the principal components analysis. To insure that any one climate variable did not account for 209 the majority of the variation in a given food web metric, we ran a set of preliminary univariate regressions 210 with each climate variable, latitude, PC1, and PC 2 as predictor variables and food web metrics as 211 response variables. We ranked models using the Akaike Information Criterion (AIC), and used the AIC 212 score to select the best fitting model(s) among the candidate set (Burnham &Anderson, 2002). If any 213 single climate variable had a ΔAIC < 2 when compared to the first primary component axis (PC1), we 214 selected that climate variable for further consideration in the regression models described below. If no 215 single variable distinguished itself as a better fit (ΔAIC < 2) than PC1, PC1 was selected for further 216 consideration in the regression models described below. PC1 was the climate variable used in regression 217 models for all but one case. 218 219

Pitcher scale 221
We used linear mixed effects models (function "lme" in package "nlme" in R v.2.11.1) to assess the 222 influence of predictor variables on food web structure for the 780 individual pitchers (the pitcher scale of 223 analysis). We used food web metrics as response variables, site as a random effect, and PC1, pitcher 224 volume, nitrogen deposition, mosquito abundance, and bacterial abundance as fixed effects for the 225 pitcher-scale analyses. 226 We built a set of candidate models for each response variable that included a null model (i.e. 227 random intercept only), global model (with random intercept and all predictor variables entered), 228 univariate models for each predictor variable, and all subsets of variables that had a P-value < 0.1 for the 229 slope coefficient in both the global and univariate models (model structure is given in appendix S4). We 230 ranked models using the Akaike Information Criterion (AIC), and used the AIC score to select the best fitting model(s) among the candidate set (Burnham &Anderson, 2002). We calculated the variance 232 explained (R 2 ) by the fixed effects in this mixed-effects model using Xu"s (2003) method (see appendix 233 S5 for calculation). 234

Site scale 235
Variation in species richness and compositional turnover are greater within sites than across sites in the 236 Sarracenia food web (Buckley et al., 2010). Therefore, to measure variation in food web structure within 237 sites, we calculated the coefficient of variation (CV) of each web metric in the 20 webs at each of our 39 238 sites (site scale). We term this measure structural turnover and use it in the same sense as compositional 239 turnover (i.e. β-diversity). High structural turnover means that when moving from one web to the next we 240 are likely to encounter different network structure; high structural turnover results in a high CV at that 241 site. Low structural turnover (measured as a low CV) means that web structure is similar from pitcher to 242 pitcher within a single site. 243 We used linear models (lm in R v.2.11.1) to assess the influence of predictor variables on 244 structural turnover at each of our 39 sites. We regressed structural turnover (the CV of each food web 245 metric for the 20 pitchers at each of the 39 sites) on PC1, and the CVs of pitcher volume, mosquito 246 abundance, and bacterial abundance. The model for the food web metric Top included the climate 247 variable CV of precipitation instead of PC1 based on the climate variable model selection (See Climate 248 Data). We used the CV of pitcher volume, mosquito abundance, and bacterial abundance for this analysis 249 because we were interested in how pitcher-to-pitcher variation of predictor variables within each site was 250 correlated with pitcher-to-pitcher variation in food web metrics across all pitchers within each site. 251 We built a set of candidate models for each response variable that included a global model, 252 univariate models for each predictor variable, and all subsets of variables that had a P-value < 0.1 for the 253 slope coefficient in both the global and univariate models. We ranked models using the Akaike 254 Information Criterion (AIC), and used this score to select the best fitting model(s) among the candidate set (Burnham &Anderson, 2002). We calculated the adjusted R 2 to determine the proportion of variance 256 explained by each model. 257

Predictor variables 258
The correlations among predictor variables had correlation coefficients < 0.4. All variables were 259 transformed into standard deviation units (positive values indicate observations that were greater than the 260 mean and negative values were less than the mean) for the pitcher-scale analyses except for PC1. N 261 deposition was also transformed to standard deviation units for the site-scale analyses. 262

Principal components analysis of pitcher-scale variation in food web metrics 265
The first two principal components explained 70% of the variation in network structure for the set of 780 266 Sarracenia food webs. The first principal axis (PC1) explained 45% of the variation and was related to 267 complexity and chain length. This axis was negatively correlated with metrics related to complexity (e.g. 268 connectance, species richness, links per species) and chain length (e.g. mean trophic level, chain length) 269 (Table 2). Percentage of omnivores and intermediate species were also negatively correlated with PC1 270 (Table 2). Webs with negative scores on PC1 were species-rich and contained many omnivores, which 271 increased chain length, linkage density, and connectance (Fig 3a). PC1 was positively correlated with the 272 percentage of top species in a web, the percentage of detritivores in the web, variation in the number of 273 consumers and links per taxon, and the mean path length across the network (Table 2). An example of a 274 web with a high positive PC1 score contains only detritivores (e.g. bacteria), which are all top species in 275 this context (Fig. 3b). The second principal axis (PC2) explained 25% of the variation and was related to 276 trophic redundancy and variation in prey and predator strategies. Webs with a positive score for PC2 (Fig.  277 3c) tended to have more species, and these species were intermediate detritivore species (i.e. bacteria) and negative side of this axis (Fig. 3d) had fewer species and the species that dropped out were bacteria. 280

Principal components analysis of climate variables and latitude 281
The first two principal components explained 92% of the climatic variation across our 39 sites. The first 282 principal axis (PC1) explained 78% of the variation. Sites with high scores on this axis were at high 283 latitudes, had short growing seasons with low mean temperatures and precipitation, and had high annual 284 variation in both of these variables (Table 3). Sites representative of these conditions were located in the 285 northern US and Canada (Fig. 1). Sites with low scores are located at low latitudes and experience higher 286 mean temperatures and precipitation, but lower variability in both variables (Table 3). These sites are 287 located in the south-eastern US. PC2 explained 14% the climatic variation across our 39 sites. Sites with 288 positive scores had a high mean diurnal temperature range. 289 290 Resource availability, food chain length, and food web complexity 291 The hypothesized positive relationships between resource availability measured as bacterial abundance 292 and measures of complexity and chain length were not observed. Bacterial abundance showed no 293 relationship with any of the food web metrics at either the pitcher or site scale and was absent from all 294 best-fit models (Tables 4 and 5). 295 the site scale. However, chain length increased with ecosystem size only at the pitcher scale and showed 303 no relationship with variation in pitcher volume at the site scale. 304 Nitrogen deposition, food chain length, and food web complexity 305 Nitrogen deposition, which was predicted to have a negative effect on complexity and chain length 306 metrics, showed no relationship with food web metrics and was not a significant predictor in any of the 307 best-fitting models. 308

Predator-prey interactions, food chain length, and food web complexity 309
Mosquito abundance was not correlated with food web structure at the pitcher scale. However, variation 310 in mosquito abundance at the site scale was positively correlated with structural turnover in two measures 311 of chain lengthmean trophic level and mean chain length (Table 4; Fig 4 a, b). Although the 312 hypothesized relationship between mosquito abundance and chain length was observed at the site scale, 313 food web complexity did not increase with mosquito abundance at the pitcher scale and variation in 314 mosquito abundance within site did not increase structural turnover in complexity metrics. 315 Biogeographic correlates of food chain length, and food web complexity 316 PC1, which is positively correlated with latitude (Table 3) was present in all of the best-fitting models 317 that explained more than 5% of the variation in food web structure (Tables 4, 5). At the pitcher scale, PC1 318 was positively correlated with two measures of food web complexitylinkage density and species 319 richness. Species richness was shown to follow the same patterns and increase with latitude in previous 320 analyses of these data (Buckley et al., 2003;2010). Structural turnover at the site scale showed a 321 consistent negative relationship with PC1 for more than half of the food web metrics measured (Table 5; 322  Overall, network structure of the Sarracenia food web was only weakly influenced by all 324 predictor variables at the pitcher scale across 780 webs. The best-fit models left a large portion (> 95%) of the variance in food web structure unexplained at the pitcher scale (Table 4). At the site scale, predictor 326 variables explained more (8% -35%) variance in structural turnover in food webs (Table 5). Resource availability showed no relationship with food web structure at the pitcher or site scales. 339 One possibility why the predicted relationship was not observed is that bacterial abundance is not an Although mosquito abundance had no effect on food web structure at the pitcher scale, increased 347 variation in mosquito abundance between pitchers was positively correlated with structural turnover in 348 chain length and mean trophic level within sites. Mosquito larval density varied from 0 to over 11 larvae  providing nitrogen deposition estimates for our sites using their AURAMS model. Support for this 401 research was provided by NSF grants 0083617 to TEM, AME, and NJG, and 0541680 to AME.   Table 3. Factor loadings for the first two principal components axis (PC1 and PC2) describing climatic 584 variation across our 39 sites. PC1 explained 78% of the variation and PC2 explained 14%. PC1 describes 585 latitudinal variation in temperature and precipitation. Sites with high scores on this axis were at high 586 latitudes, had short growing seasons with low mean temperatures and precipitation, and had high annual 587 variation in both of these variables. Sites representative of these conditions were located in the northern 588 US and Canada (Fig. 1). Sites with low scores are located at low latitudes and experience higher mean 589 temperatures and precipitation, but lower variability in both variables. Climate variables are defined in 590 Appendix S1.     Table 1 for web metric definitions.) 625 calculated for 780 Sarracenia food webs. Component 1 is related to complexity and chain length and 626 component 2 is related to trophic redundancy and variation in prey and predator strategies. These two 627 components explain 70% of the variation in Sarracenia food web structure. Inlay, four networks (a,b,c,d) 628 representing the extremes of each axis. For these four food webs, arrows are drawn from their position in 629 the PCA plot to the food web. White nodes represent the resource (dead prey items), grey nodes represent 630 bacteria, and black nodes represent consumers. 631