Study sites and flower community sampling
The study was conducted on Yulong Snow Mountain (27°00′ N, 100°10′ E) in the Himalayan-Hengduan Mountains region, southwestern China. This region is recognized as one of the world’s biodiversity hotspots [36], and Yulong Mountain hosts a high diversity of sub-alpine and alpine flowering plants with more than 2815 species [37]. The vegetation type is characterized by pine forests at low elevations, Abies-Rhododendron forests at mid-elevations, and dwarf Rhododendron forests at high elevations. The region is known for its warm and rainy monsoon season, from May to October, and a colder and drier period marked by occasional snow storms from November to April. The average temperature and relative humidity in the two growing seasons from June to September (2019 and 2020) were 13.5 °C and 85.2%, respectively, with an average precipitation of 11.2 mm in 2019; precipitation data for 2020 was lacking. We selected five different sub-alpine and alpine meadow communities ranging from 2709 to 3896 m a.s.l., spaced approximately 200–500 m a.s.l. apart on the eastern slope of the mountain (Fig. 1, Table S1). The three lowest sites are located within the Lijiang Forest Biodiversity National Observation and Research Station, Kunming Institute of Botany, Chinese Academy of Sciences (CAS). All except the lowest elevation community are extensively grazed by cows, horses, and yaks at lower elevations during the growing season, whereas yaks are the sole grazers at high elevations. All data (including flower buds for quantifying pollen and ovule production and flowers for various floral trait measurements) were collected from June to August over three consecutive flowering seasons (2019 to 2021).
Map of the study area showing the five sub-alpine/alpine meadow communities along an elevation gradient from 2709 m a.s.l. (site 1) to 3894 m a.s.l. (site 5) on Yulong Mountain, southwestern China. Google Satellite Image exported from QGIS software
We collected data for 84 species belonging to 57 genera and 23 families across the five different communities. Since some species occur at multiple elevations, we recorded data regarding each species, considering each community independently with a total of N = 24 (site 1), N = 40 (site 2), N = 22 (site 3), N = 22 (site 4), and N = 22 (site 5). In this way, we treated each population separately for species with populations in multiple communities (i.e., spanning elevation). Our sampling represents the majority of the flowering plant species within the study area. Given that we expected pollinator availability along the elevation gradient to impact plant reproductive traits, we focused our sampling regime on insect-pollinated herbaceous species. Dominant families include Asteraceae (14 spp.), Lamiaceae (9 spp.), Ranunculaceae (7 spp.), Campanulaceae (6 spp.), and Fabaceae (5 spp.).
All plant species studied and collected in the study area are not listed as protected or endangered. All fieldwork and collections were permitted by the Lijiang Forest Ecosystem National Observation and Research Station and the Special Foundation for National Science and Technology Basic Research Program of China (2021FY100200). They followed the guidelines and legislation of the Kunming Institute of Botany, Chinese Academy of Science, Yunnan Government, and the Government of China, as well as the rules of the Convention on the Trade in Endangered Species of Wild Fauna and Flora (https://www.cites.org/). Voucher specimens of all the species (from Anemone rivularis: SN10078 to Viola biflora var. rockiana: SN10207 from the Yulong Snow Mountain; Table S2) were identified by Shristhi Nepal and deposited in the herbarium of the Kunming Institute of Botany (KUN), Chinese Academy of Sciences, Kunming, China.
Quantification of pollen and ovule production
To quantify pollen and ovule production per flower for each species, we collected at least ten mature flower buds, close to opening but anthers not dehisced, each from a different individual. All plant species sampled in this study have hermaphroditic flowers; therefore, we collected pollen samples and ovule data from the same flowers. All flower buds from each species were preserved in separate vials in 70% ethanol. We then randomly selected five of the collected flower buds per species and took one stamen per bud, each placed in a separate 2 mL centrifuge tube containing 1 mL distilled water. Each anther was gently ground using a glass rod to release the pollen into suspension. 20 µL of each suspension was placed on a Haemocytometer (Neubauer improved Haemocytometer, MARIENFELD, Tiefe-depth pro-founder 0.10 mm), and the number of pollen grains was counted under a light microscope at 10 × magnification [38]. The number of pollen grains in the 20 µL suspension was multiplied by the amount of suspension prepared (1 mL). Total pollen production per flower was calculated by multiplying the number of pollen grains in one anther by the number of anthers per flower. Finally, for each population, we calculated the mean ± SE for pollen counts across the five bud samples per population.
The ovule number per flower was counted from the same buds as anthers by dissecting the ovary and counting all ovules under a stereo-microscope [38]. We estimated the mean ± SE of ovule production by averaging across the five bud samples per population. In species of Asteraceae, we considered all the pollen and ovules present within one flower head as the number of pollen and ovules per flower following Arroyo et al. [3]. The ratio of pollen to ovule production was then calculated by dividing the number of pollen grains per flower by the number of ovules per flower. We also calculated the total pollen number, ovule number, and P/O per individual by multiplying each value by the average number of flowers per individual for each species (see below).
Quantification of floral traits
We measured inflorescence height and counted flower production (number of flowers per individual) for 15–30 individuals per species directly in the field. For each individual, we measured 2–3 fresh, fully opened flowers for a series of flower morphological traits, including floral display area (the product of the vertical and horizontal distance between the two tips of the corolla), tube depth (length from the base of the tube to the tube opening), stamen height (length from the base of the stamen and the tip), stigma height (length from the base of the stigma to the tip), and stigma-stamen separation (distance between the tip of stigma and the anther tip). We measured floral traits using either (1) digital measurement of photographs of the flowers (front and side views) with the digital imaging software ImageJ 1.38e (http://rsbweb.nih.gov/ij/) [39] or (2) direct measurement of freshly collected flowers in the lab using a digital calliper with 0.01 mm precision. Using the same two methods, we also recorded a series of floral categorical (i.e., qualitative) traits such as flower shape (open/fused), flower cluster (solitary/clustered), flower symmetry (actinomorphic or zygomorphic), and pollen presentation (open/enclosed).
Phylogeny construction
We constructed a phylogenetic tree using the function phylo.maker from the R package “V.Phylomaker” to determine the phylogenetic relationships among the studied plant species. We considered the GBOTB.extended mega phylogeny in Newick format as the backbone, default option “scenario 3” [40], where the tip for a new genus is bound to the upper 1/3 of the family branch (i.e., the branch between the family root node and the basal node). Polytomies in the final phylogenetic tree were resolved randomly using the function multi2di, available in the R package “ape” [41].
Data analysis
Based on Shapiro-Wilk’s normality test, all data were Log10 transformed to meet the assumptions of normality for all statistical analyses. To visualize the phenotypic trait differences across the different communities, we applied a non-metric multidimensional scaling (NMDS) using the function metaMDS from the Vegan package [42] to each species. The analyses were based on a Bray-Curtis distance matrix and were run for a maximum of 100 iterations. We considered the reduced-dimension representations of our data to be acceptable if NMDS stress scores were ≤ 0.2 [43, 44]. We used the package ggplot2 to generate NMDS plots with confidence ellipses for each community [45].
We ran a correlation analysis (CA, Table S3) and principal component analysis (PCA) to investigate the correlation between phenotypic traits across the different communities. To further determine whether variation in pollen and ovule number and P/O significantly co-vary with different trait expression of morphological floral traits and elevations, we performed both a non-corrected and a phylogenetically-corrected analysis using Ordinary Least Square (OLS) regression and Phylogenetic Generalized Least Square (PGLS) regression, respectively. The main aim of comparing OLS and PGLS regression was to test the importance of phylogenetic relationships among the studied species and traits. PGLS were estimated using Pagel’s lambda (λ) transformation and a Brownian model (λ = 1) in the r package CAPER 3.1.3 [46]. In our model, we first used pollen number per flower as the response variable and all traits individually as the source of variation. Secondly, we used the ovule number per flower as the response variable and all traits individually as the source of variation. Third, we set all the quantitative traits separately as the response variable and elevation (i.e., community) as the source of variation. We used p-values and AIC values from the analysed models to compare the non-corrected and phylogenetically corrected models [3, 47, 48]. To test the effects of the categorical traits on pollen number and ovule number per flower, we performed Phylogenetic ANOVA and ANOVA, respectively.
We estimated the phylogenetic signal for all the traits to determine if the variation in trait expression among species correlates with their phylogenetic relationships. We used Pagel’s λ and Blomberg’s K to estimate phylogenetic signals, i.e., non-random trait evolution under a Brownian motion model [49, 50]. As Pagel’s λ and Blomberg’s K differ in their methods for testing phylogenetic signals, we use both estimates to ensure that our interpretation of the patterns found in our plant community is accurate [51]. The objective was to detect possible inconsistent results between co-flowering species across different communities. Pagel’s λ is typically used for community studies and Blomberg’s K for closely related species [49, 50, 52]. For Pagel’s λ, 0 indicates no phylogenetic signal, while 1 indicates a strong phylogenetic signal. Similarly, for Blomberg’s K, 0 indicates little or no phylogenetic signal, and 1 indicates a strong phylogenetic signal. Pagel’s λ and Blomberg’s K were calculated in the R package GEIGER [33]. Their significance was tested against a null distribution generated by 1000 random permutations of the tips of the phylogeny using the Picante package [53]. The overall data analysis was conducted in R version 4.0.2 [54].
In addition to floral shape and presentation, we found a significant positive correlation between flower size and pollen and ovule number per flower. In general, in our communities, larger flowers tended to have higher male and female gamete production than smaller flowers. These results are consistent with a widely reported phenomenon; significant positive correlations between flower size and pollen number, and ovule number per flower [30, 58, 59]. In contrast, we found a significant negative relationship between flower production per individual and pollen number, and ovule number per flower. Collectively, these results suggest that to maximize fertilization success, plants in our sub-alpine and alpine communities may be allocating resources to a few but large flowers with high gamete production [60], although this relationship remains to be experimentally tested.
We found that pollen production and P/O per flower increased with increasing elevation only after accounting for phylogeny but did not affect the production of ovules per flower with or without phylogeny included in the analysis. This result is similar to previous studies, e.g., Cunha et al. [48], where increased pollen production across increasing latitude is more significant with phylogeny, which may indicate an evolutionary response to an unpredictable stochastic pollination environment. Given that pollen number, and P/O per flower, varied with elevation before and after taking phylogeny into account, variation in pollen production may be more indicative of a response to the environment rather than an evolutionary adaptation. Further, our findings that increased pollen production, and hence P/O, but no concomitant changes in ovule number with the increasing elevation may indicate that pollinator dependence and pollination efficiency have a greater effect on the evolution of pollen production compared to ovule production [61]. The pollinator guild of the Himalayan-Hengduan Mountains region is dominated by hymenopteran visitors at higher elevations [55, 62], whereas in other regions, bee abundance and richness decrease with increasing elevation and are replaced by flies as the dominant visitors [63]. Many of the herbaceous flowering plant species in our study sites are insect-pollinated [55, 62], and as such, pollinator dependence and pollination efficiency may play a large role in the trends for pollen production and P/O found in our study, although this source of variation remains to be empirically tested in the field.
Our study was limited to the sub-alpine and alpine communities across one elevation gradient (2709 m a.s.l. to 3896 m a.s.l.) without covering different environments across multiple mountain ranges. As such, it may explain why we did not find a strong significant relation between pollen production and elevation and a lack of support for concomitant increases in ovule production. This last finding, in particular, is in stark contrast to the ovule bet-hedging hypothesis by Arroyo et al. [3], which, according to Burd et al. [2], predicts the unpredictability of the pollinator environment may select for an increase in ovule number. They hypothesized that plants that package more ovules per flower might take advantage of rare pollination events through higher stigmatic pollen deposition, enabling them to produce more seeds than plants with fewer ovules per flower [28, 48]. Future studies should include a larger environmental gradient to more fully address this hypothesis by comparing plant traits within communities from lower elevation regions (sub-tropical/temperate) to those at high (sub-) alpine sites.
Finally, we found the number of flowers per individual had a strong negative relationship with elevation, as did the distance between the stigma and stamen tips. These findings support previous evidence that plants produce fewer flowers at high elevations (e.g. alpine regions) and shorten the distance between the stigma and stamen tips to ensure successful fertilization through self-pollination [3, 48]. Two findings from our study, an increased pollen number and P/O per flower in higher alpine species indicating a potentially higher outcrossing rate and a closer distance between stigma and stamen tips, may increase the potential for selfing. However, since we have not conducted bagging/hand-pollination to test the breeding systems for all the plants in this study, we cannot draw conclusive results to show the trends of breeding system change along the elevation. Here, at high-elevation sites like those found in our study, plants might reduce investment in flower production by producing structural components such as head stalks, involucres, and receptacles. In addition, they may have fewer flowers with higher pollen production [24], longer flower longevity [25] and increased stigma receptivity [71]. Although not empirically tested in this study, these strategies can directly contribute to successful fertilization and seed production and should be investigated further. We should also consider that in this region, cold-adapted social bumblebees are dominant pollinator groups in high elevations [62], and bumblebees play a key role in shaping plant reproductive strategy. Indeed, high-elevation populations of bumblebee-pollinated Incarvillea mairei also have high cross-pollination rates compared with lower-elevation populations [72].
Our research on sub-alpine and alpine meadow species suggests that high-elevation insect-pollinated communities produce more pollen per flower, resulting in higher P/O yet overall reductions in inflorescence height, flower production, tube height, and the distance between stigma tip to stamen tip. Furthermore, flower traits such as the floral display size and tube depth are highly correlated with pollen and ovule production but not with elevation, suggesting a plausible mechanism driving the pollination efficiency hypothesis. Ours is the first study to investigate pollen and ovule production, and P/O, as a function of elevation and floral traits for the majority of herbaceous species on Hengduan Mountains region, southwestern China. Until the present study, pollen and ovule production among different species spanning communities has received little attention. Until now, most pollen and ovule production studies are limited to species-level [69, 76, 77] or family-level [58, 59, 78] comparisons. However, it is important to note that our study did not assess plant fitness, which we expect to be relatively high in the sub-alpine compared to alpine plant communities via increased visitation and pollen export [14, 78, 79]. Future experimental studies in other mountain communities should include bagged and hand-pollination experiments to test additional factors governing patterns of plant reproductive character evolution along elevational gradients.
