Temperature plays a pivotal role in defining the distribution of species and the fitness of individuals within species’ ranges. Phenotypic plasticity can allow individuals to cope with varying environmental conditions, including rapid climate change. Populations at range edges experience more variable conditions than core populations and thus are hypothesized to exhibit higher thermal plasticity. However, as the strength of plasticity often varies between individuals, it can also differ among local populations at range edges. We studied the extent of and variation in thermal plasticity for several traits within and between populations of the perennial herb Plantago lanceolata L. (Plantaginaceae) at its northern range edge. We sampled seeds from nine sites within a 50 x 50 km region and grew them under three temperature regimes in a greenhouse. We measured traits related to size, flowering, pathogen responses, and inflorescence pigmentation. We expected to find higher plasticity in traits less strongly connected to fitness, and that differences between individuals would outweigh differences between populations in underpinning this variation in plasticity. Our results show thermal plasticity in leaf size and abundance, flowering probability and abundance, and pigmentation. Notably, we also found increased pathogen symptoms and higher infection rates of one of two viruses screened, highlighting the potential for changes in pathogen sensitivity and exposure under climate change. Importantly, in all traits but flower abundance, more variation in plasticity was attributable to differences within populations than between populations. Although this contribution was small in magnitude compared to thermal effects on traits, the higher intra- versus interpopulation variation in plasticity suggests that differences between individuals provide most of the variation in thermal plasticity, which may be driven by small-scale variations in habitat conditions, highlighting the need for conservation strategies that consider microhabitat variation to support short-term adaptive responses to thermal variability.

Plant material

We collected seeds of Plantago lanceolata from nine local populations (habitat patches consisting of meadows or pastures) in Åland during August 2019 (Fig. 1; in Hällfors et al. in prep and Table S1 in ‘Supplementary Information.pdf’). In the Åland Islands, P. lanceolata occurs in a network of ca. 4000 meadows (Ojanen et al., 2013). We chose the seed collection sites for this study from nine populations in three different parts of mainland Åland to represent (1) a variety of areas across the Åland main island, (2) P. lanceolata habitat area extent, and (3) connectivity to other nearby local populations of P. lanceolata (Fig. 1 in Hällfors et al. in prep).

All meadows in Åland where P. lanceolata occurs are annually surveyed for the presence of the butterfly species Melitaea cinxia, the larvae of which use it as their food plant, and for the fungal pathogen Podosphaera plantaginis that infects P. lanceolata (Ojanen et al., 2013). During the first surveys at the beginning of the 1990s, the cover of P. lanceolata and Veronica spicata (the other host species of M. cinxia larvae) were also measured and their exact locations recorded (Ojanen et al., 2013). As a proxy for P. lanceolata population size, ground cover (m²) of the plant has been estimated within each patch (Ojanen et al., 2013). Habitat area is considered the extent to which the two host plant species of M. cinxia occurs in a patch. This measure of habitat area can be used as a proxy for fragmentation and population size and thus a rough proxy for genetic diversity (González et al., 2020). Connectivity between patches is evaluated annually and can be used as a proxy for gene flow among populations (Hanski, 1999).

As P. lanceolata was the most common or the only host species in the seed collection sites for this study (Ojanen et al., 2013), the habitat area corresponds well with the extent of the focal species. Based on the connectivity metrics and the species’ dispersal ability (with both pollination and seed dispersal mainly occurring within local populations), we consider each sampling site to be independent from the others and treat them as spatially unstructured variables in our subsequent analyses. However, to test that closely located sites with high connectivity values within regions are not more similar to each other than expected, we included a sensitivity analysis testing for consistency of results by leaving out one of the sites within the region in turn (see Statistical analyses below and results in ‘Supplementary Information.pdf’ Tables S2 and S3).

To obtain a representative sample of the diversity within each population, we collected seeds from 17–50 individuals in each of the nine chosen sites (i.e., per population; individuals sampled based on availability, following seed-collecting guidelines; ENSCONET, 2009; Table S1 in ‘Supplementary Information.pdf’). We collected seeds from each mother individual separately to allow the seeds of different mother plants to be sown across all treatments to measure within-population variation in plasticity. We deposited voucher specimens from each sampling site at the herbarium of the Finnish Museum of Natural History (H sensu Thiers, 2016). We left the seeds to dry and ripen at room temperature for a few weeks before the start of the experiment. 

Experimental set-up and thermal conditions

