The effect of plant–plant interactions as a key biotic process mediating the spatial variation of phenotypes in a Pinus sylvestris forest

In this work, we highlighted the importance of the phenotypic structure of forest in regulating inter-tree competition with scattered individuals showing larger growth than close neighbours, with lower growth rates. Plant interactions are among the fundamental processes shaping the structure and functioning of ecosystems as they modulate competitive dynamics. However, the connection between the response of individual growth to neighbours and to environmental conditions and the mechanisms determining interactions in monospecific stands remain poorly understood. Here, we followed a phenotypic-based approach to disentangle the effect of plant size, neighbourhood interactions and microhabitat effects on Pinus sylvestris growth and traits, as well as their spatial variation of growth. We mapped all adult trees (1002 pines) in a 2 ha stand and measured their height, DBH and crown projection. For each individual, we assessed its growth and a competition index in relation to the closest neighbours. Soil chemical and physical properties and ground cover were also measured in a grid within the stand. We analysed the effects of tree size, neighbour competition and microhabitat variation on tree growth with a linear model. We also used spatial mark-correlation functions to explore the spatial dependence of tree age, secondary growth and phenotypic traits. Our results showed that trees with close neighbours displayed lower growth rates, whilst individuals with larger growths appeared scattered throughout the stand. Moreover, we found that growth depended on competition, tree height and crown area while tree growth poorly correlated with age or microhabitat conditions. Our findings highlight the importance of forest structure, in regulating inter-tree competition and growth in a Mediterranean pure stand and they provide insight into the causes and consequences of intraspecific variation in this system.


