Structured and unstructured intraspecific propagule trait variation across environmental gradients in a widespread mangrove

Abstract Increasing studies have shown the importance of intraspecific trait variation (ITV) on ecological processes. However, the patterns and sources of ITV are still unclear, especially in the propagules of coastal vegetation. Here, we measured six hypocotyl traits for 66 genealogies of Kandelia obovata from 26 sites and analyzed how ITV in these traits was distributed across geography and genealogy through variance partitioning. We further constructed mixed models and structural equation models to disentangle the effects of climatic, oceanic, and maternal factors on ITV. Results showed that size‐related traits decreased along increasing latitudinal gradients, which was mainly driven by positive regulation of temperature on these traits. By contrast, ITV of shape trait was unstructured along latitudinal gradients and did not show any dependence among environmental variables. These findings indicate that propagule size mainly varied between populations, whereas propagule shape mainly varied between individuals. Our study may provide useful insights into the ITV in propagule from different functional dimensions and on a broad scale, which may facilitate mangrove protection in light of ITV.


| INTRODUC TI ON
Despite being recognized as a foundation for the theory of evolution by natural selection, the importance of intraspecific trait variation (ITV) has been neglected over time in ecology (Bolnick et al., 2011).Recently, there has been a resurgence of ecological interest in ITV, stimulated by the proliferating studies that underline the tremendous effects of ITV on community assembly and ecosystem functioning (Cardou et al., 2022;Violle et al., 2012).
However, the patterns and sources of ITV in themselves are still unclear (Cope et al., 2022).For example, ITV can be structured at various spatial scales in relation to different drivers.Betweenpopulation ITV tends to be shaped by large-scale environmental gradients, whereas within-population ITV is more likely the consequence of heritable differences and/or plastic responses to the local environment (Albert et al., 2010;Martin et al., 2017).
Alternatively, ITV can be unstructured if the trait is mainly determined by stochastic processes or constrained by fitness tradeoffs (Kendall & Fox, 2003;Moran et al., 2016).Elucidating ITV distribution across multiple scales is yet crucial for understanding and predicting ecological responses to global changes (Cochrane et al., 2015;Moran et al., 2016).
Trait-based plant studies to date have largely focused on vegetative traits (e.g., leaf and root functional traits), with less attention placed on seed traits despite the significance of seed traits for plant regeneration and long-term persistence (Saatkamp et al., 2019).
Knowledge of seed traits is particularly urgent for mangroves, which once covered over 200,000 km 2 of sheltered tropical and subtropical coastlines but nowadays has been disappearing worldwide for decades at an extremely high rate (Duke et al., 2007;Friess et al., 2019;Giri et al., 2015).Due to the property of low species richness and redundancy, mangrove degradation is always followed by pronounced losses of ecological multifunctionality (Arifanti et al., 2022;Donato et al., 2011;Feller et al., 2010).Such situation emphasizes mangrove recruitment and thereby calls for characterizing seed/propagule traits that represent key dimensions of the "regeneration niche" (sensu Grubb, 1977) in mangroves (Feller et al., 2010;Peterson & Bell, 2012).
Many mangroves (e.g., Rhizophoraceae and Avicenniaceae) have evolved a special reproductive strategy of viviparous seeds (germinating precociously while attached to the maternal plant), probably reflecting an adaptation to the salty and flooding intertidal environments (Feller et al., 2010).Due to the lack of seed dormancy, mangrove forests destroyed by severe disturbances such as hurricanes and tsunamis may not have local seed reserves, necessitating recolonization through long-distance seed dispersal from relatively undisturbed populations (Nettel & Dodd, 2007).This is particularly the case for Rhizophoraceae forests, as Rhizophoraceae species lack the capacity of resprouting from damaged trees as Avicenniaceae species (Baldwin et al., 2001).
Therefore, dispersal/retention is a key dimension of propagule functions where focal traits should be targeted for mangroves (Van der Stocken et al., 2019).Whether dispersed or retained, establishment is a prerequisite for any propagule to contribute to the regeneration, representing another critical dimension of seed trait (Krauss et al., 2008).
Most ITV studies in mangroves, though few in absolute number, have investigated propagule traits (e.g., size and weight) in relation to establishment, rarely considering the dispersal/retention dimension (Saenger & West, 2018;Yang et al., 2020;Zhu et al., 2021).
These studies showed that temperature and/or precipitation may have given rise to structured ITV in mangroves across geographic gradients (Saenger & West, 2018;Yang et al., 2020;Zhu et al., 2021), a pattern commonly reported in terrestrial forests (Kumordzi et al., 2019).Nevertheless, as mangroves are coastal vegetation, their trait variation may also be shaped by oceanic factors including salinity and tidal currents (Richards et al., 2021;Sousa et al., 2007).Additionally, maternal plants can also affect propagule traits, through both genetically fixed differences and environmentally induced transgenerational plasticity (Alam et al., 2018;Cochrane et al., 2015).However, the relative importance of these abiotic and biotic factors in shaping the distribution of ITV across multiple scales, and whether the patterns differ between the propagule function dimensions (establishment vs. dispersal/retention) are still poorly understood.
Here we use Kandelia obovata as the model species to investigate the issues.Kandelia obovata (Rhizophoraceae) is the most cold-tolerant true-mangrove species and has a wide latitudinal distribution along the southeast coast of China (Sheue et al., 2003;Wang et al., 2011), providing an ideal system for studying the ITV in mangroves.Using a stratified sampling design across a 9° latitudinal gradient, we analyzed the structure of intraspecific variability for five hypocotyl functional traits in relation to propagule dispersal/retention and establishment.The stratified sampling design and contrasting environmental conditions allowed us to address the following questions: (i) How is ITV structured spatially (between populations, maternal trees, hypocotyls)?(ii) What is the major driver (climatic, oceanic, or maternal factors) shaping the distribution of ITV? Based on the results from previous studies (Saenger & West, 2018;Zhu et al., 2021), the variability of traits on the establishment axis is expected to be higher between populations than within populations and predominantly shaped by abiotic factors.By contrast, the trait on the dispersal/retention axis, likely more reflecting a trade-off between post-disturbance recolonization and local recruitment (Sousa et al., 2007;Van der Stocken et al., 2019), may be less structured or even unstructured with regard to particular environmental gradients.

