The study was conducted at two field trials in Monroe (33°51’ N, 88°17’ W) and Pontotoc (34°08’ N, 88°59’ W) counties, in northeastern Mississippi, US. These sites are located about 70 km apart and are located on the Gulf Coastal Plain of the US. The climate of the region is humid subtropical characterized by mild winters and long, hot summers. These sites have a 30 year (1991–2020) mean annual rainfall of about 1,490 mm, mean maximum and minimum temperatures of 22.5°C and 10.3°C, respectively and a mean annual temperature of 16.4°C (NOAA, 2023). Table 1 contains specifics on site elevation, depth to the water table and soil composition. Specifically, the Monroe site is categorized as an “alluvial terrace” located approximately 500 m from the Buttahatchee River with fine sandy loam soils and a water table within 2 m of the surface (Natural Resources Conservation Service, 2023). The Pontotoc site is characterized as an “upland” site with an elevation about 65 m higher than the Monroe site and a water table over 2 m below the surface. Soils are categorized as silt loams (Table 1; Natural Resources Conservation Service, 2023). Volumetric soil moisture contents were measured at each field site with 30 cm long time domain reflectometry probes (Campbell Scientific Inc. Logan UT, USA) attached to a CR1000 datalogger (Campbell Scientific Inc.) that recorded data half-hourly. Relative soil moisture content in the top 30 cm was calculated as soil moisture divided by the highest soil moisture in the sample period (representing field capacity) and converted to a percentage. Figure 1 presents relative soil moisture contents for both the Monroe (alluvial terrace) and Pontotoc (upland) sites and shows that they did not differ significantly from one another during the growing season with the alluvial terrace site being near field capacity during the dormant season (Figure 1).

The alluvial terrace was established and planted in April 2018, and the previous use of this site was peanut production. The upland site was established and planted in April 2019 and was previously used for annual row crop production. Both sites were tilled and subsoiled along planting rows that were 2.7 m apart to a depth of 45.7 cm and treated with pre-emergent herbicide before planting. One hundred and 108 unique genotypes were planted at the alluvial terrace and upland sites, respectively, (Table 2). Six taxa including improved eastern cottonwoods (D × D) as well as hybrid poplars were used including P. deltoides × P. maximowiczii (D × M), P. deltoides × P. nigra (D × N), P. deltoides × P. trichocarpa (D × T), P. trichocarpa × P. deltoides (T × D), P. trichocarpa × P. maximowiczii (T × M) and P. deltoides × P. nigra × P. maximowiczii, (D × N × M) which was grouped with D × N genotypes since only one genotype was represented (Supplementary Table S1). Dormant and unrooted hybrid poplar cuttings (22.9 cm long) were obtained from GreenWood Resources (Portland, OR, USA) and kept in a walk-in freezer until planting. Dormant, unrooted, eastern cottonwood cuttings (40.6 cm long) were obtained from the Department of Forestry Cottonwood Cutting Orchard at Mississippi State University. Before planting, hybrid poplar cuttings were soaked in insecticide (Admire® Pro, Bayer Corp., Whippany, NJ, USA) for 3–4 days, and eastern cottonwood cuttings were also soaked in water and insecticide for 18 h, wrapped in plastic film, and placed in the freezer prior to planting. Cuttings were planted in subsoiled rows at a 1.8 m spacing and two ramets of the same genotype were planted adjacent to one another in each block. Taxa were contained together in subplots within blocks, when possible, given the number of genotypes per taxa at each site, to provide consistent competition around trees.

For the first growing season, sites were maintained with minimal competition from herbaceous vegetation through hand weeding combined with tilling and chemical herbicide (Select, Valent USA LLC., Walnut Creek CA, USA) application. Additional fertilization and irrigation were not added to the site. If trees were being impacted by cottonwood leaf beetle (Chrysomela scripta Fabricius) larval feeding, sites were aerially sprayed with insecticide (Admire® Pro, Bayer Corp., Whippany, NJ, USA). Tree heights were measured at the end of the first growing season using a height pole.