Introduction
Understanding tree growth variation within a forest is critical for ecologists because it determines the spatial distribution of biomass across the whole forest ecosystem and the services it provides (Burkhart and Tomé 2012). Unmanaged natural forests normally exhibit a great variability in tree sizes and functional characteristics, reflecting profound differences in tree growth as a result of the combined effects of biotic and abiotic factors which act at fine spatial scales (Gómez-Aparicio et al. 2011;Kunstler et al. 2012;Forrester et al. 2013;Fraver et al. 2014;Calama et al. 2019). Intrinsic features of each tree such as its age, architecture and access to the forest canopy (Wyckoff and Clark 2004;Madrigal-González and Zavala 2014;Chi et al. 2015) are major determinants that, together with extrinsic factors such as competition with neighbours, influence tree growth (Contreras et al. 2011). Abiotic heterogeneity (like soil nutrients, microclimatic conditions or ground covers) also plays a critical role in tree growth (Naithani et al. 2014) and influences its spatial variation within and among forests (Toledo et al. 2011;Khairil et al. 2014;Prada & Stevenson 2016). Consequently, tree growth is a multi-faceted biological process that can depend simultaneously on several interacting factors. Thus, size variability in plants is the result of a deterministic (genetically designed) growth pattern that relates individual growth rate to current plant size, and, stochastic patterns, due to microsite heterogeneity and neighbourhood effects (Bonan 1988). Examining the links between individual tree growth, tree size, neighbourhood interactions and microsite environmental conditions may thus help to understand the processes underlying how individual tree biomass is accumulated and how biomass is spatially structured in forests (John et al. 2007;Gupta and Pinno 2018).
Competition takes place at the very fine spatial scales where plant-plant interactions for resources occur (Lortie et al. 2004;Tilman 2004). Empirical evidence shows that niche differences are critical for coexistence and, consequently, a high degree of phenotypic dissimilarity among coexisting species is the norm. This implies a competition release at small spatial scales (Forrester et al. 2017;Gusmán et al. 2018;Bastias et al. 2020) and the maintenance of rich plant communities (Götzenberger et al. 2012). These phenotypic differences may also be critical to reduce the conspecific competition in the case of monospecific or species-poor stands (MacArthur and Levins 1967;Weiner 1990;Escudero et al. 2021). In fact, intraspecific phenotypic variation can be larger than heterospecific variation (Siefert et al. 2015;Bastias et al. 2017;Benavides et al. 2019b), and many recent works advocate for its explicit consideration to explain plant coexistence and community assembly (Messier et al. 2010;Violle et al. 2012;Kraft et al. 2015;Benavides et al. 2019a).
Understanding the importance of intraspecific variation on community assembly requires the so-called plant´s eye perspective (sensu Murrell et al. 2001), i.e. overcoming approaches based on species average simplifications and explicitly considering the variation among co-occurring individuals. This recognizes the role played by each individual in a community, shifting the focus from the inter-to the intraspecies phenotypic variation perspective (Cadotte et al. 2011).
Since phenotypic variation among individuals, regardless of the species, affects the functional structure of the whole assemblage (Violle et al. 2012;Hart et al. 2016) and the interactions occur between individual trees (Arroyo et al. 2015), a forest stand can be seen as a collection of interacting phenotypes located in specific points in the space (Aarssen 1983). The spatial occurrence of individuals is mediated by ecological processes (Ackerly and Cornwell 2007), which in turn affects tree performance according to individual-specific environmental conditions (Lasky et al. 2013;Michelaki et al. 2019) and neighbourhood interactions with conspecifics (He and Biswas 2019) affecting differentially each tree life stage from dispersal to reproduction. The relative importance of abiotic and biotic factors on tree growth and coexistence has been debated for decades, but there is a consensus that these mechanisms are not mutually exclusive. On the one hand, there is ample evidence indicating that conspecific neighbours in monospecific stands could result in strong competition for resources (Weiner 1990;Fichtner et al. 2017;Chen et al. 2017). On the other hand, intraspecific variation seems critical to respond to fine-scale edaphic variation (Cornwell and Ackerly 2009;Wright and Sutton-Grier 2012;Fajardo and Siefert 2018).
Here, we assessed the relative importance of plant size (i.e. ontogeny), biotic (i.e. neighbours' competition, finescale biotic heterogeneity) and abiotic factors (i.e. soil properties) on individual tree growth. We mapped 1002 trees in a well-preserved Mediterranean Pinus sylvestris L. pure stand found at its highest altitude in Central Spain located about three degrees latitude from the southernmost limit of the species distribution. In addition, we investigated the spatial variability of some phenotypic traits such as tree height, crown size and secondary growth, robust indicators of individual plant performance. More specifically, we explored the effect of biotic and abiotic environmental factors on tree growth, taking into consideration spatial variation of growth and phenotypic traits in the forest stand at contrasting spatial scales. We assumed that phenotypically dissimilar individuals minimize intraspecific competition (see Escudero et al. 2021). Accordingly, we hypothesize that neighbour interactions would optimize spatial phenotype distribution, increasing local variation on resource acquisition-related traits and we expected that, if competition is the leading process in this monospecific forest, secondary growth of each individual tree will be influenced by the phenotypic characteristics of their neighbouring trees resulting in a structured pattern of the tree phenotypes and growth in space. Namely, it is expected that the closer the neighbouring trees, the smaller do the individual growth and size occur; and additionally, the greater the differences among neighbours, the greater is individual growth.

Study area and phenotypic measurements at tree level
The study was conducted in a forest stand in Pingarrón (40 N 48′ 46"; 3 W 57′ 14"), at 1900 m. a. s. l., located in Sierra de Guadarrama National Park, Spain ( Fig. 1a-b). The climate is typical mountainous Mediterranean, characterized by wet, cold winters and warm, dry summers. Specifically, the annual mean temperature is 6.5 °C and the annual precipitation is 1328 mm (Puerto de Navacerrada station at 1894 m. a.s.l., Spanish National Meteorological Agency, AEMET: http:// www. aemet. es/ es/ porta da). The site is a natural pure stand of Pinus sylvestris L. in the north-facing slope very Fig. 1 Site description. a Location of the study site (orange circle) in Central Spain, Pingarrón, is located at Sierra de Guadarrama National Park (black square) at 1900 m a.s.l. b Picture of the study standphoto: B. Carvalho; c spatial layout of surveyed pines within the close to the timberline with an uneven structure and no signs of recent management. The bedrock is mainly composed of gneiss, and soils are acid and relatively homogeneous, being predominantly humic cambisols and leptosols (Fernández-González 1991). The understory vegetation is dominated by shrubs, such as Juniperus communis L. subsp alpina (Neilr.) Čelak, Cytisus oromediterraneus (G. López & C.E. Jarvis) Rivas Mart., Adenocarpus complicatus (L.) J. Gay. and mostly by grasses such as Avenella flexuosa subsp. iberica (Rivas Mart.) Valdés & H. Scholz.
From June to September 2017, we mapped and measured all adult pines (diameter at breast height (DBH) greater than 7.5 cm) over a 2 ha area, totalling 1002 trees (Fig. 1c). Each tree was mapped with a sub-metric GPS (GeoXT, Trimber® GeoExplorer® 2008 Series). Moreover, we measured two perpendicular crown diameters to assess crown projection area assuming an elliptical area, DBH and total tree height. We also extracted one wood core at breast height (1.3 m) using increment borers (5 mm diameter; Pressler Haglöf).

Tree-ring analysis
Wood cores were air dried, glued, and polished using a series of sandpaper until all tree rings' boundaries were clearly visible. Wood cores were scanned (1200 dpi) and then we visually determined the age (years) of all trees by counting the rings (Stokes and Smiley 1968) (Fig. 1d). When cores did not include the pith, we estimated the age by fitting a template of concentric circles with known radii to the curve of the innermost rings, which allowed the estimation of the missing radius length and transforming it into the number of missing rings (Norton et al. 1987). We measured the width of the tree rings (RW, Figure S1) to the nearest 0.01 mm using the CooRecorder software (version 9.1; Cybis Elekronik & DataAB, Saltsjöbaden, Sweden). Cross-dating accuracy was repeatedly checked using COFECHA (Holmes 1983). As the last ring (i.e. 2017) was not completely formed when the sampling took place, it was removed for further analyses. Finally, we selected the data from the last 5 years (to include the youngest adults) and estimated secondary growth as the basal area increment (BAI) using the formula BAI = πr 2 -(π (r-r 2012-2016 ) 2 ). This metric standardizes growth, accounting for decreasing ring width with increasing stem size, which is considered a better estimate of tree growth than individual ring width (Biondi and Qeadan 2008). Note that r is the total tree radius (cm) and r 2012-2016 corresponds to the radius (cm) of the last 5 years of tree-ring formation. In total, tree rings were accurately measured in 702 pines. In subsequent analyses, we employed the logarithm of BAI as a measurement of tree growth.

Biotic and abiotic environmental characterization
The study area was divided into 48 20 × 20 m plots, where we collected five soil samples (soil cores of 10 cm depth and 5 cm of diameter) per plot, one at the plot centre and the others at four distant points randomly selected in four directions (NE, NW, SE and SW) from the centre. Each plot, in turn, was subdivided into four subplots (10 × 10 m) where the percentage cover of leaf litter, shrubby vegetation, and ground stoniness were visually estimated always by the same field technician (Table 1).
Soil samples were air dried and sieved through a 2 mm mesh before analysing in the laboratory Nutrilab at the Universidad Rey Juan Carlos (Spain) (https:// nutri lab-urjc. es/). Five soil physical and chemical properties were measured: pH, electrical conductivity (µS/cm), total organic carbon (C; %), total nitrogen (N, mg/g) and available potassium (K, mg/g). We selected these soil variables because they are known to be important determinants for plant growth (Hasanuzzaman et al. 2018). The soil pH and electrical conductivity were measured in a soil and water suspension at a mass:volume ratio of 1:3 using a pH meter (GLP 21; Crison, Barcelona, Spain) and a conductivity meter (GLP 31; Crison, Barcelona, Spain), respectively. Soil organic C was determined by colorimetry after oxidation with a mixture of potassium dichromate and sulphuric acid (Yeomans and Bremmer 1988). Total N was determined with a SKALAR + + San Analyzer (Skalar, Breda, The Netherlands) in the laboratory after digestion with sulphuric acid and Kjedahl's catalyst (Anderson and Ingram 1994). Potassium (K) was measured with the same analyser system after shaking the soil samples with distilled water (1:5 ratio) for 1 h. Finally, we averaged the five soil properties obtained in each plot. The intensity of intraspecific competition experienced by each tree was evaluated using a competition index (CI). This index assumes a circular neighbourhood centred on the focal tree, thereby it defines the focal tree's competitive neighbourhood. For each i focal tree, CI was calculated where d i is the diameter of the focal tree, d j is the diameter of each of the N(r) neighbouring trees within a r distance, and dist i-j is the distance between the focal tree and each neighbouring tree (Lorimer 1983). The threshold distance r for the calculation of CI was r = 12 m. This value represents the distance that maximized the negative correlation between BAI and CI within the range 1-16 m ( Figure S2).

Building surface soil variables by kriging
Maps for each soil property and ground cover (leaf litter, shrubby vegetation and stoniness) were generated for the 2 ha area by kriging. We first transformed the soil and cover data with a Box-Cox transformation to meet normality assumption (Asar et al. 2017), using function boxcoxnc of 'AID' R package (Graves 2019). We then fitted trend-surface regressions and computed empirical variograms from the residuals. We fitted variogram models from the empirical variograms and obtained spatial predictions for 20 × 20 m blocks, using block kriging. The trend was then added back to the kriged means, and the values were back-transformed to the original scale (Table S1). The back transformation was conducted using Ecfun function in 'AID' R package (Dag et al. 2019). Geostatistical modelling was carried out with gstat R functions (Pebesma 2004), accessed throughout the R package 'automap' (Hiemstra et al. 2009).

Data analyses
As a preliminary analysis, we computed minimum and maximum values and the first and third quantiles for each variable looking for outliers, and based on them we removed seven individuals from our database. We performed pairwise Pearson´s correlations for the plant size traits (crown projection area, DBH and tree height), soil properties (pH, conductivity, C, N, K), and ground cover (vegetation, stoniness, leaf litter) variables for 695 individuals. We also checked if age was correlated with tree growth and the other plant size variables, and paired correlation analysis showed that there is no collinearity between age and tree growth (Table S2). The correlation analysis and Bonferroni correction for multiple testing were conducted using function corr.test in the 'psych' R package (Revele 2018).
We fitted least-square regressions to assess the effects of competition (CI), ground cover and soil composition on tree growth. Age, tree height and crown projection area were also included to account for allometric effects on tree growth. As CI was computed using a neighbourhood radius of 12 m, only trees located 12 m away from the border of the plot (i.e., 393 individuals) were considered in the regression. Prior to model fitting, all predictor variables were standardized, i.e. we calculated the mean and standard deviation (sd), then scaled by subtracting the mean and divided it by the standard deviation. We calculated the variance inflation factor (VIF) for all predictors before including them in the regression ( Figure S3). DBH was greatly correlated with crown area and age; we therefore removed it from further analyses. After model fitting, we checked that residuals of the model fulfilled the conditions of normality (p > 0.05) and homoscedasticity. We also checked the absence of residual spatial autocorrelation (Moran I statistic standard deviate = -1.32, p-value = 0.90; Dormann et al. 2007). Moran's I test was computed using function moran.test in the 'spdep' R package (Bivand et al. 2019). We used commonality analysis to assess the variance accounted for each individual variable and by groups of variables relative to plant size (crown projection area, tree height) and microhabitat (soil properties and ground cover variables) using commonalityCoefficients function in the 'yhat' R package (Nimon et al. 2008).
The spatial structure of functional variables such as tree growth, tree height, crown projection area and tree age were evaluated with the mark-correlation function (Stoyan and Penttinen 2000) using all 695 individuals. For a quantitative variable m (i.e. the "mark"), which varies throughout the points of a spatial point pattern, the mark-correlation function is defined as k mm (r) = c mm (r) tional mean of the product of the marks (i.e. E[m i m j ]) of all point pairs (i, j) separated by a distance r, and is the mean of m (Illian et al. 2008). It measures the spatial dependence of the marks (Baddeley et al. 2015). Values of k mm (r) larger or smaller than 1 indicate, respectively, positive and negative spatial dependence of mark values (Baddeley et al. 2015). k mm (r) function was estimated up to a distance of r = 37.38 m, with 0.07 m steps, using the isotropic correction to account for edge effects (Baddeley et al. 2019). Inference about the mark correlation function was obtained computing simulation envelopes from a null model of random labelling, i.e. randomly permuting the marks among the points (Wiegand and Moloney 2014). For the test, we computed 199 simulations and we used envelope and markcorr in the 'spatstat' package (Baddeley et al. 2019). All statistical analyses were conducted in R (version 4.1.3, R Core Team 2022).

Results
The magnitude of plant size (crown projection area, DBH and tree height) showed substantial variability, with CV varying between 37 and 93% (Table 1) and average crown projection of 21.8 m 2 ; DBH 24.5 cm and tree height 8.39 m. Some soil chemical properties (N and K) and conductivity showed high variability as well, with CV between 43 and 72%. N showed average values 3.15 mg/g, while K was on average ranging between 0.06 and 0.019 mg/g, and conductivity appears with the highest average 57.1 µS/cm. Similarly, ground cover variables showed CV over 40% with average ranging from 17.3% for stoniness, 33.6% for vegetation cover and 51.1% for leaf litter. In contrast, C and pH showed low CV, with 18% and average value 13.7% for C and 12% and 4.16 average value for pH.
We found that tree growth was strongly and negatively affected by the competition index (CI hereafter), which explained 11.1% of tree growth total variance (Table 2; Figs. S4, S5). Age also showed a significant negative effect on tree growth and explained 4.1% of total variation (Table 2; Figs. S4, S5). Tree height was positively related and explained 7%, whereas crown projection area, with a positive effect on growth, explained a very small fraction of 1.72%. Microhabitat conditions (soil properties and ground cover) explained only a fraction of 4.6% of the total growth variation with a positive effect of pH, conductivity and stoniness (Table 2; Fig. S4). However, our variance partition analysis also showed that CI, combined with crown projection area, contributed 9.9% and CI combined with tree height and crown projection area with 8.7% of the explained variance in tree growth (Fig. S5). Also, age combined with tree height, crown projection area and CI explained 4.43% of the total variation in tree growth (Fig. S5).
We also found that the correlation of tree growth values of closer tree pairs was lower than the expected by chance (negative dependence). However, for r > 10 m, k mm (r) values fluctuated around 1, indicating that at larger distances there was no spatial dependence of tree growth. The same pattern was observed for crown projection area (Fig. 2) and height, although for r > 15 m, k mm (r) values for tree height were higher than the expected by chance (Fig. 2). Mark correlation analysis showed that tree age in the entire area was randomly distributed at the forest stand scale (Fig. 2) with k mm (r) values fluctuating around 1.

Discussion
With this individual-based approach, we intended to disentangle the relative effects of biotic and abiotic microscale heterogeneity on tree growth in the presence of inter-tree competition, and to explore the spatial pattern of phenotypic traits and growth at fine scale. The relative importance of competition on tree growth has been subject to analysis for a long time (Goldberg and Novoplansky 1997;Woodall et al. 2003;Rozas 2014;Li et al. 2015;Zhang et al. 2016;Lustosa Junior et al. 2019). As expected, we found a negative effect of conspecific competition (CI) on growth in terms of basal area increment. Neighbouring trees compete for the same local resources, limiting the individual availability that eventually may reduce

Fig. 2
Spatial pattern of age and growth. Graphs on the right (b, d, f and h) show the estimation of the functional mark correlation function ( k mm (r)) calculated for tree age, basal area increment (log-transformed BAI values), tree height and crown projection area of 695 pines. The grey zone represents the envelopes obtained from 199 random labelling simulations. We considered random patterns when the observed function remained within the envelopes for all distances r (in metres), while we considered that the function deviates from a random pattern when the line remains outside the envelopes. Dashed ellipses indicate higher or lower values than expected by chance for pairs of individuals. Graphs on the left are maps of spatial distribution of the  Das et al. (2012) captured local competitive effects, also describing a steep decline in growth rates for four coniferous species in response to competition among neighbourhood trees in stands within the coniferous forests. Although to perform our model there was a substantial reduction in the number of individuals, it has been possible to keep a very considerable sample n. This was needed to achieve a very robust result of competition between individuals, allowing us to see clearly the mechanism of population structure related to resource partitioning resources.
We also found that its intrinsic tree height was the most important variable explaining the secondary growth of each tree, which agreed with expectations for structure and dynamics of forests based on individual-level allometric scaling relationships that guide trees' resources use, space filling, and growth (Enquist 2002;Enquist et al. 1999). We observed that together with tree height and crown projection area, neighbourhood competition explained an important fraction of the variation of tree growth (Fig. S5). Our results also showed that older trees grew more slowly than younger trees. In our work, the age effect on tree growth was negative, contrary to the positive influence of tree size, suggesting an ontogenetic reduction in tree growth with age (Stoll et al. 1994). After accounting for neighbourhood effects and tree size, we observed a weak influence of underlying local resources on tree growth, which only explained a small fraction of its variation. All the soil properties included in the regression are known for determining plant productivity (Almendro-Candel et al 2018;Hasanuzzaman et al. 2018); however, only pH, conductivity and, to a lesser extent, stoniness affected tree growth. Soil pH is a good indicator of the chemical status of the soil and informs of plant growth because soil pH influences nutrients uptake and consequently the tree growth. Low values for pH can make some nutrients unavailable to plants (Jensen 2010). Moreover, soil electrical conductivity is an indirect indicator of the amount of water and water-soluble nutrients available for plant uptake. The regression shows that the growth of Scots pines is affected negatively when the conductivity values decrease. It is worth noting that our sampling grain may be too coarse to detect variation in the other environmental factors considered, which would explain their relatively small effect on tree growth. Other environmental factors not assessed in this study such as soil moisture (Baker et al. 2003;Calama et al. 2019) and micro-topography (John et al. 2007) could also play an important role in the spatial pattern of tree growth.
The spatial variation of phenotypic characteristics and growth was mediated mainly by competition. We found a negative spatial dependence of growth at fine spatial scales (up to 10 m) and independence at larger distances (Fig. 2), which suggested a clear incidence of inter-tree competition. Similarly, in a study conducted by Fraver et al. (2014) with Picea abies in a Swedish boreal forest, they found that tree clustering locally intensified competition and reduced tree growth. Parallel to this result, we did not find in our study any spatial dependence of tree age in the population. This means that age was not a factor shaping individuals' spatial distribution in our stand, reflecting the natural dynamics in unmanaged stands, i.e. the spatial distribution of these phenotypes is not age dependent. It is worth mentioning, however, the tendency for older pines, which are also the largest, to be separated among them more than 20 m (Fig. 2a). We found that growth and crown projection area for any pair of pines separated up to 10 m were less correlated than expected by chance, suggesting the relevant negative effect of competition, after discarding individuals' age as we mentioned above. A similar pattern was also observed in the study of Getzin et al. (2008) in a monoculture of Pseudotsuga menziesii, suggesting that the crown projection area is very sensitive for the detection of competition in stands. Also, our results showed that the higher individuals tend to be more separated among them than expected by chance. This suggests that taller individuals exert asymmetric competition to their closest neighbours and require a suppressed neighbourhood of at least 15 m radius to establish themselves as dominant trees. Although plants are sessile, they can modify their growth pattern to minimize light interception and the use of soil water by neighbours. Some studies have already shown that P. sylvestris has an inefficient water transport system to deal with drought events, which are very common in the Mediterranean region (Vilalta and Piñol 2002, Eilmann and Rigling 2012), the southernmost region of distribution of Scots pine, where it is located at the limit of its water capacity to survive.
Based on this, tree growth of our study trees can be effectively limited by drought and asymmetric competition may favour taller phenotypes. This phenomenon appears to result from competition for light and water, which is 'one-sided,' with larger plants shading smaller plants and taking up more soil water, while smaller plants exerting almost no effect on the available resource to their larger neighbours (Ding et al 2019;Freckleton and Watkinson, 2001;Picard 2019). Steckel et al. (2020) showed that larger Scots pine trees were significantly more able to maintain growth levels during drought than small Scots pine trees. This result suggests that strong resource competition can lead to neighbouring trees with less variable trait distribution. Our findings indicate competition as an important mechanism in this system demonstrating that competition continues to influence forest processes and functional trait space in an old-growth system. These results hold our hypothesis posing that competition is key to understand stand structure and tree spatial layout, and support the community assembly theory that emphasizes the prevalence of inter-tree interactions at fine scales as drivers of coexistence (Law 1999).

Conclusions
The individual responses within a species using a point pattern analysis indicated that individuals with very close neighbours (< 10 m) grew less, suggesting that competition is a key biotic process mediating the spatial structure of the population. Moreover, secondary growth was positively correlated with tree height, an important trait related with tree competitive ability (Rijkers et al. 2000;Valladares and Niinemets 2007), and negative with the neighbourhood competition. Our findings demonstrate that individuals of different age have their size and growth affected as a result of competition with their neighbours and highlight the importance of the spatial arrangements of trees, i.e. forest structure, in regulating inter-tree competition and growth in this population of P. sylvestris at the rear edge of its geographic distribution and near the altitudinal limit of the species, which actually was the treeline in this mountain range.