Pollinators are declining globally, with climate change implicated as an important driver. Climate change can induce phenological shifts and reduce floral resources for pollinators, but little is known about its effects on floral attractiveness and how this might cascade to affect pollinators, pollination functions and plant fitness. We used an in situ long-term drought experiment to investigate multiple impacts of reduced precipitation in a natural Mediterranean shrubland, a habitat where climate change is predicted to increase the frequency and intensity of droughts. Focusing on three insect-pollinated plant species that provide abundant rewards and support a diversity of pollinators (Cistus albidus, Salvia rosmarinus and Thymus vulgaris), we investigated the effects of drought on a suite of floral traits including nectar production and floral scent. We also measured the impact of reduced rainfall on pollinator visits, fruit set and germination in S. rosmarinus and C. albidus. Drought altered floral emissions of all three plant species qualitatively, and reduced nectar production in T. vulgaris only. Apis mellifera and Bombus gr. terrestris visited more flowers in control plots than drought plots, while small wild bees visited more flowers in drought plots than control plots. Pollinator species richness did not differ significantly between treatments. Fruit set and seed set in S. rosmarinus and C. albidus did not differ significantly between control and drought plots, but seeds from drought plots had slower germination for S. rosmarinus and marginally lower germination success in C. albidus.

Synthesis. Overall, we found limited but consistent impacts of a moderate experimental drought on floral phenotype, plant reproduction and pollinator visits. Increased aridity under climate change is predicted to be stronger than the level assessed in the present study. Drought impacts will likely be stronger and this could profoundly affect the structure and functioning of plant-pollinator networks in Mediterranean ecosystems.

1. Study site: CLIMED long-term drought experiment

All field data were collected in February-June 2018. We used a subset of established plots that were part of the CLIMED (CLImate change effects on MEDiterranean biodiversity) long-term drought experiment situated at Massif de l’Étoile in Marseille, France (43° 22’ N, 5° 25’ E). This site has a typical woody shrub community dominated by three species: Quercus coccifera Linnaeus, 1753 (Fagaceae; anemophilous and a resource of very limited use to pollinators in the region; Ropars et al., 2020a), Salvia rosmarinus Spenn., 1835 (Lamiaceae; previously Rosmarinus officinalis; Drew et al., 2017), and Cistus albidus Linnaeus, 1753 (Cistaceae; Montès et al., 2008). Local cumulative precipitation between January and May 2018 (the flowering period surveyed) reached 291 mm, while the average precipitation between January and May for the period 2008-2018 was 205 mm (Marseille-Marignane meteorological station; www.infoclimat.fr). The site is equipped with 46 metallic control and 46 4 × 4 m rain-exclusion shelters established in October 2011, spaced by 1 to 30 m (Santonja et al., 2017). Plot locations were chosen randomly at the time of establishment of the long-term experiment, and were assigned at random to control or drought treatment (Montès et al., 2008). Gutters from rain-exclusion shelters in drought plots were designed to exclude up to 30 % and excluded on average (± SE) 12 ± 2% of precipitation between 2011 and 2018 at the centre of the plots; the intercepted water was carried away from the site with a pipe system. In control plots, the upside-down gutters intercepted a very small fraction of precipitation and rainfall reached the ground (Montès et al., 2008; Santonja et al., 2017). This water deficit attempts to mimic the mean predicted changes during the dry season in the Mediterranean area by the end of this century except in winter when rainfall is expected to increase (Giorgi & Lionello, 2008: averages for 2071-2100 relative to 1961-1990: December to February +0 to +10 %, March to May -10 to -20 %, June to August -20 to -30 %, September to November -0 to -10 %; Mariotti et al., 2015: averages for 2071-2098 relative to 1980-2005: December to February -0.1 to +0.2 mm/day, June to August -0.1 to -0.3 mm/day). The moderate but chronic experimental water deficit induced by the CLIMED experiment can alter plant physiology: carbon assimilation was reduced in C. albidus, and transpiration was reduced in C. albidus and S. rosmarinus but water use efficiency was not significantly changed in 2014 (Rodriguez-Ramirez, 2017). Between January and May 2018, permanent soil moisture probes (TDR100, Campbell Scientific Inc., Logan, Utah) measured soil moisture at 10, 20 and 40 cm in two control and two drought plots. For clarity we use the term drought to refer to the drought treatment in our study.

We selected 10 control plots and 10 drought plots out of the 92 plots, based on: (i) where Thymus vulgaris Linnaeus, 1753 (Lamiaceae) was present (four plots for each treatment only) because it is an important resource for pollinators (Ropars et al., 2020a); and (ii) a high and similar percentage cover of C. albidus and S. rosmarinus. The chosen control and drought plots were homogeneously distributed throughout the site. We measured the percentage cover of each species in selected plots twice (February and June 2018). The percentage cover of S. rosmarinus, C. albidus and Q. coccifera and T. vulgaris was 21, 19, 15 and 0.5 % on average respectively in the 20 plots selected, and the community composition did not differ significantly between treatments throughout the long-term experiment. Despite such low diversity, this plant community is natural, and is representative of the site and of the type of dense, closed vegetation plant communities found in the region in areas where wildfires are ancient (> 10 years; Pimont et al., 2018). Thymus vulgaris, C. albidus and S. rosmarinus are all perennial, entomogamous shrub species; T. vulgaris is gynodioecious and obligate entomogamous (dichogamous; Arnan et al., 2014), while S. rosmarinus and C. albidus are self-compatible but with limited self-pollination (Hammer & Junghanns, 2020; Blasco & Mateu, 1995). A fourth shrub species, Ulex parviflorus Pourr., 1788, was also present but very rare (0.3 % percentage cover) with very few flowers during the study period, and other flowering species were even rarer. We did not observe any insect visit to these very rare species and hence excluded them from our study.

2. Floral traits involved in pollinator attraction

2.1. Floral scent sampling and GC-MS analysis

