Stability analysis of stem solidness, grain yield, and grain protein concentration in spring wheat

Abstract: The wheat stem sawfly, Cephus cinctus Norton (Hymenoptera: Cephidae), is a major pest of wheat (Triticum aestivum L.) in the northern Great Plains, where it is a constant threat in Montana and is resurging in Alberta, Saskatchewan, and North Dakota. Adoption of solid-stemmed cultivars is an important management tool for wheat growers; however, the inconsistent pith expression first noted with the release of ‘Rescue’ has been repeatedly observed in modern cultivars such as ‘Lillian’ in Canada. Given the extensive hectares planted to solid-stemmed wheat cultivars during an outbreak, the identification of cultivars that display stable stem solidness, grain yield, and grain protein concentration across a wide range of environments where stem sawfly infestations occur is desirable. We assessed spring wheat plant responses in eight solid-stemmed and two hollow-stemmed genotypes grown across diverse environments using multiple statistical models. Study sites included southern Alberta and Saskatchewan, Montana, and North Dakota. Most models agreed that the genotypes ‘Choteau’, ‘BW925’, and ‘Mott’ consistently displayed high and stable stem solidness concomitant with high grain yield. ‘Choteau’ and ‘BW925’ also consistently met or exceeded the desired threshold of a 3.75/5 pith rating (averaged from the lower four stem internodes) for optimum resistance, whereas ‘Mott’ developed optimal pith at a specific (early) phenological stage when resistance to wheat stem sawfly infestation is critical. Exploring the stability of stem solidness identified ideal genotypes that would enhance germplasm development efforts, which exemplifies how this approach can facilitate the selection, production, and adoption of solid-stemmed wheat cultivars in regions prone to wheat stem sawfly attack.


