Ecological Archives E096-099-A2
James T. Cronin, Ganesh P. Bhattarai, Warwick J. Allen, and Laura A. Meyerson. 2015. Biogeography of a plant invasion: plant–herbivore interactions. Ecology 96:1115–1127. http://dx.doi.org/10.1890/14-1091.1
Appendix B. Detailed census procedures and statistical methods.
Census procedures
Herbivore load and damage. Because of the diversity of arthropods associated with P. australis, a variety of methods were used to assess herbivore abundance and damage. For each P. australis stand, we walked a straight line from the edge to the interior (transect) and at every 2 m, the 3 stems nearest the line were inspected for the presence of herbivores or their traces (galls, dead apical meristems, leaf mines, chewing damage). A total of 50 stems were inspected and this was repeated for 3 transects that were widely spaced apart. We focused on the presence or damage from 3 main feeding guilds, internal stem feeders (gall formers and borers), chewing herbivores, and sucking herbivores (specifically, the mealy plum aphid H. pruni).
For internal feeders, we focused on their incidence of occurrence at the end of the summer census (July-August, 2011 and 2012) when gall formers and boring damage were visible. Most internal feeders can be identified by changes to the external appearance of the stem caused by feeding damage (Lipara spp., Trypticus sp., Tetramesa phragmitis). Species identity and composition of the internal-feeder guild was determined/confirmed by dissecting a random subset of stems with damage consistent with internal feeders (up to 50), preserving insects in alcohol and identifying them in the laboratory as necessary. Other internal feeders such as root gallers and leaf miners are either very difficult to quantify (the former) or are uncommon (the latter) and therefore were ignored.
For stems that possessed chewing damage, we collected and photographed all damaged leaves from each stem. Using ImageJ (Rasband 2012), the amount of leaf tissue consumed per stem (cm²) was computed (consumed). For comparison, we also photographed and determined the total leaf areas for each of 10 randomly chosen stems (leaf area). The mean proportion of leaf-area consumed per stem was determined as mean consumed/mean leaf area. Because leaves were green and persisted for most of the growing season (J. T. Cronin, unpublished data), a late-season census represented the cumulative chewing damage for the year.
Finally, for the aphids, we estimated their abundance per stem using the following method. Using the first 25 aphid-infested stems encountered along each transect, we counted the number of leaves with aphids. One infested leaf per stem was chosen at random and the aphids were counted. Mean aphids per stem was then estimated as the proportion of stems infested × mean number of infested leaves per infested stem × mean number of aphids/infested leaf. Because aphid colonies were relatively ephemeral and reached peak abundance in the spring at southern latitudes and in late summer at northern latitudes, we summed aphid densities from the 2 census periods within a year as our overall density estimate for the season. For number of aphids per stem, we only had data for 2012.
Plant nutritional quality, defensive trait levels and fitness. For each of the 3 census periods, the top 3 leaves were collected from 10 stems per P. australis patch and used for assessment of stem nutritional condition (percent carbon [%C], percent nitrogen [%N], carbon:nitrogen [CN ratio]), and attributes putatively associated with plant resistance (leaf toughness and total phenolics). The 3 leaves per stem were combined and placed together in a ziploc bag with silica gel dessicant packs. Packs were replaced as needed until the leaves were dry. In the laboratory, leaves were lyophilized and ground to a fine powder. P. australis nutritional condition was determined using an elemental analyzer at one of two independent laboratories, either LSU AgCenter Soil Testing & Plant Analysis Lab (https://www.lsuagcenter.com/en/our_offices/departments/SPESS/ServiceLabs/soil_testing_lab/) or Brown University Environmental Chemistry Facilities (http://www.brown.edu/Research/Evchem/facilities/). %N in plant tissues is often positively correlated with herbivore oviposition preference and/or performance (Price 1991, Cook and Denno 1994, Rossi and Stiling 1998, Salgado and Pennings 2005). The effects of carbon on nutritional condition are less straightforward in that increased %C can cause a reduction in %N or water content, or an increase in starch, leaf toughness or various defenses (see Hamilton et al. 2004).
Leaf toughness is often considered a resistance trait that affects the ability to consume plant tissues (e.g., Denno et al. 1990, Agrawal 2005, Salgado and Pennings 2005). Little work has been done on P. australis defensive chemistry, but grasses often produce secondary defense chemicals such as phenolics that can affect herbivore preference and performance (e.g., Gershenzon 1994, Abrahamson and Weis 1997). Leaf toughness was assessed in the field using a penetrometer (Itin Scale Co., Brooklyn, New York) which measures the force (in kg) necessary to puncture tissue (Salgado and Pennings 2005). Total leaf phenolics (nM/g of dried leaf tissue) were estimated using a modified version of the Folin-Ciocalteu method (Waterman and Moles 1994, Haynes et al. 2007, Hakes and Cronin 2012) that requires less dried leaf material. This modification involved using 96-well plates and a plate reader instead of cuvettes. The full protocol is available by contacting LAM.
Statistical Methods
For both the NA native-invasive comparison and the range of EU genotypes comparison, the data were analyzed with separate repeated-measures mixed-model ANOVAs (SAS 9.3, Proc Mixed). Within each comparison, we focused on the following dependent variables: %C, %N, CN ratio, leaf toughness (kg), total phenolics (nM/g of dried leaf tissue), proportion of total leaf area consumed, proportion of stems occupied by a gall former or borer, and number of aphids per stem. Census date (late summer 2011, early summer 2012, late summer 2012) for the plant nutritional and defense traits, and year (late summer 2011, late summer 2012) for the herbivory traits were treated as the repeated measure. Genotype (native or invasive) or Range (EU or NA) were treated as fixed factors and latitude, patch size (m²) and stem density (number per m²) were used as covariates. We included a quadratic term (latitude2) to evaluate whether the relationship between the trait and latitude was nonlinear and a genotype (or range) by latitude2 interaction to assess whether genotypes or the European genotypes from different ranges exhibited different nonlinear responses to latitude. For the herbivory traits, we also included all of the nutritional and defense traits as covariates in the model. Data were analyzed using SAS 9.3 Proc Mixed and maximum likelihood was used as the estimation method for the covariance parameters (SAS Institute, Cary, North Carolina). Finally, for aphid herbivory, i.e., the total number of aphids per stem, the statistical model was the same as above except for the absence of a repeated measure (year).
For each dependent variable, we used Akaike's Information Criteria corrected for finite sample size (AICc) to select the most informative model (Burnham and Anderson 2010). Our full model included the main effects and covariates listed above and only first-order interactions. First, AICc was used with the full model to choose the best covariance structure among repeated measures in the mixed-model design (Barnett et al. 2010). All methods available in SAS in the Repeated Statement were considered and the variance structure producing the smallest AICc value was used. Based on a combination of Cook's and likelihood distances (SAS Model option INFLUENCE), we assessed whether individual patches (unit of replication) had undue influence on parameter estimates and model outcomes (Christensen et al. 1992, Schabenberger 2004). These highly influential cases were removed from the model but in no case did their removal qualitatively change the conclusions of the model. Undue influence was evaluated in the full and best fit models.
Next, we employed a backward selection procedure to reduce the total number of variables in the model (Burnham and Anderson 2010). The interaction term or quadratic term with the lowest value of the F statistic was removed from the model and if the AICc score was reduced by ≥ 2, the term was omitted from the model. This procedure was continued for all factors excluding census date/year, genotype and latitude until there were ≤ 8 variables in the model or when the removal of any single variable did not achieve an AICc reduction of < 2. We considered eight parameters to be an upper limit for computing a manageable number of possible parameter combinations (see below). As part of our criteria, interaction terms were only allowed in the model if the two main effects making up the interaction were also retained in the model. It was possible that a main effect was retained in the model only because the interaction term contributed in a significant way (based on AICc scores) to the model fit.
Based on the variables in the reduced model, candidate models were constructed using all possible combinations of the variables. Restrictions to the possible combinations of variables included the requirement that interaction terms could only be present in the model if their main effects were also present in the model, and the repeated measures was included in all models. Without this underlying structure to the model, the design would be pseudoreplicated.
Candidate models were ranked by AICc from lowest to highest value and AICs with a Δi value (= AICci - AICcmin) of ≤ 2 were deemed to have substantial support (Burnham and Anderson 2010). We also report the AICc weights (wi) which indicate the weight of evidence (as a proportion) in favor of model i being the best model given the set of candidate models. As the Proc MIXED procedure does not report goodness-of-fit for the models, we emphasize effect sizes of the factors in the model (i.e., proportional differences in least-squares means or slopes in relationships). Because census date/year was included in every model, we report the associated P value for that factor as an indicator of whether there was any important temporal variation in the dependent variable.
To help normalize data distributions and homogenize variances among categories (i.e., census dates/years, genotypes), stem density, total phenolics, and all herbivory variables were ln-transformed. Even with the transformations, the distributions of herbivory had an excess of zeros. Otherwise, the data were approximately normally distributed. For the best model(s), we repeated the analyses with zeros removed. In almost all cases, the significances of the different factors in the model were not changed in any qualitative way. The exception was the analysis of aphids per stem in the NA-EU comparison of the European genotypes. However, the estimated effect sizes were very similar with and without the zeros. Because these zeros represent real data, we left them in the model.
Finally, one of the main comparisons that we were interested in was whether genotypes (NA-natives vs. NA-invasives) or ranges (European native genotypes in EU or NA) were different and exhibited parallel responses to the latitudinal gradient. Using the best model for each dependent variable, we re-ran the analysis with genotype/range and latitude removed from the model and generated model residuals. We then plotted the residuals with respect to latitude for each genotype/range. This approach allowed us to more clearly reveal the pure effects of genotype/range and latitude on the dependent variable after taking into account the effects of all other variables. We used a similar approach when graphing the relationship between other covariates (e.g., %C) and the dependent variable.
Latitudinal differences in the distributions of genotypes and ranges
An issue that we had to contend with regarding our data set was that the latitudinal ranges of the native and invasive genotypes in North America, and the European genotypes in Europe and North America were not exactly the same. In North America, the native populations included in the study ranged between 36.5° and 48.4° latitude, whereas populations of the European invasive genotypes ranged from 29.7° to 48.1° latitude. Populations of the European genotypes in the native range had a latitudinal range of 37.2° to 59.3°. To confirm that the differences in ranges were not the source of any differences between native and invasive genotypes in NA (or ranges of the European genotypes), we conducted an additional set of analyses in which we controlled for latitudinal range. For the comparison of native and invasive genotypes in North America, we used only populations between 36° and 48° latitude (± 1°). For the comparison of the European genotypes in their native and invasive range, we used populations between 37° and 48° latitude (± 1°). We analyzed these data using the same repeated-measures mixed-model ANOVAs as outlined above. Because we found no evidence of nonlinear relationships between our dependent variables and latitude in the previous analyses that spanned the entire latitudinal range of the study, we did not include any quadratic terms in the model. Also, because of the reduced statistical power associated with fewer number of replicate populations, we reduced the number of predictor variables to (1) those variables of primary interest in this study (i.e., genotype [or range], latitude, and the genotype [or range]*latitude interaction), (2) the repeated measure (season or year), and (3) variables that were common to the AICc best models in the previous analyses (genotype [or range]*year [or season], latitude*year [or season], and stem density; see Appendix C).
Literature cited
Abrahamson, W. G., and A. E. Weis. 1997. The Evolutionary Ecology of a Tritrophic-Level Interaction: Goldenrod, the Stem Gallmaker and its Natural Enemies. Princeton University Press, Princeton, New Jersey, USA.
Agrawal, A. A. 2005. Natural selection on common milkweed (Asclepias syriaca) by a community of specialized insect herbivores. Evolutionary Ecology Research 7:651-667.
Barnett, A. G., N. Koper, A. J. Dobson, F. Schmiegelow, and M. Manseau. 2010. Using information criteria to select the correct variance–covariance structure for longitudinal data in ecology. Methods in Ecology and Evolution 1:15-24.
Burnham, K. P., and D. R. Anderson. 2010. Model selection and multimodel inference: a practical information-theoretic approach, Second edition. Springer-Verlag, New York, New York, USA.
Christensen, R., L. M. Pearson, and W. Johnson. 1992. Case-deletion diagnostics for mixed models. Technometrics 34:38-45.
Cook, A. G., and R. F. Denno. 1994. Planthopper/plant interactions: feeding behavior, plant nutrition, plant defense, and host-plant specialization. Chapman and Hall, New York, New York, USA.
Denno, R. F., S. Larsson, and K. L. Olmstead. 1990. Role of enemy-free space and plant quality in host-plant selection by willow beetles. Ecology 71:124-137.
Gershenzon, J. 1994. Metabolic costs of terpenoid accumulation in higher plants. Journal of Chemical Ecology 20:1281-1328.
Hakes, A. S., and J. T. Cronin. 2012. Successional changes in plant resistance and tolerance to herbivory. Ecology 93:1059-1070.
Hamilton, J. G., A. R. Zangerl, M. R. Berenbaum, J. Pippen, M. Aldea, and E. H. Delucia. 2004. Insect herbivory in an intact forest understory under experimental CO2 enrichment. Oecologia 138:566-573.
Haynes, K. J., F. P. Dillemuth, B. J. Anderson, A. S. Hakes, H. B. Jackson, S. E. Jackson, and J. T. Cronin. 2007. Landscape context outweighs local habitat quality in its effects on herbivore dispersal and distribution. Oecologia 151:431-441.
Lambertini, C., I. A. Mendelssohn, M. H. G. Gustafsson, B. Olesen, T. Riis, B. K. Sorrell, and H. Brix. 2012. Tracing the origin of Gulf Coast Phragmites (Poaceae): a story of long-distance dispersal and hybridization. American Journal of Botany 99:538-551.
Meyerson, L. A., and J. T. Cronin. 2013. Evidence for multiple introductions of Phragmites australis to North America: detection of a new non-native haplotype. Biological Invasions 12:2605-2608.
Price, P. W. 1991. The plant vigor hypothesis and herbivore attack. Oikos 62:244-251.
Rasband, W. S. 2012. ImageJ, U. S. National Institutes of Health v. 1.47t, http://imagej.nih.gov/ij/, Bethesda, Maryland.
Rossi, A. M., and P. Stiling. 1998. The interactions of plant clone and abiotic factors on a gall-making midge. Oecologia 116:170-176.
Salgado, C. S., and S. C. Pennings. 2005. Latitudinal variation in palatability of salt-marsh plants: are differences constitutive? Ecology 86:1571-1579.
Saltonstall, K. 2002. Cryptic invasion by a non-native genotype of the common reed, Phragmites australis, into North America. Proceedings of the National Academy of Sciences of the United States of America 99:2445-2449.
Schabenberger, O. 2004. Mixed model influence diagnostics. Paper 189-29. Cary, North Carolina, USA.
Waterman, P. G., and S. Moles. 1994. Analysis of plant metabolites. Methods in ecology series. Blackwell Scientific Publications, Oxford, United Kingdom.