We conducted the experiment in the Viikki Plant Growth Facilities at the University of Helsinki during September 2019 – February 2020, in six greenhouse compartments with three temperature treatments each replicated once. The temperatures chosen for the treatments were based on average June and July temperatures in Åland, 1959–2018 (Fig. S1 in ‘Supplementary Information.pdf’). The Cold treatment day temperature was set at 17°C with night temperature 8°C; while the Mean temperature treatment was 20°C during daytime and 11°C at night; and the Warm treatment was set at 23°C during daytime and 14°C at night. High-pressure sodium lamps created light conditions and photoperiod were consistent across all treatments (18 h light and 6 h dark), with day and night temperatures following the photoperiod.

From each of the nine sampled populations, we randomly chose 12 mother individuals to provide 12 seeds for the experiment, with the expectation that this would result in at least 6 germinated seedlings from 10 mother individuals (one for each replicate, two per treatment; a total of 1296 seeds). Seeds were sown in mid-September 2019 during 4 consecutive days (September 17–20). We placed each seed individually in germination tray cells (160 per tray) filled with a potting mix (2:1:1 of sowing mix (Kekkilä WHS R8017), sand (Weber 1–2 mm), and perlite). After sowing the seeds, we sprinkled the soil surface with a 1:1 mixture of vermiculite and sand. We placed the germination trays on moisture-retaining matting on the greenhouse benches. They were watered automatically from below every 3 days for 3 minutes initially, but this was changed to 3 minutes every second day from September 24. This watering regime kept the soil moist in all temperature treatments.  

One month after sowing, we chose the experimental individuals by stratified random sampling from the entire pool of germinated seedlings. We repotted the chosen plants into 0.7 L pots containing soil of 4:2:1:1 parts of potting compost (Kekkilä WHS R8030), perlite, fine sand, and coarse sand. We allocated one seedling for each of the ten randomly selected successfully germinated mother individuals to each treatment-by-replicate combination. When fewer than two individuals per mother plant germinated in a treatment, we randomly chose a germinated seedling from another mother plant for which more than two individuals had germinated or from one that had not been allocated to the group of ten mother individuals. Because of uneven germination, each treatment ended up containing representatives of 7-12 mother plants, as sometimes there was zero or only one individual available from one mother plant. Therefore, the other replicate was randomly assigned a leftover individual from another mother. If there were not enough germinated seedlings from another mother individual, either, the replicate in question ended up with fewer individuals than the others. This was the case for the mean and warm treatments, where the germination rate was lower than in the cold treatment.

Altogether 517 plants were included in the experiment, with 82–92 individuals in each replicate (90 and 92 in the cold replicates; 85 and 86 in the mean replicates, and 82 in both warm replicates). We placed the plants on two benches in each replicate. They were initially watered automatically for three minutes every second day, but on November 14, we changed the watering regime to three minutes every three days, as the decline in sunshine towards the winter reduced evapotranspiration. On December 13, we changed the watering regime to 3 minutes on Mondays, Tuesdays, and Fridays, to ensure the plants would stay sufficiently moist until the next watering event. For the light spectrum of the High-Pressure Sodium lamps used in the greenhouses, see Guo et al. (2016).To avoid the effects of differing local conditions on the benches (for example, higher light and temperature directly underneath lamps or being closer to the windows where the temperature may differ slightly), we rotated the plants systematically 3 or 4 times a month. In addition to light from lamps, there was ambient daylight from windows. 

On November 22 and January 13, we fertilized the plants using Kekkilä Professional Superex NPK 12-5-27: 2kg/1000 L = 20g/10 L bucket, with 5 x 10 L buckets per bench. We submerged the pots in the fertilizer for 30 minutes, after which the benches were drained. On the same days, we applied nematodes to the soil, using 1 teaspoon of nematodes in 3 L of water. We used two species of entomopathogenic nematodes: Steinernema feltidae to control sciarid flies and Steinernema carpocapsae to control other flies. We applied them after the fertilization and automatic watering. We re-fertilized four additional times, once a week from January 30 to February 19, to induce flowering to be able to collect inflorescences for conducting pigmentation and reflection measurements at the end of the experiment. During this same period, we increased watering so that the plants were watered from the top using 5 L of water per treatment, in addition to the watering regime from underneath.

Motivation for measuring the chosen traits