During approximately July and early August of the first growing season at each site, leaf gas exchange measurements were made on an attached, fully expanded and developed leaf near the top of the tree using a LiCOR 6400 portable photosynthesis system (calibrated in 2015; LI-COR Biosciences Inc., Lincoln, NE, USA). One individual per genotype (generally the larger of the two) in three blocks at each site was measured for leaf gas exchange during midmorning on sunny days after dew had evaporated, but before significant midday leaf stress. Generally, one replicate block was measured per day and blocks were chosen that represented a range of study site conditions within the larger study sites. Light levels in the leaf chamber were 1,500 μmol m⁻² s⁻¹, CO₂ concentrations were 400 ppm and relative humidity was controlled to be between 40 and 60%. Measurements were made between about 8:30 and 12:30 local time, leaf temperatures ranged from about 27 to 37°C and leaf VPD ranged from about 0.8 to 3.8 kPa. Net assimilation rates (A_net; μmol m⁻² s⁻¹) and stomatal conductance (g_s; mol m⁻² s⁻¹) were monitored and allowed to equilibrate for approximately 2 minutes before measurements were taken with the ratio of the two representing intrinsic water use efficiency (iWUE).

After gas exchange measurements, leaves were placed in labeled, zip top bags and into a cooler to be transported back to the laboratory. Leaves were kept in a refrigerator until their fresh areas could be measured using a LI-COR 3100 leaf area scanner. Leaves were then placed in envelopes and oven dried for approximately 1 week at 65°C, then their dry weights were obtained. The ratio of leaf fresh area and dry mass was used to calculate leaf dry mass per unit area (LMA; g m⁻²) which was used to calculate net assimilation on a per unit leaf mass basis (A_mass). Dry leaves were then separated from the petiole and midrib and ground to pass through a 2 mm sieve. Ground leaf material was encapsulated in tin capsules, and sent to the University of California, Davis Stable Isotope Facility (Davis, CA, USA) for total sample carbon (C), and total sample nitrogen (N), and stable isotopes of carbon and nitrogen (δ¹³C, δ¹⁵N) using elemental combustion analysis (PDZ Europa ANCA-GSL elemental analyzer; Sercon Ltd., Cheshire, UK) and inductively coupled plasma optical emission spectrometry (PDZ Europa 20–20 isotope mass spectrometer; Sercon Ltd., Cheshire, UK) respectively. Leaf nitrogen content was calculated on a per unit area basis (N_area; g N m⁻²) by multiplying leaf N concentration by LMA. Photosynthetic nitrogen use efficiency (PNUE; μmol g⁻¹ N s⁻¹) was calculated by dividing leaf photosynthetic rates by nitrogen content on a per unit area basis.

Carbon isotope analysis provides an indirect measure of stomatal diffusion, but more importantly, of water use efficiency due to the inherent discrimination of the RuBisCO enzyme against ¹³C (Farquhar et al., 1989) and potentially provides a temporally and spatially integrated measure of intrinsic water use efficiency. Nitrogen isotope analysis can be more difficult to interpret compared with δ¹³C due to the multitude of soil nitrogen transformations and soil microbial species performing these chemical reactions and each exhibiting differing discrimination against ¹⁵N (Craine et al., 2015). However, in general, δ¹⁵N can distinguish between sources of N (soil or atmospheric; Craine et al. (2015)) as well as potentially indicating N availability and if a site is nearing N saturation (Craine et al., 2009).

Because gas exchange data were measured throughout the day and on multiple days, net assimilation and stomatal conductance data were scaled to a common temperature before further analysis (Willis et al., 2023). Using these corrected data, intrinsic water use efficiency was calculated as the ratio of net assimilation and stomatal conductance. Significant differences in physiological parameters across sites were accessed using ANOVA (aov) in R version 4.2.2 (R Core Team, 2022). Graphs were produced and p-values, correlation coefficients, and data quartiles were calculated using Sigmaplot version 13 (Systat Software Inc. San Jose, California, USA). Spearman rank correlations (ρ) and p-values were calculated in R using the “corr.test” function on genotypes that were present at both sites to evaluate ranking stability between sites (40 genotypes; Supplementary Table S1). Weak correlation between parameters (high p-value, low ρ) indicate a large difference in genotype ranks between sites and suggests the presence of a G × E interaction and traits with high plasticity. To determine total rank differences, genotypes were ranked from high to low values for each physiological parameter and absolute differences between sites for each genotype were calculated. These absolute differences were summed for each genotype to determine their total rank difference. Principal component analysis was performed for each site in R using the “prcomp” function with the datasets scaled and zero centered.