| Species, study sites, and sample collection
Kandelia obovata was first reported by Sheue et al. (2003), which was mainly different from the relative species K. candel (L.) Druce in leaf shape and cold tolerance.The two species are well-differentiated sets of geographical populations separated by the South China Sea.In China, K. obovata ranges from Hainan, Guangxi, Guangdong, Fujian, Zhejiang and Taiwan.Mature hypocotyls of K. obovata were collected during 2020-2021 from 26 sites along the coastline of southern China, spanning from 19°37′ N to 28°41′ N in latitude and from 108°05′ E to 121°24′ E in longitude (Figure 1).This geographic range provides wide climatic gradients in annual average precipitation (1102-2373 mm, mean 1578 mm), annual mean temperature (14.7-24.7°C,mean 21.1°C), and surface seawater salinity (13‰-35‰, mean 21.8‰).We sampled 1-6 maternal trees from each site (66 trees in total), depending on the population size.To ensure genetic independence, the sampled trees were spaced at least 30 m apart according to Geng et al. (2008).Each maternal tree was thus considered as a single genealogy.We randomly collected 30 hypocotyls from each tree for the study.

| Trait measurement for hypocotyls and maternal trees
We measured traits including fresh weight (FW), fresh length (FL), maximum transverse diameter (TD max ), and minimum transverse diameter (TD min ) for each hypocotyl.We also calculated the shape index of hypocotyl as the ratio of maximum transverse diameter to the minimum (RTD).This index is important for propagule floating because the slender morphology propagule disperses fast (De Ryck et al., 2012).We measured height (H) and diameter at breast height facets of how maternal trees can affect the progeny.We thus used AGB as a synthesis indicator to predict potential effects of maternal performance on hypocotyl traits.