Vegetative traits like the size of leaves and the number of leaves can tell us about the amount of resources that the plant is able to allocate towards biomass production. Vegetative traits have been shown to be rather plastic in P. lanceolata, even to the degree that this plasticity masks differences in local adaptation across large geographic extents (Villellas et al., 2021). Thus, we expect to see high plasticity in these traits for these high-latitude populations.

By measuring reproductive traits like flowering and number of flowers, we could assess fitness-related plasticity, i.e., to what extent is the reproductive output is affected by differing environmental conditions. P. lanceolata has been shown to express phenotypic plasticity in inflorescence reflectance and color (Lacey & Herr, 2005). The reflectance of inflorescences affects their internal temperature when there is incoming solar radiation (Anderson et al., 2013). Flower pigmentation may adjust reflectance and flower temperature since dark flowers absorb more solar energy than light, reflective flowers (Jewell et al., 1994; McKee & Richards, 1998). Pigmentation would be affected by, among other things, flavonoid and hydroxycinnamic acid (HCA) derivatives content, and less pigmented inflorescences should reflect more blue light.

Individuals and populations of P. lanceolata vary in their susceptibility to various fungal and viral pathogens (Norberg et al., 2023; Sallinen et al., 2023; Susi & Laine, 2015). Symptoms may be caused by a variety of pathogens, including fungi, but two of the most common viruses within the metapopulation on Åland include Plantago lanceolata latent virus (PlLV) and Plantago latent caulimovirus (PlCaV) (Susi et al., 2019). The prevalence of these two viruses varies across the local populations in Åland. It is not known whether temperature affects the susceptibility of the plants or the virulence of the pathogens, and thus whether viruses can differentially affect the populations and their symptoms when grown in different thermal conditions (Trebicki, 2020).

Trait measurements

To test thermal plasticity, we measured vegetative traits (leaf size and number), floral traits (flowering probability and abundance), floral bract pigmentation - content of both flavonoids and hydroxycinnamic acid (HCA) derivatives (including phenylethanoids), and floral blue-light reflection as a proxy of pigmentation, and the proportion of leaves showing symptoms of pathogen infection. We also tested for the presence of two common viruses: Plantago lanceolata latent virus (PlLV; Susi et al., 2017) and Plantago latent caulimovirus (PLCaV) and co-infections by the two viruses (Susi et al., 2019). The measured trait values across populations are given in Fig. S4 and separately for each population in Fig. S5 in ‘Supplementary Information.pdf’. See data ‘P.lanceolata_plasticity_data.csv’.

We took leaf size measurements during three consecutive days, starting 45 days after sowing began. We measured the length and width in cm of the leaf that, based on visual inspection, appeared to be the longest. We counted the number of leaves on each individual plant 12 weeks from sowing. Simultaneously, we assessed the degree to which plants showed symptoms of vertical pathogen infection (pathogens that likely had arrived with the seeds) by separately counting the number of leaves with typical signs of pathogen infection symptoms. We noted the number of leaves showing different kinds of pathogen symptoms (red, yellow, curly, necrosis, undefined). We used the proportion of leaves showing any symptoms as a pathogen response in our analyses. The information on different categories of symptoms is available in the accompanying data (‘P.lanceolata_plasticity_data.csv’). To measure infection rate by two common pathogens, we also took leaf samples from each individual for virus detection by cutting a 1 cm² piece from a leaf of each individual and placing the sample in a microtube. We froze and kept the samples in -20°C until DNA extraction. We extracted the DNA following Lodhi et al. (1994). See below for a description of how viruses PlLV and PlCaV were detected. The pathogen responses thus measure how the load of pathogens introduced in the experiment through the seeds affects infection rate and symptoms in plants grown at different temperatures. While these measurements allow us to capture the outcome of pathogen effects in the different thermal treatments, this does not allow us to distinguish any differences in resistance among the plants themselves or whether the pathogens differed in their virulence across the different temperatures tested.

19 weeks into the experiment, we began collecting inflorescences. We collected only fresh mature inflorescences, just before the stigma, petals, and anthers emerged from between the sepals. We collected inflorescences 1–3 times per week, based on their availability. By week 22, we had collected inflorescences from 72, 155, and 130 individuals in the Cold, Mean, and Warm treatments, respectively. These were used to measure floral bract pigmentation (flavonoids and HCA derivative content) and floral blue light reflection (see below for description of how this was conducted). For those individuals from which several inflorescences were collected, we averaged all subsequent measurements per individual prior to statistical analyses. After inflorescence collection for reflectance and pigmentation analyses was completed, we counted the total number of flowers by tallying the remaining inflorescences and previously cut flowering stalks. We summarized flower presence by truncating flower numbers >1 to 1. 