Although growth was higher across taxa at the alluvial terrace site, many leaf physiological and biochemical metrics were higher at the upland site including A_net, g_s, LMA, leaf N concentrations and δ¹³C (signifying more moisture stress and/or more water use efficiency; Table 3). The alluvial terrace site also had higher δ¹⁵N across taxa and higher leaf-level iWUE in D × D and D × M individuals (Table 3).

Physiological and biochemical parameters to predict productivity in terms of yearly height growth were site specific with no common parameters across sites (Figure 2). At the alluvial terrace site, height growth was positively correlated with leaf-level iWUE, and negatively correlated with A_net, g_s and PNUE (less so for D × M than other taxa; Figure 3). At the upland site, height growth was negatively correlated with leaf level iWUE and δ¹³C and positively correlated with g_s (Figures 2, 3). When separating out D × D and hybrid taxa that were planted at both field sites, no parameters were significantly correlated with height growth for D × D genotypes at the upland site while height growth was negatively correlated with A_net and positively correlated with δ¹⁵N at the alluvial terrace site (Figure 4). For hybrid poplar taxa, height growth at the alluvial terrace site was negatively correlated with water use efficiency in terms of δ¹³C and negatively correlated with N concentrations, while being positively correlated with δ¹⁵N and negatively correlated with A_net at the upland site (Figure 4).

Several physiological parameters exhibited consistent relationships across all planted Populus taxa at both the alluvial terrace and upland sites including positive relationships between A_net and g_s as well as PNUE, negative relationships between g_s and iWUE and δ¹³C (Figure 2), and negative correlations between leaf N concentration and LMA (Figure 5) and well as PNUE (Figure 2). Particularly, δ¹³C and iWUE were positively correlated at both sites (Figure 5). Stomatal conductance was also positively correlated with PNUE across taxa and sites suggesting a consistent relationship between water use and nitrogen use efficiency (Figure 5). LMA was negatively correlated with δ¹⁵N at both sites across all taxa (Figure 5) and δ¹⁵N and δ¹³C were also negatively correlated with one another (Figure 2). However, δ¹⁵N exhibited opposing relationships with leaf N concentration being positively correlated at the upland site, but negatively correlated at the alluvial site (Figure 5). For the 40 common genotypes (D × D and hybrids) planted at both sites, physiological relationships including positive correlations between A_net and g_s as well as PNUE and negative relationships between g_s and iWUE were evident (Figure 4). PNUE and g_s were also positively correlated across sites and Populus types suggesting a broad relationship between water use and PNUE. For D × D genotypes, leaf nitrogen parameters exhibited more similar relationships across sites with leaf N concentration being positively correlated with A_net, and δ¹³C and negatively correlated with LMA and PNUE (Figure 4). For Populus hybrids, fewer consistent relationships were evident across sites except for iWUE which was positively correlated with δ¹³C and negatively correlated with δ¹⁵N at both sites (Figure 4). Nevertheless, P. deltoides tended to have more significant correlations among leaf traits particularly with regard to iWUE at the alluvial terrace site, while hybrid poplars exhibited more trait correlation particularly for leaf nitrogen traits at the upland site (Figure 4).

For principal component analysis, the first two PCs explained about 59 and 53% of the variation in data for the alluvial terrace and upland field sites, respectively, (Figure 6). At both sites, iWUE and δ¹³C were positively correlated and δ¹⁵N and PNUE were positively correlated with one another. At the alluvial terrace site, assimilation rates were positively correlated with g_s and negatively correlated with LMA, while at the upland site, assimilation rates were positively correlated with leaf N_conc and PNUE and negatively correlated with LMA and leaf CN ratios (Figure 6). At the alluvial terrace site, yearly height growth was positively correlated with iWUE and negatively correlated with PNUE and δ¹⁵N. At the upland site, height growth was positively correlated with g_s and negatively correlated with δ¹³C and iWUE (Figure 6).

