*Article* **Establishment of Regional Phytoremediation Buffer Systems for Ecological Restoration in the Great Lakes Basin, USA. I. Genotype** × **Environment Interactions**

**Ronald S. Zalesny, Jr. 1,\* , Andrej Pilipovi´c <sup>2</sup> , Elizabeth R. Rogers 1,3 , Joel G. Burken <sup>4</sup> , Richard A. Hallett <sup>5</sup> , Chung-Ho Lin <sup>3</sup> , Bernard G. McMahon <sup>6</sup> , Neil D. Nelson <sup>6</sup> , Adam H. Wiese <sup>1</sup> , Edmund O. Bauer <sup>1</sup> , Larry Buechel <sup>7</sup> , Brent S. DeBauche <sup>3</sup> , Mike Peterson <sup>7</sup> , Ray Seegers <sup>8</sup> and Ryan A. Vinhal 1,3**


**Abstract:** Poplar remediation systems are ideal for reducing runoff, cleaning groundwater, and delivering ecosystem services to the North American Great Lakes and globally. We used phytorecurrent selection (PRS) to establish sixteen phytoremediation buffer systems (phyto buffers) (buffer groups: 2017 × 6; 2018 × 5; 2019 × 5) throughout the Lake Superior and Lake Michigan watersheds comprised of twelve PRS-selected clones each year. We tested for differences in genotypes, environments, and their interactions for health, height, diameter, and volume from ages one to four years. All trees had optimal health. Mean first-, second-, and third-year volume ranged from <sup>71</sup> <sup>±</sup> <sup>26</sup> to <sup>132</sup> <sup>±</sup> 39 cm<sup>3</sup> ; 1440 <sup>±</sup> 575 to 5765 <sup>±</sup> 1132 cm<sup>3</sup> ; and 8826 <sup>±</sup> 2646 to 10,530 <sup>±</sup> 2110 cm<sup>3</sup> , respectively. Fourth-year mean annual increment of 2017 buffer group trees ranged from 1.1 ± 0.7 to 7.8 <sup>±</sup> 0.5 Mg ha−<sup>1</sup> yr−<sup>1</sup> . We identified generalist varieties with superior establishment across a broad range of buffers ('DM114', 'NC14106', '99038022', '99059016') and specialist clones uniquely adapted to local soil and climate conditions ('7300502', 'DN5', 'DN34', 'DN177', 'NM2', 'NM5', 'NM6'). Using generalists and specialists enhances the potential for phytoremediation best management practices that are geographically robust, being regionally designed yet globally relevant.

**Keywords:** ecosystem services; multi-environmental trials (MET); phenotypic plasticity; phyto buffers; phyto-recurrent selection; phytotechnologies; poplars; *Populus*

## **1. Introduction**

Ninety-five percent of the United States' surface freshwater and 20% of the world's freshwater reserve are contained within the Great Lakes Basin [1,2]. The Basin provides ecosystem services (including clean drinking water) to over 34 million people in the United States and Canada [3–5], and contributes substantially to the United States' economy, with a gross regional product (GRP) estimated at nearly 4.1 trillion USD [6]. Anthropogenic

**Citation:** Zalesny, R.S., Jr.; Pilipovi´c, A.; Rogers, E.R.; Burken, J.G.; Hallett, R.A.; Lin, C.-H.; McMahon, B.G.; Nelson, N.D.; Wiese, A.H.; Bauer, E.O.; et al. Establishment of Regional Phytoremediation Buffer Systems for Ecological Restoration in the Great Lakes Basin, USA. I. Genotype × Environment Interactions. *Forests* **2021**, *12*, 430. https://doi.org/ 10.3390/f12040430

Academic Editors: Francisca C. Aguiar

Received: 6 March 2021 Accepted: 30 March 2021 Published: 2 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

activities in the region, combined with population growth and land use changes, induce a vast number of stressors such as toxic pollution, increases in invasive species populations, and climate change [7], all of which can disturb terrestrial water cycles [8] and impact water quality of the Great Lakes and their watersheds [9,10]. Non-point sources of pollution such as landfills and similar sites contribute to watershed contamination by runoff and leakage. Although pollutant levels within landfills generally decrease over time through chemical and biological alteration and degradation [11–13], treatment of leachate and wastewaters can help mitigate soil and water contamination [14]. Sustainable, long-term restoration practices are needed to preserve and enhance ecosystem services in the Great Lakes Basin. In the early 2010s, over 1.5 billion USD were invested in restoration projects to halt or mitigate pollution of the Great Lakes ecosystem [15]. Ecological restoration techniques and technologies are more economically viable and sustainable long-term solutions than off-site treatment methods for mitigating contamination.

Ecological restoration may help reverse land degradation, improve the resilience of biodiversity, and deliver important ecosystem services [16,17]. Phytotechnologies are an ideal solution for preventing water contamination, and they include biological recovery activities that have been classified into four categories: remediation, reclamation, restoration, and rehabilitation [18]. The most common phytotechnology, phytoremediation, involves the use of plants for remediation and prevention of water contamination [19]. Often, phytoremediation is accomplished by installing vegetative covers, riparian buffers, and/or hydraulic control systems. Depending upon their type and chemical properties, contaminants are accumulated, immobilized, metabolized or volatilized [20–22].

Fast-growing trees such as *Populus* L. species and their hybrids (hereafter referred to as poplars) have been studied and used in plantation-based silvicultural systems known as short rotation woody crop (SRWC) systems for over 100 years [23,24]. The primary focus for SRWCs has been on the production of wood and wood products. More recently, poplar production systems have been implemented to provide a variety of ecosystem services [25–27], including their application in phytotechnologies for ecological clean up and restoration [28]. Phytoremediation research starting in the 1990s showed that poplars were great candidates for remediation systems because of their vigorous growth, easy vegetative propagation, well-developed root systems, and high productivity on marginal and liability lands [27,29,30]. The results of such studies indicate the potential of poplars for phytoremediation of various inorganic and organic contaminants [20,31–35]. Another factor that adds to the value of poplars in phytoremediation is that they exhibit a broad range of genetic diversity given their ability to undergo spontaneous and controlled intra- and inter-species hybridization, which thereby creates a high number of simple and complex hybrids [36].

Phytoremediation projects are often installed on sites and in soils that are not ideal for plant growth. Variability in environmental conditions affects evolutionary processes of populations or individuals within species, resulting in broadening or narrowing of their ecological niches, defined as generalists and specialists, respectively [37,38]. Poplars are ideal candidates for tree improvement by hybridization and breeding [39–41], with further aims to fully exploit poplar genomic and genetic potential and maximize wood biomass productivity from fast-growing trees worldwide. The ability of genotypes to produce different phenotypes in distinct environmental conditions is defined as phenotypic plasticity, and includes variation in the morphology, physiology, behavior, and life history of organisms [42]. This can lead to unpredictability in performance that can complicate selection of proper genotypes for phytoremediation projects. The genotype by environment (G × E) interaction, in which a phenotype is a function of the genotype, the environment, and the differential phenotypic response of genotypes to site-specific edaphic and climatic conditions, is a determining factor of clonal site performance [43]. Such interaction exists when comparative performances of genotypes vary according to local site conditions, with superiority of a genotype in a certain environment converted to inferiority in another environment [44]. Genotypic variability across a variety of environments can lead to unpredictability in performance that can complicate selection of proper genotypes for commercial use. Therefore, in tree breeding, multi-environmental trials (MET) are implemented to evaluate the degree and pattern of G × E interactions, as well as to test the robustness of genotypes in different environments [45,46].

