How much of morphological variability in pollen from genus Rubus L. might be explained by climate variability

The aim of this study was to assess the influence of 19 climatic factors during flowering periods and taxonomic proximity on the morphological features of the pollen from the genus Rubus L., which comprises numerous species, often with small ranges of natural occurrence.. It was hypothesized that the pollen morphology would be driven more by the effects of taxonomic proximity than by climatic variables, due to the conservatism of the pollen features, connected with a shared evolutionary history. The analyses revealed that climatic variability can explain an additional 2.5% to 14.0% of pollen morphology. The majority of the modelled pollen features were not correlated with the bioclimatic factors studied, except for the P/E ratio, which was positively correlated, and E, which was negatively correlated with PC3. However, most of the variability was explained by random effects connected with the taxonomic affiliation of the studied species to the genus Rubus L., which is very difficult in taxonomic terms. The study, therefore, showed how much additional interspecific variability in pollen morphology might be explained by the climatic variability of the species distributions.


INTRODUCTION
Environmental factors which affect plant growth can be classified as abiotic factors and biotic factors. The abiotic factors that affect plant characteristics include topography, soil, and climatic factors (light, temperature, moisture etc.). They are the non-living components of the environment. They also affect plant adaptation (Eyduran et al. 2015;Unlukara 2019;Marsic et al. 2019).
Pollen production and pollen morphology may be strongly constrained by environmental factors, including climate (Charlesworth et al. 1987;Murcia 1990;Delph et al. 1997;Walther et al. 2002;Rao et al. 2019). Nevertheless, knowledge of their impact on the morphological features of pollen is scarce (Ejsmond et al. 2011(Ejsmond et al. , 2015. The most important climatic factors affecting the growth and development of pollen are air temperature and humidity (Delph et al. 1997;Harder and Aizen 2010;Zinn et al. 2010;Ejsmond et al. 2011Ejsmond et al. , 2015Hinojosa et al. 2018). Furthermore, pollen grains are sensitive to abiotic stresses such as high temperature (Paupière 2014). Pollen reaction to heat stress during the flowering period was also observed by Prasad et al. (2011) and Omidi et al. (2014).
Several studies have confirmed that plants have to choose between the quantity and size of pollen grains produced (Mione and Anderson 1992;Vonhof and Harder 1995;Cruden 2000;Sarkissian and Harder 2001;Ashman et al. 2004;Yang and Guo 2004;Knight et al. 2005). However, empirical support of this compromise does not explain which ecological or functional factors determine the optimal combination of pollen size and quantity produced by a plant growing under given conditions (Ashman et al. 2004;Ejsmond et al. 2011Ejsmond et al. , 2015. Therefore, there is a need for a new trend in palyno-climatic research, the results of which may be helpful in solving the abovementioned research problems. Ejsmond et al. (2011) proved that pollen production is closely connected to environmental temperature through the optimization of the number and size of the pollen grains produced. In the opinion of these authors, temperature does not significantly affect pollen shape. According to Ejsmond et al. (2015), pollen size increases with temperature. Indeed, it is likely that the intensity of the pollen competition on stigma increases the optimal temperature of the flowering period, which in turn is expected to promote large pollen grains.
The genus Rubus L. is species rich and includes from 750 to more than 1000 species distributed worldwide (Weber 1995); The genus includes 108 species in Poland (Kosiński et al. 2018). Genus Rubus is highly complex and is one of the most taxonomically challenging genera of flowering plants (Robertson 1974;Ling-Ti 1983;Richards et al. 1996) and circumscription of the species is complicated by hybridization, polyploidy, agamospermy, and the lack of a universal species concept (Gustafsson 1943;Weber 1996;Zieliński 2004;Zieliński et al. 2004). The very large and growing number of Rubus species are resulted from the small and local geographic distribution of their natural occurrence. A recent species concept for European Rubus agamosperms, only allows as species those biotypes whose distribution exceeds an area of 50 km in diameter (Weber 1996).
In this study, for the first time, 11 morphological features of the pollen from 57 Rubus species were tested for their correlation with 19 climatic factors during flowering periods. The aim was to assess the influence of bioclimatic variables and taxonomic proximity on the morphological features of the pollen in the genus Rubus, which comprises numerous species, often with small ranges of natural occurrence. It was hypothesized that the pollen morphology would be driven more by of taxonomic proximity than by climatic variables, due to the conservatism of the pollen features, connected with a shared evolutionary history.

STUDIED TAXA
The study was conducted on 57 Rubus species, representing four out of five subgenera, all three sections and 22 series found in Poland, including six endemic species (R. capitulatus Utsch, R. chaerophylloides Sprib., R. ostroviensis Sprib., R. posnaniensis Sprib., R. seebergensis Pfuhl ex Sprib. and R. spribillei Sprib.) ( Table 1). The taxonomic classification of the studied taxa followed Zieliński (2004). The verification of the taxa was made by Prof. Jerzy Zieliński (from the Institute of Dendrology, the Polish Academy of Sciences in Kórnik), an outstanding taxonomist and specialist of the genus Rubus.

Pollen sampling and preparation
Several randomly-selected inf lorescences (f lowers) were collected from 57 natural bramble localities in Poland (Table 1). The plant material was stored in the herbarium of the Department of Forest Botany, Poznan University of Life Sciences (PZNF).
Acetolysis was carried out on the pollen grains according to the method used by Erdtman (1952Erdtman ( , 1960. The grains were mixed with the acetolysis solution, which consisted of nine parts acetic anhydrite and one part concentrated sulphuric acid. The mixture was then heated to boiling point and kept in a water bath for 2-3 min. The samples were centrifuged in the acetolysis mixture, washed with acetic acid and centrifuged again. The pollen grain samples were then mixed with 96% alcohol and centrifuged 4 times, with the processed grains subsequently divided into two groups. One half of the processed sample was immersed in an alcohol-based solution of glycerin for LM, while the other was placed in 96% ethyl alcohol in preparation for scanning electron microscopy (SEM). The SEM observations were made using a Zeiss Evo 40 and the LM measurements of the acetolysed pollen grains were taken using a Biolar 2308 microscope at a magnification of 640x. The pollen grains were immersed in glycerin jelly and measured using an ocular eyepiece with a scale. The measurement results were converted into micrometers by multiplying each measurement by two. Each sample consisted of 30 mature, randomly selected, properly developed pollen grains. In total, 1710 pollen grains were measured. This

Features analyzed
The pollen grains were analyzed for 11 quantitative characters: length of the polar axis (P) and equatorial diameter (E), length of the ectoaperture (Le), thickness of the exine along the polar axis and equatorial diameter (Exp and Exe), distance between apices of two ectocolpi (d) and P/E, Le/P, Exp/P, Exe/E, d/E (apocolpium index P.A.I) ratios.

Climatic data
In the analysis, 19 bioclimatic variables were used (Table 2), developed for species distribution models BIO-CLIM (Booth 2018;Booth et al. 2014)SDM is one of the most active areas of global ecology. Three books published in 2009, 2011 and 2017 have reviewed SDM, and the closely related areas of ecological niche modelling and habitat suitability modelling. All three books provide excellent introductions to these topics, but give very little information on the role that BIOCLIM played in laying the foundation for these research areas. Understanding the history of BIOCLIM is vital because it was the first package to implement the basic SDM process in an easy-to-use integrated system. It provided what are still the most commonly used set of 19 bioclimatic variables and contributed to the development of the interpolation routines used to prepare the most commonly used source of bioclimatic data (WorldClim. These variables were obtained from the WorldClim 1.4 database (Hijmans et al. 2005) using raster::getData() function in 2.5' resolution (~5 km in the study area). Due to intercorrelations between the variables, Principal Component Analysis (PCA) was used to use the main gradients of the bioclimatic variables in the reduced space (Fig. 1 a-c). Prior to PCA, the variables were scaled and centered to avoid artifacts connected with differences in ranges and units. PCA was performed using vegan::rda() function. Analysis of the inertia shared by particular principal components (screeplot, Fig. 1 d) revealed that the PC1-PC3 axes explained more variance than the null model of random variance distribution (broken stick model). Thus, all of them were used in further analyses. PC1 described the transition between the wetter and colder parts of the study area (a positive correlation with mean annual temperature) and the warmer and drier parts (a negative correlation with precipitation in both cold and warm periods). Therefore, PC1 was seen as representing aridity and temperature gradient. PC2 increased as the temperature of the driest and coldest month rose, as well as the precipitation in the coldest quarter, but decreased with the increasing temperature range, isothermality, precipitation in the wettest month and precipitation seasonality. Thus, PC2 was considered to represent aridity and climate variability. PC3 described seasonal variation in the climate, showing a positive correlation with temperature isothermality and a negative correlation with seasonality. Therefore, PC3 represented continentality -the seasonal variation gradient from a maritime to a continental climate (high PC3 values indicated a low seasonality of temperature). These three axes of PCA explained 53.96%, 22.24% and 15.97% of the variability, respectively. A lack of variance inflation in the models was also ensured by calculating the variance inflation factors.

Data analysis
To assess the studied relationships which acknowledged intraspecific variability and species-specific Precipitation of Driest Quarter mm bio18 Precipitation of Warmest Quarter mm bio19 Precipitation of Coldest Quarter mm effects, we used linear mixed-effects models (LMM) would be used. For data aggregated at species level, full models were used (eq. 1) Y = β 0 + ∑ n i=1 β i X i + u p:q:r:s + u q:r:s + u r:s + u s + ε p,q,r,s,j (1) u p:q:r:s~N (0,σ 2 ) u q:r:s~N (0,σ 2 ) u r:s~N (0,σ 2 ) u s~N (0,σ 2 ) ε p,q,r,s,j~N (0,σ 2 ) where Y = dependent variable (pollen morphological feature), X = predictors (particular climate PCA axes and their interactions), u p:q:r:s = random effects connected with series (nested in subsection, section and subgenus), u q:r:s = random effects connected with subsection (nested in section and subgenus), u r:s = random effects connected with section (nested in subgenus), u s = random effects connected with subgenus, and ε p,q,r,s,,j = residual error of particular samples.
A mixed-effects models lmerTest package (Bates et al. 2015;Kuznetsova et al. 2017) was developed. Firstly, a model compromising all three main bioclimatic components (PC1-PC3) was developed and then reduced according to Akaike Information Criterium corrected for small samples (AICc) using MuMIn::dredge() function (Bartoń 2017). From the list of candidate models the best fit was chosen, according to AICc and Akaike weights. In the case that the best fit was a null model (interceptonly), the second best model was used as the final model. Metrics (AICc and Akaike weights) were also provided for the null and full models in order to show the increase in information in the models. Information was provided on the amount of variance explained by the fixed effects using only the marginal coefficient of determination (R 2 m ) and by both the random and fixed effects using the conditional coefficient of determination (R 2 c ), calculated following Nakagawa and Schielzeth (2013), using the MuMIn::r.squaredGLMM() function (Bartoń 2017).
We assumed random effects connected with taxonomic nestedness as a proxy of interspecific variability, to compare with effect size of climatic variables. This made it possible to acknowledge phylogenetic non-independence in the data using random effects in the models. We ignored p-values as a measure of statistical significance since these can be biased by sample size or not connected with biologically meaningful effects (Wasserstein and Lazar 2016).

RESULTS
The two most important climatic factors -temperature and humidity -were analyzed based on 19 biocli-matic variables ( Table 2). Analyses of the mixed-effects models revealed that in the case of the majority of the modeled pollen features, the null model had the lowest AIC, which means that these traits were not correlated with the bioclimatic variables studied (Table 3). However, it was found that the P/E ratio was positively correlated and E was negatively correlated with PC3 (Table  4). Nevertheless, the estimates indicated low effect sizes. In the case of these morphological features of the pollen, the climatic data explained 14.0% and 2.8% of the variability, respectively. In contrast, the taxonomic affiliation, included in the models as a random effect, explained from 2.5% to 75.2% of the variability (Table 3).
Among the taxonomic random effects, the most important were these of the subgenus and series, while the least were those of the section (Table 4).

DISCUSSION
The results confirmed the hypothesis that pollen morphology is driven more by taxonomic proximity effects than by climatic variables, due to pollen feature conservatism associated with shared evolutionary history. The reproductive parts of plants (pollen grains and seeds) characters are more conservative and constant than their vegetative ones (Cruden 1977(Cruden , 2009). Therefore, pollen grains of related taxa usually have similar morphology, as it is the case also with pollen grains of the majority of the genera of the Rosaceae family (e.g. Crataegus, Malus, Rosa, Rubus, and Spiraea). They have isopolar monads, are generally medium-sized (rarely  Table 4).
small-sized), with tricolporate or tricolpate pollen grains and mostly striate exine ornamentation (Nazeri 2008;Polyakova and Gataulina 2008;Wrońska-Pilarek and Jagodziński 2011;Wrońska-Pilarek et al. 2013, 2019Lechowicz et al. 2020). Moreover, intra-generic studies indicated that pollen grains were so similar that usually it was possible to distinguish only a few sections or series and from a few to several individual species, and most often there were groups of taxa with similar pollen characteristics. In this research, among the random effects describing the phylogenetic relatedness of the species studied, the most important ranks were those of subgenus and series, while the least were those of section (Table 4). In contrast, the latest palynological study on Rubus (Lechowicz et al. 2020) revealed a low agreement between pollen morphological differentiation and taxonomic division. Pollen traits were most useful at the species level. In the case of the subgenus and series, it was observed that species belonging to these taxa did not generally form separate groups. Other genera of the Rosaceae family (e.g. Spiraea, Rosa, and Crataegus) showed a greater correlation between pollen morphology and infrageneric taxonomic classification (Wrońska-Pilarek and Jagodziński 2011; Wrońska-Pilarek et al. 2013, 2019. In Rubus this may be the result of apomixes, that is the replacement of the normal sexual reproduction by asexual reproduction, without fertilization, which could reduce natural variability (Weber 1996;Zieliński 2004). In older papers, Bell (1959) and Aizen and Raffaele (1998) indicated differences in pollen size due to fluc-tuations in the temperature occurring under the influence of different climatic conditions. According to Déri's (2011) modern studies on Cydonia oblonga, the pollen size and shape of this species were dependent on different climatic factors, such as temperature and humidity. Ejsmond et al. (2011), based on theoretical findings, showed a general trend for plants in environments with higher temperatures and potential evapotranspiration (PET) to produce less numerous and larger pollen grains (pollen with larger values of P and E) and to exhibit a slight change in grain shape, which may be more difficult to detect. Dainese and Sitzia (2013) as well as Maiti and Rodriguez (2015) and Azzazy (2016) also confirmed this opinion. Ejsmond et al. (2011) claimed that temperature does not significantly affect pollen shape. The research presented here did not fully confirm these results. Our study shows that pollen shape (P/E ratio) was positively correlated and equatorial diameter (E) was negatively correlated with PC3 (Table 4). This means that the less seasonal variability of the climate (higher PC3 values), the shorter the equatorial diameter of the pollen (E) and the larger the value of the P/E ratio , that is, the larger the share of elongated pollen grains. Lawrence and Campbell (1999) sampled 57 Rubus taxa including 20 species of subgenus Rubus, one to seven species from other 11 subgenera. Their genetic analyzes indicated that species from this genus were generally consistent with biogeography and ploidy, but traditionally important morphological characters, such as stem armature and leaf type, appeared to have a limited phylogenetic value in Rubus. This confirms SD -standard deviations of random effects: u p:q:r:s -random effects connected with series (nested in subsection, section and subgenus); u q:r:s -random effects connected with subsection (nested in section and subgenus); u r:s -random effects connected with section (nested in subgenus); u s -random effects connected with subgenus; ε p,q,r,s,,j, -residual error of particular samples.
the results of our earlier palynological studies (Lechowicz et al. 2020), which also showed that the morphological features of the pollen had considerable but limited impact on the taxonomy of genus Rubus. Lawrence and Campbell (1999) proved that ITS sequences were most informative among subgenera, and variability was low between closely related Rubus species. They distinguished three large clades in the genus Rubus. The first one contained all the sampled species of nine of the 12 studied subgenera, including subgenera Cylactis, Anoplobatus and Idaeobatus analyzed in this paper. The second clade included extreme Southern Hemisphere species of subgenera Comaropsis and Lampobatus, and the third consisted of subgenus Rubus (from which came 54 of the 57 examined species) and R. alpinus of subgenus Lampobatus. Such research results seem to confirm the hypothesis presented in this study that in Rubus pollen, the "impulse" caused by taxonomic and genetic factors is stronger than the influence of climatic factors. The cited authors showed the compatibility of the studied species with biogeography, which would indicate the great importance of the ranges of natural occurrence of the individual blackberry species for their diagnosis. The cited studies also indicate that the impact of geographical factors associated with climate factors on pollen morphology could perhaps be greater than demonstrated in this paper. It cannot be excluded that the results obtained in this study were influenced by the fact that the pollen grains (pollen samples) were collected from one natural site of a given blackberry species.
In contrast to palyno-climatic studies, research on the relationships between leaf traits and climatic data have been conducted more often. Wright et al. (2017) analyzed leaf data for 7670 plant species, along with climatic data from 682 sites worldwide. The authors provided a fully quantitative explanation for the latitudinal gradient in leaf size, with implications for plant ecology and physiology, vegetation modelling, and paleobotany. Large-leaved species predominate in wet, hot, sunny environments; small-leaved species typify hot, sunny environments only in arid conditions; small leaves are also found in high latitudes and at high elevations. By modelling the balance of leaf energy inputs and outputs, they showed that daytime and nighttime leaf-to-air temperature differences were key to geographic gradients in leaf size. Midolo et al. (2019) performed a global metaanalysis of leaf traits in 109 plant species located in four continents and demonstrated that there were common cross-species patterns of intraspecific leaf trait variation across elevational gradients worldwide. Irrespective of whether such variation is genetically determined via local adaptation or attributed to phenotypic plasticity, the leaf trait patterns quantified here suggest that plant species are adapted to living in a range of temperature conditions. The comprehensive studies cited above showed clear relationships between the morphological characteristics of leaves (e.g. leaf size) and climatic factors. It is therefore probable that in the case of a much larger pollen grain sample of different species from different geographical locations, the relationships between the pollen features and climatic factors would be much more pronounced.

CONCLUSIONS
The study revealed than climatic variability can explain an additional 2.5% to 14.0% of pollen morphology. However, most of the variability was explained by random effects connected with the taxonomic affiliation of the studied species to the genus Rubus L., which is very difficult in taxonomic terms.
Although the study analysed data concerning the interactions between the study site climate and the taxonomic affinity of the species on the pollen morphology, the results are biased due to the lack of species-specific replications. For this reason, phylogenic differences in the pollen morphology may have been masked by site specific effects. Despite this, as previously no attempts had been made to differentiate these two effects in the case of apomictic genus comprised of species often with small geographic ranges, it is assumed that these results might be a preliminary finding in the further exploration of such relationships. Further studies are required in order to determine whether knowing the biometric features of the pollen grains can significantly improve predictions of the impact of climate change on plant populations and to reconstruct past environmental conditions.