For the 40 genotypes planted at both field sites, rankings among genotypes changed the most, and G × E interactions were the most evident, in terms of physiological parameters, in particular A_net, nitrogen parameters, and to a lesser extent, water use parameters (Table 4). For A_net, D × M genotypes were the most likely to have high rates at the upland site but low rates at the alluvial terrace site, while the opposite was true for D × D which tended to have high rates at the alluvial terrace site as opposed to the upland site (Figure 7). For leaf N concentrations, D × D genotypes were more likely to be in the extremes with high leaf N concentration at one site and low leaf N for the same genotype at the other field site (Figure 7). Most other physiological traits exhibited a mix of taxa with opposing ranking at one site vs. another (Figure 7). Leaf isotope ratios (δ¹³C and δ¹⁵N) and leaf anatomical parameters (LMA) tended to have the lowest plasticity and low G × E interaction as they were correlated across sites (Table 4 and Figure 7). Finally, height growth was consistent across sites suggesting minimal site effects on height rankings among genotypes (Figure 7).

Among the genotypes that were common at both field sites, certain genotypes were more likely to exhibit trait plasticity across sites when all measured parameters were taken into account while others exhibited fairly low plasticity (Figure 8). Taxa did not factor into trait plasticity with each taxa containing genotypes that exhibited high and low plasticity across sites. For D × D genotypes, ‘83–2’, ‘S7C2’, ‘S7C4’ and ‘27–5’ exhibited high trait plasticity while genotypes ‘122’, ‘95A-6’, ‘6–1’ and ‘6–5’ exhibited relatively low plasticity (Figure 8). For D × T and T × D genotypes ‘5077’ and ‘9552’ (both T × D) exhibited high plasticity while ‘10029’ and ‘7903’ (both D × T) exhibited relatively low plasticity (Figure 8). Finally, for D × M genotypes ‘9252’, ‘13724’, ‘8019’ and ‘13738’ exhibited relatively high plasticity while ‘9225’ exhibited low plasticity between sites (Figure 8). Overall, trait plasticity was not correlated with height growth (p = 0.52; data not shown) with both low and high plasticity genotypes of each taxa capable of high growth.

Year one height growth was significantly greater at the alluvial terrace site for all taxa groups compared with the upland site and site and taxonomic differences meant that no single measured physiological parameter predicted growth across sites and taxa. Across all taxa at the alluvial terrace site, a tradeoff was evident in terms of either using water or nitrogen efficiently with a positive correlation between height growth and iWUE but a negative correlation between height growth and stomatal conductance as well as PNUE. Several other studies found a tradeoff between water and nitrogen use efficiency (Sage and Pearcy, 1987; DeLucia and Schlesinger, 1991; Hikosaka et al., 1998; Chen et al., 2005; Broeckx et al., 2014), since restricting CO₂ intake and exhibiting a high water use efficiency should make photosynthetic enzymes, and the nitrogen used to make them, less efficiently utilized. Other studies have also reported a positive relationship between growth and water use efficiency in Populus (Zhang et al., 2004; Monclus et al., 2009; King et al., 2013; Renninger et al., 2021b; Renninger et al., 2023). At the upland site, the opposite was true with height growth being negatively correlated with both iWUE and δ¹³C and positively correlated with stomatal conductance which was seen in other studies as well (Monclus et al., 2009; Renninger et al., 2021a) and fits with the fast/slow hypothesis (Reich, 2014) of high resource acquisition leading to higher growth rates. These data as well as the results from other research suggest that whether growth is positively or negatively related with water use efficiency is related to site characteristics and potentially water availability. However, it is surprising that at the alluvial terrace site with a seemingly more abundant water supply, height growth would be positively associated with water use efficiency, while height growth was negatively correlated with water use efficiency parameters at the drier upland site. This could suggest that access and acquisition to water leads to greater height growth at the upland site, while at the alluvial terrace site, efficient use of water at the leaf level may allow for more overall leaf area and therefore greater growth. This is corroborated by the results of the PCA analysis wherein height growth was positively correlated with iWUE at the alluvial terrace site but was inversely correlated with iWUE at the upland site. Other studies found that leaf area was correlated with productivity in Populus (Pellis et al., 2004; Marron et al., 2005; Monclus et al., 2005; Renninger et al., 2021b). Likewise, the large growth rates seen at the alluvial terrace site may have led to water resource limitations, or at least limitations in resource acquisition compared with the rate of growth occurring.