Measuring floral blue light reflection and floral bract pigmentation proxies

Right after we collected the inflorescence, we captured images in the near-infrared (NIR) region using an Olympus E-M1 camera and manual exposure with a Sigma 30mm F2.8 DN lens to visualize reflectance in the floral bracts. We converted the camera with DSLR AstroTEC to full spectrum using quartz glass to replace the sensor's built-in UVIR-cut filter. Lacey and Herr (2005) have reported that temperature influences the reflectance of P. lanceolata inflorescences from 725-850 nm in the NIR region, so we additionally used a 780 nm Heliopan RG780 ES 52 filter, which attenuates wavelengths below 780 nm. As a light source, we used a lamp with an incandescent bulb. Comparison of the NIR images allowed any differences in absorbance characteristics of flowers to be discerned. We processed the images using RawDigger: Research edition (Version 1.4.2). We evaluated blue channel sensor counts by selecting two 300x300 pixel samples from each flower at the tip and the widest section of the flower. We used a white Teflon slab with 95% reflectance as a white reference for each image. We used the following equation to calculate percentage reflection in the NIR region, which was subsequently normalized across samples:

*((sensor counts for flower) / (sensor counts for Teflon slab / 0.95)) 100

The inflorescences were stored for about 18 months in dry and dark conditions before having the opportunity to analyse phenolic compounds from their bracts. The extraction and analysis of phenolic compounds from the bracts of P. lanceolata followed the protocol of Stiles et al. (2007) with some modifications detailed here. Approximately 10 mg DW of bracts (previously dried at room temperature) from each sampled flower was soaked in 1 ml of 100% MeOH (CH₄O, ≥99.9%, Mr = 32.04 g·mol ⁻¹, Sigma-Aldrich, Schnelldorf, Germany) in an Eppendorf tube. This sample was shaken for 0.5 min using a vortex mixer and ultrasonicated (K-51E, Kraintek Czech s.r.o., Hradec Králové, Czech Republic) for 30 min at room temperature, then macerated for 24h at -22°C. The liquid part was collected using a pipette, and the remaining plant material was re-extracted in 1 ml of aqueous 50% (v/v) MeOH solution, but otherwise repeating the same protocol. The liquid fractions were combined (in 1:1 ratio, 0.5 ml each) and filtered through a 0.2 µm PTFE syringe filter.

The extracts of the bracts were subsequently analysed using the HPLC-DAD system (Agilent 1260 Infinity II, Agilent Technologies, Santa Clara, CA, USA) equipped with a Hypersil Gold chromatographic column (C18, 50 mm × 2.1 mm, 1.9 μm, Thermo Scientific, Waltham, NJ, USA). 3 µL of extract was injected into the system, the phenolic compounds were then separated in gradient of two mobile phases: A - water:COOH (995:5, v/v) B – acetonitrile:COOH (995:5, v/v) (C₂H₃N, ≥99.9%, Mr = 41.05 g·mol ⁻¹, Sigma-Aldrich, Schnelldorf, Germany; CH₂O₂, Mr = 46.03 g·mol ⁻¹, Sigma-Aldrich, Schnelldorf, Germany). Specification of the mobile phase gradient: t = 0-20 min A: 100 % -> 90 %; t = 20-35 min A: 90 % -> 85 %; t = 35-45 min A: 85 % -> 70 %; t = 45-50 min A: 70 % -> 0 %; t = 50-60 min A: 0 %; t = 60-70 min A: 0 % -> 100 % (+7 min. equilibration time). The flow was set as constant during the whole analysis (0.2 mL min^-1), likewise the temperature of the column (25°C). Absorption spectra of detected compounds were recorded in the range from 190-750 nm, chromatograms at 270, 314, 360, 440, 530, and 680 nm. Compounds were tentatively classified as flavonoids and hydroxycinnamic acid derivatives based on their absorption spectra and retention behaviour. For the relative quantification of phenolic compound content (mainly of flavonoids and hydroxycinnamic acid derivatives), the area under those peaks detected at 270 nm was used. The adjusted peak areas of compounds were adjusted by sample dry weight (DW) prior to further data processing (providing a measure of peak area mAU.s/ d.w. mg).