We randomly selected up to 14 plant individuals per species in each treatment (control vs. drought) with a maximum of two (four for T. vulgaris) plants in the same plot. A few samples were lost during laboratory analysis, hence final sample sizes were 23 (control: 11; drought: 12) for S. rosmarinus, 22 (control: 11; drought: 11) for C. albidus, and 19 (control: 6 female, 6 hermaphroditic; drought: 5 female, 2 hermaphroditic) for T. vulgaris. Branches of the selected flowering plants bearing around 30-50, 2-3 or 100-400 flowers [1st-3rd quantiles] for S. rosmarinus, C. albidus and T. vulgaris respectively, were enclosed in a Nalophan bag (NA CAL, 30 cm × 30 cm, thickness 25 µm, volume ~ 2L; ETS Charles Frères, Saint-Étienne, France) connected to a pumping system maintaining a 1000 mL/min and a 200mL/min inlet and outlet air flows, respectively, provided by pumps (DC 12V, NMP850KNDC, KNF Neuberger SAS, France) powered by batteries (RS Pro 5Ah, 12V, RS Components SAS, France) and controlled by debit-metres (F65-SV1 Porter, Bronkhorst, France). Inlet air was first purified with activated charcoal (untreated, Mesh 4-8, Sigma Aldrich, USA) to limit the amount of volatiles from ambient air. Second, excess of humidity was removed using drierite (W.A. Hammond Drierite^TM Indicating Absorbents Mesh size 8, USA). Finally, ozone was filtered out through a fiberglass filter disk impregnated with sodium thiosulfate (Na₂S₂O₃) following Pollmann et al. (2005) to limit oxidation of plant volatile organic compounds (VOCs). Air flow was first stabilized for 15 min (the time required to entirely renew the air inside the 2L-bags). VOCs were then adsorbed on a cartridge placed at the bag outlet for 10 min for S. rosmarinus and T. vulgaris, and 15 min for C. albidus. This protocol optimizes the signal-to-threshold ratio without exceeding the breakthrough volume of each VOC in the conditions of our experiment, which would distort the estimated relative composition of chemical profiles (Ormeño et al., 2007). The cartridges were made of glass tubes (Gerstel OD 6 mm for TDS2/3, RIC SAS, Lyon, France) filled with 0.120 g Carbotrap® adsorbent (matrix Carbotrap® B, 20-40 mesh, Sigma-Aldrich, France) then 0.050 g Tenax® Porous Polymer Adsorbent (matrix Tenax® GR, 20-35 mesh, Sigma-Aldrich) separated by glass wool and maintained in the tube by a fixing screen (Gerstel for TDS 2 ID 4.0 mm, RIC SAS, France) at the entrance side and glass wool at the exit side.

To discriminate VOCs emitted by plants from possible environmental contamination, ambient air was sampled after every five plant samples using the same protocol. VOCs from four leaf-only plant samples per plant species were also measured to investigate which VOCs contribute most to floral scent versus leaf scent, enclosing branches of comparable size than the inflorescences of floral samples in collection bags. Sample cartridges were stored in a cooler immediately after collection, and transferred to a freezer at -20 °C as soon as possible. Prior to sampling, all cartridges had been cleaned in a Thermal Adsorbent Regenerator (RTA EcoLogicSense RG1301002, TERA Environment SARL, France) at 300 °C for 4 h. To reduce environmental variation from flowering phenology in scent emissions, each species was sampled over three days maximum around the flowering peak, during sunny weather and between 10:00 and 15:00. Throughout sampling, temperature and humidity were recorded inside and outside plant bags with data loggers (OM-EL-USB2-LCD, Omega Engineering Limited, UK). Plant parts inside bags were cut after sampling, dried in an oven at 50°C, and weighed after mass had stabilized (3-5 days, depending on plant species).