Introduction
The wheat stem sawfly (WSS) Cephus cinctus Norton (Hymenoptera: Cephidae), has been a serious pest of wheat (Triticum aestivum L.) in the northern Great Plains of North America since the late 19th century (Comstock 1889). In Canada, severe infestations of WSS began to occur in the southern Prairies of Canada shortly after the cultivation of wheat commenced in this region (Criddle 1923). In the United States, the areas historically prone to WSS are eastern and northern Montana, western North Dakota, northern South Dakota, and western Minnesota . Infestations of WSS have dropped in Canada after a resurgence period from 1999 to 2010 . The factors responsible for this decline may include climate, natural enemies of WSS, the adoption of solid-stemmed cultivars, and diverse cropping systems, such as crop rotation with non-host species. However, these factors may not be long-lasting, as the WSS appears to be resurging in Canada (Meers 2020) and remains a constant economically devastating insect pest in Montana. The WSS is also resurging in North and South Dakota, and has expanded outside of the traditional distribution area, where it has now caused damage in Colorado, Kansas, Wyoming, and Nebraska (Lesieur et al. 2016).
A comprehensive review of WSS biology and management practices have been reported by Beres et al. (2011). Briefly, in our study areas, WSS emerges from the infested stubble of the previous host crop such as wheat or barley (Hordeum vulgare L.), usually from late-May to mid-July. After mating, the female deposits an egg in a suitable host plant by puncturing the stem using a specialized saw-like ovipositor. Within 5-7 d, the egg will hatch and the larva begins feeding upon the tissues within the culm of the stem (Ainslie 1929). Larval feeding injury destroys the vascular bundles, interferes with water and nutrients reaching the developing wheat head, and consequently reduces photosynthesis, especially under lower light intensity (Macedo et al. 2007). Reduced plant moisture and light penetrating through the stem at senescence serve as environmental cues, whereby the larva moves towards the base of the plant, girdles the stem by chewing a v-shaped groove ("cutting") around it, fills that region with frass, and then encases itself in a hibernaculum to overwinter. In the field and under laboratory conditions, WSS infested plants have 5%-30% lower kernel weights compared with uninfected plants (Wallace et al. 1973;Delaney et al. 2010). The "girdled" or "cut" stems readily topple over when exposed to windy or rainy conditions. The toppled stems are often difficult to recover with a harvester, thus increasing attainable yield losses by another ∼15%, for a total loss of around 30%-35% Bekkerman and Weaver 2018) when WSS pressure is high. Based on modern commodity prices, stem cutting alone that is close to 50% could result in economic losses in excess of $400 million annually across the area of distribution now prone to attack (Beres et al. , 2017. An important strategy to prevent yield losses from WSS is to grow solid-stemmed wheat cultivars that develop an abundance of pith within the culm of the stem, creating a mechanical restriction for the larva. These cultivars can affect the larval survivorship, health, and fitness of C. cinctus (Beres et al. 2017). Pith increases larval mortality either by physically crushing the egg (Holmes and Peterson 1960), or by acting as a mechanical barrier that restricts larval movement inside the internodes, thereby reducing the boring activity of larvae (O'Keeffe et al. 1960;Cárcamo et al. 2005). Pith may also cause antibiosis by impairing the ability of WSS larvae to thrive inside the stem (Cook et al. 2018). The common solid-stemmed spring and winter wheat cultivars grown in North America derive their solidness from the Portuguese landrace line, 'S-615' (Kemp 1934). The underlying genetics conferring solidness across these cultivars are rather complex because of multiple recessive genes (Cook et al. 2004;Lanning et al. 2006). Therefore, cultivars often express pith (stem solidness) inconsistently (Hayat et al. 1995); for example, greater susceptibility to stem-cutting was noted in the first solid-stemmed cultivar 'Rescue' (S-615 derived) in 1949 at Regina, SK (Platt and Farstad 1949). More recently, Nilsen et al. (2017) suggested that the level of solidness in 'Lillian' was inadequate under economic infestations of WSS in Canada. Beres et al. (2017) developed an artificial neural network using precipitation-related variables to predict when 'Lillian' was most susceptible to damage by WSS, which could result in stem cutting levels >50%, depending on precipitation levels during the critical period of 1 June to 5 July.
The efficacy of WSS resistance is greatly influenced by both genotype (G) and growing environment, and (or) their interactions (GE) . Genetic mapping confirmed that stem-solidness in wheat cultivars is controlled predominantly by the Qss.msub-3BL locus on chromosome 3B introduced from the Portuguese landrace (McNeal et al. 1959;Varella et al. 2019) and other minor genes located on chromosomes 3A, 3B, 3D, 5A, 5B, and 5D (Larson and MacDonald 1959;Cook et al. 2004Cook et al. , 2019Houshmand et al. 2007;Nilsen et al. 2017;Appels et al. 2018;Varella et al. 2019). The recently released hexaploid wheat reference genome reported a gene called SSt1 (Appels et al. 2018), which is responsible for stem solidness controlled by the Qss.msub-3BL.b allele. The strong expression of stem solidness that has been observed in 'Choteau' over other cultivars was influenced by the presence of both Qss.msub-3BL and Qss.msub-3DL (Lanning et al. 2006). In addition to genetics, environmental factors, such as photoperiod, light intensity, and precipitation all greatly influenced pith development (Beres et al. 2017). For example, intense sunlight optimizes pith expression, whereas shading or cloudy conditions limit it (Eckroth and McNeal 1953;Holmes 1984).
The inconsistent nature of pith development in solidstemmed wheat cultivars is a concern for growers in WSS prone areas. Producers often express a reluctance to adopt solid-stemmed cultivars over concern of a purported yield penalty (Beres et al. 2017). In some cases, however, solid-stemmed common and durum wheat can be agronomically superior to hollow-stemmed cultivars, particularly in environments with high WSS pressure (Beres et al. , 2009. For example, 'Lillian' , the solid-stemmed wheat cultivar released during the last outbreak in Canada (1999Canada ( -2007, occupied one-third of the wheat hectares in Saskatchewan and 20% of the prairie-wide wheat hectares (Beres et al. 2017). Thus, the problem of inconsistent pith development can result in a widespread failure to mitigate WSS damage if additional agronomic practices are not adopted (Nilsen et al. 2017). Advances in cultivar development to manage WSS would benefit from an improved understanding of the pith stability across regions prone to infestations in Canada and the USA. One strategy would be to characterize cultivars for pith expression using genotype by environment interaction (GEI) analyses, commonly used for yield stability assessments and used on a range of crops, including triticale and wheat (Navabi et al. 2006;Goyal et al. 2011).
Several regression-based and other statistical methods have been proposed and applied to determine the yield stability and magnitude of GEI (Finlay and Wilkinson 1963;Piepho 1999;Goyal et al. 2011;Mohammadi and Amri 2013). The classical stability analysis, such as Finlay and Wilkinson (1963), used a coefficient of regression βi value to measure the stability of a cultivar in multi-environments. A cultivar with βi < 1.0 is categorized as having above average stability and is particularly adapted to low performing environments; βi > 1.0 is classified as below average stability and is specifically adapted to high performance environments, and βi = 1.0 has average stability and is well-or poorlyadapted to all environments depending on having a high or low mean performance. Another approach suggested by Eberhart and Russell (1966) uses both the regression coefficient and deviation from linear regression S 2 d values to identify stable genotypes, which can be considered stable if βi = 1.0 and S 2 d = 0.0. These regression methods are best to predict stability performance across large environments using linear or nonlinear data.
A grouping methodology previously described by Francis and Kannenberg (1978) and later adapted to agronomy studies (Gan et al. 2009;Beres et al. 2010;May et al. 2010) has also been used to explore treatment responses over environments. More recently, linearbilinear models, including the additive main effects and multiplicative interaction (AMMI) model and the Sites Regression (SREG) model have been used in applied research to dissect GEI and stability parameters in wheat and other crops (Crossa 1990;Piepho 1999;Burgueno et al. 2000). The AMMI model combines ANOVA (with additive parameters) and principal components analysis (PCA) (with multiplicative parameters) into a single analysis (Zobel et al. 1988;Gauch and Zobel 1997;Gauch 2006), and the SREG model accounts for both the genotype main effects and GE effects (Yan et al. 2007). The genotype plus genotype by environment (GGE) biplots of two principal components (PC) generated either from a SREG model or an AMMI model are useful techniques to visually identify similarity and dissimilarity among genotypes or environments (Yan 2001). The first two PC axes generated in a GGE biplot have been used to display two-way data and allow visualization of the interrelationship of environments, genotypes, and their interactions. The GGE plots also allow the identification of highly discriminating and representative test environments (Yan et al. 2007). The GGE studies declare a stable genotype should have large primary effects (high mean yield) and near zero secondary effects, whereas, ideal environments (sites) are those with high power to discriminate and will have large primary effects and small secondary effects (Burgueno et al. 2000). Our objective was to explore variations in responses using an array of statistical models for stem solidness, grain yield, and grain protein concentration of solid-and hollow-stemmed wheat genotypes when grown in WSS regions in the northern Great Plains.

Materials and Methods
Experimental sites, design, and plant materials Field experiments were conducted from 2011 to 2013 at six study locations in Canada and the USA where wheat production is traditionally affected by WSS. The locations were Lethbridge (49°41′N, 112°45′W) and near Bow Island (49°44′N, 111°20′W), Alberta; Swift Current (50°17′N, 107°79′W) Saskatchewan, in Canada; Amsterdam (45°73′N, 111°33′W) and Loma (45°95′N, 110°48′W), Montana and Hettinger (46°00′N, 102°38′W), North Dakota, in the United States. The soil at the Lethbridge site is an Orthic Dark Brown Chernozem clay loam soil (Typic Boroll) with approximately 3.0% organic matter content and a pH of 7.5. Bow Island is a Brown Chernozem loam soil (Aridic Boroll) with around 2.0% organic matter content and a pH of 6.0. Swift Current is an Orthic Brown Chernozemic silt loam (Aridic Haploboroll) with 3.0% organic matter content and a pH of 6.5. Amsterdam has silt clay loam soil (Typic Haplustolls) with 3.5% organic matter content and a pH of 6.6-7.3, and Loma is similar with clay to silt loam, 3.0% organic matter, and a pH of 6.5-7.0. Hettinger is a Shambo-loam soil with 3.6% organic matter content and a pH of 6.1-7.0. Commercially available spring wheat genotypes and breeding lines were selected from both Canada and the USA. Genotypes conferring the stem solidness trait were all derivatives of the historical landrace solid-stemmed line, 'S-615'. The genotypes (10) were arranged in a randomized complete block design (RCBD) with at least four replicates per site (Table 1).
Sites within provincial or state boundaries were considered a unique agroecosystem based on the relative similarity of soil and weather-related variables, as well as host plant responses to WSS.

Experimental procedures and data collection
Experimental units were approximately 90 cm wide (four-rows, spaced 23 cm apart) and 3-5 m long, depending on the site collaborator's access to land and plot equipment. Plots at all locations were usually seeded directly into standing stubble of the previous crop, which was typically canola, durum, or spring wheat (depending on site-year), using a zero-tillage plot drill configured with disc or single-shoot hoe-type openers. Plots were sown at a density of approximately 300 seeds·m −2 in late-April to early-May in each year ( Table 1). The experimental area was treated with glyphosate (RoundUp®, Monsanto, and St. Louis, MO) a few days prior to seeding at a rate of 900 g a.i.·ha −1 using a motorized sprayer calibrated to deliver a carrier volume of 45 L water·ha −1 at 275 kPa pressure.
All plots were maintained using local best management practices and received amendments of macroand micro-nutrients based on residual soil levels and plant nutrient requirements. In-crop herbicides were chosen based on the weed spectrum present at each site-year and applied in early-June at label rates. No insecticides were used at any site during the study period. To manage foliar and root diseases, seeds were treated with fungicides (Dividend XL RTA, Syngenta Crop Protection Canada) prior to sowing.
Weather data (precipitation, temperature, and light intensity) were collected at each site using Hobo a Cultivar registered as a solid-stemmed cultivar but, now classified as hollow-stemmed spring wheat.
Pendant temperature and light loggers (Onset Computer Corporation, Bourne, MA). The pendants were mounted on wooden stakes configured with a pigtail style stake at the top and secured in a horizontal position above the canopy. A pendant was placed in each replicate range between two plots near the centre of the range, initialized with a 30-min logging interval and installed when plant growth reached the two-leaf stage (Zadoks 12 GS) (Zadoks et al. 1974). All weather data are summarized in Fig. 1 and Table 2.
To ensure an adequate estimate of stem solidness ), a 0.50 m section of plant row was collected in late-July or early-August  from the two random locations in each plot to determine stem diameter and degree of stem solidness in the culm of the main stem. Each main stem was then split lengthwise starting from the crown to the neck, and was visually assessed at each internode for pith development. Ratings were as follows: 1-hollow-stem, no pith development; 2-some degree of pith development, may appear "cotton-like"; 3-moderate pith development with large hollow tunnel in the stem, or, a large cavity at a particular point in the internode; 4-continuous pith with the size of hollow equivalent   to a pencil lead, or, some cavitation has occurred at a particular point in the internode; and 5-solid-stem (DePauw and Read 1982). The average of four bottom to top internodes were used to estimate mean solidness rating using 30-50 randomly selected plants from each plot. In Montana, due to severe WSS infestations-related stem cutting damage, the treatments were rated for "% stems cut" and corresponding samples were assessed for "% infestation". Plots were harvested at crop maturity using a Wintersteiger Expert (Wintersteiger AG, Salt Lake City, UT) plot combine equipped with a straight cut header, pickup reel, and crop lifters. Grain yield was calculated from harvested samples representing the entire plot and corrected to 14% seed moisture. Post-harvest measurements also included seed weight and grain bulk density. To assess a quality parameter, grain protein concentration was determined from whole grain using near-infrared reflectance spectroscopy technology (Foss Decater GrainSpec, Foss Food Technology Inc., Eden Prairie, MN).

Data analyses
All agronomic and grain protein concentration data from 17 site-years were analyzed using the MIXED procedure of SAS (version 9.4) software. The data from the Amsterdam 2013 plot were lost due to hail damage. Prior to the MIXED procedure, homogeneity of error variances was tested using the UNIVARIATE procedure of SAS; any outlier observations detected using box plots, quartile-quartile (Q-Q) plots, and histogram were removed before a combined analysis over years and environments. Normality assumptions using the Shapiro-Wilk test were also tested for the categorical data "pith expression" as multiple categories were used and values were generally not extreme (Cochran 1954). For analyses within site-year, replicate was considered random, while genotype effects were considered fixed and significance was declared at P ≤ 0.05. For analyses within agroecozones (year and provincial or state zone), we considered each agro-ecozones and year as one environment, environment and agro-ecozones, their interaction considered random, while genotypes effects were considered fixed. Significance was declared at P ≤ 0.05. The results by environment indicated similar treatment response patterns among environments; therefore, a combined analysis was performed with replicate, years, environments, and their interactions considered random effects and genotype effects treated as fixed effects and significance declared at P ≤ 0.05. Genotypic least square means were computed using the LSMEANS statement in SAS. If fixed effects were significant, mean comparison tests were performed using Tukey's honestly significant difference (HSD) values at P ≤ 0.05. The COVTEST test statement in PROC mixed was used to compute tests of significance for random effects. Response variable least square means generated for each environment (agro-ecozone and year) were used to create a Pearson's correlation coefficient matrix of stem solidness, grain yield, and grain protein concentrations data using the CORR procedure of SAS.
Stability analyses for stem solidness, grain yield, and grain protein concentration were computed using multiple statistical procedures. A grouping methodology developed by Francis and Kannenberg (1978) was used to further explore treatment responses over environments. In this method, the mean and coefficient of variation (CV) were estimated for each treatment, and then means were plotted against the CV for each treatment. The overall mean of the treatments and CVs were used to categorize the biplot ordination area into four quadrats/ categories: Group I: High mean, low variability (optimal); Group II: High mean, high variability; Group III: Low mean, high variability (poor); and Group IV: Low mean, low variability (very poor).
The regression models originally introduced by Finlay and Wilkinson (1963) and Eberhart and Russell (1966), using two stages of analysis as described by Piepho (1999), were used to further determine genotype stability. For this analysis, we considered each agro-ecozone and year as one environment. An AMMI model using the Agricolae R Package (R Development Core Team 2019) was used to determine the magnitude of the main and treatment interaction effects generated by additive main effects and multiplicative interactions (AMMI). In this model, each agro-ecozone and year was considered as one environment, the effects of genotype (G), environment, and genotype by environment (G × E) interactions were considered as random. In this case, the best linear unbiased prediction (BLUP) of G and G × E effects was calculated (Piepho et al. 2008). The components of genotypic variance (σ 2 g), variance of G × E interaction (σ 2 i), and residual (σ 2 e) were estimated by the method of restricted maximum likelihood. The genotype by environment interaction was partitioned into two principal components effects (IPCA1 and IPCA2). Stable genotypes across site-years were identified by analyzing the contribution of the variation to the total sums of squares. The AMMI 1 biplots were created using the mean of a parameter with its respective IPCA1 (PC1) genotype scores tested across environments (Figs. 5A, 5C, and 5E) Abamu and Alluri 1998;Torres and Henry 2018). The AMMI 2 biplots, PC1 vs. PC2, were drawn to understand the differentiating environment and responsive genotypes (Figs. 5B, 5D, and 5F) ). The ranking of genotypes was conducted using both AMMI stability values (ASV) and yield stability index values (YSI) (Sabaghnia et al. 2008;Farshadfar et al. 2011). The ASV value is the distance from zero in a two-dimensional scatter graph of IPCAI (interaction principal component analysis axis 1) scores against IPCA2 scores, whereas yield stability index (YSI) is calculated by ranking the mean grain yield of genotypes (RY) across environments and rank of AMMI stability (rASV) value. The YSI incorporates both mean and stability in a single criterion (YSI = rASV + RY). A low value of YSI identifies stable genotypes with a high mean yield environment.
Finally, a SREG model was also performed where the genotypic least square means generated in the mixed model analysis of individual agro-ecozone and year (environment) were used for analysis with SAS procedures developed by CIMMYT and the GGEBiplotGUI R package (R Development Core Team 2019). The genotype main effect (G) plus genotype × environment biplot, commonly known as a GGE biplot (Yan 2001), was generated from the SREG model in the GGEBiplotGUI package with the first two principal components (PCs). These two principal components, PC1 and PC2, referred to as the primary and secondary effects, were derived from subjecting double-centered data to singular value decomposition (Yan 2001). Following the graphical method, a polygon was formed on the biplot by connecting markers for the most responsive cultivars. The biplot was then divided into sectors by drawing perpendicular lines from the origin of the biplot to the sides of the polygon. The perpendicular lines divide the polygon in sections (quadrants) and facilitate comparison of the two vertex genotypes connected by its respective polygon side (Yan 2001).

Environmental conditions
A wide range of environmental conditions was encountered over the 17 site-years of the study. The average monthly temperature, light intensity, and total precipitation pattern is summarized in Fig. 1 and Table 2. In general, the Swift Current, SK, site received higher precipitation while Amsterdam and Loma in Montana had less precipitation compared with other sites in all years (Table 2). Trends in light intensity and temperature were similar at all sites, with the exception of Hettinger, ND, and Loma, MT. While mean temperatures for other sites peaked in July, Hettinger did not record the highest mean temperature until August in 2011 and 2012 (Fig. 1). Light intensity levels peaked in either June or July for all sites, and often correlated to growing season precipitation. For example, the lowest light intensity for the entire study period was observed in Loma, MT, 2013, which was also characterized by high precipitation accumulation. A similar relationship was observed at Bow Island, AB, in 2012; however, low light intensity at Amsterdam was incongruent with precipitation for the same period (Fig. 1).

Stem solidness
Wheat genotypes differed among agro-ecozones and years for stem solidness at each internode and when averaged across the four lower main stem internodes (Table 3, Fig. 2). As expected, stem solidness was greater in solid-stemmed genotypes across all internodes, but responses varied. The highest overall stem solidness was observed in 'BW925' and 'Choteau' (Table 3). In contrast, the lowest solidness expression was observed in the hollow-stemmed genotypes 'Infinity', 'Reeder', and 'AAC Bailey' and the solid-stemmed landrace line 'S-615', which is the source of stem solidness for all registered cultivars in North America. 'Mott' displayed higher stem solidness similar to 'BW925' and 'Choteau' in
Montana. Low pith expression was observed in 'AAC Bailey' and 'S-615' in all agro-ecozones, particularly in North Dakota (Fig. 2). Except for Montana, notable WSS infestations or stem cutting damage were not observed in any site-year. In Montana, greater WSS infestations were observed in 2011 and moderate WSS infestations were observed in 2012 and 2013. Among the spring wheat genotypes, and in all years, the lowest WSS infestation and incidence of stem cutting damage were observed in 'Mott', 'BW925', and 'AC Eatonia' (Table 4).
Stability of stem solidness among genotypes across environments and years is an important consideration; therefore, different models were explored to determine if methods employed displayed alternative outcomes. The results from the Francis and Kannenberg method for constructing mean vs. CV biplots are summarized in Fig. 3A, which indicate that 'BW925' and 'Choteau' consistently displayed high pith expression (averaged over the lower four stem internodes), i.e., Group I (Fig. 3A). In Alberta and Saskatchewan, 'AC Eatonia' and 'Fortuna' also displayed high solidness and low variability across site-years; whereas, in Montana and North Dakota, 'Mott' displayed high stem solidness similar to 'BW925' and 'Choteau' (Supplementary Fig. S1 1 ).
The Finlay-Wilkinson method defines a well-adapted stable cultivar as one with high mean yield (trait mean) and a regression coefficient that is within one deviation of the standard error of the mean regression coefficient. Our results for the Finlay-Wilkinson test indicated the mean regression coefficient was 0.36, and 'BW925', 'Choteau', 'Mott', 'AC Eatonia', and 'Fortuna' possessed the highest mean solidness and fell within one deviation of the standard error of the mean regression coefficient (Fig. 4A). Among these lines, 'BW925' and 'Choteau' were classified as the most stable genotypes, since these genotypes had the highest overall and consistent stem solidness. 'Reeder', 'Infinity', and 'AAC Bailey' were classified as the most unstable cultivars as they displayed low and variable stem solidness.
The results of Eberhart-Russell stability analysis for mean stem solidness rating are summarized into two graphs (Figs. 4B and 4C). The results obtained by this model for solidness were similar to Finlay-Wilkinson, whereby the regression coefficient was regressed upon the cultivar's mean pith expression. The differences between these two methods are shown in Fig. 4C, where the deviation from the regression coefficient were regressed upon the cultivar mean value to account for cultivar response patterns that may not be linear.
Using this method, 'BW925', 'Mott', 'Fortuna', and 'Choteau' were more stable for solidness expression. The highest and most stable response was for 'BW925', falling within one deviation of the regression coefficient (Fig. 4C).
Additive main effects and multiplicative interaction (AMMI) analysis indicates all components; genotype, environment, and genotype by environment were highly significant for mean solidness (averaged across the lower four stem internodes) (Supplementary Table S1 1 ). The portioning of the total sum of square (TSS) indicates that the genotype effect was a predominant source of variation, which explained 64.5%, followed by the environment (11.8%) and GEI (11.3%). The first three interactions with environment PC (IPCA1, IPCA2, and IPCA3) were highly significant and explained the majority of variation. IPCA1 explained 42% of the variability related to GEI and 19% of the interaction of degrees of freedom. Likewise, the second IPCA2 accounted for 29% of the  variability and IPCA3 explained 17.6% of the variability in GEI. Jointly, they accounted for almost 89% of GEI variability (Supplementary Table S1 1 ). Genotypes with IPCA1 scores near zero had little interaction with environment, while genotypes with very high IPCA1 had considerable interaction with environments. The AMMI biplot 1 analysis indicates genotypes 'Choteau' and 'BW925' with high mean and positive high PC interaction were more responsive for the most environment. Conversely, the genotype 'Reeder', 'Infinity', and 'AAC Bailey' with low mean and negative PC scores indicate they were undesirable genotypes for solid-stemmed traits (Fig. 5A). 'Choteau' and 'BW925' displayed a higher mean yield and large PCI scores but were less stable than 'AC Eatonia', 'S-615', 'Fortuna', and 'Mott'. While the latter-mentioned genotypes displayed inferior grain yield, they did have small interaction effects, indicating that these genotypes were less influenced by the environment (Fig. 5A).
The environments with minor G × E interactions for stem solidness were AB 013, MT 013, and ND 011, whereas G × E interactions were notable at MT 011, MT 012, ND 012, and SK 012. All other environments were generally intermediate with respect to the magnitude of a G × E interaction. The maximum stem solidness expression was recorded in ND 012, MT 012, and MT 013, Fig. 4. Comparisons of three different models to estimate stability of stem solidness, grain yield, and grain protein concentration across environments. Finlay-Wilkinson yield stability model used for (A) stem solidness rating (average over whole stem), (D) grain yield (Mg·ha −1 ), and (G) grain protein concentration (g·kg −1 ) of spring wheat genotypes grown across the agro-ecozones in North America in 2011-2013. Eberhart-Russell stability model used for (B) stem solidness rating, (E) grain yield (Mg·ha −1 ), and (H) grain protein concentration (g·kg −1 ) for spring wheat genotypes grown across the agro-ecozones in North America in 2011-2013. Eberhart-Russell stability model used for (C) deviation from regression vs. genotype mean stem solidness, (F) grain yield (Mg·ha −1 ), and (I) grain protein concentration (g·kg −1 ) as the parameter for spring wheat genotypes grown across the agro-ecozones in North America in 2011-2013. MSE, mean coefficient of regression means standard errors obtained by Finlay Wilkinson; MCRE, mean coefficient of regression means standard errors obtained by Eberhart-Russell; MDRSE, mean deviation from regression stability parameters mean std errors; AMPR, average mean pith rating for all wheat genotypes; AMGY, average mean grain yield for all wheat genotypes; AMGP, average mean grain protein for all genotypes.
indicating that these were favorable environments for high pith expression (Fig. 5A). The factors leading to high pith expression in these environments were likely high light intensity and temperatures with relatively low precipitation (Table 2 and Fig. 1).
The GGE biplot results from SREG model are summarized in Fig. 6, which shows the association between stability parameter PC2 and mean stem solidness (PC1) (Fig. 6A). The portion of sums of squares accounted by each PC1 and PC2 were 49.5% and 22.5%, respectively. The biplot accounted for 71% of GE. The biplots identify the best performing genotypes in each environment and group of environments. Those genotypes that form the polygon were the most responsive and representative of the best or poorest performing at some and all environments (locations). Perpendicular lines divided polygons into six sections. The proximity to the origin of the vertices between PC1 and PC2 indicated the discriminating power of the sites. The scores of the genotypes furthest from the origin were connected to form a polygon (Fig. 6A), with all other genotypes contained within the polygon, indicating which genotypes "won where" based on their association with the site scores (Yan et al. 2007). For example, 'Mott' and 'AC Eatonia' were the most responsive in the Montana agroecozone in all years (Fig. 6A). Except for Montana and Alberta in 2011, 'BW925' and 'Choteau' were most stable for the stem solidness trait. According to Yan et al. (2007), in the SREG-based GGE biplot, the best genotype is the one with large PC1 scores (high mean) and near zero PC2 scores (high stability). In this study, 'S-615', 'BW925', and 'Choteau' were the best as indicated by values close to zero on the Y-axis, but 'S-615' displayed lower pith expression that did not meet the optimum threshold criteria to be considered optimal. 'Choteau' and 'BW925' were considered the best performing lines for stem solidness as they had large PC1 values and close to zero PC2 values. The GGE biplots can be also a useful tool for the identification of ideal test environments (Yan et al. 2007). The best selection environments are those with small (absolute) PC2 scores (more representative of the overall environment) and large PC1 scores (more power to discriminate genotypes in term of the genotype main effect). Based on these criteria, the conditions observed in AB (2012) and SK (2012 would be considered ideal environments for characterizing pith expression.
The Eberhart-Russell analysis results for yield are summarized into two graphs (Figs. 4E and 4F). The results obtained from the first approach were similar to the Finlay-Wilkinson's result, where the regression coefficient was regressed with the cultivar mean yield. The only differences in results between these two methods for yield are illustrated in Fig. 4F; the deviation from regression coefficient was regressed upon cultivar mean yield to account for the cultivar response pattern that may not be linear. Our results indicate that 'Mott' had the highest yield and yield stability followed by 'Reeder', 'Choteau', 'BW925', and 'Lillian' (Fig. 4F).
Similar to stem solidness, the AMMI results for genotype, environment, and genotype by environment highly influenced grain yield (Supplementary Table S3 1 ). In contrast to stem solidness, the environment explained the largest variability (83.3%) to the total sum of squares, GEI (9.0%), and genotype (3.1%), respectively. This indicates there was a low cross over interaction effect. The first two principal components for interaction (IPCA1 an IPCA2) were significant and explained the most variability. IPCA1 and IPCA2 explained 36.8% and 21.8% of the variability relating to GEI, respectively. These PC components jointly accounted for 58.4% of the variability in GEI (Fig. 5D).
The AMMI 1 biplot indicates 'Reeder' and 'Infinity' had higher grain yields and high PC1 scores and are adapted to high-yielding environments but displayed lower stability (Fig. 5C). Conversely, 'Choteau', 'Mott', and 'AC Eatonia' displayed high grain yield and close to PC1 zero score, indicating that these cultivars were more stable and well-adapted to most growing environments. The most stable growing environments for yield were SK 012, SK 013 and AB 011, as they had a relatively small G × E interaction, whereas ND 011, ND 012, ND 013, MT 011, MT 012, and MT 013 were most variable for yield and had a corresponding high contribution to G × E interaction (Fig. 6C). The maximum yield was recorded in ND 012, AB 012, and AB 013, indicating that these were favorable environments to obtain high grain yield (Fig. 5C).
The GGE biplot results for grain yield produced PC1 and PC2 values of 53.9% and 19.06%, respectively (Fig. 6B), which accounted for 73.1% of GE. GGE biplots for grain yield (Fig. 6B) shows that except for the North Dakota agro-ecozone, the solid-stemmed wheat genotypes 'BW925', 'Choteau', and 'AC Eatonia' were more responsive and the best cultivars for higher grain yield in all agro-ecozones. 'Mott' and 'Choteau' were most stable with yield values close to zero for PC2 scores, and 'Choteau' was considered the winning and most stable cultivar.
In contrast to yield, the CV vs. mean biplot for grain protein indicated the solid-stemmed cultivars 'Lillian' and 'AC Eatonia' consistently produced higher grain protein with low CVs across the agro-ecozones (Table 3 and Fig. 3C, Supplementary Fig. S3 1 ). In Montana, 'Reeder' and 'AAC Bailey' had higher grain protein similar to 'AC Eatonia'. The Finlay-Wilkinson analysis for protein resulted in a mean regression coefficient of 1.88, and 'AC Eatonia' and 'Lillian' produced high mean grain protein concentrations and had regression coefficients within one standard error deviation (Fig. 4G). Therefore, both cultivars were considered stable and with the highest overall protein concentration. The Eberhardt-Russell analysis for protein concentration (Figs. 4H and 4I) shows that the cultivars 'AC Eatonia', 'Lillian', and 'Mott' produced higher and stable grain protein concentration across the agro-ecozones. 'Lillian' and 'AC Eatonia' were the top-performing cultivars for grain protein concentration.
The AMMI analysis for protein indicated that only genotype and environment were highly significant for grain protein (Supplementary Table S5 1 ). Like yield, environment contributed the largest variability (73.0%), followed by genotype (2.25%) and GEI (3.44%) to the total sum of squares (Fig. 5F). The AMMI 1 biplot shows 'Choteau' had high PC1 score and moderate protein concentration and were more unstable. The cultivar 'AAC Bailey' had PC1 score close to zero and similar protein concentration and had more stable and consistent protein concentration for all environments. The most stable growing environment for grain protein concentration were MT 011, MT 013, and AB 011 as they had a small G × E interaction, whereas conditions at SK 013, SK 011, and ND 013 produced unstable grain protein value and had a high contribution to G × E interaction (Fig. 6C).
The AMMI 2 biplot indicates genotype 'BW925' was more responsible for protein concentration for six growing environments (MT 011, MT 012, MT 013, AB 013, SK 011, and SK 013) although 'BW925' had lower protein concentration than other genotypes. The genotypes 'AAC Bailey' and 'AC Eatonia' were best and most responsive to protein concentration in all years in North Dakota and two growing years in Alberta and Saskatchewan agro-ecozones (Fig. 5F). The ASV (AMMI stability value and ranking) values were low for 'AAC Bailey', while 'Lillian' and 'AC Eatonia' had the highest protein stability index (Supplementary Table S6 1 ). The GGE biplots results for grain protein in the PC1 and PC2 were 35.4% and 27.84%, respectively (Fig. 5F), which jointly accounted for <70% of the GE. The GGE biplots shows that 'BW925' had low but more consistent grain protein concentration in seven growing environments (MT 011, MT 012, MT 013, AB 013, SK 011, and SK 013). Genotypes 'AAC Bailey' and 'AC Eatonia' were more responsive for grain protein in five growing environments (ND 011, ND 012, ND 013, SK 012, and AB 012 (Fig. 6C). Although 'BW925', 'Mott', and 'S-615' accumulated consistent grain protein, it was not at a higher level like 'AC Eatonia' and 'Lillian'. No correlations were observed between stem solidness and grain yield or grain protein concentration and stem solidness. As expected, a negative and significant correlation was observed between grain yield and grain protein concentration (r = −0.38, P < 0.001) (Supplementary Table S7 1 ).

Discussion
Regression, AMMI, and GGE biplot analyses are tools that provide a reasonable estimation of stability and performance of cultivars, particularly for yield and other complex traits in cereals (Zobel et al. 1988;Yan 2001;Gauch 2006). In this study, we used multiple statistical models to explore the stability and variability of stem solidness, grain yield, and protein concentration of solid-stemmed spring wheat genotypes that are used to manage WSS infestations. The focus of our study was to examine the stability of pith expression, which is the main pillar of an integrated pest management strategy for C. cinctus. Adoption of solid-stemmed cultivars by producers will only occur if the genetic package also optimizes agronomic performance with respect to grain yield and quality; therefore, these traits were also included in our analysis. Wallace et al. (1973) suggested that a threshold rating of 3.75 was required to achieve effective resistance against WSS. The problem is that many of the cultivars grown today fall below this threshold in some or all environmental conditions. Thus, a better understanding of the stability of pith expression is required to help minimize the risk of losses during an outbreak .
Using the regression-based model and the GGE biplot, the Canadian genotype 'BW925' and the US cultivar 'Choteau' consistently expressed high and stable pith greater than a rating of 3.75 across the agro-ecozones. Both genotypes derive their stem solidness trait from 'S-615', which was used to develop the first North American commercial solid-stemmed wheat cultivar, 'Rescue' (Platt et al. 1948). In general, genotypes that display low ASV are considered more stable, while those with high values are less stable (Sabaghnia et al. 2008). Stability does not warrant selection since a consistently low-yielding genotype is not normally desirable. A yield stability index that integrates high mean yield and stability across environments into a single index would facilitate rapid selection of stable varieties (Zobel et al. 1988;Purchase et al. 2000). In our study, 'BW925' and 'Choteau' are considered the most desirable and stable genotypes in terms of pith expression as they had the lowest (YSI) ranking index combined with both high mean pith expression and stability along with large PC1 scores and near zero PC2 scores (Figs. 5A and 6A). Some breeders use the lowest YSI index to assess genetic stability for yield and related traits (Sabaghnia et al. 2008;Farshadfar et al. 2011). A number of studies have shown that the best adapted and stable genotype is the one with large PC1 scores and near zero PC2 scores (high stability) (Yan et al. 2007;Beres et al. 2008;Goyal et al. 2011). Furthermore, our results of AMMI 2 and GGE biplots confirm that 'BW925' and 'Choteau' had stable and optimum pith expression even during unfavorable, wet, and cloudy weather conditions (i.e., North Dakota 2013), which usually inhibits pith development in solid-stemmed wheat cultivars (Beres et al. 2017). Environmental factors such as high precipitation, cloudy conditions, and reduced photoperiod tend to impede pith expression (Holmes 1984;Beres et al. 2017). Some research such as DePauw and Read (1982) and Beres et al. (2012) reported that environmental factors were more important for pith expression than fertility management, as levels of nitrogen that would cause excess canopy development and shading, and subsequently inhibit pith expression, exceeded crop nutrient requirements. As such, our results suggest that 'BW925' and 'Choteau' are the most stable in their expression of stem solidness across WSS affected environments. The other solid-stemmed genotypes such as 'AC Eatonia', 'S-615', and 'Mott' also had stable pith across environments, but they did not meet the required threshold rating of 3.75/5 for pith development (Wallace et al. 1973).
Interestingly, our results showed that 'Mott' expressed the highest pith in internodes 1 and 2 (>3.75 ratings) compared with later-developed internodes, which would make it more solid and protective from WSS in the early season. Growing 'Mott', therefore, may offer optimal protection from WSS when an outbreak occurs earlier in the season. In general, in Montana, WSS infestation can occur as early as late-May to mid-June when spring wheat stem elongation has been initiated (one or two internode stages) (Fulbright et al. 2017). Thus, early stages of spring wheat seem to be more vulnerable to WSS compared with later stages due to WSS ovipositional preferences (Holmes and Peterson 1960). The synchronization of WSS emergence and plant development by modifying the planting dates may also reduce WSS damage (Morrill and Kushnak 1999;Morrill et al. 2001). Our results also indicate that 'Mott' may be more vulnerable as the season progresses, as it did not produce acceptable pith in the third or fourth internodes. Therefore, plant physiological responses like this may allow for the partitioning of greater resources to form grain if pith development subsides. 'Mott' generally displayed high yield with optimal stem solidness occurring only in the early part of the growing season. It has been observed that some solid-stemmed cultivars lose solidness by the time they reach maturation. This interaction with phenology suggests breeding programs could develop genetics tailored for a specific agroecozone where solidness is required at a specific growth stage and then terminated in favour of partitioning resources to yield and quality. In contrast, other cultivars such as 'Choteau' display uniform pith expression through all stages of reproductive plant development. These phenomena suggest that traditional selection approaches for WSS resistance may need to be modified and based on regional-specific requirements for when optimal pith expression is critical (Varella et al. 2016).
The consistent expression of stem solidness in 'BW925' and 'Choteau' and variable pith expression in other genotypes in our study indicates the genotypic effect may be the predominant and largest source of variation for stem solidness. This conclusion agrees with the findings of Nilsen et al. (2017), who reported that the interaction of minor genes with the major SSt1 gene located in 3B contributed 78%-92% variability for pith stability in some lines of Canadian durum and spring wheat populations. The major QTL Qss.msub-3BL explained at least 76% of the variation for stem-solidness in a winter wheat mapping population and contains multiple alleles conferring varying levels of stem solidness (Cook et al. 2004). The strong expression of stem-solidness in 'Choteau' over other cultivars was influenced by presence of both Qss.msub-3BL and Qss.msub-3DL localized in chromosomes B and D (Lanning et al. 2006). The underlying genetics for stemsolidness in the 'S-615' source are quite complex, and may be controlled by the action of a major gene coupled with four or more additional recessive genes (Larson and MacDonald 1959;Clarke et al. 2002). The phenotypic expression of stem solidness in hexaploid wheat results from quantitative factors with an additive nature rather than dominance or a recessive kind of gene action (Bainsla et al. 2020). Due to the heterogeneity in SSt1 in 'S-615' (Beres et al. 2013) and genetic suppression effects in some wheat backgrounds, any 'S-615' derived cultivars can suffer from inconsistent pith expression (Larson and MacDonald 1959). However, in some genotypes, stem solidness is predominately controlled by major QTLs, which confer high heritability (h 2 = 0.83) to this trait (Cook et al. 2018). This would support the concept that direct selection of stem solidness is possible to mitigate stem cutting damage by WSS.
To ensure widespread adoption of solid-stemmed cultivars, they also must possess acceptable or superior agronomic attributes that are consistently displayed across an array of WSS prone environments. Regressionbased analyses indicates that 'Mott', 'Choteau', 'Lillian', and 'BW925' all consistently expressed high and stable grain yield. Finlay and Wilkinson (1963) and Eberhart and Russell (1966) declared that the most stable genotypes should have high mean grain yield with a coefficient of regression βi close to 1 and the deviation from linear regression should be zero. These methods have been used to measure yield stability in a number of cereal crops including soft white spring wheat and triticale (Beres et al. 2008;Goyal et al. 2011). Furthermore, 'Choteau', 'Mott', and 'AC Eatonia' had low ASV values and high mean grain yield with near zero PC1 (Fig. 5), and it was confirmed that these genotypes had stable yield in most growing environments. The ASV and YSI values are often used in many cereals crops to determine yield stability and crossover interaction while selecting the best genotypes for the target areas (Gauch and Zobel 1997). We also observed high yield performance in the Canadian line 'BW925' and the US cultivar 'Mott', particularly in Montana where moderate to heavy WSS pressure was observed. This suggests 'BW925' and 'Mott' seem to possess greater yield potential and would be well-adapted in a region where WSS infestation levels are consistently high. These results were comparable to those from Beres et al. (2007), who reported that some solid-stemmed common and durum wheat could be agronomically superior and even produce higher yield under sawfly pressure. Conversely, 'Choteau', 'Lillian', and 'AC Eatonia' would be well-adapted where WSS pressure is low to moderate.
Our results for grain protein concentration differed somewhat depending on the model used, but generally 'AC Eatonia' and 'Lillian' were the only genotypes that consistently produced higher grain protein concentration across the range of environments. The findings for higher grain protein in 'Lillian' were not unexpected, as 'Lillian' carries the high protein gene Gpc-B1 from its parent '90B07-AU2B' and was the first solid-stemmed Canadian cultivar with high yield and protein concentration similar to hollow-stemmed cultivars (DePauw et al. , 2007. In contrast, the AMMI model and GG biplots indicated that 'Lillian' displayed high and consistent protein concentration in some environments such as two site-years in North Dakota (i.e., ND 013) that experienced very wet and humid weather conditions (Table 2). In general, good water availability increases grain yield but reduces grain protein concentration (Welch 2011). Several studies, such as Flagella et al. (2010), reported that water deficit or drought conditions throughout the growing season drastically increase grain protein concentration in durum wheat. In wet, warm soil conditions, nitrogen may be depleted from soils later in the growing season either by leaching or denitrification. Lack of adequate nitrogen during the grain filling stage will reduce protein accumulation (Lotfollahi and Malakouti 1999). The hollow-stemmed cultivar 'AAC Bailey' and solid-stemmed cultivars 'AC Eatonia' and 'BW925' displayed stable grain protein concentration in most growing environments. However, the grain protein concentration in 'BW925' was at times less than the ideal level 135 g·kg −1 . A common issue with solid-stemmed bread wheat genotypes in Canada is suboptimal quality characteristics for meeting the parameters to be included in the designated list for the Canada Western Red Spring (CWRS) premium export class. When a new market class was established by the Canadian Grain Commission, Canada Northern Hard Red, all solidstemmed hexaploidy cultivars were moved down to this lower quality class, leaving growers of CWRS without a WSS tolerant variety. The negative correlation observed in most genotypes between grain yield and grain protein was expected, as it is prevalent in wheat and other crops (Terman et al. 1969;Oury and Godin 2007;Bogard et al. 2010).

Conclusions
This study was conducted to explore the stability of pith expression, grain yield, and grain protein content of solid-stemmed spring wheat cultivars that are used to manage WSS infestations in the northern Great Plains of North America. To our knowledge, this is the first comprehensive report to explore G × E interactions and genetic stability using an array of statistical models for pith expression of wheat genotypes grown across large geographical agro-ecozones. Most models identified the genotypes 'Choteau', 'BW925', 'AC Eatonia', and 'Mott' as stable with respect to pith expression across all environments. Cultivars 'Choteau', 'Mott', 'Lillian', and 'BW925' had comparable yield with the hollowstemmed cultivar 'Reeder'. Only 'AC Eatonia' and 'Lillian' produced consistently higher grain protein concentrations. Based on a minimum desired pith level (|>3.75) for stem-boring tolerance to WSS, 'Choteau', 'BW925', and 'Mott' would possess both optimal stability and pith expression across the agro-ecozones concomitant with higher yield potential; however, these genotypes have relatively lower grain protein compared with other spring wheat genotypes. For protein premium and quality parameters to meet high quality bread-making market class designations, the challenge remains for breeders to develop new solid-stemmed cultivars that optimize pith, quality, and high yield traits. Our results also suggest that regression-based and other stability models are robust techniques for breeders or agronomists who want to measure stability of stem solidness parameters for whole agro-ecozones. Moreover, the AMMI and GGE biplots are unique as they are useful to identify the stable and best genotypes for specific and whole environmental regions. These approaches helped us to identify 'BW925' and 'Choteau' as lines that would offer stability across a range of environments, and 'Mott' as a line that will develop optimal pith at a specific phenological stage when resistance to WSS infestation is critical. Our results will help facilitate the decision-making process for the selection, production, and adoption of solid-stemmed wheat cultivars in regions prone to WSS infestations.