Detecting viruses

We detected PlLV and PlCaV using specific primers (for PlLV: PiLVi2_forward_1 5′-GTGTTTAACAATGAAGTGAGCC-3′C-3′ and PiLVi2_reverse_4 5′-AATCCATCCACACATCCAATC-3′; and for PlCaV forward primer 5′-AGGAGATGCCCATACTTTACC-3′ and reverse primer ′-GACTTGCCAGAACCTGATTTAC-3′) as in Sallinen et al. (2020). We ran PCR reactions to detect viruses in a final volume of 10 μL containing 1 μL of DNA, and GoTaq Green® polymerase 5x Mastermix (Promega Corporation, USA) according to the manufacturer’s instructions. We subjected samples to initial denaturation at 95 °C for 2 min, following 35 cycles of denaturation at 95 °C for 40 s, annealing at 60 °C for 40 s, and extension at 72 °C for 1 min with a final extension step of 72 °C for 5 min. The PCR product sizes were 117 bp for PlLV and 100 bp for PlCaV. We resolved the amplicons on a 1.2% agarose gel and visualized them using the Gel Doc XR System (Bio-Rad Laboratories, Inc., USA).

Statistical analyses

Our main aim was to assess, for the different traits, 1) the degree of thermal plasticity and 2) the origin of the majority of the variation in this plasticity, i.e., whether most variation in plasticity could be assigned to the population level or to the mother individual level. A mixed modelling approach allowed questions 1) and 2) to be assessed simultaneously, by choosing temperature treatment as a fixed effect and population and mother individual as nested random effects (Arnold et al., 2019). Here, the fixed effect describes the overall average response to the temperature treatment of the focal population, while the random effects describe how much variation is partitioned among populations and mother individuals, i.e., the contribution of variation in response attributable to interpopulation and intrapopulation differences, respectively. Using this approach, we were able to answer the main research questions concerning the contribution to variation in the measured traits by different local populations and mother individuals as opposed to the average value of all populations. Simultaneously, we accounted for non-independence between plants grown in the same greenhouse chamber by also using the greenhouse chamber as a random effect (Replicate). We evaluated the independence of the populations within each of the three main sampling areas by their connectivity values (Fig. 1 in Hällfors et al. in prep). For populations with relatively high connectivity values and situated in close proximity to each other (the region containing populations G, E, and F), we conducted a sensitivity analysis to test the consistency of the results. We left out one of the highly connected populations at a time and compared the model estimates, posterior probabilities, and standard deviations (Tables S2 and S3 in ‘Supplementary Information.pdf’). 

We used the brms package (version 2.15.0; Bürkner, 2021) in the R environment (R version 4.2.2; R Core team, 2023) to fit generalised and linear mixed-effect models (GLMMs and LMMs) in a Bayesian framework. We modelled leaf size, total flavonoid content, HCA content and blue light reflectance with a Gaussian error distribution, number of leaves with a negative binomial, flowering probability, PLCaV infection and PlLV infection with a Bernoulli error distribution, number of flowers and number of infections with a zero-inflated Poisson distribution, and the proportion of leaves with pathogen responses with a binomial error distribution. We estimated leaf size as the product of leaf length and leaf width to obtain a rough estimate in cm². These values were divided by 100 and square-root transformed, providing a better posterior predictive fit of the resulting model. See analysis code ‘P.lanceolata_plasticity_brms_analysis.R’.

For all models, we used informative prior distributions, which were needed to avoid computational issues. The informative prior distributions were justified for use here because the magnitude of the effects was expected to be small for all response variables. For the fixed effects and the random effects, we used Gaussian distributions with a mean of zero for both and a standard deviation of 2 and 1, respectively. However, for the models on the HCA derivative and flavonoid content, we set the standard deviation to 100 for the fixed effects and 50 for the random effects to account for the larger scale of these response variables. For the intercept, we used Student’s t as defined by the default settings in brms (Bürkner, 2021). For all models, we ran four chains with 4000 iterations, including a warm-up of 2000 iterations and a thinning rate of 1. Thus, we obtained 8000 posterior samples for each model.