Across sites and taxa, broadscale relationships and coordination among physiological parameters were evident including well established relationships between A_net and g_s as well as PNUE, negative relationships between g_s and iWUE and δ¹³C and negative correlations between leaf N concentration and LMA and well as PNUE. The negative correlation between LMA and leaf N concentration is in line with the worldwide leaf economic spectrum in which “quick-return” leaves exhibit high leaf N and low LMA with the opposite being true of “slow-return” leaves (Wright et al., 2004). Other broadly significant relationships were the positive correlation between leaf gas exchange-estimated iWUE and δ¹³C, an estimate of water use efficiency based on carbon isotope ratios which has been seen in other studies as well (Soolanayakanahally et al., 2009; Bogeat-Triboulot et al., 2019). This suggests that the instantaneous measures from single leaves on the canopy provided a broad representation of the water use efficiency of the individual as is captured by leaf isotope data.

Another broad-scale relationship is the positive correlation between stomatal conductance and PNUE likely explained by the aforementioned mechanism in which higher stomatal conductance allows leaf nitrogen in photosynthetic enzymes to be used more efficiently (Field et al., 1983). Additionally, a positive correlation between leaf N concentration and δ¹³C seen across taxa at the alluvial terrace site, in D × D at both sites and in hybrids at the upland site suggests the same scenario in which greater leaf nitrogen and photosynthetic enzymes allow water to be used more efficiently (Brueck, 2008). Finally, a negative correlation between LMA and δ¹⁵N was observed at both sites. This relationship is more difficult to interpret as a wide variety of soil characteristics and processes affect δ¹⁵N of nutrients taken up by plants (Craine et al., 2015). However, N saturated sites have been shown to have higher δ¹⁵N values (Craine et al., 2009) and more microsite N availability within sites may lead to “quick return” leaves with lower LMA, particularly given the negative relationship between leaf N and LMA. This explanation is corroborated by the positive relationship between leaf N concentration and δ¹⁵N at the upland site, however the opposite relationship is seen at the alluvial terrace site. Nitrogen isotope ratios at the alluvial terrace site may be further complicated by the loss of the lighter isotopes of N due to denitrification (Blackmer and Bremner, 1977) which potentially explains why that site has more positive δ¹⁵N values than the upland site. Therefore, issues of greater denitrification and potential nitrate leaching (Ostrom et al., 1998) due to the greater soil saturation at the alluvial terrace site may be affecting the isotopic ratio of nitrate and ammonium and affecting its relationship with leaf N concentration. Overall, these differing relationships help determine what may be site-specific trait interactions and which are more broadly held across environmental conditions.

For genotype × environment interactions, leaf gas exchange, particularly A_net and PNUE, as well as leaf N parameters exhibited the most plasticity between sites, while tree height growth and leaf anatomical and isotopic parameters were most similar in terms of displaying similar rankings among genotypes between sites. This means that site effects should not impact studies that aim to identify the most productive genotypes in field trials at least for similar geographic conditions. Similar findings were reported by Marron and Ceulemans (2006) and Nelson et al. (2018), but differ from findings from Renninger et al. (2021a) who found significant G × E interactions in productivity parameters across study sites, although sites were in different physiographic regions and exhibited different disease pressures from Septoria stem canker (Sphaerulina musiva). Other studies also report G × E interactions in terms of productivity (Marron et al., 2006; Marron et al., 2007; Pliura et al., 2007; Rönnberg-Wästljung et al., 2022) with some genotypes exhibiting more stability in productivity across sites than others (Alvarez et al., 2020). Plasticity in physiological functioning has been reported by other researchers who found changes in gene expression as a result of water stress conditions (Vigar et al., 2016) and plasticity in sylleptic branching (Wu and Stettler, 1998; Cooke et al., 2005; Marron et al., 2006), as well as in xylem hydraulic traits (Plavcová and Hacke, 2012; Garavillon-Tournayre et al., 2018; Lemaire et al., 2021) and leaf physiological traits (Garavillon-Tournayre et al., 2018). For leaf metrics, isotope (δ¹³C and δ¹⁵N) and anatomical (LMA) parameters had the lowest plasticity with genotype rankings remaining similar across sites. Renninger et al. (2021a) also found that LMA, δ¹³C and δ¹⁵N exhibited similar rankings among genotypes across different field sites and Himes et al. (2021) and Voltas et al. (2006) found that δ¹³C exhibited low plasticity across irrigation treatments. However, Guet et al. (2015) found that LMA and δ¹³C exhibited significant phenotypic plasticity across sites in P. nigra genotypes. Nevertheless, these leaf traits may represent robust parameters for broad phenotypic comparison since genotypes maintained their ranking for these parameters regardless of site.