Pollinator Proboscis Length Plays a Key Role in Floral Integration of Honeysuckle Flowers (Lonicera spp.)

Pollinator-mediated selection is supposed to influence floral integration. However, the potential pathway through which pollinators drive floral integration needs further investigations. We propose that pollinator proboscis length may play a key role in the evolution of floral integration. We first assessed the divergence of floral traits in 11 Lonicera species. Further, we detected the influence of pollinator proboscis length and eight floral traits on floral integration. We then used phylogenetic structural equation models (PSEMs) to illustrate the pathway through which pollinators drive the divergence of floral integration. Results of PCA indicated that species significantly differed in floral traits. Floral integration increased along with corolla tube length, stigma height, lip length, and the main pollinators’ proboscis length. PSEMs revealed a potential pathway by which pollinator proboscis length directly selected on corolla tube length and stigma height, while lip length co-varied with stigma height. Compared to species with short corolla tubes, long-tube flowers may experience more intense pollinator-mediated selection due to more specialized pollination systems and thus reduce variation in the floral traits. Along elongation of corolla tube and stigma height, the covariation of other relevant traits might help to maintain pollination success. The direct and indirect pollinator-mediation selection collectively enhances floral integration.


Introduction
Phenotypic integration refers to the coordinated variation of morphological traits within functional modules and results from the simultaneous occurrence of historical, physiological, developmental, and adaptive processes [1][2][3][4]. When phenotypic integration is the consequence of natural selection, acting on the functioning of those modules, it is called functional integration [2,5,6]. Disentangling the factors that influence phenotypic integration can enhance our understanding of the evolution of organs with convoluted modules.
Flowers are an ideal model system for studying the evolution of phenotypic integration. On the one hand, modules within a flower may be highly correlated due to genetic or developmental correlations [7,8]. On the other hand, flowers are functional structures, consisting of different modular units that play roles on different functions, such as pollinator attractiveness [9,10], pollinator accessibility [11], and pollen transfer (pollination efficiency) [12]. Pollinator-mediated selection may help to enhance floral integration [3,13,14]. However, natural selection usually acts on individual floral traits or on a

The Relationship between Floral Integration and Floral Traits and Pollinators' Proboscis Length
The magnitude of floral integration varied over 7.47-fold range; the value for L. webbiana (INT = 9.346) was the minimum, while L. tragophylla exhibited the maximum value (INT = 33.724; Table 1). Results of PGLS (phylogenetic generalized least squares) revealed that floral integration was positively related to corolla tube length (PGLS: model = OU, b = 0.417 ± 0.166, t = 2.513, p = 0.0332; Figure 2a Table S2.  Table S1). The composite pollinator proboscis length of the pollinator assemblage (a weighted proboscis length by considering all of the pollinators for each species, PLa) for L. maackii was the minimum, while L. tragophylla exhibited the longest value (Table S1). Results of PGLS indicated that floral integration was positively related to PLa (PGLS: model = PL, b = 0.486 ± 0.119, t = 4.075, p = 0.0028; Figure 2e).  The proboscis length of different pollinators of the study Lonicera species varied over a 16.8-fold range (from an average proboscis length of 2.45 mm in small carpenter bees to 41.15 mm in hawk moths; Table S1). The composite pollinator proboscis length of the pollinator assemblage (a weighted proboscis length by considering all of the pollinators for each species, PLa) for L. maackii was the minimum, while L. tragophylla exhibited the longest value (Table S1). Results of PGLS indicated that floral integration was positively related to PLa (PGLS: model = PL, b = 0.486 ± 0.119, t = 4.075, p = 0.0028; Figure 2e).

Discussion
In this study, the pollinator proboscis was found to influence the magnitude of floral integration in Lonicera plants. By PGLS, we found significant positive correlations between floral integration and pollinator proboscis length and four floral traits (corolla tube length, stigma height, upper lip length, and lower lip length). The best PSEM indicated that pollinator proboscis directly influenced the length of corolla tube (a trait that determine the accessibility to resources for pollinators) and stigma height (a trait with efficiency function), while upper lip length (a trait with attractiveness function) and lower lip length (attractiveness function) had a tight relationship with stigma height, which helped to maintain and to enhance the integration of honeysuckle flowers. The results, thus, illustrate a potential pathway by which pollinators may drive floral integration among closely related species.