However, sometimes phenotypic plasticity can hinder expression of G × E interactions. For example, testing of hybrid poplar clones for growth performance and robustness at different sites in the Midwestern USA showed low G × E interaction and high plasticity of tested genotypes, indicating their suitability for growing on a wider spectrum of habitats [46,47]. In poplar research, G × E interactions and definitions of generalists vs specialists have been studied for various traits ranging from biomass production [48–50] and rooting ability [51], to physiological traits related to productivity [36] and drought tolerance [52], to wood properties [53,54]. Often, G × E interactions are also determined by the genomic groups of the hybrids, where components of different *Populus* species (e.g., *P. deltoides* Bartr. ex. Marsh, *P. maximowiczii* A. Henry, *P. nigra* L., *P. trichocarpa* Torr. et Gray) significantly contribute to the expression of specific traits by a given genotype, and are reflected in rooting ability and adaptability to certain climate, latitude or soil-water conditions [48,49,55,56]. Therefore, clonal testing to maximize understanding of G × E interactions is paramount in establishing successful poplar plantings for a broad range of end uses and ecosystem services, including environmental objectives such as phytoremediation. For example, optimizing G × E interactions may be most appropriate for small-scale applications with site-specific needs, while minimum G × E interactions are needed to deploy superior and robust clones at a justifiable cost for large-scale commercial systems [46].