| Access to climate and tide data
We extracted the data of each site from the WorldClim dataset (Hijmans et al., 2005) for mean annual temperature (MAT) and mean annual precipitation (MAP), and from Global Tidal Forecasting Service (http:// globa l-tide.nmdis.org.cn/ ) for tide and tidal current datasets.Specifically, the data on the highest tide, the lowest tide, the fastest tidal current, and the slowest tidal current were obtained from the nearest tide station.Since the annual average of the tide data is not available, we used the data of 4 months (January, April, July, and October) to represent the whole year.We calculated the annual highest tide, annual lowest tide, annual fastest tide current, and annual slowest tide current.We also measured sea surface salinity at each site using the Portable Refractometer (Apu, China).

| Data analysis
Bivariate analysis with ordinary least squares linear regression and quadratic regression were used to quantify how hypocotyl trait values varied with latitudinal gradients and biotic/abiotic factors.
Because of high correlations among most hypocotyl traits (r > .36,p < .001),we performed a principal component analysis with multiple traits using the "princomp" function in R 4.1.3(R Core Team, 2022) and used the two first principal component axes to represent the hypocotyl traits.
To evaluate how environmental factors, maternal plants, or inherent factors explained variation in hypocotyl traits, we used a nested analysis of variance coupled with variance partitioning techniques (Martin et al., 2017).We carried out linear mixed model for the first axis (PC1) and RTD by using the "lme" function in "nlme" R package (Pinheiro et al., 2022).In each model, all nested levels (i.e., site > genealogy > within [individual]) were entered as sequential random effects and the intercept was the only estimated fixed effect.We then used the "varcomp" function in the "ape" R package (Paradis et al., 2004) to calculate the variance components associated with each nested level.
To quantify how hypocotyl traits were affected by climatic factors, oceanic factors, and maternal performance, we implemented linear mixed models using the "lme" function in the R package "nlme."The fixed-effect terms included the climatic, oceanic, and maternal variables.To account for additional variation potentially caused by some missing site-specific effects (e.g., other environmental factors), and that caused by other maternal effects uncaptured by AGB, we treated sampling site and tree genealogy as random factors.All variables were standardized before the modeling such that each variable had a mean of zero and a standard deviation of one.To reduce the adverse influence of multicollinearity, we removed multicollinear variables until the variance inflation factors of all variables in the model were less than three (Ouyang et al., 2019).Both primary and quadratic mixed models were considered, and only the better-fitted model was shown (based on the Akaike information criteria).We calculated the variance inflation factor using the R package "car" (Fox & Monette, 2019).The pseudo-R 2 was calculated using the function "r.squaredGLMM" in the R package "MuMIn" (Bartoń, 2022) After standardizing all variables, multicollinear variables were removed based on variance inflation factor.We first considered a full model that included all variables and all reasonable pathways.Nonsignificant pathways were then sequentially removed, unless the pathways were biologically informative.The removing and adding of pathways were repeated until both p χ2-test ≥ .05(i.e., no significant difference between model predictions and the observed data) and root mean square error of approximation <.08 were reached (Wu et al., 2022).The structural equation modeling was performed using the "lavaan" R package (Rosseel, 2012).

| Patterns of ITV in hypocotyl along latitudinal gradients
The first two principal component axes (associated eigenvalues >1) accounted for 83.6% of the total variation among all these traits, where the PC1 mainly reflected hypocotyl mass and size (including fresh weight, fresh length, maximum transverse diameter, and minimum transverse diameter), and the second axis was primarily determined by the shape index (RTD; Figure 2).The PC1 decreased broadly with increasing latitudinal gradients (R 2 = .42,p < .001),yet statistically quadratic regression had a better explanation for their relationship (R 2 = .53,p < .001; Figure 3a).Intriguingly, RTD did not show any significant correlation with latitude (p = .6;Figure 3b).

| Sources of ITV in hypocotyl
Variance partitioning showed that the majority of variation in PC1 (63.75%) existed at the site level, with much less variation explained by genealogy (7.62%) or by individual (28.63%; Figure 4a).In contrast, the variability of hypocotyl shape did not appear to be structured either spatially or genealogically, as shown by the small part of RTD variation attributed to site (15.10%) or genealogy (15.87%).
Variability within individuals represented a large part of the variance in RTD (69.03%; Figure 4e).The linear mixed model indicated that 58% of the variation in PC1 could be explained by the predictive variables (quadratic mixed model).Specifically, the quadratic term of MAT (i.e., MAT 2 ) had the strongest negative effect on PC1 (standardized coefficient = −0.55 ± 0.12, p < .001),while MAT had a significant positive effect on PC1 (standardized coefficient = 0.38 ± 0.16, p < .001).
It is urgent to study the ITV in mangrove propagule at the background of global changes.We analyzed the relationship between biotic and abiotic factors and hypocotyl traits on a large scale for a typical viviparous mangrove species, Kandelia obovata.We found that hypocotyl size-related traits were mainly structured between populations, a pattern mainly shaped by climate.However, hypocotyl shape trait was almost unstructured either spatially or genealogically, with its variance predominantly attributed to the within-individual level.This pattern may suggest that variation in hypocotyl shape is independent of environmental variables.Our study provides insights into the ITV mechanism of hypocotyl traits in mangrove plants.

| Structured trait variation in mangrove hypocotyl
Previous studies have showed the presence of ITV in propagule size, especially for the genera Rhizophora, Bruguiera, Kandeli, and Ceriops (Tomlinson, 2016).Our results showed that variation in PC1 was significantly negatively correlated with latitude, which is similar to the results of a previous study (Yang et al., 2020).Nevertheless, an inverse relationship has been reported in a black mangrove (Avicennia marina), where small propagules were found in the northern (low latitude) populations of Australia yet larger propagules in the southern populations (Saenger & West, 2018).The opposite patterns of structured ITV are possibly due to the different strategies for improving fitness between true viviparous and cryptic viviparous mangrove species.Moreover, genetic analyses also suggested that differences evident in populations of A. marina at the different sites of its range are explained as evolutionary adaptations in response to the local physical, climatic, or biological characteristics (Maguire et al., 2000).Therefore, ITV in PC1 was influenced by the local biotic and abiotic factors at different sites.
We found that climate (especially temperature) was the main driver of intraspecific variation in PC1.In general, temperature and precipitation have a negative relation with latitude (Hijmans et al., 2005), while oceanic factors such as salinity and tides have little to do with latitude.It is straightforward and intuitive that more materials can be accumulated in hypocotyl under higher temperature.Precipitation had little effect on PC1, probably because mangrove plants have adapted to tidal environments and are able to obtain fresh water from seawater (Parida & Jha, 2010).
Although oceanic factors were not found to be the main factors affecting the variation in hypocotyl traits, salinity did have a significant effect on PC1 variation.Salinity has long been considered an important factor limiting the growth and distribution of mangroves (Chen & Ye, 2014;Richards et al., 2021), and negative correlations between the lengths and weights of mangrove propagules and sea salinity have been reported (Alam et al., 2018).In addition, salt stress may also be a factor that shapes the evolution of vivipary in mangrove plants.The hypocotyl of viviparous seeds contains a large amount of nutrients and water, which are necessary for the growth of newly colonized seedlings, thus avoiding the vulnerable life stage (Joshi, 1933;Zhou et al., 2016).When salinity is too high, maternal tree will invest more energy to against salt stress (Alam et al., 2018;Parida et al., 2004), which may reduce the investment on the hypocotyl and eventually lead to variation in its traits.Tide is an important factor in mangrove growth and hypocotyl dispersal process (Clarke et al., 2001;Duke et al., 1998;Van der Stocken et al., 2015;Zhang et al., 2021).However, the effect of tides on the variation in hypocotyl traits was marginal in our study.A possible reason is that previous studies have only considered the influence of tidal effects, rather than the combination of multi-biotic and abiotic effects.

| Unstructured trait variation in mangrove hypocotyl
It is interesting that the shape index RTD did not vary with latitude and was much less affected by abiotic and biotic factors than sizerelated traits.Few previous studies have focused on the variation in propagule shape at large scales.This lack may be due to the limited distribution areas of other true viviparous mangroves, such as Bruguiera sexangular, B. gymnorrhiza, and Rhizophora stylosa (Wu et al., 2018).Nevertheless, if global warming continues and these thermophilic mangroves expand to higher latitudes, we may also expect unstructured RTD along latitude gradients for their hypocotyls, because they share very similar reproductive strategies as K. obovata.By contrast, cryptic viviparous mangroves (e.g.Avicennia marina, Aegiceras corniculatum) may not show such a pattern, because their propagules are much smaller and even have no elongated hypocotyl (Zhang et al., 2021).
Many studies have emphasized the importance of maternal effects on the variation in offspring traits (Alam et al., 2018;Cochrane et al., 2015;Galloway, 2005), but the strength of maternal effect on ITV was minor in this study.As we only recorded the height and DBH of maternal plants, which reflects the size and age may have less to do with morphology than with yield (Hangelbroek & Santamaria, 2004).Rather, heredity may be the major contributor to maternal effects.In general, genetic diversity is low in natural mangrove populations (Hodel et al., 2018;Richards et al., 2021), but high epigenetic diversity may maintain the high plasticity of traits in mangrove (Mounger et al., 2021).This may explain why the ITV in RTD was mainly attributed to the within-individual level.Additionally, when traits are related to plant fitness trade-offs, trait variation may be independent among environmental variables (Bonte et al., 2012).For example, variation in reproductive traits of conifer tree was mainly explained at the scale of individual trees in the western United States (Vasey et al., 2022).This hypothesis needs to be further studied in coastal vegetation.
Due to the lack of dormant seed periods, true mangrove species do not have seed banks and are difficult to regenerate after anthropogenic disturbances or natural disasters (Nettel & Dodd, 2007).Therefore, the collection of germplasm resources High ITV represents high trait diversity and plant may also have higher adaptive capacity (Moran et al., 2016).But if RTD determines some important ecological function (e.g., dispersal), then the low ITV may not allow the mangroves adapt to global change quickly, which in turn affects the mangrove's natural regeneration and persistence.Thus, we recommend that the ITV in hypocotyl should be considered in germplasm collection and restoration so as to facilitate the regeneration of mangrove populations.Future studies that combine biotic and abiotic factors will also help us to have a comprehensive understanding of the trait variation in mangroves and to better protect mangroves.

(
DBH) for maternal trees, and transformed them into aboveground biomass (AGB) using an allometric model of tropical wood tree: AGB = 0.0673 × (ρ × DBH 2 × H) 0.976 (Chave et al., 2014), where wood density ρ is set as 0.57 g/cm 3 (Jiang, 2014).The AGB is a comprehensive outcome of the interaction between the genotype and the environment experienced by maternal trees, thereby capturing multiple F I G U R E 1 Sampling sites of Kandelia obovata hypocotyl on the south coast of China.The color gradient indicates the salinity (‰) of surface seawater.Circle size is proportional to sampling size of maternal trees (N).
, to represent the variance explained by the fixed effect in the linear mixed model.The effect sizes of fixed factors were measured by the regression coefficients in the linear mixed model.Structural equation modeling was used to disentangle direct and indirect effects of all predictive factors on hypocotyl traits.
).In contrast, all these predictive factors together explained only 12% of the F I G U R E 2 Principal component analysis of hypocotyl traits.FW, fresh weight; FL, fresh length; TD max , maximum transverse diameter; TD min , minimum transverse diameter; RTD, ratio of TD max to TD min .Patterns of hypocotyl intraspecific traits variations along latitudinal gradients.(a) Relationship between PC1 and latitude.(b) Relationship between RTD and latitude.Blue line indicates ordinary least squares linear regression and red line indicates quadratic regression.Solid line indicates statistical significance (p < .05)and dashed line indicates statistical non-significance (p > .05).Filled area means standard error.PC1, the first principal component.RTD, ratio of TD max to TD min .y = −1.24x+ 29.59 p < .001,R² = .42y = −0.29x²+ 13.03x −142.76 p < .001,R² = .53variation in the shape index RTD, though significant effects of MAP (standardized coefficient = −0.30± 0.07) and AGB (standardized coefficient = −0.15± 0.07) were indicated by the model (Figure 4f).Bivariate analyses also showed a significant correlation between MAP or AGB and RTD, yet with very low explanation of the model (R 2 = .06and R 2 = .02,respectively; Figure4g,h).