For error distributions with an additional parameter to the location parameter, we considered a version of the model where the additional parameter was allowed to vary by treatment temperature. For response variables with a Gaussian error distribution (leaf size, flavonoid and HCA content, blue light reflectance), the additional parameter was the residual standard deviation (σ). For the number of flowers and number of infections modelled with a zero-inflated Poisson distribution, the additional parameter was the zero-inflation probability.  For the number of leaves, modelled using a negative binomial distribution, the additional parameter allowed to vary was shape. We compared the basic models to those models including this varying error distribution using the loo_compare function within the brms package, which compares the predictive performance of the models with leave-one-out cross-validation (Vehtari et al., 2017). Because of computational issues, the comparison of the number of leaves was conducted using the widely applicable information criterion (WAIC; Vehtari et al., 2017). If the additional model did not perform better (uncertainty intervals overlapped), we used the basic model for interpretation and result presentation.

We evaluated model convergence by investigating whether the Rhat values were <1.05 and the bulk effective sample size and tail effective sample size were each >400. We also visually inspected plots of posterior predictive checks using 10 posterior samples to identify potential discrepancies between the observed and predicted data.

In addition to interpreting the results based on mean posterior probability distribution and credible intervals, we used the hypothesis function in the brms package to interpret the models vis-à-vis the effect of temperature treatment. This function computes an evidence ratio for a one-sided hypothesis. In other words, we asked what the posterior probability is of the response being bigger (or in the case of phenolic absorbance where the direction of the effect was negative across treatments, smaller) than zero: a) for the mean temperature, b) between the mean and the warm temperature, and c) in the warm temperature.

We did a correlation analysis using the ggpairs function from the GGally package (Schloerke et al., 2024) on all measured traits, both across and within treatments, to evaluate the degree to which, e.g., high values in certain traits correlate with high values in others.

To estimate the effect on variance in thermal plasticity resulting from intrapopulation and interpopulation differences, we compared the estimated standard deviation attributable to the random effect of population and mother individual.

References

Anderson, E. R., Lovin, M. E., Richter, S. J., & Lacey, E. P. (2013). Multiple Plantago species (Plantaginaceae) modify floral reflectance and color in response to thermal change. American Journal of Botany, 100(12), 2485–2493. https://doi.org/10.3732/ajb.1300180

Arnold, P. A., Kruuk, L. E. B., & Nicotra, A. B. (2019). How to analyse plant phenotypic plasticity in response to a changing climate. New Phytologist, 222(3), 1235–1241. https://doi.org/10.1111/nph.15656

Bürkner, P.-C. (2021). Bayesian Item Response Modeling in R with brms and Stan. Journal of Statistical Software, 100(5), 1–54. https://doi.org/10.18637/jss.v100.i05

ENSCONET (2009) ENSCONET Seed Collecting Manual for Wild Species. Accessed in 2025 at: https://brahmsonline.kew.org/Content/Projects/msbp/resources/Training/ENSCONET_Collecting_protocol_English.pdf

González, A. V., Gómez-Silva, V., Ramírez, M. J., & Fontúrbel, F. E. (2020). Meta-analysis of the differential effects of habitat fragmentation and degradation on plant genetic diversity. Conservation Biology, 34(3), 711–720. https://doi.org/10.1111/cobi.13422

Guo, X., Hao, X., Zheng, J. M., Little, C., & Khosla, S. (2016). Response of greenhouse mini-cucumber to different vertical spectra of LED lighting under overhead high pressure sodium and plasma lighting. Acta Horticulturae, 1134, 87–94. https://doi.org/10.17660/ActaHortic.2016.1134.12

Hanski, I. (1999). Habitat Connectivity, Habitat Continuity, and Metapopulations in Dynamic Landscapes. Oikos, 87(2), 209–219. https://doi.org/10.2307/3546736

Jewell, J., McKEE, J., & Richards, A. J. (1994). The keel colour polymorphism in Lotus corniculatus L.: Differences in internal flower temperatures. New Phytologist, 128(2), 363–368. https://doi.org/10.1111/j.1469-8137.1994.tb04020.x

Jousimo, J., Tack, A. J. M., Ovaskainen, O., Mononen, T., Susi, H., Tollenaere, C., & Laine, A.-L. (2014). Ecological and evolutionary effects of fragmentation on infectious disease dynamics. Science, 344(6189), 1289–1293. https://doi.org/10.1126/science.1253621

Lacey, E. P., & Herr, D. (2005). Phenotypic plasticity, parental effects, and parental care in plants? I. An examination of spike reflectance in Plantago lanceolata (Plantaginaceae). American Journal of Botany, 92(6), 920–930. https://doi.org/10.3732/ajb.92.6.920