The current study is a component of a regional phytoremediation testing network originally funded by the Great Lakes Restoration Initiative (GLRI) [2] to select and test the ability of new poplar genotypes to reduce surface runoff and prevent groundwater contamination at landfills and similar sites in the Great Lakes Basin. We implemented phyto-recurrent selection, a stepwise greenhouse- and field-based approach that combines crop and tree improvement strategies to evaluate, identify and select superior-performing clones based on multiple testing cycles [32,57]. Results of the initial screening of candidate clones in greenhouse experiments, previously published by Rogers et al. [58], were used to develop the MET testing network consisting of sixteen phytoremediation buffer systems (i.e., phyto buffers) at ten sites located in the Lake Superior (i.e., Michigan's Upper Peninsula) and Lake Michigan (i.e., eastern Wisconsin) watersheds. Our objective was to test for differences in genotypes, environments, and their interactions for health, height, diameter, and volume during early field establishment (i.e., from one to four years after planting). Our results identify poplar clones with maximum phytoremediation potential that can be established across a wide range of environmental conditions. These data are useful for clonal selection to maximize phytoremediation potential in future phytotechnologies, regardless of specific site conditions or genotypes deployed.

#### **2. Materials and Methods**

#### *2.1. Site Description*

A regional phytotechnologies network consisting of sixteen phytoremediation buffer systems (i.e., phyto buffers) was established in 2017 (×6 phyto buffers), 2018 (×5), and 2019 (×5) in the Lake Superior watershed of the Upper Peninsula of Michigan, USA and the Lake Michigan watershed of eastern Wisconsin, USA. Multiple phyto buffers were installed at some locations, resulting in ten field testing sites throughout the network (Figure 1). The sites ranged in latitude from 46.7840 to 42.8382 ◦N and in longitude from −89.1291 to −86.5976 ◦W, which is consistent with poplar productivity supplysheds in the region [47,59–61]. Site-related climate properties are shown in Table 1. Twentyyear historical monthly averages for precipitation, temperature, and drought were determined across each growing season (April to October) for the time period 2000 to 2020 and summed (precipitation) and averaged (temperature, drought) to obtain final annual

values. Based on the nearest weather station to each site, precipitation (P, mm) along with maximum (Tmax, ◦C) and minimum (Tmin, ◦C) air temperatures were obtained from the National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (https://www.ncdc.noaa.gov/cdo-web/, accessed on 20 January 2021). Average temperature (i.e., Tavg = ((Tmax + Tmin)/2), ◦C) and the difference between maximum and minimum temperatures (i.e., Tdiff = Tmax − Tmin, ◦C) were calculated for each site. Daily growing degree days (GDDday) were calculated as the average temperature minus a base temperature of 10 ◦C (i.e., GDDday = Tavg − 10 ◦C) and summed across each growing season. Average annual growing degree days for each site (GDDannual) were estimated by averaging the seasonal GDD values from 2000 to 2020. The United States Drought Monitor (https://droughtmonitor.unl.edu/, accessed on 20 January 2021) was accessed to obtain drought index scores according to percent area within each county belonging to four drought index categories: D0 (abnormally dry), D1 (moderate drought), D2 (severe drought), D3 (extreme drought), and D4 (exceptional drought). Categories D3 and D4 were negligible and, therefore, not reported in Table 1. Buffer-specific soil properties were acquired from the USDA Natural Resources Conservation Service (NRCS) Web Soil Survey (https://websoilsurvey.sc.egov.usda.gov/, accessed on 20 January 2021) and are provided in Table 2. ues. Based on the nearest weather station to each site, precipitation (P, mm) along with maximum (Tmax, °C) and minimum (Tmin, °C) air temperatures were obtained from the National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (https://www.ncdc.noaa.gov/cdo-web/ accessed on 20 January 2021). Average temperature (i.e., Tavg = ((Tmax + Tmin)/2), °C) and the difference between maximum and minimum temperatures (i.e., Tdiff = Tmax−Tmin, °C) were calculated for each site. Daily growing degree days (GDDday) were calculated as the average temperature minus a base temperature of 10 °C (i.e., GDDday = Tavg−10 °C) and summed across each growing season. Average annual growing degree days for each site (GDDannual) were estimated by averaging the seasonal GDD values from 2000 to 2020. The United States Drought Monitor (https://droughtmonitor.unl.edu/ accessed on 20 January 2021) was accessed to obtain drought index scores according to percent area within each county belonging to four drought index categories: D0 (abnormally dry), D1 (moderate drought), D2 (severe drought), D3 (extreme drought), and D4 (exceptional drought). Categories D3 and D4 were negligible and, therefore, not reported in Table 1. Buffer-specific soil properties were acquired from the USDA Natural Resources Conservation Service (NRCS) Web Soil Survey (https://websoilsurvey.sc.egov.usda.gov/ accessed on 20 January 2021) and are provided in Table 2.

region [47,59–61]. Site-related climate properties are shown in Table 1. Twenty-year historical monthly averages for precipitation, temperature, and drought were determined across each growing season (April to October) for the time period 2000 to 2020 and summed (precipitation) and averaged (temperature, drought) to obtain final annual val-

*Forests* **2021**, *12*, x FOR PEER REVIEW 4 of 25

**Figure 1.** Regional phytotechnologies network consisting of sixteen phytoremediation buffer systems (i.e., phyto buffers) established in 2017 (×6 phyto buffers), 2018 (×5), and 2019 (×5) in the Lake Superior watershed of the Upper Peninsula of Michigan, USA and the Lake Michigan watershed of eastern Wisconsin, USA. **Figure 1.** Regional phytotechnologies network consisting of sixteen phytoremediation buffer systems (i.e., phyto buffers) established in 2017 (×6 phyto buffers), 2018 (×5), and 2019 (×5) in the Lake Superior watershed of the Upper Peninsula of Michigan, USA and the Lake Michigan watershed of eastern Wisconsin, USA.

**Table 1.** Precipitation, temperature, growing degree days, and drought indices of ten field testing sites in a regional phyto technologies network consisting of sixteen phytoremediation buffer systems (i.e., phyto buffers) established from 2017 to 2019 in the Lake Superior watershed of the Upper Peninsula of Michigan, USA and the Lake Michigan watershed of eastern Wisconsin, USA.


Climate and drought data are means ± one standard error across each growing season (April to October) from 2000 to 2020. Climate source: National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (https://www.ncdc.noaa.gov/cdo-web/, accessed on 20 January 2021). Drought source: United States Drought Monitor (https://droughtmonitor.unl.edu/, accessed on 20 January 2021); percent of area within the county in each category (D0 to D2).

#### *2.2. Clone Selection*

Phyto-recurrent selection was conducted through a polycyclic evaluation process to choose superior genotypes for out planting [32,58,62]. For each phyto buffer, soils were collected in the field, brought to the USDA Forest Service, Institute for Applied Ecosystem Studies in Rhinelander, Wisconsin, USA (Figure 1), and processed (i.e., dried and sifted to remove large debris) for use as soil treatments during three greenhouse selection cycles [58]. The number of clones tested in each cycle decreased from 140 (cycle 1) to 60 (cycle 2) to 15 (cycle 3) based on multiplicative rank summation indices incorporating tree survival, growth, physiology, and health [63]. Twelve clones were selected for cycle 4 field testing, resulting in three buffer groups related to phyto buffer establishment in 2017, 2018, and 2019. Phyto buffers within buffer group years had the same twelve clones, but the twelve clones differed among the three buffer groups based on cycle 3 phyto-recurrent selection

results. In the current study, separate analyses were conducted for each buffer group. Each set of buffer group clones consisted of genotypes that have: (1) been commonly used for commercial and/or research purposes in the region (i.e., Common clones), (2) a rich history of testing but are still at the experimental stage (i.e., Experimental), and (3) been bred, tested, and selected at the University of Minnesota Duluth, Natural Resources Research Institute (NRRI) and show promise for broad-ranging applications (i.e., NRRI) [46,56]. Specific clones and their genomic groups are listed in Table 3.

**Table 2.** Soil properties of sixteen phytoremediation buffer systems (i.e., phyto buffers) comprising a regional phytotechnologies network established from 2017 to 2019 in the Lake Superior watershed of the Upper Peninsula of Michigan, USA and the Lake Michigan watershed of eastern Wisconsin, USA.


Source: USDA Natural Resources Conservation Service (NRCS) Web Soil Survey (https://websoilsurvey.sc.egov.usda.gov/, accessed on 20 January 2021). <sup>a</sup> Phyto buffers: BC: Bellevue (Central); BE: Bellevue (East); BW: Bellevue (West); CE: Caledonia (East); CW: Caledonia (West); EE: Escanaba (East); EW: Escanaba (West); MA: Manitowoc; ME: Menomonee Falls (East); MW: Menomonee Falls (West); MQ: Marquette; MU: Munising; ON: Ontonagon (North); OS: Ontonagon (South); SL: Slinger; WH: Whitelaw. <sup>b</sup> Drainage classes: MWD: moderately well drained; PD: poorly drained; SED: somewhat excessively drained; SPD: somewhat poorly drained; WD: well drained. <sup>c</sup> Textures: L: loam; S: sand; SCL: sandy clay loam; SiCL: silty clay loam; SL: sandy loam. <sup>d</sup> na: not available.

**Table 3.** Genomic groups, clones, buffer groups (i.e., years of planting), and clone groups for *Populus* genotypes tested in a regional phyto technologies network of sixteen phytoremediation buffer systems (i.e., phyto buffers) established from 2017 to 2019 in the Lake Superior watershed of the Upper Peninsula of Michigan, USA and the Lake Michigan watershed of eastern Wisconsin, USA.


<sup>a</sup> Species authorities: *P. deltoides* Bartr. Ex Marsh; *P. nigra* L.; *P. maximowiczii* A. Henry NRRI = promising genotypes bred, tested, and selected at the University of Minnesota Duluth, Natural Resources Research Institute (NRRI) for broad-ranging applications [46,56]; <sup>b</sup> Experimental = genotypes with a rich history of testing but that are still at the experimental stage; Common = genotypes commonly used for commercial and/or research purposes in the region.

#### *2.3. Phyto Buffer Establishment and Experimental Design*

Dormant, unrooted hardwood cuttings (measuring 25.4 cm in length, and containing at least one primary bud within 2.54 cm from the top of each cutting) were processed from one-year-old whips collected during the dormant season (i.e., January through March) from: the USDA Forest Service Hugo Sauer Nursery (Rhinelander, WI, USA); University of Minnesota NRRI Clonal Orchard at the North Central Research and Outreach Center Nursery (Grand Rapids, Minnesota, USA); Iowa State University Clonal Orchard at the Iowa Department of Natural Resources State Nursery (Ames, IA, USA); and Michigan State University Clonal Orchard at the Tree Research Center (Lansing, MI, USA).

Processed cuttings were stored in polyethylene bags at 5◦ C and, prior to planting, soaked in water to a height of 16.93 cm for 48 h in a dark room at 21 ◦C during May and June in 2017, 2018, and 2019 (Table S1). Prior to planting, rocks and other large obstructions were removed and the soil was tilled to a depth of 30 cm. Thus, at planting, there was no competition from weeds and/or grasses. To reduce potential competition effects over time, subsequent site maintenance throughout the study period included: (1) tilling to a depth of 30 cm, (2) hand weeding to a minimum diameter of 0.61 m around each individual tree, and (3) removing rocks and other obstructions. For each phyto buffer, a minimum of one maintenance cycle per month throughout each growing season was conducted.

At each of the six (2017) or five (2018, 2019) phyto buffers, trees were planted in a randomized complete block design (RCBD) with eight blocks and twelve clones at a spacing of 2.44 m <sup>×</sup> 2.44 m (i.e., 1680 trees ha−<sup>1</sup> ). Due to field space constraints, four blocks were planted at Slinger, Wisconsin. Clones were arranged in randomized complete blocks to minimize effects of potential environmental gradients, especially those related to runoff and soil physico-chemical properties. Two border rows were established on the perimeter of each phyto buffer to reduce potential border effects [64,65]. All phyto buffers were fenced using 2.3 m tall Trident extra strength deer fencing (Trident Enterprises, Waynesboro, PA, USA) to eliminate potential impacts from white-tailed deer (*Odocoileus virginianus* Zimmerman) browse. At the end of the first growing season, tree survival across all phyto buffers and clones was 97.1 %, with 95.5%, 96.2%, and 99.6 % survival for the 2017, 2018, and 2019 buffer groups, respectively. All trees that died were replanted with the same genotype to ensure full stocking of 1680 trees ha−<sup>1</sup> . However, the replanted trees were not included in the analyses below.

#### *2.4. Field Measurements*

At the end of each growing season, tree height to the nearest 0.1 m was measured from the ground to the apical bud. Diameter was measured at 10 cm above the soil surface for one- and two-year-old trees and at breast height (i.e., DBH at 1.37 m) for three-year-old trees. All diameter estimates were determined to the nearest 0.1 cm. Tree volume (V) was calculated from height (H) and diameter (D; including one- and two-year diameter and DBH) using the model: V = D<sup>2</sup> <sup>×</sup> H [66]. Starting in 2020, trees of the 2017 buffer group were too tall to continue height measurements at the requisite precision noted above. As a result, 2020 DBH measurements collected from 2017 buffer group trees were used to estimate mean annual increment (MAI; Mg ha−<sup>1</sup> yr−<sup>1</sup> ). In this particular case, individualtree biomass was estimated using the general model: Biomass = 10a0 <sup>×</sup> DBHa1 while applying genomic-group specific coefficients from Headlee and Zalesny [67] (*P. deltoides* 'D': a<sup>0</sup> = -0.65, a<sup>1</sup> = 2.01; *P. deltoides* × *P. nigra* 'DN': a<sup>0</sup> = −1.02, a<sup>1</sup> = 2.36; *P. deltoides* × *P. maximowiczii* 'DM': a<sup>0</sup> = −1.03, a<sup>1</sup> = 2.33; *P. nigra* × *P. maximowiczii* 'NM': a<sup>0</sup> = −0.50, a<sup>1</sup> = 1.94). Individual-tree biomass estimates were then scaled to MAI using standard metric conversion factors and stocking of 1680 trees ha−<sup>1</sup> at four years after planting.

#### *2.5. Health Assessments*

Based on a modified methodology of Rogers et al. [58], individual tree health parameters were scored using a five-category qualitative scale ranging from 1 to 5, where 1 = optimal health, 2 = good health, 3 = moderate health, 4 = poor health, and 5 = dead (i.e.,

health score was inversely related to health). Two researchers scored each parameter to promote consistency in ratings. There were six parameters: (1) vigor, (2) defoliation, (3) leaf discoloration, (4) chlorosis, (5) leaf scorch, and (6) leaf spots. A multiplicative weighted summation index was used to calculate final health index values, with vigor receiving a coefficient of 0.25 and all other parameters having a health coefficient equal to 0.15. Health assessments were not conducted in 2020.

#### *2.6. Data Analysis*

Health (of all buffer groups) and MAI (of the 2017 buffer group) data were subjected to analyses of variance (ANOVA) and analyses of means (ANOM) using SAS® (PROC GLM; PROC ANOM; SAS INSTITUTE, INC., Cary, NC, USA) assuming a two-way factorial design including six (2017) or five (2018, 2019) buffers, twelve clones, and their interactions. Fisher's Least Significant Difference (LSD) was used to identify significant differences among least-squares means for main effects and interactions at *P* < 0.05.

Height and volume (of all buffer groups) and diameter (excluding 2020 diameter of 2017 buffer group trees) data were subjected to analyses of variance (ANOVA) and analyses of means (ANOM) using SAS® (PROC MIXED; PROC ANOM; SAS INSTITUTE, INC., Cary, NC, USA) assuming a three-way, repeated measures factorial design including six (2017) or five (2018, 2019) buffers, twelve clones, three (2017, 2018) or two (2019) ages, and their interactions. The ages (representing tree growth after each growing season) were analyzed as the repeated measure. To account for pseudo-replication over time, six different covariance structures (i.e., vc, cs, ar (1), toep, ante (1), un) were tested in PROC MIXED to determine which one provided the best model fit based on the lowest BIC scores. Using these covariance structures, ANOVA were conducted in PROC MIXED for all traits, and multiple comparisons analyses were conducted to identify significant differences among least-squares means for main effects and interactions as noted above.

#### **3. Results**

#### *3.1. Health*

There were significant differences in health of the 2017 buffer group trees (HEALTH2017) among buffer and clone main effects for the first, second, and third year of growth (*P2017,2018,2019* < 0.0001), yet the buffer × clone interaction governed health for all three years (*P<sup>2017</sup>* < 0.0001; *P<sup>2018</sup>* < 0.0001; *P<sup>2019</sup>* = 0.0152) (Table S2). HEALTH<sup>2017</sup> of trees measured in 2017 (i.e., HEALTH2017(2017)) ranged from 1.0 ± 0.0 ('NM2' and 'NM5' at Menomonee Falls (East); 'NM2' at Whitelaw; most healthy] to 1.6 ± 0.1 ('DM114' at Caledonia (East); least healthy), with an overall mean of 1.1 ± 0.1 (Figure 2). Thus, all trees were of optimal health (i.e., health index ranging from 1 to 2). The healthiest trees were grown at Whitelaw, which had 14.1% better HEALTH2017(2017) than at Caledonia (East), the buffer with trees exhibiting the poorest health. The range in health scores was broader for clones, with 'NM5' having 37.0% healthier trees than 'DM114'. Common clones had the healthiest and Experimental clones the least-healthiest trees, with NRRI genotypes being intermediate for HEALTH2017(2017). Seven buffer × clone interactions resulted in HEALTH2017(2017) values that were significantly greater than the overall mean ['DM114', 'DN177', and 'NC14106' at Caledonia (East); 'DM114' and 'NC14106' at Menomonee Falls (East); 'DM114' at Menomonee Falls (West); 'DM114' at Slinger]. With the exception of 'DM114' and 'DN177', all clones were generalists for HEALTH2017(2017), with health index scores not varying by more than 0.3 for any buffer × clone combinations. Trees grown at Bellevue (West) and Whitelaw were 38.4% healthier than those at Menomonee Falls (East) and Caledonia (East) for 'DM114'. Similarly, for 'DN177', trees at Whitelaw were 32.2% healthier than at Caledonia (East) (Figure 2). Trends in second- (HEALTH2017(2018)) and third-year (HEALTH2017(2019)) health of the 2017 buffer group trees were similar to HEALTH2017(2017) (Figures S1 and S2).

years (*P2017* < 0.0001; *P2018* < 0.0001; *P2019* = 0.0152) (Table S2). HEALTH2017 of trees measured in 2017 (i.e., HEALTH2017(2017)) ranged from 1.0 ± 0.0 ('NM2' and 'NM5' at Menomonee Falls (East); 'NM2' at Whitelaw; most healthy] to 1.6 ± 0.1 ('DM114' at Caledonia (East); least healthy), with an overall mean of 1.1 ± 0.1 (Figure 2). Thus, all trees were of optimal health (i.e., health index ranging from 1 to 2). The healthiest trees were grown at Whitelaw, which had 14.1% better HEALTH2017(2017) than at Caledonia (East), the buffer with trees exhibiting the poorest health. The range in health scores was broader for clones, with 'NM5' having 37.0% healthier trees than 'DM114'. Common clones had the healthiest and Experimental clones the least-healthiest trees, with NRRI genotypes being intermediate for HEALTH2017(2017). Seven buffer × clone interactions resulted in HEALTH2017(2017) values that were significantly greater than the overall mean ['DM114', 'DN177', and 'NC14106' at Caledonia (East); 'DM114' and 'NC14106' at Menomonee Falls (East); 'DM114' at Menomonee Falls (West); 'DM114' at Slinger]. With the exception of 'DM114' and 'DN177', all clones were generalists for HEALTH2017(2017), with health index scores not varying by more than 0.3 for any buffer × clone combinations. Trees grown at Bellevue (West) and Whitelaw were 38.4% healthier than those at Menomonee Falls (East) and Caledonia (East) for 'DM114'. Similarly, for 'DN177', trees at Whitelaw were 32.2% healthier than at Caledonia (East) (Figure 2). Trends in second- (HEALTH2017(2018)) and third-year (HEALTH2017(2019)) health of the 2017 buffer group trees were similar to HEALTH2017(2017) (Figures S1 and S2).

**Figure 2.** Tree health (± one standard error) determined after the 2017 growing season of twelve poplar clones tested in six phytoremediation buffer systems (i.e., phyto buffers) (**A**–**F**) established in 2017 (i.e., the 2017 buffer group) in the Lake Michigan watershed of eastern Wisconsin, USA. The dashed line represents the overall mean, and asterisks indicate means different than the overall mean at *P* < 0.05. Bars with different letters across all buffer × clone combinations are different at **Figure 2.** Tree health (± one standard error) determined after the 2017 growing season of twelve poplar clones tested in six phytoremediation buffer systems (i.e., phyto buffers) (**A**–**F**) established in 2017 (i.e., the 2017 buffer group) in the Lake Michigan watershed of eastern Wisconsin, USA. The dashed line represents the overall mean, and asterisks indicate means different than the overall mean at *P* < 0.05. Bars with different letters across all buffer × clone combinations are different at *P* < 0.05. See Materials and Methods for complete tree health definitions (1 = optimal health, 2 = good health, 3 = moderate health, 4 = poor health, and 5 = dead).

Differences among buffer and clone main effects were significant for first- (HEALTH2018(2018)) and second-year (HEALTH2018(2019)] health of the 2018 buffer group trees (*P<sup>2018</sup>* = 0.0009 for Clone; *P2018,2019* < 0.0001 for all other main effects), yet the buffer × clone interaction governed health for both years (*P<sup>2018</sup>* = 0.0029; *P<sup>2019</sup>* < 0.0001) (Table S2). HEALTH2018(2018) ranged from 1.0 ± 0.0 ('9732-11', '9732-24', 'DM114', 'DN2', 'NM2', 'NM5', 'NM6' at Bellevue (East); '9732-11', 'DN2' at Bellevue (Central); most healthy) to 1.5 ± 0.1 ('DM114' at Manitowoc; least healthy), with an overall mean of 1.1 ± 0.1 (Figure 3). Similar to HEALTH2017(2017) above, all trees were of optimal health. The healthiest trees were grown at Marquette and Bellevue (East), which had 25.5 % better HEALTH2018(2018) than at Manitowoc that exhibited the poorest health. The range in health scores was narrower for clones, with 'NM5' having 12.9 % healthier trees than 'DM114' that had the poorest health. Common and NRRI clones had the best overall health scores that were similar to one another yet better than Experimental genotypes. Seven buffer × clone interactions resulted in HEALTH2018(2018) values there were significantly greater (i.e., of poorer health) than the overall mean, with all occurring at Manitowoc: '9732-11', '9732-24', '9732-31', '9732-36', '7300502', 'DM114', 'DN2', and 'DN5'. With the exception of '7300502' and 'DM114', all clones were classified as generalists for HEALTH2018(2018). Trees grown at

health, 4 = poor health, and 5 = dead).

Bellevue (East) and Marquette were 52.5% healthier than those at Manitowoc for '7300502', while those at Bellevue (East), Bellevue (Central), and Marquette were 48.8% healthier than at Manitowoc (Figure 3). Trends in second-year health of the 2018 buffer group trees (HEALTH2018(2019)) were similar to HEALTH2018(2018) (Figure S3). clones were classified as generalists for HEALTH2018(2018). Trees grown at Bellevue (East) and Marquette were 52.5% healthier than those at Manitowoc for '7300502', while those at Bellevue (East), Bellevue (Central), and Marquette were 48.8% healthier than at Manitowoc (Figure 3). Trends in second-year health of the 2018 buffer group trees (HEALTH2018(2019)) were similar to HEALTH2018(2018) (Figure S3).

Differences among buffer and clone main effects were significant for first- (HEALTH2018(2018)) and second-year (HEALTH2018(2019)] health of the 2018 buffer group trees (*P2018* = 0.0009 for Clone; *P2018,2019* < 0.0001 for all other main effects), yet the buffer × clone interaction governed health for both years (*P2018* = 0.0029; *P2019* < 0.0001) (Table S2). HEALTH2018(2018) ranged from 1.0 ± 0.0 ('9732-11', '9732-24', 'DM114', 'DN2', 'NM2', 'NM5', 'NM6' at Bellevue (East); '9732-11', 'DN2' at Bellevue (Central); most healthy) to 1.5 ± 0.1 ('DM114' at Manitowoc; least healthy), with an overall mean of 1.1 ± 0.1 (Figure 3). Similar to HEALTH2017(2017) above, all trees were of optimal health. The healthiest trees were grown at Marquette and Bellevue (East), which had 25.5 % better HEALTH2018(2018) than at Manitowoc that exhibited the poorest health. The range in health scores was narrower for clones, with 'NM5' having 12.9 % healthier trees than 'DM114' that had the poorest health. Common and NRRI clones had the best overall health scores that were similar to one another yet better than Experimental genotypes. Seven buffer × clone interactions resulted in HEALTH2018(2018) values there were significantly greater (i.e., of poorer health) than the overall mean, with all occurring at Manitowoc: '9732-11', '9732-24', '9732-31', '9732-36', '7300502', 'DM114', 'DN2', and 'DN5'. With the exception of '7300502' and 'DM114', all

*Forests* **2021**, *12*, x FOR PEER REVIEW 10 of 25

*P* < 0.05. See Materials and Methods for complete tree health definitions (1 = optimal health, 2 = good health, 3 = moderate

**Figure 3.** Tree health (±one standard error) determined after the 2018 growing season of twelve poplar clones tested in five phytoremediation buffer systems (i.e., phyto buffers) (**A**–**E**) established in 2018 (i.e., the 2018 buffer group) in the Lake Superior watershed of the Upper Peninsula of Michigan, USA and the Lake Michigan watershed of eastern Wisconsin, USA. The dashed line represents the overall mean, and asterisks indicate means different than the overall mean at *P* < 0.05. Bars with different letters across all buffer × clone combinations are different at *P* < 0.05. See Materials and Methods for complete tree health definitions (1 = optimal health, 2 = good health, 3 = moderate health, 4 = poor health, and 5 = dead).

Differences among buffer and clone main effects were significant for first-year health of the 2019 buffer group trees (HEALTH2019(2019)) (*PBuffer* < 0.0001; *PClone* = 0.0159), yet the buffer × clone interaction governed health (*P* = 0.0036) (Table S2). HEALTH2019(2019) ranged from 1.0 ± 0.0 ['9732-24', 'DM114', 'DN2', 'NM2', 'NM6' at Munising; '9732-31' at Ontonagon (North); '99038022', '9732-11', 'DM114', 'DN2', 'NM6' at Ontonagon (South); most healthy] to 1.3 ± 0.0 ['DM114', 'DN2' Escanaba (West); least healthy], with an overall mean of 1.1 ± 0.0 (Figure 4). All trees were of optimal health. The healthiest trees were grown at Ontonagon (South), which had 13.7% better HEALTH2019(2019) than at Escanaba (West) where trees had the poorest health. The range in health scores was narrower for clones, with 7.2% separating the healthiest trees of 'NM2' from the least healthy trees of 'DN2', 'DN34', and 'DN177'. Similar to HEALTH2017(2017), Common clones had the healthiest and Experimental clones the least-healthiest trees, with NRRI genotypes being intermediate

for HEALTH2019(2019). Two buffer × clone interactions resulted in HEALTH2019(2019) values there were significantly greater than the overall mean ('DM114'and 'DN2' at Escanaba (West)). No clones varied by more than 0.3 health index points for any buffer × clone combinations, thus indicating that all had generalist health performance. The largest variation in clonal responses to specific buffers was where trees were 21.7% healthier at all buffers relative to Escanaba (West) for 'DM114', and where trees were 26.3% healthier at Munising and Ontonagon (South) than Escanaba (West) for 'DN2' (Figure 4). *Forests* **2021**, *12*, x FOR PEER REVIEW 12 of 25

**Figure 4.** Tree health (±one standard error) determined after the 2019 growing season of twelve poplar clones tested in five phytoremediation buffer systems (i.e., phyto buffers) (**A**–**E**) established in 2019 (i.e., the 2019 buffer group) in the Lake Superior watershed of the Upper Peninsula of Michigan, USA. The dashed line represents the overall mean, and asterisks indicate means different than the overall mean at *P* < 0.05. Bars with different letters across all buffer × clone combinations are different at *P* < 0.05. See Materials and Methods for complete tree health definitions (1 = optimal health, 2 = good health, 3 = moderate health, 4 = poor health, and 5 = dead). **Figure 4.** Tree health (±one standard error) determined after the 2019 growing season of twelve poplar clones tested in five phytoremediation buffer systems (i.e., phyto buffers) (**A**–**E**) established in 2019 (i.e., the 2019 buffer group) in the Lake Superior watershed of the Upper Peninsula of Michigan, USA. The dashed line represents the overall mean, and asterisks indicate means different than the overall mean at *P* < 0.05. Bars with different letters across all buffer × clone combinations are different at *P* < 0.05. See Materials and Methods for complete tree health definitions (1 = optimal health, 2 = good health, 3 = moderate health, 4 = poor health, and 5 = dead).

#### *3.2. Biomass and Growth 3.2. Biomass and Growth*

(East), and Whitelaw.

Differences for buffer and clone main effects were significant for fourth-year mean annual increment of the 2017 buffer group trees (MAI2017(2020)) (*P* < 0.0001), yet the buffer × clone interaction governed MAI2017(2020) (*P* < 0.0001) (Table S2). MAI2017(2020) ranged from 1.10 ± 0.73 ('7300502' at Whitelaw) to 7.67 ± 0.5 Mg ha−1 yr−1 ('NM5' at Menomonee Falls (East)), with an overall mean of 3.20 ± 0.51 Mg ha−1 yr−1 (Figure 5). The largest trees were grown at Menomonee Falls (West), which had 55.4% greater MAI2017(2020) than at Whitelaw, the buffer with the smallest trees. The range in MAI2017(2020) was broader for clones, with 'NM5' exhibiting 69.8% more biomass than 'NC14106' that had the smallest trees of any geno-Differences for buffer and clone main effects were significant for fourth-year mean annual increment of the 2017 buffer group trees (MAI2017(2020)) (*P* < 0.0001), yet the buffer × clone interaction governed MAI2017(2020) (*P* < 0.0001) (Table S2). MAI2017(2020) ranged from 1.10 <sup>±</sup> 0.73 ('7300502' at Whitelaw) to 7.67 <sup>±</sup> 0.5 Mg ha−<sup>1</sup> yr−<sup>1</sup> ('NM5' at Menomonee Falls (East)), with an overall mean of 3.20 <sup>±</sup> 0.51 Mg ha−<sup>1</sup> yr−<sup>1</sup> (Figure 5). The largest trees were grown at Menomonee Falls (West), which had 55.4% greater MAI2017(2020) than at Whitelaw, the buffer with the smallest trees. The range in MAI2017(2020) was broader for clones, with

type. While MAI2017(2020) varied across genotypes, trends across clone groups were non-

ues that were significantly greater than the overall mean, with the most notable being 'NM5' outperforming the mean at four of the six buffers: Caledonia (East), Menomonee Falls (East), Menomonee Falls (West), and Slinger. In contrast, 'NC14106' had significantly less MAI2017(2020) than the overall mean at three buffers: Bellevue (West), Menomonee Falls

'NM5' exhibiting 69.8% more biomass than 'NC14106' that had the smallest trees of any genotype. While MAI2017(2020) varied across genotypes, trends across clone groups were non-existent, with Common, Experimental, and NRRI clones performing similarly across buffer × clone combinations. Many buffer × clone interactions resulted in MAI2017(2020) values that were significantly greater than the overall mean, with the most notable being 'NM5' outperforming the mean at four of the six buffers: Caledonia (East), Menomonee Falls (East), Menomonee Falls (West), and Slinger. In contrast, 'NC14106' had significantly less MAI2017(2020) than the overall mean at three buffers: Bellevue (West), Menomonee Falls (East), and Whitelaw. *Forests* **2021**, *12*, x FOR PEER REVIEW 13 of 25

**Figure 5.** Mean annual increment (Mg ha<sup>−</sup>1 yr−1) (±one standard error) determined after the 2020 growing season of twelve poplar clones tested in six phytoremediation buffer systems (i.e., phyto buffers) (**A**–**F**) established in 2017 (i.e., the 2017 buffer group) in the Lake Michigan watershed of eastern Wisconsin, USA. The dashed line represents the overall mean, and asterisks indicate means different than the overall mean at *P* < 0.05. Bars with different letters across all buffer × clone combinations are different at *P* < 0.05. **Figure 5.** Mean annual increment (Mg ha−<sup>1</sup> yr−<sup>1</sup> ) (±one standard error) determined after the 2020 growing season of twelve poplar clones tested in six phytoremediation buffer systems (i.e., phyto buffers) (**A**–**F**) established in 2017 (i.e., the 2017 buffer group) in the Lake Michigan watershed of eastern Wisconsin, USA. The dashed line represents the overall mean, and asterisks indicate means different than the overall mean at *P* < 0.05. Bars with different letters across all buffer × clone combinations are different at *P* < 0.05.

In addition to these changes in magnitude, there were distinct changes in MAI2017(2020) ranks that defined the genotypes as generalists or specialists. In particular, clones exhibited generalist MAI2017(2020), with four exceptions (Table 4). First, '99059016' had stable performance across five of the six buffers, with Whitelaw (i.e., the buffer with its lowest rank of eleventh) having 60.9% less MAI2017(2020) than Bellevue (West), where '99059016' ranked fifth. Second, '7300502' had broad variation in MAI2017(2020) across all six buffers, with a 75.3% reduction in MAI2017(2020) at Whitelaw (rank = 12) versus Bellevue (West) (rank = 3). Third, 'DN177' was a specialist because of its high MAI2017(2020) at Slinger (where it ranked In addition to these changes in magnitude, there were distinct changes in MAI2017(2020) ranks that defined the genotypes as generalists or specialists. In particular, clones exhibited generalist MAI2017(2020), with four exceptions (Table 4). First, '99059016' had stable performance across five of the six buffers, with Whitelaw (i.e., the buffer with its lowest rank of eleventh) having 60.9% less MAI2017(2020) than Bellevue (West), where '99059016' ranked fifth. Second, '7300502' had broad variation in MAI2017(2020) across all six buffers, with a 75.3% reduction in MAI2017(2020) at Whitelaw (rank = 12) versus Bellevue (West) (rank = 3). Third, 'DN177' was a specialist because of its high MAI2017(2020) at Slinger (where it ranked

fourth), which was 36.7% greater than Menomonee Falls (East) (rank = 9). Fourth, despite higher ranking (rank = 7) at Whitelaw relative to other buffers for 'DN34', MAI2017(2020) at

response (Table 4). However, this classification for 'DN34' should be interpreted with caution, especially given that 'DN34' ranked tenth at four buffers and eleventh at one buffer,

indicating stable performance across buffers.

fourth), which was 36.7% greater than Menomonee Falls (East) (rank = 9). Fourth, despite higher ranking (rank = 7) at Whitelaw relative to other buffers for 'DN34', MAI2017(2020) at Whitelaw was 60.6% less than that of Caledonia (East) (rank = 11), resulting in its specialist response (Table 4). However, this classification for 'DN34' should be interpreted with caution, especially given that 'DN34' ranked tenth at four buffers and eleventh at one buffer, indicating stable performance across buffers.

**Table 4.** Clonal rank for mean annual increment (MAI2017(2020)) measured after the fourth growing season of twelve poplar clones tested in six phytoremediation buffer systems (i.e., phyto buffers) established in 2017 (i.e., the 2017 buffer group) in the Lake Michigan watershed of eastern Wisconsin, USA.


<sup>a</sup> BW: Bellevue (West); CE: Caledonia (East); ME: Menomonee Falls (East); MW: Menomonee Falls (West); SL: Slinger; WH: Whitelaw. <sup>b</sup> Generalist = clone exhibiting stable MAI2017(2020) across phyto buffers (i.e., minimal rank changes); Specialist = clone exhibiting exceptional MAI2017(2020) at one or more phyto buffers relative to the other buffers (i.e., broad rank changes).

The buffer × clone × year interaction was significant for height (*P* = 0.0483), diameter (*P* = 0.0018), and volume (*P* = 0.0001) of the 2017 buffer group trees (Table S3). VOLUME<sup>2017</sup> of trees measured in 2017 (VOLUME2017(2017)) ranged from 3.1 ± 54.4 ('7300502' at Whitelaw) to 459.6 <sup>±</sup> 22.6 cm<sup>3</sup> ('7300502' at Slinger), with an overall mean of 132.0 <sup>±</sup> 38.7 cm<sup>3</sup> , while VOLUME2017(2018) (i.e., 2017 buffer group trees measured in 2018) ranged from 503.5 <sup>±</sup> 1619.6 ('7300502' at Whitelaw) to 13,027.0 <sup>±</sup> 1402.6 cm<sup>3</sup> ('NM2' at Slinger), with an overall mean of 5765.2 <sup>±</sup> 1132.3 cm<sup>3</sup> (Table 5). VOLUME2017(2019) ranged from 1668.7 <sup>±</sup> 3018.1 ('7300502' at Whitelaw) to 26,652.0 <sup>±</sup> 1848.0 cm<sup>3</sup> ('NM5' at Menomonee Falls (East)), with an overall mean of 10,530.6 <sup>±</sup> 2109.5 cm<sup>3</sup> (Table 5). Across all buffer × clone × year combinations, volume increased 43.7-fold from the first year to the second year after planting, and then 1.8-fold from the second year to the third year. After the first and second growing seasons, the largest trees were grown at Slinger, which had 1301 and 268% greater volume than the buffer with the smallest trees (Whitelaw), respectively. For the third growing season, trees at Menomonee Falls (East) were largest, with VOLUME2017(2019) being 338% greater than Whitelaw, which had the smallest trees.

The range in volume was narrower for clones, with 'DN5' having 204% greater volume than the least productive clone ('99059016') for VOLUME2017(2017), 'NM5' being 86% greater than '7300502' for VOLUME2017(2018), and 'NM5' being 114% greater than 'NC14106' for VOLUME2017(2019). With the exception of the first growing season where NRRI clones exhibited the least overall volume followed by Experimental and Common (most volume) genotypes, Common clones had the largest and Experimental clones the smallest trees, with NRRI genotypes being intermediate for VOLUME2017(2018) and VOLUME2017(2019). Trends in height and diameter of the 2017 buffer group trees were similar to volume (Tables S4 and S5). Furthermore, in addition to the broad variability in the magnitude of differences among buffer × clone × year combinations, the frequency and magnitude of changes in rank within and across years defined genotypes as generalists (i.e., minimal rank changes)

or specialists (i.e., broad variation resulting in ≥5 rank changes for at least one buffer × clone × year pair) (Table S6). In particular, the NRRI clones ('99038022', '99059016', '9732-36') were high-level generalists characterized by nearly universal stability in ranks across buffers within measurement years, less than three substantial (i.e., >5 ranks) rank changes across all three-way combinations, and moderate to high rank stability over time. Clones 'DM114' and 'NC14106' were also generalists, exhibiting moderate rank stability within years and consistent ranks over time. The remaining genotypes were specialists. Clone '7300502' had the most variability in early ranks of all clones, followed by high rank stability in later years that were not consistent over time. For example, '7300502' ranked first at Slinger in 2017 but then twelfth at this buffer in 2018 and 2019. Clones 'NM2', 'NM5', and 'NM6' were high-level specialists characterized by broad variability in ranks across buffers within measurement years, frequent substantial rank changes across all three-way combinations, and moderate stability over time. Similarly, 'DN5', 'DN34', and 'DN177' were specialists, albeit with more moderate rank variation than the 'NMx' genotypes (Table S6).

**Table 5.** Volume (cm<sup>3</sup> ) (±one standard error) of twelve poplar clones tested in six phytoremediation buffer systems (i.e., phyto buffers) established in 2017 (i.e., the 2017 buffer group) in the Lake Michigan watershed of eastern Wisconsin, USA. Trees were measured following the 2017, 2018, and 2019 growing seasons. Volume values with different letters within a clone column across measurement years are different at *P* < 0.05.


<sup>a</sup> BW: Bellevue (West); CE: Caledonia (East); ME: Menomonee Falls (East); MW: Menomonee Falls (West); SL: Slinger; WH: Whitelaw.

The buffer × clone × year interaction was significant for height, diameter, and volume (*P* < 0.0001) of the 2018 buffer group trees (Table S3). VOLUME2018(2018) ranged from 12.2 <sup>±</sup> 38.4 ('DN5' at Marquette) to 185.7 <sup>±</sup> 25.3 cm<sup>3</sup> ('NM2' at Manitowoc), with an overall mean of 71.2 <sup>±</sup> 26.3 cm<sup>3</sup> , while VOLUME2018(2019) ranged from 287.8 ± 1518.7 ('DN5' at Marquette) to 11,085.0 <sup>±</sup> 930.0 cm<sup>3</sup> ('NM2' at Manitowoc), with an overall mean of 3418.4 <sup>±</sup> 1035.2 cm<sup>3</sup> (Table 6). VOLUME2018(2020) ranged from 261.0 ± 3882.2 ('DN5' at Marquette) to 27,220.0 <sup>±</sup> 2377.4 cm<sup>3</sup> ('NM5' at Manitowoc), with an overall mean of 8826.0 <sup>±</sup> 2646.2 cm<sup>3</sup> (Table 6). Across all buffer × clone × year combinations, volume increased 48-fold from the first year to the second year after planting, and then 2.6-fold from the second year to the third year. After the first growing season, the largest trees were grown at Caledonia (West), which had 529% greater volume than Marquette, the buffer with the smallest trees. During the second and third growing seasons, the largest trees were from Manitowoc, which had 1079 and 1744% greater volume than the buffer with the smallest trees (Marquette), respectively. As with volume for the 2017 buffer group trees, there was less variability among clones than buffers, with '9732-31' having 156% greater volume than the least productive clone ('7300502') for VOLUME2018(2018), 'NM5' being 163% greater than '7300502' for VOLUME2018(2018), and 'NM5' being 143% greater than 'DM114' for VOLUME2018(2018). During the first growing season, NRRI clones had the greatest overall volume, while Common clones exhibited the largest trees at two and three years after planting. For all three years, Experimental clones had the least volume. Trends in height and diameter of the 2018 buffer group trees were similar to volume (Tables S7 and S8). Moreover, as with 2017 buffer group trees, changes in magnitude and rank of 2018 buffer group clones across buffer × year combinations defined them as generalists or specialists (defined as for Table S6 above) (Table S9). Similar to 2017, the NRRI clones ('9732-11', '9732-24', '9732-31', '9732-36') were high-level generalists with universal stability in ranks, few rank changes greater than five ranks, and high rank stability over time. Clone 'DM114' was also a generalist, having only two substantial rank changes and nearly identical ranks from 2017 to 2019. All other genotypes were specialists. As in 2017, clone '7300502' had a high level of early rank variability and low to moderate rank consistency over time; clones 'NM2', 'NM5', and 'NM6' were high-level specialists with broad rank variability, numerous substantial rank changes, and moderate stability as trees aged; 'DN2', 'DN5', and 'DN34' were consistent specialists with moderate levels of rank changes and age-dependent stability (Table S9).

The buffer × clone × year interaction was significant for height (*P* = 0.0079), diameter (*P* < 0.0001), and volume (*P* < 0.0001) of the 2019 buffer group trees (Table S3). VOLUME2019(2019) ranged from 8.6 <sup>±</sup> 26.9 ('DN177' at Ontonagon (North)) to 396. <sup>±</sup> 26.7 cm<sup>3</sup> ('99038022' at Escanaba (West)), with an overall mean of 88.7 <sup>±</sup> 26.7 cm<sup>3</sup> , while VOLUME2019(2020) ranged from 92.6 ± 612.7 ('NM5' at Ontonagon (North)) to 8909.6 <sup>±</sup> 573.1 cm<sup>3</sup> ('NM2' at Escanaba (West)), with an overall mean of 1440.4 <sup>±</sup> 575.1 cm<sup>3</sup> (Table 7). Across all buffer × clone × year combinations, volume increased 16.2-fold from the first year to the second year after planting. After the first and second growing seasons, the largest trees were grown at Escanaba (West), which had 1008 and 1066% greater volume than the buffer with the smallest trees (Ontonagon (North)), respectively. For VOLUME2019(2019), '99038022' had 101% greater volume than the least productive clone ('DN177'), while 'NM2' was 164% greater than '9732-11' for VOLUME2019(2020). NRRI clones exhibited the greatest first-year volume, followed by Common and Experimental genotypes. For VOLUME2019(2020), ranks of NRRI, Experimental, and Common clones changed, with Common genotypes having the most VOLUME2019(2020) followed by Experimental and NRRI (least volume) clones. Trends in height and diameter of the 2019 buffer group trees were similar to volume (Tables S10 and S11). Furthermore, 2019 clones were classified as generalists or specialists (defined as for Table S6 above) (Table S12).



<sup>a</sup> BC: Bellevue (Central); BE: Bellevue (East); CW: Caledonia (West); MA: Manitowoc; MQ: Marquette.

As in previous buffer groups, 'DM114' was a high-level generalist. With the exception of '99038022' that was a high-level generalist with nearly universal rank stability, only two substantial rank changes across all three-way combinations, and moderate stability over time, the NRRI clones were specialists for the 2019 buffer group. In particular, '9732-11', '9732-24', '9732-31', and '9732-36' were moderate- to high-level specialists characterized by broad rank variability, substantial rank changes, and moderate stability over time. Similarly, the performance of 'DN2' as a generalist in the 2019 buffer group was different relative to its specialist volume in the 2018 buffer group. Across 2019 phyto buffers, 'DN2' exhibited moderate rank stability within years and consistent ranks over time. All remaining clones were moderate- ('DN34', 'DN177') and high-level ('NM2', 'NM5', 'NM6') specialists with trends for changes in magnitude and rank similar to their performance in 2017 and 2018 buffer groups (Table S12).



<sup>a</sup> EE: Escanaba (East); EW: Escanaba (West); MU: Munising; ON: Ontonagon (North); OS: Ontonagon (South).

#### **4. Discussion and Conclusions**