Samples were analysed one to 20 days after sampling. VOCs were thermodesorbed (cool trap and flash heating -50 to 250 °C at 12°C/s for 10 min; TDS 3 Gerstel equipped with an autosampler TDS A Gerstel). They were analysed with a gas chromatograph coupled with a quadrupole low-resolution mass spectrometer in solvent vent and CIS splitless mode (GC 6890N; MS 5973N; Agilent) equipped with a HP-5MS non-polar capillary column (5 % phenyl-methylsiloxane; length 30 m; internal diameter 0.25 mm; film thickness 0.25 μm; Agilent 19091S-433). The temperature gradient applied to the column was 40°C for 5 min, then up to 245 °C at 3 °C/min and maintained for 2 min (total run time 75.33 min). The carrier gas was helium at 7.1 psi and 1 mL/min. Mass spectra were recorded in the scan EMV mode (EM voltage 1295 eV and scanned from m/z 40 to 400, with one scan every 0.004 min.

Chromatograms were analysed with MZmine2 (version 2.18.1 developed for gas chromatography; Pluskal et al., 2010) in a 11-steps batch. Briefly, the baseline of each chromatogram was adjusted to 0, then values of m/z in each scan extracted and attributed to a peak. Peak heights (sum of m/z), areas and retention times, as well as mass spectrum at peak maximum were then exported from MZmine2 and imported into R (R Core Team, 2020; version 3.6.3) for peak identification. Retention indexes of each peak were calculated via a linear approximation built on a C5-C20 n-alkane series injected externally (Van Den Dool & Kratz, 1963). The retention indexes of 21 external standards were also verified with this method and checked against literature (Adams, 2007). The similarity between each peak’s mass spectrum and reference mass spectra was calculated using the R function ‘SpectrumSimilarity’ (R library ‘OrgMassSpecR v0.5-3’; Stein & Scott, 1994). The reference mass spectra were that of the 21 standards, and of two libraries converted in JPS format: the Adams 2007 library (Adams, 2007), and the NIST11 (NIST, 2011) library, was used when all similarity hits from the Adams 2007 library were lower than 0.7. Only similarity to reference molecules with a retention index within ± 15 of the peak’s calculated retention index was calculated. Identification was processed sequentially starting with the most common VOCs and retention indexes were adjusted locally (+/- 10) at each step based on these new identifications. Peaks with highest similarity < 0.6 were discarded. Peaks with the same identity were then aligned for each species. True absence of a VOC was verified, and areas not integrated with the first MZmine2 round were manually added, similarity calculated as above, with peak discarded if similarity < 0.6 and using a smaller retention index tolerance of ± 5 (0.2 min on RT). Only VOCs previously reported as known plant volatiles in Pherobase (El-Sayed, 2019) or in a comprehensive review (Knudsen et al., 2006). Most of them had also previously been reported in oil extracts or as plant volatiles from the three study species (Katerinopoulos et al., 2005; Ormeño et al., 2007; Satyal et al., 2016).

Following Campbell et al. (2019), we removed ambient air contaminants, by selecting only VOCs whose areas exceeded three times that of ambient air samples, by performing two-sample t-tests on , where A is the area of a peak (MZmine2 integration), Q and q are the inlet and outlet flows, respectively, and t is the sampling time, and using the type of sample (air vs. flower sample) as a factor, and with a false-discovery rate of 5 % to control for multiple comparisons. We also removed VOCs quantified in fewer than three samples of each species, because they added too much variance. We then removed the contribution of air samples to the retained VOCs in plant samples (each plant sample matched with the air sample taken closest in time) and calculated emission rates (in µg.h^-1.g_DM^-1): (Sabillón & Cremades, 2001), with k the response coefficient calculated for each chemical family based on external calibration of pure standards (see below), and m_DM the total dry mass of the inflorescence branch inside the sampling bag (flowers and leaves). Emission rates were then standardized by temperature (measured inside the bag during air flow stabilization and sampling; Ormeño et al., 2007; Sabillón & Cremades, 2001): where T is the temperature inside the sampling bag (in °K).

2.2. MZMine2 batch parameters

1. Raw data import

Select all raw chromatogram files associated with one species

2. Filter scans

Filter selected “Round resampling”

“Sum duplicate intensities” = True

“Remove zero intensity m/z peaks” = True

3. Crop filter

Retention time:

“min” = 4.5 min (remove solvents and water eluted in the first min)

“max” = 64.0 min (after Docosane ~ 311 da; upper limit for molecule volatility at ambient temperature is ~ 300 da).

m/z:

“min” = 40 (lower limit of mass spectrometer acquisition)

“max” = 315 (mass of Docosane molecular peak)

4. Baseline correction

Correction method: “RollingBall baseline corrector”

“wm (number of scans)” = 200

“ws (number of scans) = 5

5. Mass detection

mass detector selected: “Centroid”

“Noise level” = 2000 for 1^st round, 1000 for 2^nd round

6. Chromatogram builder

“min time span (min)” = 0.05

“min height” = (same as 5.)

“m/z tolerance” = 0.5 (absolute) / 0.001 (ppm)

7. Smoothing

“Filter width” = 15

8. Deconvolution

Algorithm selected : “Local minimum search”

“Chromatographic threshold” = 0.50

“Search minimum in RT range (min)” = 0.04

“Minimum relative height” = 0.001

“Minimum absolute height” = 2000

“Min ratio of peak top/edge” = 1.2

“Peak duration range (min)”: “min” = 0.05; “max” = 2.0

9. Peak merging

“m/z tolerance” = 500 (absolute) / 5.0 (ppm)

“RT tolerance window (number of scans)” = 15

“Use original raw data file” = False

“Use detected peaks only” = False

“Cumulative computing mode (TIC)” = True

10. Join Aligned (GC module)

Our aim here was to align as few peaks as possible, so as to have one line per peak in the final matrix

“m/z tolerance” = 0.5 (absolute) / 5.0 (ppm)

“Weight for m/z” = 0.8

“Retention time tolerance” = 0.001

“Weight for RT” = 0.2

“Minimum score” = 0.7

“Use RT recalibration” = False

“Use detected m/z only” = False

“RT tolerance post-recalibration” = 0.4

“Export dendrogram as TXT” = False

11. Export CSV

Export Peak RT, RTstart, RTend, Height, Area

2.3. External calibrations and calculation of response coefficients k

Three different mixtures of pure standard chemicals (Sigma-Aldrich) were made using cyclohexane as solvent. Cartridges were loaded with 1 µL of the dilutions, except for Mixture 3 in which they were loaded with 1 µL of the solid mixture and 1 µL of the E-Caryophyllene solution. Three samples of each dilution were analysed for each mixture, in the same conditions as plant samples (see Main text). Two linear regressions were calculated: at low injected masses, and high injected masses (DIL6 and DIL5B for Mixtures 1 and 2, and DIL4 and DIL3 for Mixture 3). This is because for high injected masses the column is overloaded and this leads to a smaller than expected peak area. An intercept was used in linear regressions, because areas at low injected masses were indistinguishable no matter the dilution, showing that the quantification threshold (= intercept) was higher than lowest dilutions tested. In the conditions of the analysis it was impossible to properly quantify exact mass of VOCs below that threshold.

2.4. Nectar production

For each plant individual on which we sampled floral scent we also selected five flowers (only 2.5 on average for C. albidus) to measure nectar production (as well as flower size and colour, see section 2.5.). We measured nectar standing crop (otherwise referred to as nectar production for simplicity) in each flower using 0.5 µL microcapillary tubes (Hirschmann Laborgeraete, Germany), placed in contact with nectaries. Flowers chosen were about to open for S. rosmarinus (any prior pollinator visit would have been prevented by the closed corolla) and sampled before 9:00 for C. albidus and T. vulgaris, to reduce variability due to flower ageing (3-4 days for S. rosmarinus and T. vulgaris and 1-2 days for C. albidus; authors pers. obs. and Flo et al., 2018) and nectar removal by pollinators (Flacher et al., 2020). Nectar volume was calculated (in μL/flower) by measuring the amount of liquid in a microcapillary tube with a digital calliper (Digit-Cal MK IV, Brown & Sharpe, USA). Sugar concentration was measured immediately afterwards using a hand-held refractometer (Eclipse 0-30°Brix, Bellingham & Stanley Ltd, UK) calibrated by the manufacturer. The whole volume of nectar collected was deposited on the measuring surface of the refractometer and diluted with 1 µL of pure water (measured with a 1 µL microcapillary tube) to increase the total volume and measure reproducibility. The sugar concentration (in µg/µL) and total sugar content (in µg) were calculated a posteriori using the conversion table from ° Brix to g/L by Kearns and Inouye (2003) and by applying the dilution factor and a temperature correction after the manufacturer’s correction table (Flacher et al., 2020).

2.5. Measure of flower colour and size, and statistical analysis

Flower colour. On the same flowers as used for nectar production, we recorded reflectance spectra from 300 to 700 nm at 45° using an optic fibre (FCR-7UV200-2 45°; Avantes, Apeldoorn, the Netherlands) connected to a Jaz portable spectrometer equipped with a Pulsed xenon light source (Ocean Optics, USA). Reflectance spectra were calibrated with a WS-2 white reflecting standard (Avantes, The Netherlands). Measurements were taken as such: for S. rosmarinus, in the middle of the lower petal; for C. albidus, in the middle of both the magenta and the yellow parts of two petals; and for T. vulgaris, on all the petals simultaneously, due to their very small size. Each reflectance spectrum was corrected by an offset calculated as the mean reflectance between 300 and 370 nm (the three species do not reflect in UV wavelengths) to correct for the bias due to the unavoidable light exposition variation during measurement. We calculated the position of each flower in the Apis mellifera hexagon colour space (Chittka, 1992), using the spectral sensitivity curves of its three photoreceptors (ultraviolet, blue and green; normalized each by their area under the curve) and the standard D65 as the illuminant spectrum: the position in the hexagon was calculated as the Euclidian distance between the stimulus and the background colour at the centre of the hexagon. For the background, we used the spectral reflectance of leaves of each species, which are the dominant flower background in the study site: we used the mean reflectance over five spectra for each species and each treatment, disregarding the natural variation in background colour. The hexagon representation illustrates how colours appear to a given visual system (Chittka, 1992).

Flower size. On each flower, we measured the distance between the base of the corolla tube end and the apex of the upper petal for S. rosmarinus. For C. albidus, we measured the distance between the centre of the flower and the apex of a petal, as well as the radius of the yellow centre on two petals per flower, and we calculated the corresponding means over the two values for each flower. For T. vulgaris, we measured both the length between the base of the corolla tube and the apex of the upper petals, as well as the largest diameter of the open petals.

2.6. Flowering phenology and number of flowers and other floral traits

When measuring plant-pollinator interactions, it is necessary to control for the flowering phenology of each species and the number of flowers at each sampling date, because pollinators are naturally attracted to patches with more flowers (Vrdoljak et al., 2016) and the flowering phenology of each species also influences the visitor community. To assess flowering phenology, the total number of flowers per plant was measured each week from March 9^th to May 23^rd, 2018 (12 weeks total), on two marked plant individuals per species (S. rosmarinus, C. albidus, and T. vulgaris) in each plot (20 individuals per treatment, except T. vulgaris: control: 11, drought: 8). On each individual, all flowers present in a volume delimited by the maximum height of the plant and a marked 40 × 40 cm quadrat were counted or estimated. For individuals smaller than the quadrat, their percent cover in the quadrat was measured, and all their flowers were counted, to obtain a number of flowers per m² for all plant individuals. When S. rosmarinus plants had many flowers around the flowering peak, flowers were counted by groups of ten, visually estimating the area covered by ten flowers.. For each individual, we calculated the date of its flowering peak as the week with the maximum number of flowers counted, and the total number of flowers as the sum of all flowers counted over all sampling dates. We also calculated a proxy for the flowering period as the number of weeks with at least half as many flowers as the maximum number of flowers recorded.

3. Plant-pollinator interactions

To quantify plant-pollinator interactions, pollinator visits were observed once a week for 12 weeks on the same days as flowering phenology, by monitoring for 5 min in each of the 20 plots: in each plot we chose the three patches (~1.5 m²) with the highest flower density to increase the chance of observing interactions (Ropars et al., 2020b). This represented a total of 20 hours of observation. Each patch was assigned to one observer, and all observed flowers of each flowering species were counted using the same method as for flowering phenology (see section 2.6.) to control for plot attractiveness to pollinators (patches with higher flower density are often more attractive and receive more visits; Feldman, 2006; Vrdoljak et al., 2016). Similarly to the plant individual level, the total number of observed flowers per plot and the flowering phenology at the plot level were not significantly different between control and drought plots (data not shown).

We defined a single visit when the observed insect touched the reproductive parts of a flower. We recorded the number of visits, the functional group of the visiting insect, and the plant species visited (S. rosmarinus, C. albidus or T. vulgaris). Functional groups recorded were: (1) Bombus gr. terrestris, (2) the managed honeybee Apis mellifera, (3) small wild bees (smaller than A. mellifera), (4) large wild bees (larger than A. mellifera), (5) Coleoptera, and (6) Diptera. Bombus gr. terrestris is by far the most abundant Bombus group species in the region and the other very rare species are functionally similar (Ropars et al., 2020a). No wild Apis mellifera colony was recorded close to our study site and we assumed that 100 % of A. mellifera we observed were managed (see also Herrera, 2020). For clarity we refer to Apis mellifera as the managed honeybee in this study. The six pollinator groups defined above are easily identified in situ. Visits by non-Apoidea Hymenoptera and by Lepidoptera were extremely rare and therefore not considered. Observations were carried out under optimal weather conditions, i.e., on sunny days with temperatures above 15° C and without strong wind, between 10:00 and 16:00 from 9^th March to 5^th April and from 8:00 to 14:00 from 13^th April to 23^rd May. Observations were recorded in the 20 plots (grouped in four blocks) successively in a block-randomized order at each sampling date. Contrary to our expectations and because T. vulgaris was rare (but always selected in the observed most flower-dense patches when flowering), insect visits to T. vulgaris were rare. In particular, all 85 recorded visits by A. mellifera on T. vulgaris took place in a single plot on a single day. We therefore discarded these visits from our analyses.

Immediately after the 5 min observation round in each plot, visiting insects were caught with a sweep net during the following 5 min, adding 30 s to that duration for each insect caught to account for insect handling time and disturbance to the patch. Insects were then frozen overnight and prepared for identification in the laboratory. They were pinned and dried prior to identification by professional taxonomists: David Genoud for Andrenidae, Anthophorinii, Colletes spp. and Halictidae, Matthieu Aubert for Megachilidae, Ceratinii and Hylaeus spp., Eric Dufrêne for Nomada and Mellecta spp., Christophe Lauriaut for Bombyliidae, Gabriel Nève for Syrphidae and other Diptera, and Jean-Yves Meunier and Jean-Pierre Hebrard for Coleoptera.

4. Plant reproduction

4.1. Fruit and seed set

As an estimate of plant female fitness we measured the fruit set of each marked plant individual of C. albidus and S. rosmarinus (two per plot). As a proxy for fruit set, we counted the total number of fruits per m² with the same method as for the total number of flowers on May 31^st for C. albidus and on four successive weeks (May 3^rd to May 23^rd) for S. rosmarinus. This was because first ripe fruits dehisced before the last ones were fully formed in S. rosmarinus; the maximum count was used as a proxy for fruit number. We also collected up to 40 and 13 fruits per individual for S. rosmarinus and C. albidus, respectively on May 31^st. This was the maximum number of intact (non-parasitized) fruits in C. albidus. We counted the number of seeds per fruit, and we calculated the mean and variance in number of seeds per fruit for each plant individual. Finally, we measured the mean seed mass, by weighing all collected seeds of a plant individual together (Ohaus® Discovery semi-micro analytical balance; resolution 0.1 mg).

4.2. Seed germination experiment

Seeds collected in May 2018 were conserved indoors at room temperature in a dry, dark environment. They were sown in Petri dishes on May 16^th 2019 in a greenhouse (controlled temperature 25 ± 2 °C; natural light) in the John Krebs Field station (University of Oxford, Oxford, UK). Both Salvia rosmarinus and Cistus albidus are obligate seeders (they regenerate after fire from seeds only; Céspedes et al., 2012), and C. albidus is pyrophile: fire is necessary to break physical dormancy in hard-coated seeds but has limited impact on germination of soft-coated S. rosmarinus seeds (Salvador & Lloret, 1995; Reyes & Trabaud, 2009; Chamorro et al., 2017). For this reason, half of the seeds of each species were flash-heated at 100 °C for 5 min in a drying oven, which was found to result in the highest germination rate of both species (Moreira et al., 2010). For each treatment (control vs. heated seeds) and each plant individual with enough seeds collected, up to 25 seeds were placed in a petri dish (diameter: 90 mm, height: 14 mm) on a filter paper (Whatman^TM, Grade 2) and watered with 1 mL of distilled water so that the whole filter paper was humid. In C. albidus, three boxes had between 14 and 22 seeds, and in S. rosmarinus 17 boxes had between 11 and 24 seeds, due to seed shortage. Petri dishes were sealed with Parafilm® sealing film. Every two days, every dish was watered with 1 mL of distilled water and the seal replaced. Every ten days, viable and germinated seeds were transferred in a new box with a new filter paper to avoid mould development. Rotted seeds were counted as unviable. Every two or three days from May 18^th to July 17^th, the number of seeds that had germinated in each box was recorded. The status of the seedlings (dead or alive) and their stage (two visible cotyledons) was also recorded.

We calculated seed viability as the number of viable seeds over the total number of seeds in each dish. The corrected germination rate and the survival rate were calculated as the total number of seeds that had germinated over the total number of viable seeds per dish, and the total number of seedlings that reached the cotyledon stage over the total number of seeds that had germinated in each dish and had enough time to develop (hence, seedlings that died during the experiment vs. seedlings that developed to the cotyledon stage). The average time to germination and the average time from germination to cotyledon stage were calculated as the mean time from the start of the experiment to the day germination was observed in germinated seeds in each dish, and the mean time from observation of germination to observation of cotyledon stage in each dish.

For each variable and each species independently, we analysed the effect of the maternal treatment (plants in control vs. drought conditions) and of the seed treatment (control vs. heated) using GLMERs. These two effects as well as the interaction between the two were implemented as fixed effects, and the identity of the mother plant was implemented as a random effect. A binomial distribution was chosen for rates (link function ‘logit’, function ‘glmer’, R library ‘lme4’; Bates et al., 2015), and the response variable was implemented as a two-column matrix (viability: number of viable seeds vs. unviable seeds; corrected germination rate: number of germinated seeds vs. number of viable seeds that did not germinate; survival rate: number of seedlings having survived to cotyledon stage vs. number of dead seedlings). A linear model was used for the time to germination and the time from germination to cotyledon stage. The significance of fixed effects and interactions was estimated through a step-wise regressive type-II comparison of models, based on χ² tests.

5. Statistical analyses

All statistics were performed using R version 3.6.3 (R Core Team 2020). We performed multivariate analyses of floral scent and pollinator community compositions (detailed below). Otherwise (unless specified), we used linear (LMMs) or generalised linear (GLMMs) mixed models with a poisson or a negative binomial error distribution function to account for data structure and overdispersion (frequent in count data), specified based on the distribution of model residuals. Random effects were the plot (repeated measures in the same plot), and the plant nested within the plot for plant traits (repeated measures on the same plant). The fixed effect was the treatment (control versus drought). For tests on plant traits involving T. vulgaris, the sex (female versus hermaphroditic) was also added as a cofactor to account for sexual dimorphism (Thompson et al., 2002; Arnan et al., 2014). The absence of residual heteroscedasticity and overdispersion was verified in the best model (functions ‘plot(simulateResiduals())’ and ‘testDispersion’; R library ‘DHARMa’; Hartig, 2019). The significance of fixed effects was estimated through a stepwise regressive type-II model comparison with an ANOVA based on χ² tests.

5.1. Floral scent

First, we analysed how drought altered floral scent chemical profiles of the three species independently. We performed a multivariate analysis on ⁴√-transformed emission rates of all quantified volatile organic compounds (VOCs). Such transformation reduces variance heterogeneity across VOCs spanning several orders of magnitude in natural plant chemical profiles (Hervé et al., 2018). Typical metabolomic datasets often comprise many more molecules than samples, and these molecules or explanatory variables are often strongly inter-correlated, notably due to shared metabolic pathways (Junker et al., 2017). The canonical powered partial least squares discriminant analysis overcomes these biases by combining classification and regression and is commonly used for the analysis of chemical datasets (function ‘cppls’ with parameters ‘centre’ and ‘scale’ to true, R library ‘pls’; Mevik et al., 2019; Indahl et al., 2009; Hervé et al., 2018). We implemented the treatment as factor and three components (or latent variables; further components always captured less than 5 % of the variance). The three resulting synthetic components hence reflect the maximized covariance between the drought response variable and the VOC predictor variables. To remove noise, i.e. VOCs with little contribution to covariance and group differentiation, we performed one round of variable selection, retaining only VOCs contributing most to group separation based on their relative projection to the axis separating group barycenters in the 3-dimensional space, and up to 2/3 of total axis contribution. To assess drought impact, we then performed a cross-validation significance test on the matrix of selected VOCs (Westerhuis et al., 2008), using the treatment as factor and three components (function ‘MVA.test’, R library ‘RVAideMemoire’; Hervé, 2019; model: “PLS-DA”, cmv=TRUE, and otherwise defaults parameters). In short and for each species separately, all samples were randomly divided into three sets of equal length and respecting the proportions of control and drought samples in the full dataset. The first was set aside as the test group (outer loop). The second and third sets were used for model training and validation (inner loop). The first set was then used to test the consensus model. Inner and outer loops were repeated in 999 permutations. The function outputs the classification error rate and a P-value computed using the Benjamini and Hochberg (1995) correction for multiple testing. For T. vulgaris flower dimorphism is also likely to be reflected in floral scent, so we added the sex (female or hermaphroditic) as an additional response in both the multivariate analysis and the cross-validation test (parameter ‘Y.add’).

Second, we analysed how drought affected the total emission rate (all VOCs summed per sample), and the total emission rates per chemical class for each plant species using two-sample t-tests on ⁴√-transformed data. Similarly, we looked at how drought affected the diversity of VOCs quantified in emissions, using t-tests on the number of VOCs quantified in each sample.

Finally, we analysed the impact of the sample type (flowers and leaves vs. leaves only) on the chemical profile to identify typically floral versus leafy VOCs in our experimental conditions. Similar to the analysis of drought impact on floral scent, we used a multivariate analysis built on the ‘cppls’ function, and one round of VOC selection. Because of the small sample size of leaf-only samples (N = 4), we could not perform a cross-validation test on the multivariate analysis. Instead, we performed a constrained correspondence analysis on the reduced matrix including selected VOCs only and a permutation test with 999 permutations (functions ‘cca’ and ‘anova.cca’, R library ‘vegan’; Oksanen et al., 2019). We also performed analyses on the total emission rate and diversity of VOCs emitted as for the drought impact on floral scent.

5.2. Other floral traits

Nectar. We calculated the proportion of flowers with nectar present as the number of flowers with nectar volume >0 over the total number of flowers sampled per individual and analysed it with a GLM with a binomial distribution (link function ‘logit’; function ‘glm’; R library ‘stat’; R Core Team, 2020). The response variable was implemented as a two-column matrix containing the number of flowers with and without nectar (Phillips et al., 2018). Using flowers with nectar only, we then analysed the drought impact on nectar volume (*100 and rounded, which corresponded to an estimated measure error of 2 %) and sugar content (µg, rounded) using GLMMs with poisson and negative binomial distributions, respectively (functions ‘glmer’ and ‘glmer.nb’, R library ‘lme4’; Bates et al., 2015).

Flower colour. We analysed the impact of drought on the polar coordinates r (radiant, spectral purity) and θ (angle, primary colour; Kantsa et al., 2017) using independent LMMs for each species, and for each colour for C. albidus (petal tip: magenta, flower centre: yellow; function ‘lmer’, R library ‘lme4’; Bates et al., 2015).

Flower size. Finally, we analysed drought impacts on flower size metrics for each species using LMMs (function ‘lmer’).

Flowering phenology and number of flowers. We analysed the impact of drought on flowering peak and the duration of the flowering period using two-sample Wilcoxon tests, because of the very small variance due to the per-week count. To analyse the total number of flowers counted throughout the season per plant per m², we used a GLMM with a negative binomial distribution.

5.3. Pollinator observations

We first analysed the drought impact on the total number of visits per plot throughout the season for each plant species using a GLMM with a poisson distribution (function ‘glmer’) and the treatment in interaction with the plant species (S. rosmarinus or C. albidus; N = 20 for each species, 10 per treatment) as fixed effects. Then, we analysed the impact of drought on the weekly number of visits by each pollinator functional group in each plot and to each plant species, using a two-step Hurdle modelling approach (Geslin et al., 2020). The first test worked on the presence or absence of visits of each group, assessing the attractiveness of a plot (how likely is each plant species in a plot likely to be visited by each functional group). We used a GLMM with a binomial distribution (function ‘glmer’), implementing the treatment in interaction with the plant species (S. rosmarinus versus C. albidus) and the functional group (A. mellifera, bumble bee, large wild bees, small wild bees, Coleoptera and Diptera) as fixed effects, along with the scaled observation duration and the scaled counted number of flowers as co-variables to control for differences in the number of flowers per plot. The second test worked on the number of visits, assessing the interaction strength between each plant species and each visitor functional group (how intensely are each plant species expected to be visited?). We selected data with visits only (removing all zeros), and we used the same fixed and random effects as in the first test, but with a negative binomial distribution (function ‘glmer.nb’). Since the interaction term was significant in this second test, we performed a post hoc comparison of means within functional groups and plant species and between treatments (function ‘emmeans’: ‘specs = pairwise ~ treatment | pollinator functional group * plant species’; R library ‘emmeans’; Lenth, 2019).

Finally, we analysed the drought impact on the pollinator community composition at a species level. We calculated the abundance and species richness of the total number of insect visitors caught in each plot through the season. Species identified to genus only were not counted in species richness (except one single specimen of Glyphipterix sp.), and such data was also used to estimate sampling completeness. Drought impact on the abundance and species richness were analysed using t-tests. We also performed a constrained correspondence analysis to test whether drought affected community composition, using a permutation test with 999 permutations and a Fisher test (functions ‘cca’ and ‘anova.cca’, R library ‘vegan’; Oksanen et al., 2019) on the same species used for species richness, but removing singletons.

5.4. Plant reproduction

The drought impact on the number of fruits per m² was tested with a GLMM with a negative binomial error distribution for each species. For each species, the impact of drought on the mean number of seeds per fruit, the variance in seed number per fruit (multiplied by ten and rounded), and the mean seed mass was tested using a LMM, a GLMM with a negative binomial error distribution, and a LMM, respectively.

References

• Adams, R. P. (2007). Identification of essential oil components by gas chromatography/ mass spectrometry, 4^th Edition. Allured Publ., Carol Stream, IL.
• Annamma Abraham, A., Verghese, A., & Muthangi, S. (2018). Role of colour and volatile in foraging behaviour of honeybee Apis cerana on Jacquemontia pentanthos. Journal of Asia-Pacific Entomology, 21, 1122-1128.
• Arnan, X., Escolà, A., Rodrigo, A., & Bosch, J. (2014). Female reproductive success in gynodioecious Thymus vulgaris: pollen versus nutrient limitation and pollinator foraging behaviour. Botanical Journal of the Linnean Society, 175, 395-408. http://dx.doi.org/10.1111/boj.12173.
• Bates, D., Maechler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67, 1-48.
• Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society. Series B, 57, 289-300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.
• Bhowmik, B., Lakare, S., Sen, A. & Bhadra, K. (2016). Olfactory stimulation of Apis cerana indica towards different doses of volatile constituents: SEM and EAG approaches. Journal of Asia-Pacific Entomology, 19, 847-859.
• Blasco, S., & Mateu, I. (1995). Flowering and fruiting phenology and breeding system of Cistus albidus L. Acta Botanica Gallica, 142, 245-251. https://doi.org/10.1080/12538078.1995.10515713.
• Campbell, D. R., Sosenski, P., & Raguso, R. A. (2019). Phenotypic plasticity of floral volatiles in response to increasing drought stress. Annals of Botany, 123, 601-610.
• Céspedes, B., Torres, I., Luna, B., Pérez, B., & Moreno, J. M. (2012). Soil seed bank, fire season, and temporal patterns of germination in a seeder-dominated Mediterranean shrubland. Plant Ecology, 213, 383-393.
• Ceuppens, B., Ameye, M., Langenhove, H. V., Roldan-Ruiz, I., & Smagghe, G. (2015). Characterization of volatiles in strawberry varieties ‘Elsanta’ and ‘Sonata’ and their effect on bumblebee flower visiting. Arthropod-Plant Interactions, 9, 281-287.
• Chamorro, D., Parra, A., & Moreno, J. (2016). Reproductive output, seed anatomy and germination under water stress in the seeder Cistus ladanifer subjected to experimental drought. Environmental and Experimental Botany, 123, 59-67.
• Chamorro, D., Luna, B., Ourcival, J.-M., Kavgacı, A., Sirca, C., Mouillot, F., Arianoutsou, M., & Moreno, J. M. (2017). Germination sensitivity to water stress in four shrubby species across the Mediterranean Basin. Plant Biology, 19, 23-31.
• Chamorro, D., & Moreno J. M. (2019). Effects of water stress and smoke on germination of Mediterranean shrubs with hard or soft coat seeds. Plant Ecology, 220, 511-521.
• Chao, A., Ma, K. H., Hsieh, T. C., & Chiu, C.-H. (2016). SpadeR: Species-Richness Prediction and Diversity Estimation with R. R package version 0.1.1. https://CRAN.R-project.org/package=SpadeR.
• Chittka, L. (1992). The colour hexagon: a chromaticity diagram based on photoreceptor excitations as a generalized representation of colour opponency. Journal of Comparative Physiology A, 170, 533-543.
• Dicke, M., & Baldwin, I. T. (2010). The evolutionary context for herbivore-induced plant volatiles: beyond the 'cry for help'. Trends in Plant Science, 15, 167-175.
• Dötterl, S., & Vereecken, N. J. (2010). The chemical ecology and evolution of bee-flower interactions: a review and perspectives. Canadian Journal of Zoology, 88, 668-697.
• Drew, B. T., Gonzalez-Gallegos, J., Xiang, C.-L., Kriebel, R., Drummond, C. P., Walker, J. B., & Sytsma, K. J. (2017). Salvia united: the greatest good for the greatest number. Taxon, 66, 133-145. https://doi.org/10.12705/661.7.
• El-Sayed, A. M. (2019). The Pherobase: Database of Pheromones and Semiochemicals. http://www.pherobase.com.
• Feldman, S. T. (2006). Pollinator aggregative and functional responses to flower density: does pollinator response to patches of plants accelerate at low-densities?. Oikos, 115, 128-140. https://doi.org/10.1111/j.2006.0030-1299.14493.x.
• Flacher, F., Raynaud, X., Hansart, A., Geslin, B., Motard, E., Verstraet, S., Bataille, M., & Dajoz, I. (2020). Below-ground competition alters attractiveness of an insect-pollinated plant to pollinators. AoB PLANTS, 12, https://doi.org/10.1093/aobpla/plaa022.
• Geslin, B., Gachet, S., Deschamps-Cottin, M., Flacher, F., Ignace, B., Knoploch, C., Meineri, É., Robles, C., Ropars, L., Schurr, L., & Le Féon, V. (2020). Bee hotels host a high abundance of exotic bees in an urban context. Acta Oecologica, 105, 103556. https://doi.org/10.1016/j.actao.2020.103556.
• Giorgi, F., & Lionello, P. (2008). Climate change projections for the Mediterranean region. Global and Planetary Change, 63, 90-104. https://doi.org/10.1016/j.gloplacha.2007.09.005.
• Hammer, M., & Junghanns, W. (2020). Rosmarinus officinalis L.: Rosemary. In J., Novak & W.D., Blüthner (Eds). Medicinal, Aromatic and Stimulant Plants (pp. 501-521). Springer International Publishing. https://doi.org/10.1007/978-3-030-38792-1_15.
• Hartig, F. (2019). DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.2.4. https://CRAN.R-project.org/package=DHARMa.
• Herrera, C. M. (2020). Gradual replacement of wild bees by honeybees in flowers of the Mediterranean Basin over the last 50 years. Proceedings of the Royal Society B: Biological Sciences, 287, 20192657. https://doi.org/10.1098/rspb.2019.2657.
• Hervé, M. R., Nicolè, F., &. Lê Cao, K.-A. (2018). Multivariate analysis of multiple datasets: a practical guide for chemical ecology. Journal of Chemical Ecology, 44, 215-234. https://doi.org/10.1007/s10886-018-0932-6.
• Hervé, M. (2019). RVAideMemoire: Testing and Plotting Procedures for Biostatistics. R package version 0.9-73. https://CRAN.R-project.org/package=RVAideMemoire.
• Indahl, U. G., Liland, K. H., & Næs, T. (2009). Canonical partial least squares – a unified PLS approach to classification and regression problems. Journal of Chemometrics, 23, 495-504. https://doi.org/10.1002/cem.1243.
• Junker, R. R., Kuppler, J., Amo, L., Blande, J. D., Borges, R. M., van Dam, N. M., Dicke, M., Dötterl, S., Ehlers, B. K., Etl, F., Gershenzon, J., Glinwood, R., Gols, R., Groot, A. T., Heil, M., Hoffmeister, M., Holopainen, J. K., Jarau, S., John, L., Kessler, A., Knudsen, J. T., Kost, C., Larue-Kontic, A.-A. C., Leonhardt, S. D., Lucas-Barbosa, D., Majetic, C. J., Menzel, F., Parachnowitsch, A. L., Pasquet, R. S., Poelman, E. H., Raguso, R. A., Ruther, J., Schiestl, F. P., Schmitt, T., Tholl, D., Unsicker, S. B., Verhulst, N., Visser, M. E., Weldegergis, B. T., & Köllner, T. G. (2017). Covariation and phenotypic integration in chemical communication displays: biosynthetic constraints and eco-evolutionary implications. New Phytologist, 220, 739-749. http://dx.doi.org/10.1111/nph.14505.
• Kantsa, A., Raguso, R. A., Dyer, A. G., Sgardelis, S. P., Olesen, J. M., & Petanidou, T. (2017). Community-wide integration of floral colour and scent in a Mediterranean scrubland. Nature Ecology & Evolution, 1, 1502-1510.
• Katerinopoulos, H. E., Pagona, G., Afratis, A., Stratigakis, N., & Roditakis, N. (2005). Composition and insect attracting activity of the essential oil of Rosmarinus officinalis. Journal of Chemical Ecology, 31, 111-122.
• Kearns, C. A., & Inouye, D. W. (1993). Techniques for pollination biologists. Niwot, CO: University Press of Colorado.
• Klatt, B. K., Burmeister, C., Westphal, C., Tscharntke, T., & von Fragstein, M. (2013). Flower volatiles, crop varieties and bee responses. PLoS ONE, 8, 1-7.
• Knudsen, J. T., Eriksson, R., Gershenzon, J., & Stahl, B. (2006). Diversity and distribution of floral scent. Botanical Review, 72, 1-120.
• Lenth, R. (2019). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.3.5. https://CRAN.R-project.org/package=emmeans.
• Leonhardt, S. D., Baumann, A.-M., Wallace, H. M., Brooks, P., & Schmitt, T. (2014). The chemistry of an unusual seed dispersal mutualism: bees use a complex set of olfactory cues to find their partner. Animal Behaviour, 98, 41-51.
• Mevik, B.-H., Wehrens, R., & Liland, K. H. (2019). pls: Partial Least Squares and Principal Component Regression. R package version 2.7-2. https://CRAN.R-project.org/package=pls.
• Montès, N., Maestre, F. T., Ballini, C., Baldy, V., Gauquelin, T., Planquette, M., Greff, S., Dupouyet, S., & Perret J.-B. (2008). On the relative importance of the effects of selection and complementarity as drivers of diversity-productivity relationships in Mediterranean shrublands. Oikos, 117, 1345-1350. https://doi.org/10.1111/j.0030-1299.2008.16910.x.
• Moreira, B., Tormo, J., Estrelles, E., & Pausas, J. G. (2010). Disentangling the role of heat and smoke as germination cues in Mediterranean Basin flora. Annals of Botany, 105, 627-635.
• Moreno, J. M., Zuazua, E., Pérez, B., Luna, B., Velasco, A., & Resco de Dios, V. (2011). Rainfall patterns after fire differentially affect the recruitment of three Mediterranean shrubs. Biogeosciences, 8, 3721-3732.
• Oksanen, J., Blanchet, F.G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O'Hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. M., Szoecs, E., & Wagner, H. (2019). vegan: Community Ecology Package. R package version 2.5-4. https://CRAN.R-project.org/package=vegan.
• Ormeño, E., Mévy, J., Vila, B., Bousquet-Mélou, A., Greff, S., Bonin, G., & Fernandez, C. (2007). Water deficit stress induces different monoterpene and sesquiterpene emission changes in Mediterranean species. Relationship between terpene emissions and plant water potential. Chemosphere, 67, 276-284.
• Pham-Delegue, M. H., Masson, C. Etievant, P., & Azar, M. (1986). Selective olfactory choices of the honeybee among sunflower aromas: A study by combined olfactory conditioning and chemical analysis. Journal of Chemical Ecology, 12, 781-793.
• Phillips, B. B., Shaw, R. F., Holland, M. J., Fry, E. L., Bardgett, R. D., Bullock, J .M., & Osborne, J. L. (2018). Drought reduces floral resources for pollinators. Global Change Biology, 24, 3226-3235. https://doi.org/10.1111/gcb.14130.
• Pimont, F., Dupuy, J. L. & Rigolot, E. (2018). A simple model for shrub-strata-fuel dynamics in Quercus coccifera L. communities. Annals of Forest Science, 75, 44. https://doi.org/10.1007/s13595-018-0713-y.
• Pluskal, T., Castillo, S., Villar-Briones, A., & Orešič, M. (2010). MZmine 2: Modular framework for processing, visualizing, and analyzing mass spectrometry-based molecular profile data. BMC Bioinformatics, 11, 395.
• Pollmann, J., Ortega, J., & Helmig, D. (2005). Analysis of atmospheric sesquiterpenes: Sampling losses and mitigation of ozone interferences. Environmental Science & Technology, 39, 9620-9629.
• Quintana, J. R., Cruz, A., Fernández-González, F., & Moreno, J. M. (2004). Time of germination and establishment success after fire of three obligate seeders in a Mediterranean shrubland of central Spain. Journal of Biogeography, 31, 241-249.
• R Core Team. (2020). R: A Language and Environment for Statistical Computing. Vienna, Austria, https://www.R-project.org/.
• Reyes, O., & Trabaud, L. (2009). Germination behaviour of 14 Mediterranean species in relation to fire factors: smoke and heat. Plant Ecology, 202, 113-121.
• Rodriguez-Ramirez, N. (2017). Aridification of the Mediterranean climate and biotic interactions : functional consequences on plant communities of a shrubland ecosystem. [Doctoral dissertation]. Aix-Marseille University, France. http://www.theses.fr/2017AIXM0467.
• Ropars, L., Affre, L., Aubert, M., Fernandez, C., Flacher, F., Genoud, D., Guiter, F., Jaworski, C., Lair, X., Mutillod, C., Nève, G., Schurr, L., & Geslin, B. (2020a). Pollinator specific richness and their interactions with local plant species: 10 years of sampling in Mediterranean Habitats. Environmental Entomology, 49, 947-955. https://doi.org/10.1093/ee/nvaa061.
• Ropars, L., Affre, L., Schurr, L., Flacher, F., Genoud, D., Mutillod, C., & Geslin, B. (2020b). Land cover composition, local plant community composition and honeybee colony density affect wild bee species assemblages in a Mediterranean biodiversity hot-spot. Acta Oecologica, 104, 103546. https://doi.org/10.1016/j.actao.2020.103546.
• Sabillón, D., & Cremades, L. V. (2001). Diurnal and seasonal variation of monoterpene emission rates for two typical Mediterranean species (Pinus pinea and Quercus ilex) from field measurements – relationship with temperature and PAR. Atmospheric Environment, 35, 4419-4431.
• Salvador, R., & Lloret, F. (1995). Germinación en el laboratorio de varias especies arbustivas mediterráneas: efecto de la temperatura. Orsis, 10, 25-34.
• Santonja, M., Rancon, A., Fromin, N., Baldy, V., Hättenschwiler, S., Fernandez, C., Montès, N., & Mirleau, P. (2017). Plant litter diversity increases microbial abundance, fungal diversity, and carbon and nitrogen cycling in a Mediterranean shrubland. Soil Biology and Biochemistry, 111, 124-134. https://doi.org/10.1016/j.soilbio.2017.04.006.
• Satyal, P., Murray, B. L., McFeeters, R. L., & Setzer, W. N. (2016). Essential oil characterization of Thymus vulgaris from various geographical locations. Foods, 5, 70.
• Saura-Mas, S., Saperas, A., & Lloret, F. (2019). Climatic and fire determinants of early life-history stages in the Mediterranean shrub Cistus albidus. Journal of Plant Ecology, 13, 3-11.
• Siles, L., Müller. M., Cela, J., Hernández, I., Alegre, L., & Munné-Bosch, S. (2017). Marked differences in seed dormancy in two populations of the Mediterranean shrub, Cistus albidus L. Plant Ecology & Diversity, 10, 231-240.
• Stein, S. E., & Scott, D. R. (1994). Optimization and testing of mass spectral library search algorithms for compound identification. Journal of the American Society for Mass Spectrometry, 5, 859-866.
• Thompson, J. D., Rolland, A.-G., & Prugnolle, F. (2002). Genetic variation for sexual dimorphism in flower size within and between populations of gynodioecious Thymus vulgaris. Journal of Evolutionary Biology, 15, 362-372. https://doi.org/10.1046/j.1420-9101.2002.00407.x.
• Van Den Dool, H., & Kratz, P. D. (1963). Journal of Chromatography, 11, 463-471.
• Vrdoljak, S. M., Samways, M. J., & Simaika, J. P. (2016). Pollinator conservation at the local scale: flower density, diversity and community structure increase flower visiting insect activity to mixed floral stands. Journal of Insect Conservation, 20, 711-721. https://doi.org/10.1007/s10841-016-9904-8.
• Westerhuis, J. A., Hoefsloot, H. C. J., Smit, S., Vis, D. J., Smilde, A. K., van Velzen, E. J. J., van Duijnhoven, J. P. M., & van Dorsten, F. A. (2008). Assessment of PLSDA cross validation. Metabolomics, 4, 81-89. https://doi.org/10.1007/s11306-007-0099-6.

MZMine: Pluskal, T., Castillo, S., Villar-Briones, A., & Orešič, M. (2010). MZmine 2: Modular framework for processing, visualizing, and analyzing mass spectrometry-based molecular profile data. BMC Bioinformatics, 11, 395.

R Core Team. (2020). R: A Language and Environment for Statistical Computing. Vienna, Austria, https://www.R-project.org/.