Lodhi, M. A., Ye, G.-N., Weeden, N. F., & Reisch, B. I. (1994). A simple and efficient method for DNA extraction from grapevine cultivars andVitis species. Plant Molecular Biology Reporter, 12(1), 6–13. https://doi.org/10.1007/BF02668658

Mckee, J., & Richards, A. J. (1998). Effect of flower structure and flower colour on intrafloral warming and pollen germination and pollen-tube growth in winter flowering Crocus L. (Iridaceae). Botanical Journal of the Linnean Society, 128(4), 369–384. https://doi.org/10.1111/j.1095-8339.1998.tb02127.x

Norberg, A., Susi, H., Sallinen, S., Baran, P., Clark, N. J., & Laine, A.-L. (2023). Direct and indirect viral associations predict coexistence in wild plant virus communities. Current Biology, 33(9), 1665-1676.e4. https://doi.org/10.1016/j.cub.2023.03.022

Ojanen, S. P., Nieminen, M., Meyke, E., Pöyry, J., & Hanski, I. (2013). Long-term metapopulation study of the Glanville fritillary butterfly (Melitaea cinxia): Survey methods, data management, and long-term population trends. Ecology and Evolution, 3(11), 3713–3737. https://doi.org/10.1002/ece3.733

R Core Team (2024). R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing

Sallinen, S., Norberg, A., Susi, H., & Laine, A.-L. (2020). Intraspecific host variation plays a key role in virus community assembly. Nature Communications, 11(1), 5610. https://doi.org/10.1038/s41467-020-19273-z

Sallinen, S., Susi, H., Halliday, F., & Laine, A.-L. (2023). Altered within- and between-host transmission under coinfection underpin parasite co-occurrence patterns in the wild. Evolutionary Ecology, 37(1), 131–151. https://doi.org/10.1007/s10682-022-10182-9

Schloerke, B., Cook, D., Larmarange, J., Briatte, F., Marbach, M., Thoen, E., Elberg, A., Toomet, O., Crowley, J., Hofmann, H., & Wickham, H. (2024). GGally: Extension to “ggplot2” (Version 2.2.1). https://cran.r-project.org/web/packages/GGally/index.html

Stiles, E. A., Cech, N. B., Dee, S. M., & Lacey, E. P. (2007). Temperature-sensitive anthocyanin production in flowers of Plantago lanceolata. Physiologia Plantarum, 129(4), 756–765. https://doi.org/10.1111/j.1399-3054.2007.00855.x

Susi, H., Filloux, D., Frilander, M. J., Roumagnac, P., & Laine, A.-L. (2019). Diverse and variable virus communities in wild plant populations revealed by metagenomic tools. PeerJ, 7, e6140. https://doi.org/10.7717/peerj.6140

Susi, H., & Laine, A.-L. (2015). The effectiveness and costs of pathogen resistance strategies in a perennial plant. Journal of Ecology, 103(2), 303–315. https://doi.org/10.1111/1365-2745.12373

Susi, H., Laine, A.-L., Filloux, D., Kraberger, S., Farkas, K., Bernardo, P., Frilander, M. J., Martin, D. P., Varsani, A., & Roumagnac, P. (2017). Genome sequences of a capulavirus infecting Plantago lanceolata in the Åland archipelago of Finland. Archives of Virology, 162(7), 2041–2045. https://doi.org/10.1007/s00705-017-3298-0

Trebicki, P. (2020). Climate change and plant virus epidemiology. Virus Research, 286, 198059. https://doi.org/10.1016/j.virusres.2020.198059

Thiers B. (2016). Index herbariorum: a global directory of public herbaria and associated staff. New York Botanical Garden’s Virtual Herbarium.

Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432. https://doi.org/10.1007/s11222-016-9696-4

Villellas, J., Ehrlén, J., Crone, E. E., Csergő, A. M., Garcia, M. B., Laine, A.-L., Roach, D. A., Salguero-Gómez, R., Wardle, G. M., Childs, D. Z., Elderd, B. D., Finn, A., Munné-Bosch, S., Bachelot, B., Bódis, J., Bucharova, A., Caruso, C. M., Catford, J. A., Coghill, M., … Buckley, Y. M. (2021). Phenotypic plasticity masks range-wide genetic differentiation for vegetative but not reproductive traits in a short-lived plant. Ecology Letters, 24(11), 2378–2393. https://doi.org/10.1111/ele.13858