Assimilation Efficiencies and Gas Exchange Responses of Four Salix Species in Elevated CO₂ under Soil Moisture Stress and Fertilization Treatments

Abstract

Assimilation to the internal CO₂ (ACi) response curve and gas exchange parameters were quantified for four North American willows ((Salix cordata (COR), S. discolor (DIS), S. eriocephala (ERI), and S. interior (INT)) grown in a 2 × 2 factorial of atmospheric CO₂ and soil moisture treatments to see how they would respond to climate change factors. After the first year of greenhouse growth under said treatments, we saw no difference in the aboveground stem biomass between CO₂ treatments. Thus, in the second year, a second experiment on a subset of well-watered, coppiced willows was conducted in a 2 × 2 factorial of atmospheric CO₂ and soil fertilization (FERT) treatments. In both experiments, the maximum rate of carboxylation (V_cmax) significantly declined for all four species in response to elevated CO₂ (eCO₂). In response to a drought treatment (DRT), V_cmax declined, except for INT, which increased V_cmax. In both experiments, INT had the greatest V_cmax, maximum rate of electron transport (J_max), and triose phosphate utilization, followed by COR and ERI, with DIS having the lowest values. FERT resulted in a strong increase in assimilation (A) and stomatal conductance (G_wv) by 92 and 119%, respectively. G_wv is the primary driver and A is a minor driver of water use efficiency (WUE) under DRT. FERT mitigated the V_cmax and A downregulation in eCO₂, but eCO₂ did not mitigate the DRT downregulation effect. Differences between INT and the other three willows in a number of adaptive traits and responses related to drought may reflect the evolutionary origins of INT and the taxonomic group Longifoliae in the arid southwest USA and Mexico.

1. Introduction

Ecological adaptations related to carbon assimilation are important to plant fitness because they affect growth and survival in different environments and this has implications for interspecific competition. Assimilation traits vary widely among species, environments, and plant age [1,2]. Assimilation downregulation (A_dr) of the photosynthetic apparatus (e.g., chlorophyll and biochemical efficiency traits), also known as biochemical downregulation, often occurs in response to elevated CO₂ (eCO₂) and results in the reduction in maximum carboxylation efficiency (V_cmax) and maximum electron transport (J_max) among other related traits [3,4,5,6]. Studies have shown that photosynthetic apparatus downregulation generally ranges from no downregulation to complete downregulation (when assimilation (A) is the same under eCO₂ as it is under ambient levels), depending on species and study [2,7].

With climate warming, frequency and severity of drought are predicted to increase [8]. Water is a predominant factor in determining the geographic and landscape distribution of vegetation, and drought is an important factor in reducing plant growth and increasing mortality [9]. Physiologically, water stress and heat reduces plant and leaf morphology, water potential, photosynthesis, biochemical responses, and can increase photo-oxidative stress and leaf abscission, [10,11,12]. Tree physiological responses to water availability vary among and within species [10,12,13].

Salix cordata Michx. (COR), S. discolor Muhl. (DIS), S. eriocephala Michx. (ERI), and S. interior Rowlee (INT) are native to eastern and central North America and are promising as fast-growing sources of woody biomass production [14]. Salix cordata, also known as sand dune willow, is commonly associated with sand dune landscapes of the Great Lakes—St. Lawrence River of northeastern North America. Salix interior is also widespread in western North America, ranging from Mexico to Alaska and has been observed to be saline tolerant in greenhouse experiments [15,16]. Salix discolor is commonly found in wet areas on a wide variety of disturbed sites, whereas ERI and INT are more commonly associated with fast-flowing streams and rivers.

The predicted doubling of CO₂ this century and increasing drought with anticipated global warming may affect comparative physiological processes among species. There are few climate change physiological response studies for North American willow species. We could only find one study for willow response to eCO₂ and it had a water stress treatment, using a single clone of Salix sagitta [17]. For willow responses to drought, we found three studies: [18] using S. purpurea L.; [19] with S. nigra Marshall, and [20] with a larger group of willows that included S. candida Fluegge ex Willd., S. pedicellaris Pursh, S. pyrifolila Andersson, S. bebbiana Sarg., S. discolor, and S. petiolaris Sm. Our goal was to examine and compare assimilation efficiencies and gas exchange (GE) traits for four North American willow species under the interactive effect of atmospheric CO₂, soil moisture, and soil fertilization (FERT) treatments. We hypothesized that eCO₂, soil moisture stress, and FERT responses in willows would differ among species. To test the hypothesis, we (1) quantified the four willow species’ genetic variation in assimilation to internal CO₂ (ACi) response curves and GE parameters, in response to eCO₂, soil moisture stress, and FERT treatments, and (2) tested the three treatment effects on the species’ physiological interrelationships using covariate analysis.

2. Materials and Methods

This greenhouse study grew plants from dormant willow cuttings under CO₂ and soil moisture treatments over two growing seasons (2017 and 2018), with the above-ground biomass harvested after the first year, using the coppice regrowth during the second year of treatments. In the second year, a sub-set of well-watered coppice plants were given different levels of fertilizer.

2.1. Plant Material, Growing Conditions, and Treatment Delivery

Sixteen genotypes (clones), four clones for each of the four willow species, COR, DIS, ERI, and INT were used for this study (

Table 1. Species and clones of Salix used in 2017 CO₂ × soil moisture and 2018 CO₂ × fertilization experiments.

Species	Clone ID	Provenance	Lat-Long	Measured
				2017	2018
S. cordata	BIG-C4m	Big Sandy Bay, ON	44°06′ N–77°43′ W	X
	BIG-C5f	Big Sandy Bay, ON	44°06′ N–77°43′ W
	NOR-C4f	North Beach Park, ON	43°56′ N–77°31′ W	X	X
	OUT-C1m	Outlet Beach, ON	43°53′ N–77°13′ W		X
S. discolor	HAW-D4m	Hawkesbury, ON	45°36′ N–74°36′ W	X
	HAW-D5f	Hawkesbury, ON	45°36′ N–74°36′ W		X
	MON-D1f	Montmagny, QC	46°58′ N–70°33′ W	X	X
	RIC-D2m	Richmond Fen, On	45°07′ N–75°42′ W
S. eriocephala	ALL-E2m	Allumette Island, QC	45°54′ N–77°06′ W	X
	FRE-E1f	Fredericton, NB	45°57′ N–66°38′ W		X
	GRE-E1f	Green River, NB	47°34′ N–68°19′ W	X	X
	WAI-E1m	Wainfleet, ON	42°55′ N–79°20′ W
S. interior	LAF-I2f	Ottawa, ON	45°25′ N–75°41′ W	X	X
	LAF-I12m	Ottawa, ON	45°25′ N–75°41′ W
	WAI-I1m	Wainfleet, ON	42°55′ N–79°20′ W	X
	WAI-I2f	Wainfleet, ON	42°55′ N–79°20′ W		X

). Dormant stem cuttings approximately 20 cm long were collected either in the field in April 2017 and kept frozen (−8 °C), or in early May 2017 and kept cool (4 °C), in sealed plastic bags until 7 June 2017 when all the stored stem cuttings were moved to a common 4 °C cooler. Cuttings were placed in water at room temperature for two days, prior to potting on 13 June 2017.

Cuttings were rooted and grown individually in fabric root-control bags 30 cm in diameter and 23 cm deep (Smart Pot PRO 5 Gallon, High Caliper Growing Systems, Oklahoma City, OK, USA) filled with sandy soil. Four bulk soil samples that were sent for texture analysis to the University of New Brunswick nutrient analyses laboratory, showed sand, silt, and clay portions, of 87.5 ± 0.7, 7.8 ± 0.4, 4.8 ± 0.3 (mean + S.E.), respectively. Root bags were arranged on a sand-covered gravel surface, inside eight custom-made chambers, each 2.3 m wide × 4.25 m long × 2.1 m high, covered in UVA-clear 4 mil polyethylene film, with outside air delivered and circulated by positive pressure fans. Each chamber contained four randomly distributed replicates (ramets) of each genotype—thus, 64 plants per chamber (4 ramets × 16 clones), × eight chambers (a total of 32 ramets per genotype). The eight chambers were housed inside two separate greenhouse compartments (two blocks of four chambers) at the Canadian Forest Service–Atlantic Forestry Centre (CFS–AFC) in Fredericton, NB, Canada (45°52′ N, 66°31′ W).

Plants were grown under natural daylength with light levels inside the chambers approximately 65% of outside levels, determined according to [21]. Temperatures inside the chambers were controlled by venting the greenhouse chambers, and by individual air conditioning units within each chamber. Water was delivered to each rooting bag via individual spigots in a dripper system. All root bags were kept well-watered for the first two weeks prior to starting the treatments, to allow for rooting of the dormant stem cuttings. A custom-made logging and control system continuously monitored sensors inside each chamber for air temperature and humidity (model HMP 155, Vaisala, Vantaa, Finland, and model AM2302, Adafruit Industries, New York, NY, USA), soil moisture (volumetric water content VWC), temperature (model 5TM, Meter Group, Pullman, WA, USA), and CO₂ concentration sampled at plant height, mid-chamber (Model 840A, LI-COR Inc., Lincoln, NE, USA).

In each greenhouse, each chamber was randomly assigned to one of four treatment types (2 × 2 factorial of two CO₂ treatments and two irrigation levels). CO₂ treatments were ambient (aCO₂—no CO₂ added, approximately 400 ppm) and elevated (eCO₂—target of 800 ppm), with eCO₂ treatment maintained by the system opening or closing solenoid valves to control delivery of CO₂ into the outside air stream. Irrigation treatments were well-watered (WW, volumetric water content target of 20–25%) and drought (DRT, VWC 10–15%). Irrigation levels were controlled by manually powering valves to deliver water as needed to maintain treatment VWC target levels. CO₂ treatments were started on June 28, 15 days after cuttings were potted. By this time, all cuttings were leafing out. Drought treatment was allowed to progress gradually, beginning June 28.

Fertilizer (Plant Prod “Forestry Special” 20:8:20 delivered at 100 ppm N) was supplied to all plants via the dripper system on only two dates in 2017: 5 July (30 min = 1 L) and 24 July (15 min = 500 mL). Treatments ran until 31 October 2017, when above-ground portions of all plants were harvested. For the winter of 2017–2018, CO₂ enhancement was stopped, and soil in all bags was kept well watered. Temperatures were kept at a minimum, just above 0 °C, to avoid freezing greenhouse infrastructure, and blackout curtains were used to prevent solar heating.

2.2. Assimilation Efficiencies and Gas Exchange Experiment: 2017

From 12–26 September 2017 (days 76–90 of treatment), GE measurements were performed on two clones per species (Table 1. Species and clones of Salix used in 2017 CO₂ × soil moisture and 2018 CO₂ × fertilization experiments) in situ using a LI-COR 6400XT portable GE system with a CO₂ mixer and RGB (Red Green Blue) LED light source on a 2 × 3 cm leaf cuvette (LI-COR, Lincoln, Nebraska, USA). On each measurement day, approximately eight assimilation to internal CO₂ (AC_i) response curves were completed by sampling one ramet of one genotype from each of two species from one replicate chamber per treatment. This pattern was repeated for the other two species on the following day. Species pairs were changed with successive paired dates. Over four such pairs of days, this resulted in all 8 genotypes being measured in two replicate chambers for each of the four treatment combinations (64 curves). Instrument settings were set for leaf temperatures of 25–28 °C, air flow of 250–400 umol·s⁻¹, and a saturating photosynthetically active radiation (PAR) level of 1000 umol·m⁻²·s⁻¹. Relative humidity (RH) averaged 50%, and was adjusted to maintain a vapour pressure deficit (VPD) of approximately 1.75 kPa in the chamber. Gas exchange was measured 13 times at 11 different CO₂ levels in the following sequence: 400, 330, 260, 190, 120, 50, 400, 400, 600, 800, 1200, 1700, and 2200 ppm. Upon initially securing the leaf in the chamber, no measurements were taken for 15–30 min, while the leaf equilibrated to the conditions in the chamber. Equilibration was judged by stability in the rate of net photosynthesis. Each subsequent measurement had a settling time of at least 2 min after the new CO₂ was achieved. Two consecutive measurements were performed at 400 ppm following the low reading at 50 ppm, to allow for the slower acclimation following an increase in [CO₂]. Note that for ACi parameter determinations, values from all three measurements at 400 ppm were averaged to represent that CO₂ concentration level.

The ACi response curves were fitted using commercially available software (Photosyn Assistant; Dundee Scientific, Dundee, Scotland, UK). This program uses algorithms based on the biochemical model of leaf photosynthesis developed by [22], and modified by [23]. Using this model, CO₂ assimilation in the leaf can be mathematically described by saturation kinetics, which account for the change in assimilation (A) rate at increasing internal CO₂ (Ci). Through an iteration procedure, the Photosyn Assistant software provides estimates for the maximum rate of carboxylation (V_cmax), light saturation rate, maximum rate of electron transport (J_max), and phosphorus (P) limited rate of triose phosphate utilization (TPU). The model calculations are based on the assumption that one of these parameters is at a maximum and limits A. In the program, a least-squares fit is used to estimate the CO₂ compensation point.

In addition, GE parameters were examined at a common 600 ppm CO₂, the near mid-point of the two CO₂ treatments, and thus A at 600 ppm CO₂ is abbreviated as A₆₀₀. Other GE parameters examined were stomatal conductance at 600 ppm CO₂ (G_wv600) and water use efficiency at 600 ppm CO₂ (WUE₆₀₀), determined by dividing A₆₀₀ by G_wv600.

2.3. CO₂ × Fertilizer Experiment: 2018

On 23 April 2018, the blackout curtains were opened and CO₂ treatments began with the same methods and targets as in 2017. On the same day, differential watering regimes began to re-establish the two irrigation treatments. Irrigation and CO₂ treatments were continued until 10 September 2018, when all plants were harvested.

On 19 July 2018, sixteen plants from each of the four well-watered chambers were selected; two ramets from each of two clones (genotypes) per species. Paired ramets in each chamber were selected to be similar in size. These 64 plants had their root bags briefly lifted and replaced inside plastic bags, to isolate their root systems and contain water delivered to the root bags. In each chamber, one ramet of each genotype was randomly assigned to the fertilized treatment (FERT), the other to non-fertilized (NONFERT). Watering spigots were removed from all FERT bags. For the remainder of the experiment, NONFERT plants received the same watering regime as all other plants via spigots, while FERT plants were hand-watered with an equal volume of water (typically 1 L two times per week) containing fertilizer (Plant Prod “Forestry Special” 20:8:20 delivered at 100 ppm N). When water was seen to accumulate in the NONFERT plastic bags, they were all pierced at the base to allow water to drain. Water accumulation did not occur with the FERT plants.

ACi curves were measured on 10 days between 7 and 20 August 2018 (days 18–31 of fertilizer treatment) on all 64 plants in the fertilizer experiment. Species, clones, and treatments were sampled randomly on each day. Methods were as in 2017, except PAR was set to 500 umol·m⁻²·s⁻¹ and leaf temperature ranged from 24.5 to 30 °C. Although light levels were different between years, both were above the light saturation point for individual leaves. Eight bulk soil samples (one bulked from FERT and one from NONFERT pots for each chamber) were collected in early September 2018 from NONFERT and FERT treatments and sent to the University of New Brunswick nutrient laboratory for analyses. The soil analysis procedure is the same as used in [24].

2.4. Statistical Analysis

This study was established as a completely randomized block experimental design. Willow species, CO₂ level, and soil water treatments were all considered as fixed effects. The data were subjected to analyses of variance (ANOVA) using the following ANOVA model:

Y_ijklm = μ + B_i + S_j +C_k + W_l + SC_jk + SW_jl + CW_kl + SCW_jkl + e_ijklm

(1)

where Y_ijklm is the dependent seedling trait of the i^th greenhouse chamber, of the j^th willow species, of the k^th CO₂ treatment, of the l^th water treatment, of the m^th stem cutting, and μ is the overall mean. B_i is the effect of the i^th greenhouse chamber (i = 1, 2), S_j is the effect of the j^th species (j = 1, 2, 3, 4), C_k is the effect of the k^th CO₂ treatment (k = 1, 2), W_l is the effect of the l^th water treatment (l = 1, 2), SC_jk is the interaction effect of the j^th willow species and the k^th CO₂ treatment, SW_jl is the interaction effect of the j^th willow species and l^th water treatment, CW_kl is the interaction effect of the k^th CO₂ treatment and the l^th water treatment, SCW_jkl is the interaction effect of the j^th willow species, k^th CO₂ treatment and l^th water treatment, and e_ijklm is the random error component. Effects were considered statistically significant at the p = 0.05 level, although individual p values are provided so that readers can make their own interpretations. The data satisfied normality and equality of variance assumptions. The general linear model from Systat (Chicago, IL, USA) was used for analysis. Willow species source of variation, if significant was post-hoc tested using the Tukey mean separation test, was (p = 0.05). The variance component analysis in the ANOVA tables was carried out using the sum of squares as outlined by [25]. In the second experiment, the same ANOVA model was used, but water (W) was replaced with FERT (F).

Covariate analysis was used to evaluate the relationships among mean willow ACi traits and among GE traits and to test treatment effects. In these analyses, the dependent trait (i.e., A₆₀₀) is examined in relation to three sources of variation were studied: (1) covariate (i.e., G_wv600), (2) independent effect (i.e., CO₂ treatment), and (3) independent effect × covariate (i.e., CO₂ treatment × G_wv600). The analyses were done based on the following model:

Y_ij = B₀ + B_0i + B₁X_ij + B_1iX_ij + e_ij

(2)

where Y_ij is the dependent trait of the i^th species of the j^th genus treatment. B₀ and B₁ are average regression coefficients, B_0i and B_1i are the treatment specific coefficients, X_ij is the independent variable, and e_ij is the error term. Results were considered statistically significant at p = 0.10 due to fewer data points, although individual p values are provided for all traits so that readers can make their own interpretations of significance.

3. Results

3.1. 2017. CO₂ × Water Experiment

For V_cmax, species, CO₂, water effects, and species × water interaction were significant accounting for 33.3, 5.8, 4.3, and 7.9% of the total variation (

Table 2. Assimilation to internal CO₂ traits and ANOVAs from the CO₂ × moisture stress treatment experiment, 2017, including source of variation, degrees of freedom (df), mean square values (MS), variance components (VC), p-values, and coefficient of determination (R²). p-values < 0.05 are in bold print.

Source of Variation	df	Vcmax (µmol·m−2·s−1)			Jmax (µmol·m−2·s−1)			TPU (µmol·m−2·s−1)
		MS	VC (%)	p-Value	MS	VC (%)	p-Value	MS	VC (%)	p-Value
Block	1	97.1	1.0	0.310	483.0	0.6	0.427	6.80	2.5	0.094
Species	3	1119.7	33.3	<0.001	12059.5	42.2	<0.001	30.54	33.4	<0.001
CO2	1	580.6	5.8	0.016	1397.9	1.6	0.179	5.62	2.0	0.127
Water	1	431.2	4.3	0.036	1461.2	1.7	0.170	29.22	10.6	0.001
CO2 × species	3	88.0	2.6	0.422	729.1	2.6	0.414	4.26	4.6	0.155
Water × species	3	265.3	7.9	0.046	2519.2	8.8	0.027	4.43	4.8	0.142
Water × CO2	1	<0.1	<0.1	0.982	159.3	0.2	0.647	0.03	<0.1	0.909
Water × CO2 × species	3	79.4	2.4	0.468	361.0	1.3	0.697	1.90	2.1	0.491
Error	47	92.2	42.9		750.9	41.1		2.33	39.9
R2				0.571			0.589			0.601

). INT had the greatest V_cmax followed by COR and ERI, with DIS having the lowest values at 44.9, 35.0, 31.4, and 24.9 μmol m⁻² s⁻¹, respectively. All four willow species decreased V_cmax in response to eCO₂ (Figure 1A). Three species had lower V_cmax in response to DRT: COR, DIS, and ERI (Figure 1B). The significant species × water interaction resulted from INT having greater V_cmax in response to drought. For J_max, the species effect and species × water interaction were significant, accounting for 42.2 and 8.8% of the total variation (Table 2). INT had J_max of 133 μmol m⁻² s⁻¹, which was significantly greater than the other three species that were not significantly different and averaged 82 μmol m⁻² s⁻¹ (not shown). Similar to V_cmax, the species × water interaction was due to INT, which had greater J_max in response to DRT, whereas the other three species had lower J_max in response to DRT. For TPU, species and water effects were significant, accounting for 33.4 and 10.6% of the total variation (Table 2). INT had the greatest TPU with 8.3 μmol m⁻² s⁻¹ and the other three were not significantly different and averaged 5.3 μmol m⁻² s⁻¹ (not shown).

For A₆₀₀, the three main effects, species, water, and CO₂, were significant accounting for 23.3, 20.6, and 4.9% of the total variation (

Table 3. Gas exchange traits and ANOVAs from the CO₂ × moisture stress treatment experiment, 2017, including source of variation, degrees of freedom (df), mean square values (MS), variance components (VC), p-values, and coefficient of determination (R²). p-values <0.05 are in bold print.

Source of Variation	df	A600 (µmol·m−2·s−1)			Gwv600 (mmol·m−2·s−1)			WUE600
		MS	VC (%)	p-Value	MS	VC (%)	p-Value	MS	VC (%)	p-Value
Block	1	46.0	2.0	0.148	0.10	6.8	0.001	2570.9	1.1	0.126
Species	3	179.1	23.3	<0.001	0.087	9.8	0.001	5307.6	7.1	0.004
CO2	1	113.0	4.9	0.026	0.005	0.2	0.560	126.5	0.1	0.731
Water	1	474.8	20.6	<0.001	1.422	53.5	<0.001	140223.2	62.3	<0.001
CO2 × species	3	18.6	2.4	0.461	0.012	1.4	0.452	1148.3	1.5	0.365
Water × species	3	15.7	2.0	0.536	0.010	1.1	0.555	1964.8	2.6	0.151
Water × CO2	1	4.6	0.2	0.644	0.041	1.5	0.094	3310.5	1.5	0.084
Water × CO2 × species	3	8.6	1.1	0.751	0.009	1.0	0.572	1209.5	1.6	0.342
Error	47	21.3	43.4		0.014	24.7		1060.3	22.2
R2				0.566			0.753			0.778

). All species of A₆₀₀ were lower in eCO₂ and declined on average from 14.8 to 12.1 μmol m⁻² s⁻¹ or 18% (Figure 2A). All species of A₆₀₀ were lower in DRT than WW treatments with 10.7 and 16.2, μmol m⁻² s⁻¹, respectively, a decline of 33% (Figure 2B). INT had the greatest A₆₀₀ with 18.1 μmol m⁻² s⁻¹ and the other three species were not significantly different and averaged 11.9 μmol m⁻² s⁻¹, a 34% difference. For G_wv600, species and water effects were significant, accounting for 9.8 and 53.5 of the total variation (Table 3). There were no consistent eCO₂ effects on G_wv600 (Figure 3A); whereas drought had a strong effect in lowering the G_wv600 on average from 0.40 to 0.10 mmol m⁻²s⁻¹, a 75% decline (Figure 3B). INT had the greatest G_wv600 with 0.36 mmol m^{−2 s−1} and the other three were not significantly different and averaged 0.21 mmol m⁻² s⁻¹. For WUE₆₀₀, species and water effects were significant, accounting for 7.1 and 62.3% of the total variation (Table 3). There were no consistent eCO₂ effects on WUE₆₀₀ (not shown), whereas drought had a strong effect in increasing the WUE₆₀₀ from on average 49.7 to 143.4 or a 2.9 fold increase (not shown). COR had the greatest WUE₆₀₀ followed by ERI and DIS, while INT had the lowest WUE₆₀₀.

3.2. 2018. CO₂ × Fertilizer Experiment

There were no soil nutrient differences between FERT and NOFERT for organic matter, C, N, CN ratio, Ca, Na, and S (

Table 4. Soil properties (mean ± SE) during the 2018 NOFERT and FERT treatments. Treatments with different letters are significantly different using an ANOVA test, α = 0.05.

Treatment	Organic Matter (%)	Carbon (%)	Nitrogen (%)	Potassium (meq/100 g)	Calcium (meq/100 g)	Magnesium (meq/100 g)	Phosphorus (ppm)
NOFERT	0.67 ± 0.40 a	0.39± 0.02 a	0.085 ± 0.009 a	0.04 ± 0.01 b	1.21 ± 0.08 a	0.050 ± 0.004 b	6.50 ± 0.50 b
FERT	0.63 ± 0.40 a	0.37 ± 0.02 a	0.090 ± 0.009 a	0.16 ± 0.01 a	1.08 ± 0.08 a	0.068 ± 0.004 a	9.50 ± 0.50 a

and

Table 5. More soil properties (mean ± SE) during the 2018 NOFERT and FERT treatments. Treatments with different letters are significantly different using an ANOVA test, α = 0.05.

Site	Sodium (%)	Sulfur (%)	pH	C:N Ratio
NOFERT	0.11 ± 0.01 a	0.020 ± 0.007 a	6.8 ± 0.1 a	4.7 ± 0.2 a
FERT	0.11 ± 0.01 a	0.013 ± 0.007 a	6.2 ± 0.1 b	3.9 ± 0.2 a

). FERT soil nutrient traits were greater than NOFERT for K and Mg, and P. NOFERT had higher pH than FERT with 6.8 and 6.2, respectively. For V_cmax, species, CO₂, FERT, and CO₂ × FERT interaction were significant, accounting for 24.1, 1.5, 61.5, and 1.6% of the total variation (

Table 6. Assimilation to internal CO₂ traits and ANOVAs from the CO₂ × fertilizer treatment experiment, 2018, including source of variation, degrees of freedom (df), mean square values (MS), variance components (VC), p-values, and coefficient of determination (R²). p-values < 0.05 are in bold print.

Source of Variation	df	Vcmax (µmol·m−2·s−1)			Jmax (µmol·m−2·s−1)			TPU (µmol·m−2·s−1)
		MS	VC (%)	p-Value	MS	VC (%)	p-Value	MS	VC (%)	p-Value
Block	1	70.5	0.4	0.160	156.0	0.2	0.424	0.80	0.1	0.746
Species	3	1366.7	24.1	<0.001	7415.7	28.1	<0.001	19.58	8.7	0.063
CO2	1	255.5	1.5	0.009	1.0	<0.1	0.948	2.04	0.3	0.605
FERT	1	10471.1	61.5	<0.001	43202.1	54.5	<0.001	102.87	15.2	0.001
CO2 × species	3	6.5	0.1	0.904	377.1	1.4	0.208	21.36	9.4	0.048
FERT × species	3	26.9	0.5	0.514	150.8	0.6	0.599	21.10	9.3	0.050
CO2 × FERT	1	274.0	1.6	0.007	202.6	0.3	0.362	0.45	0.1	0.809
CO2 × FERT × species	3	41.5	0.7	0.321	195.1	0.7	0.492	10.71	4.7	0.248
Error	47	34.7	9.6		239.5	14.2		7.52	52.2
R2				0.904			0.858			0.478

). INT had the greatest V_cmax, followed by COR, ERI, and DIS (Figure 1C). V_cmax was consistently lower for all four species in eCO₂. FERT increased V_cmax from 29.6 to 55.2 umol m⁻²s⁻¹ for an increase of 86% (Figure 1D). The CO₂ × FERT interaction was a result of a V_cmax downregulation under NOFERT, but upregulation under FERT treatment. For J_max, species and FERT were significant, accounting for 28.1 and 54.5% of the total variation (Table 6). INT had the greatest J_max followed by COR, ERI, and DIS. There were no consistent CO₂ effects on J_max, but FERT showed a strongly increased average J_max effect of 66 and 118 umol m⁻²s⁻¹ for NOFERT and FERT treatments, respectively, or a 79% increase (not shown). For TPU, the effects of FERT, species × CO₂, and species × FERT interactions were significant, accounting for 15.2, 9.4, and 9.3% of the total variation (Table 6). The significant species × CO₂ interaction was a result of no consistent CO₂ effect on willow species (not shown). The FERT and species × FERT TPU interaction was the result of a FERT TPU greater than NOFERT for COR, ERI, and INT; however, DIS showed no TPU difference between FERT and NOFERT.

For A₆₀₀, species, CO₂, FERT, and the CO₂ × FERT interaction were significant, accounting for 20.8, 1.0, 65.0, and 1.1% of the total variation (

Table 7. Gas exchange traits and ANOVAs from the CO₂ x fertilizer treatment experiment, 2018, including source of variation, degrees of freedom (df), mean square values (MS), variance components (VC), p-values, and coefficient of determination (R²). p-values <0.05 are in bold print.

Source of Variation	df	A600 (µmol·m−2·s−1)			Gwv600 (mmol·m−2·s−1)			WUE600
		MS	VC (%)	p-Value	MS	VC (%)	p-Value	MS	VC (%)	p-Value
Block	1	3.39	0.2	0.411	0.001	0.1	0.743	321.3	0.6	0.448
Species	3	155.35	20.8	<0.001	0.064	14.3	<0.001	3292.5	19.3	0.002
CO2	1	22.81	1.0	0.036	0.004	0.3	0.452	442.2	0.9	0.374
FERT	1	1454.81	65.0	<0.001	0.588	43.5	<0.001	3373.2	6.6	0.017
CO2 × species	3	0.84	0.1	0.916	0.023	5.0	0.042	1871.2	11.0	0.025
FERT × species	3	0.91	0.1	0.906	0.015	3.4	0.125	937.4	5.5	0.178
CO2 × FERT	1	24.90	1.1	0.029	0.013	1.0	0.201	0.3	<0.1	0.983
CO2 × FERT × species	3	9.71	1.3	0.131	0.027	5.9	0.022	998.7	5.8	0.156
Error	47	4.92	10.3		0.008	26.5		548.5	50.3
R2				0.897			0.735			0.497

). INT had the greatest A₆₀₀, followed by COR, ERI, and DIS (Figure 2C). In response to eCO₂, A₆₀₀ was consistently lower under aCO₂ for all four species. FERT increased A₆₀₀ from 10.4 to 19.9 μmol m⁻² s⁻¹ or 91% (Figure 2D). The CO₂ × FERT interaction was a result of an A₆₀₀ downregulation under NOFERT, but an upregulation under the FERT treatment. For G_wv600, the species, FERT, species × CO_2, and species × FERT × CO₂ interactions were significant accounting for 14.3, 43.5, 5.0, and 5.9% of total variation, respectively (Table 7). The G_wv600 species × FERT × CO₂ interaction resulted from a varied species response to FERT and eCO₂ (Figure 3C,D). The species × CO₂ interaction is due to rank changes, with COR and ERI decreasing in G_wv600 response to eCO₂, while DIS remained the same and INT increased G_wv600 in response to eCO₂. INT and DIS seem to respond more strongly than COR and ERI to FERT (Figure 3D). Overall, this is seen in the species response: INT and DIS having the greatest G_wv600 response at 0.34 and 0.27 mmol m⁻²s⁻¹ and ERI and COR at 0.22 and 0.20 mmol m⁻²s⁻¹, respectively (Figure 3C). The overall FERT effect was 0.161 and 0.352, for NOFERT and FERT, respectively, a 2.2 fold increase.

For 2018 WUE₆₀₀, species and FERT effects and species × CO₂ interaction were significant, accounting for 19.3, 6.6, and 11.0% of the total variation (Table 7). The species × CO₂ interaction was a result of COR, DIS, and ERI increasing WUE₆₀₀ in response to eCO₂, whereas, INT decreased its WUE₆₀₀ (Figure 4A). The FERT effect was an overall decrease in WUE₆₀₀ in response to FERT from 77.3 to 62.8 or an 18.8% decrease (Figure 4B). COR had the greatest WUE₆₀₀ followed by ERI and INT, with DIS having the lowest WUE₆₀₀.

3.3. Covariate Relationships among Willow Species and Environmental Change

Covariate analysis of J_max in relation to V_cmax, testing CO₂ effect, had no significant CO₂ × V_cmax interaction (p = 0.752). There was a significant J_max to V_cmax relationship (p < 0.001) and CO₂ effect (p = 0.097) (Figure 5A). Covariate analysis of J_max in relation to V_cmax testing water effect, had no significant water × V_cmax interaction (p = 0.927) nor a water effect (p = 0.225). There was a significant J_max to V_cmax relationship (p < 0.001) (Figure 5 B). Covariate analysis of J_max in relation to V_cmax testing the FERT effect, had no significant FERT × V_cmax interaction (p = 0.675) nor a FERT effect (p = 0.114). There was a significant J_max to V_cmax relationship (p < 0.001) (Figure 5C).

Covariate analysis of A₆₀₀ in relation to G_wv600 testing CO₂ effect, had no significant CO₂ × G_wv600 interaction (p = 0.519) nor a CO₂ effect (p = 0.171). There was a significant A₆₀₀ to G_wv600 relationship (p = 0.006) (Figure 6A). Covariate analysis of A₆₀₀ in relation to G_wv600 testing water effect had a significant water × G_wv600 interaction (p = 0.057). Thus, there were two separate slopes for WW and DRT treatments (Figure 6B). Covariate analysis of A₆₀₀ in relation to G_wv600 testing FERT effect, had no significant FERT × G_wv600 interaction (p = 0.509) nor a FERT effect (p = 0.249). There was a significant A₆₀₀ to G_wv600 relationship (p = 0.004) (Figure 6C).

Covariate analysis of WUE₆₀₀ in relation to A₆₀₀ testing CO₂ effect showed no significant CO₂ × A₆₀₀ interaction (p = 0.566) nor a CO₂ effect (p = 0.492). There was a significant WUE₆₀₀ to A₆₀₀ relationship (p = 0.049, not shown). Covariate analysis of WUE₆₀₀ in relation to A₆₀₀ testing water effect, had no significant water × A₆₀₀ interaction (p = 0.705), but a significant water effect (p = 0.007). This resulted in two separate lines of the same slope (Figure 7A). Covariate analysis of WUE₆₀₀ in relation to A₆₀₀ testing FERT effect, had no significant FERT × A₆₀₀ interaction (p = 0.370) nor a FERT effect (p = 0.183). There was no significant WUE₆₀₀ to G_wv600 relationship (p = 0.615, not shown).

Covariate analysis of WUE₆₀₀ in relation to G_wv600 testing CO₂ effect, had no significant CO₂ × G_wv600 interaction (p = 0.854) nor a CO₂ effect (p = 0.903). There was a significant WUE₆₀₀ to G_wv600 relationship (p = 0.005, not shown). Covariate analysis of WUE₆₀₀ in relation to G_wv600 testing water effect, showed neither a significant water × G_wv600 interaction (p = 0.120) nor a water effect (p = 0.273). There was a significant WUE₆₀₀ to G_wv600 relationship in response to DRT (p < 0.001) (Figure 7B). Covariate analysis of WUE₆₀₀ in relation to G_wv600 testing FERT effect, showed neither a significant FERT × G_wv600 interaction (p = 0.448) nor a FERT effect (p = 0.752). There was no significant WUE₆₀₀ to G_wv600 relationship (p = 0.103, not shown).

4. Discussion

There were a number of consistent findings for biochemical efficiency and GE traits of these willow species between the two experiments (years). First, there was A_dr for all four willow species in both experiments in response to eCO₂. Second, biochemical efficiency parameters decreased for all willow species in the 2017 experiment, with reductions in V_cmax, J_max, and TPU in response to CO₂. In the second (2018) experiment, there was a reduction in V_cmax for all four willow species in response to eCO₂. However, for J_max and TPU, there was no decrease, but this was probably due to the mitigation effect from fertilization, due to an increase in sink (growth) activity. Third, the willow species effect was often the most important factor accounting for between 33 and 42% of the total variation among the three traits (V_cmax, J_max, and TPU) in the eCO₂ × DRT experiment and between 20 and 28% for the same three traits in the CO₂ × FERT experiment. Fourth, in both experiments, INT had significantly greater biochemical efficiency and GE values followed by COR and ERI, with DIS having the lowest values. DIS had the lowest biomass [26] among the four species, consistent with its low biochemical efficiency values. COR had the greatest biomass, consistent with its high biochemical efficiency values. INT had intermediate biomass results, but had among the lowest leaf dry mass. In an earlier study examining three salinity treatments including control, INT had significantly greater V_cmax and J_max and GE than DIS and ERI, which were similar [15,16].

The significant species × water interaction was due to rank changes, for example, V_cmax was lower under DRT for COR, DIS, and ERI, but 16.5% greater for INT under DRT compared to WW. Drought response has a similar physiological effect to a saline response in that they both induce moisture stress in plants and both require cellular osmotic adjustment to withstand the stress [27]. Both stress responses require additional energy from the plant, and interestingly, INT increased its biochemical efficiency by 30% in medium salinity and high salinity, and it was the only willow species of the three tested to survive high salinity [15,16]. The greater drought and saline tolerance of INT may reflect its evolutionary origins in the arid SW USA and Mexico [28,29,30], where high evapotranspiration may increase drought and soil salinity (see [31] and references therein) and result in a natural selection for increased drought and saline tolerance. The exact cause of the greater INT drought tolerance could be due to changes to the hydraulic architecture of the xylem and thus greater hydraulic conductivity [32] or as a result of greater osmotic adjustment [33,34].

Assimilation downregulation is the decrease in assimilation with time in response to eCO₂ quantified at the same measurement of CO₂ concentration [35,36]. There are generally at least three plausible causes recognized for A_dr. First, significant increases in CO₂ concentration often results in a decline in underlying photosynthetic capacities, carboxylation efficiencies, electron transport, and chlorophyll content [4,37,38], which is what we report on here. Second, the accumulation of nonstructural carbohydrates is usually associated with A_dr and is commonly interpreted as evidence of a lack of sink activity [36,39,40,41]. We observed a lack of sink activity in the present study, as there was no increase in stem dry mass in eCO₂ [26]. An increase in nonstructural carbohydrates is most likely the driver of the decline in underlying biochemical capacity. Third, it has been observed in eCO₂ conditions and on relatively infertile sites that limiting essential nutrients can cause A_dr by reducing growth [5,42].

Our preliminary data showed no CO₂ effect on willow stem dry mass after the first year of growth in 2017 (p = 0.214) [26,42] found that fertilization of the eCO₂ plots restored eCO₂ growth in forest species, and [5] found fertilization reversed A_dr. Thus, we decided to do a fertilizer × CO₂ experiment in 2018 on a subset of fully irrigated treatments to see if we could observe a CO₂ effect under FERT. First, we did see a sizable eCO₂ effect, resulting in significant increases in biochemical efficiencies and GE (A₆₀₀ and G_wv600). This change resulted in a mitigation of the biochemical downregulation found in 2017 in response to eCO₂. This was doubtless driven by the FERT effect on biomass accumulation, which resulted in a significant (p = 0.049) increase in stem dry mass and CO₂ × fertilizer interaction [26]. As in 2017, there was no stem dry mass difference between aCO₂ and eCO₂ under NOFERT in 2018, with 11.8 and 10.7 g, respectively. However, under the FERT treatment, stem dry mass was 20.0 and 25.3 g under aCO₂ and eCO₂, respectively. Thus, nutrient limitation probably caused the lack of eCO₂ effect on willow growth. Subsequent soil analysis showed very low levels of nutrients under NOFERT. Also interesting was that a soil nutrient result after the FERT treatment in 2018 showed no difference between FERT and NOFERT, particularly for soil N. Thus, everything supplied by the FERT treatment was taken up by the FERT grown willows.

As mentioned, fertilization may have mitigated the biochemical downregulation in response to eCO₂, but eCO₂ itself did not mitigate the DRT downregulation. Using the same willow species on a nutrient-poor former coal mine site reclamation study, we observed how green the willow foliage was. Willows apparently have a symbiotic relationship with diazotrophic endophytes in the stem that provide N to the plant in exchange for carbohydrates [43,44]. In the present study, we provided only two pulses of fertilizer early in the first year, assuming that endophytic bacteria would be present under both aCO₂ and eCO₂ treatments. However, it is now apparent that N fixation was either too low or that not enough other essential soil nutrients were present to allow the plants to take advantage of the eCO₂ treatment.

The J_max to V_cmax species relationship under aCO₂ and eCO₂ showed the downregulation of both traits, particularly V_cmax, resulting in two separate lines. The INT species appear at the top of the two linear relationships and at a significant separation from the other species. The J_max to V_cmax species relationship under WW and DRT showed a single positive line, similar in slope to Figure 5A, but with a treatment segregation, the DRT treatments on the lower half, except for INT-D, which is at the very top of the relationship. Again, this is counter to what the other three willow species showed, but similar to the increased INT biochemical efficiency that was found in an earlier salinity study [15]. The J_max versus V_cmax relationship by species and FERT segregated primarily by treatment, though INT still showed the greatest biochemical efficiencies among the four willows under each FERT treatment. We also saw a flattening of the J_max to V_cmax relationship by one third. This may indicate a differential N allocation between the thylakoid and soluble proteins [40] or an N reallocation between photosynthetic components [5,40] examining Picea sitchensis and Pinus taeda, respectively, under differing levels of N × CO₂, did not observe a change in the J_max:V_cmax ratio, but our soils were clearly impoverished. A review by [45] showed that these two parameters are largely species-specific, which is what we observed among the four willow species under three different treatment scenarios.

Is WUE driven by changes in A or G_wv or both? It depends on the experiment. First, the A₆₀₀ to G_wv600 species relationship under different CO₂ treatments showed a single positive line with a fair amount of variation (Figure 6A). All the eCO₂ willows were below the single diagonal line, indicating that A₆₀₀ decreased more in response to eCO₂, and G_wv600 remained largely unchanged in both experiments. However, WUE in the willows did not change in response to eCO₂. The literature is mixed on the strength of the effect of eCO₂ on WUE: Pinus radiata had lower WUE in eCO₂ [3] and both pines and spruces had lower WUE in eCO₂ [6]. In contrast, [46] found WUE was greater for Picea abies and Quercus rubra grown at eCO₂, and this was also the case for loblolly pine (Pinus taeda L.) [47]. In a study of three conifer species, Picea mariana and Larix laricina displayed a modest reduction in G_wv600 in response to eCO₂, whereas there was no response seen in Pinus banksiana [39] or P. radiata [3]. In one of the only publications on willow (S. sagitta) [17] response to CO₂ and moisture stress is suggested that eCO₂ mitigated the drought stress response in willows but not poplars. We did not find that eCO₂ mitigated the response to drought stress.

In the FERT experiment, fertilization resulted in a strong increase in A and G_wv of 92 and 119%, which accounted for 44 and 65% of the total variation, respectively. This would be in response to the greater sink demand (growth) stimulated by the FERT treatment and the need for more assimilates. Fertilization resulted in a small, but significant decrease in WUE, but this may have been affected by very large non-specific increases in A and G_wv. The A₆₀₀:Gwv₆₀₀ species slope was the same as the CO₂ species slope, but was much further extended. Nevertheless, there were no relationships of either A or G_wv with WUE, thus neither were driving the FERT WUE results from this experiment.

DRT increased the A₆₀₀:Gwv₆₀₀ species slope by a factor of three compared to the WW species slope, due mostly to the changes in G_wv induced by drought. Most species show a reduction of G_wv to control water loss during dry periods [48]. In addition, there was a four-fold reduction in G_wv600 and the average for the four willows and water effect accounted for 54% of the G_wv total variation, dwarfing all other sources of variation. In contrast, A₆₀₀ was reduced by only 50% by DRT and the water effect only accounted for 21% of the A₆₀₀ total variation. In addition, WUE in relation to G_wv by species and water treatment showed a strong negative relationship (R² = 0.917). Thus, G_wv was the primary driver and A₆₀₀ only had a minor effect on WUE when testing under DRT. While the A₆₀₀ to WUE₆₀₀ relationship within water treatments (Figure 7A) might at first glance suggest that greater A₆₀₀ results in a lower WUE, this is all driven by one species, INT. INT had lower WUE, but had greater A₆₀₀ and G_wv600, showing more evidence of a strategy of increased drought tolerance and fits with the adaptive selection pressure from its evolutionary origins in the arid SW USA and Mexico [28,29,30]. Even in the DRT experiment, INT had greater G_wv than the other three willows, pointing to it being more of a drought tolerator (keeping G_wv more open to continue to allow A) than a drought avoider (closing G_wv).

5. Conclusions

Willows are used as bioenergy crops often on marginal agricultural lands and for land restoration of severely disturbed sites; both site types are often nutrient-poor. Thus, with increasing atmospheric CO₂ it would appear that if a site is nutrient-poor and there are no stem endophytic bacteria, the addition of soil nutrients may be necessary to sequester more available CO₂. Based on physiological and morphological results, INT is the most drought and saline tolerant of the willows tested. INT and its closest relatives in the section Longifoliae are also noteworthy among willow species for their ability to easily form new stems via stem suckering from a network of shallow roots, resulting in multi-stemmed colonies [49,50]. This stem-suckering feature and the drought tolerance of INT allows more CO₂ capture to help mitigate climate change impacts. The rapid growth, high biomass yields, and ease of vegetative (clonal) propagation of INT from unrooted stem sections should be of interest for establishing biomass feedstock for bioenergy, and other resources for the bioeconomy, as well as for land reclamation activities on marginal land affected by drought or salinity.