Influence of Pollinator Proboscis Length on Floral Integration
Pollinator-mediated selection on floral integration has been reported in various flowering plants [3,4,11,15,[31][32][33]. In Lonicera, the magnitude of floral integration varied largely among species, with the two species pollinated by hawkmoths, showing higher floral integration than those pollinated by bees (Table 1). Honeysuckle flowers pollinated by hawkmoths have much longer corolla tubes than bee-pollinated flowers. The longtubed flowers may filter pollinators and form a more specialized pollination system than

Discussion
In this study, the pollinator proboscis was found to influence the magnitude of floral integration in Lonicera plants. By PGLS, we found significant positive correlations between floral integration and pollinator proboscis length and four floral traits (corolla tube length, stigma height, upper lip length, and lower lip length). The best PSEM indicated that pollinator proboscis directly influenced the length of corolla tube (a trait that determine the accessibility to resources for pollinators) and stigma height (a trait with efficiency function), while upper lip length (a trait with attractiveness function) and lower lip length (attractiveness function) had a tight relationship with stigma height, which helped to maintain and to enhance the integration of honeysuckle flowers. The results, thus, illustrate a potential pathway by which pollinators may drive floral integration among closely related species.

Influence of Pollinator Proboscis Length on Floral Integration
Pollinator-mediated selection on floral integration has been reported in various flowering plants [3,4,11,15,[31][32][33]. In Lonicera, the magnitude of floral integration varied largely among species, with the two species pollinated by hawkmoths, showing higher floral integration than those pollinated by bees (Table 1). Honeysuckle flowers pollinated by hawkmoths have much longer corolla tubes than bee-pollinated flowers. The long-tubed flowers may filter pollinators and form a more specialized pollination system than shorttubed flowers [34]. Several studies also indicate that plants with specialized pollination could lead to consistent selection on relevant floral traits and consequently, enhancing floral integration [12,35]. Based on investigations in contrasting populations of Lonicera implexa, Lázaro and Santamaría [11] showed that pollinator proboscis length could influence floral integration. Such evidence was also found in Ruellia humilis [31] and Narcissus papyraceus [16]. At a macroevolutionary scale, our results indicate that floral integration was positively correlated with the length of pollinator proboscis, confirming that pollinators could be a selection force that amplify floral integration of Lonicera.

The Potential Pathway by Which Pollinators Integrate Floral Traits
Although corolla tube length has been suggested to be evolutionarily correlated with pollinator proboscis [18,20,[36][37][38][39], it is hard to disentangle which trait (corolla tube length vs. stigma height) is the main selection target of pollinators in Lonicera, according to our current data. However, results also show that corolla tube length has a tight relationship with stigma height, suggesting it might also be a result of genetic/developmental correlation. Therefore, pollinators may exert strong selection pressures on a target floral trait, and the genetic/developmental correlations among floral traits may trigger the enhancement of floral integration [9].
Pollinators, foraging on the most fitting flowers, obtained the highest profit, and, meanwhile, they optimized the stigmatic pollen deposition on the flowers [40,41]. Floral traits, with an attractiveness function, may increase the probability of a flower to be visited, while those with accessibility function may influence pollinators' capability to access the resources of the flowers. Further, floral traits related to the pollination efficiency function may influence the pollen deposition on stigmas. A change in any of these floral traits may have an effect on the pollination success of the flower as a whole. For example, as the length of the corolla tube was selected by pollinators with different proboscis lengths, the position of stigma might, in flowers, also influence pollination success as pollinators change [33]. In this study, we found that both corolla tube length and stigma height were directly influenced by pollinators' proboscis length, increasing trait-matching between pollinators and flowers, which helps to enhance pollination success across species with different pollinator assemblages.
In addition, a honeysuckle flower with large stigma height always tends to have a long (upper and lower) corolla lip. This may help to enhance the attractiveness of the flower to its pollinators [11], but it may also be related to the position of the pollinator on the landing platform [42], which, in turn, may be influenced by how long the corolla tube and stigma are relative to the pollinator's proboscis length. The finding indicates that pollinator-mediated selection may also play an important role on the correlation pattern among floral traits [43].

Study Species and Sites
Field investigations were conducted at two sites: Taibai County, Shanxi Province (TB-SX) and Shennongjia Nature Reserve, Hubei Province (SNJ-HB) (see Table S1). The sites were sampled in 2016 and 2017. In these sites, some populations included two species (Table S1).

Measurements of Floral Traits and Estimation of Floral Integration
For each of the 11 Lonicera species, we conducted measurements of floral traits on one natural population (see Table S1 for population sizes), using a digital calliper with 0.01 mm precision (Mitutoyo, Kawasaki, Japan). Eight floral morphological traits were quantified for each of the 11 Lonicera species in the natural populations (each population included more than 100 individuals, Table S1), namely, corolla tube length (the distance between the base of the corolla and the top of floral receptacle), throat diameter (the diameter of the corolla tube opening), length of upper corolla lip (the distance between the base and the top), length of lower corolla lip (the distance from the base to the top), width of upper corolla lip (the largest distance from one side to the other), width of lower corolla lip (the largest distance between one side and the other), and the height of the stigma and anthers (the distance between the base of the corolla tube and the tip of the stigma or anthers, respectively). We measured at least 30 fully opened flowers (one flower/per plant) in the population of each species. The flowers were sampled at least 5 m away from each other to make sure they were from different individuals. Specimens of plants were deposited in the Herbarium of Wuhan University (WH).
We calculated floral integration for each studied Lonicera species based on the eight measured floral traits. Floral integration was calculated using the index of phenotypic integration (INT), which was defined as the variance of the eigenvalues (λi) of a correlation matrix through a principal component analysis (PCA) of the correlation matrix [44,45]. High variance among eigenvalues means that most traits are correlated and, thus, the first principal component (PC) accounts for most of the variation (high phenotypic integration). By contrast, low variance among eigenvalues indicates that the variation within the matrix is evenly distributed among all PCs (low phenotypic integration). PCAs were conducted through the 'PHENIX' package [46] in R to calculate INT for each species. INT values were corrected as INT = (Var (λi) − (number of traits − 1)/number of individuals per species) and expressed as a percentage of the possible maximum integration value. The possible maximum integration value equals the number of floral traits in the correlation matrix. The significance of differences between means and confidence intervals (CI) were also scored by bootstraping (n = 5000 permutations in each test).

Pollinator Observation and Measurement of Pollinators' Proboscis Length
In a previous pilot study, we recorded five insect guilds that visited the flowers of the 11 Lonicera species, namely, bumblebees, honeybees, leafcutter bees, small carpenter bees, and hawkmoths (Table S1). In the same natural populations where we measured floral traits (Table S1), we observed pollinator visitation in at least five inflorescences per species, using 20-min observation periods. A total of 30 observation periods were monitored for each of the species (10 h; Table S1). We conducted pollinator diurnal observations from 1000 h to 1500 h in three to five sunny days during the peak blooming period for each species. In addition to observation during the day, for L. japonica and L. tragophylla, pollinated by moths, we also conducted nocturnal pollinator observations from 1900 h to 2400 h [27,28,47], using a red light to reduce pollinator disturbance [11]. A visitor that encountered the anther or stigma of a flower was considered a pollinator. To better indicate the contribution of different pollinator guilds to the reproduction of each plant species, we calculated the visitation rate of each pollinator guild, namely, bumblebees, honeybees, leafcutter bees, small carpenter bees, and hawkmoths, in previously marked inflorescences. The visitation rate was defined as the total number of visits by pollinators from the same guild within an observation period, divided by the total number of open flowers in the inflorescence.
In addition, we sampled pollinators and kept them in 50-mL centrifuge tubes until the measurements of their proboscis length. All the pollinator species were identified by the Institute of Zoology at the Chinese Academy of Science. The length of pollinators proboscis was measured by a digital calliper (precision 0.01 mm) after relaxing the pollinators in a humid jar for two to three days; this allowed the tissues to soften and the proboscis to be pulled out using fine forceps. At least 20 individuals were measured for each pollinator guild, except for moths (Table S1). Morphologically specialized flowers may also be visited by multiple pollinators, rather than a single pollinator guild [48,49], and they, thus, may collectively exert selective pressures on floral traits [50]. By pollinator observations, we found pollinators with various proboscis lengths are also known to visit honeysuckle flowers. In this study, we developed a composite proboscis length by considering all of the pollinators for each Lonicera species. This composite length was calculated by the formula below: PLa is the composite proboscis length of the pollinator assemblage, n is the total number of pollinator guilds in a pollinator assemblage for each plant species, PLi is the proboscis length of the ith pollinator guild (the average proboscis length of all pollinators within the guild), VFi is the visitation rate of the ith pollinator guild, and VFn is the visitation rate of the pollinator assemblage. We used PCA to detect whether species differed among the 11 species in the 'PCAtools' package [51] of R. We performed a Tukey post hoc test by R package 'multcomp' [52] to determine whether species were statistically different among mean values for each trait after confirming that a difference between means in 11 species was supported by a one-way analysis of variance (ANOVA).

The Relationship of Floral Integration with Pollinator Proboscis Length and Floral Traits
We used PGLS to test whether floral integration was related to the floral traits and the composite pollinator proboscis length (PLa). PGLS generalizes the independent contrasts approach and can be used to incorporate a variety of models of evolutionary change [53,54]. We fitted the following evolutionary models: (1) Brownian motion (BM), whose traits evolve according to random drift; (2) Pagel's lambda (PL), which involves the rate of trait evolution that is optimized from the data; and (3) Ornstein-Uhlenbeck (OU), whose traits evolve towards an optimum (Table S2). Then, we compared model fit using AIC [55] and used the model with the best fit to estimate the relationship. In these analyses, we used the index of phenotypic integration (INT) as a response variable, and we used the composite proboscis length of the pollinator assemblage (PLa) and floral traits as predictor variables in separate models, containing one predictor variable (to avoid overparametrization). Based on the Bayesian tree of 11 Lonicera species [56], PGLS analyses were conducted in the R package 'nlme' [57].

The Potential Pathway by Which Pollinators Integrate Floral Traits
To disentangle the potential pathway by which the proboscis length of pollinators integrates floral traits, we applied phylogenetic structural equation modeling (PSEMs) by using the R package 'piecewiseSEM' [58]. The PSEMs comprised PGLS, using the gls function with the best evolutionary model in the package 'nlme', with Pagel's algorithm [59], to account for evolutionary dependence among species. By this method, we built different alternative models that combined PLa and the four floral traits significantly related to INT of the 11 Lonicera species and then compared them to find the best-fitting one (Table S3), which could help us to determine the key direct or indirect effects of PLa on variation and covariation in phenotypic traits. The four floral traits were classified into three groups according to functionality: accessibility (corolla tube length), efficiency (stigma height), and attractiveness (upper corolla lip length and lower corolla lip length), following previous studies [11,12]. In each model, PLa was thought to be a predictor variable that is directly or indirectly related to the floral traits from different functional modules. We created a total of 20 PSEMs, with every possible combination of pathways, due to all three conditions, wherein PLa was related to one, two, and three of the functional modules, respectively (Table S3). Floral traits would be correlated if they are constrained by the same genetic and/or developmental basis (e.g., anther height and stigma height) [36]. Because we could not presume the relationships between floral traits within attractiveness function (upper corolla lip length and lower corolla lip length) to be causal, they were defined as being correlated errors [58]. Shipley's test of d-separation [60,61] was used to assess the overall fit of each model, which assesses whether the model would be improved by the inclusion of identified missing paths. The d-separation test generates a Fisher's C (Fisher, 1925) test statistic, which can be used to assess overall fit of the PSEM and to calculate AIC for model selection [60,61]. The best model was defined as the model with the lowest AIC from those with p values greater than 0.05, derived from Fisher's C statistic test and with tests of d-separation, showing no statistically significant missing paths.
In this study, all data were summarized as the means ± standard errors, and all statistical tools were run in R with version 3.4.3 [62]. Except for the ANOVAs for which we applied a Bonferroni correction, the significance was considered to occur at a level of 0.05.

Conclusions
Our study illustrated a potential pathway by which pollinators integrate floral traits. Pollinator proboscis length directly influenced corolla tube length (accessibility function) and stigma height (efficiency function). In turn, covariation of stigma height with other floral traits enhances floral integration. The direct and indirect effects on floral traits might amplify the floral integration in response to a precise flower-pollinator fit. However, phylogenetic studies on the relationship between floral divergency and the pollination system, by incorporating more species, are necessary to improve our understanding of the evolution of floral integration in honeysuckle flowers.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/plants12081629/s1, Information on the study sites, pollinator information (pollinator species and guilds, pollinator proboscis length and visitation rate of each pollinator guild) and population sizes for the 11 studied Lonicera species in each natural population, Table S1. Comparison of evolutionary models in PGLS on association of floral integration (INT) with floral traits (corolla tube length, stigma height, upper lip length and lower lip length) and composite pollinator proboscis length (PLa), Table S2. Statistical results of the comparison of 20 phylogenetic structural equation models (PSEMs) based on AIC to determine the pathway that pollinators proboscis length mediate floral integration for the 11 Lonicera flowers, Table S3