F|
Sources of hypocotyl intraspecific traits variations.(a, e) Variance partitioning of PC1 and RTD across three nested levels of organization.(b, f) Summary of the linear mixed-effect modeling for multiple biotic and abiotic factors on PC1 and RTD.Data are presented as coefficients ± standard errors of the estimated effect sizes.R 2 m represents marginal R 2 and R 2 c represents conditional R 2 .*p < .05;***p < .001.(c, d, g, h) The relationships between significant biotic and abiotic factors and PC1 and RTD.Blue line indicates ordinary least squares linear regression and red line indicates quadratic regression.Filled area means standard error.SAL, sea surface salinity; MAT, mean annual temperature; MAP, mean annual precipitation; AFC, annual fastest tide current; ALT, annual lowest tide; AGB, aboveground biomass.PC1, the first principal component.RTD, ratio of TD max to TD min .Direct and indirect effects on hypocotyl traits We further used structural equation modeling to disentangle the direct and indirect effects of climatic, oceanic, and maternal factors on hypocotyl traits.Overall, these variables could explain 64% and 88% of the variation in PC1 and RTD, respectively.MAT played the strongest role in directly shaping PC1 (standardized path coefficients, β = .75).Sea surface salinity (β = −.22) was the main oceanic factor that directly affected PC1.Tide effects of annual lowest tide and annual fastest tide current on PC1 were much weaker than that of MAT.MAP had no direct effect and AGB had a weak effect on PC1 (Figure 5a,b).In contrast, all the explanatory variables, except MAP (β = −.31), had a weak direct effect on RTD (|β| < .15).Compared with the effects of climatic factors on PC1, effects of climatic factors were much smaller on RTD (Figure 5c,d).4 | DISCUSS ION Increasing studies have showed the importance of ITV on ecological processes.However, the patterns and sources of ITV in themselves are still unclear (Cope et al., 2022).There are large ITVs F I G U R E 5 Direct and indirect effects of biotic and abiotic factors on hypocotyl intraspecific trait variations.(a, c) Structural equation models showing the relationships among biotic, abiotic factors and PC1 and RTD.Blue and red arrows indicate positive and negative relationships, respectively.Solid and dashed lines indicate significant (p < .05)and non-significant relationships, respectively.Numbers near the pathway arrow indicate the standard path coefficients.R 2 represents the proportion of variance explained.(b, d) Standardized effects derived from structural equation models of PC1 and RTD.Blue and red bars indicate direct and indirect effect, respectively The relevant indicators are abbreviated as: PC1, the first principal component; RTD, ratio of TD max to TD min ; AGB, aboveground biomass; SAL, sea surface salinity; AFC, annual fastest tide current; ALT, annual lowest tide; MAT, mean annual temperature; MAP, mean annual precipitation; CFI, comparative fit index; RMSE, root mean square error of approximation.99, df = 1, p = .16,CFI = 0.99, RMSE = 0.00 χ² = 0.22, df = 1, p = .64,CFI = 1.00,RMSE = 0

| 9 of 11 YANG
et al. and artificial cultivation of true-viviparous mangrove species are particularly important, especially under global climate changes.