Winter hardiness of Miscanthus (III): Genome‐wide association and genomic prediction for overwintering ability in Miscanthus sinensis

Overwintering ability is an important selection criterion for Miscanthus breeding in temperate regions. Insufficient overwintering ability of the currently leading Miscanthus biomass cultivar, M. ×giganteus (M×g) ‘1993–1780', in regions where average annual minimum temperatures are −26.1°C (USDA hardiness zone 5) or lower poses a pressing need to develop new cultivars with superior cold tolerance. To facilitate breeding of Miscanthus, this study characterized phenotypic and genetic variation of overwintering ability in an M. sinensis germplasm panel consisting of 564 accessions, evaluated in field trials at three locations in North America and two in Asia. Genome‐wide association (GWA) and genomic prediction analyses were performed. The Korea/N China M. sinensis genetic group is a valuable gene pool for cold tolerance. The Yangtze‐Qinling, Southern Japan, and Northern Japan genetic groups were also potential sources of cold tolerance. A total of 73 marker–trait associations were detected for overwintering ability. Estimated breeding value for overwintering ability based on these 73 markers could explain 55% of the variation for first winter overwintering ability among M. sinensis. Average genomic prediction ability for overwintering ability across 50 fivefold cross‐validations was high (~0.73) after accounting for population structure. Common genomic regions for overwintering ability were detected by GWA analyses and a previous parallel QTL mapping study using three interconnected biparental F1 populations. One QTL on Miscanthus LG 8 encompassed five GWA hits and a known cold‐responsive gene, COR47. The other overwintering ability QTL on Miscanthus LG 11 contained two GWA hits and three known cold stress‐related genes, carboxylesterase 13 (CEX13), WRKY2 transcription factor, and cold shock domain (CSDP1). Miscanthus accessions collected from high latitude locations with cold winters had higher rates of overwintering, and more alleles for overwintering, than accessions collected from southern locations with mild winters.


| INTRODUCTION
The perennial C 4 grass Miscanthus is a promising bioenergy crop (Clifton-Brown, Chiang, & Hodkinson, 2008;Clifton-Brown, Stampfl, & Jones, 2004;Głowacka et al., 2014Głowacka et al., ,2015Heaton, Dohleman, & Long, 2008;Heaton, Voigt, & Long, 2004;Somerville, Youngs, Taylor, Davis, & Long, 2010). Key objectives for Miscanthus breeding are greater biomass yield and better adaptation to local production environments than the current leading biomass cultivar of M. ×giganteus (M×g). Currently, only a single triploid clone of M×g, which is derived from the interspecific hybridization between Miscanthus sacchariflorus and M. sinensis, is available to US farmers (Głowacka et al., 2015;Hodkinson & Renvoize, 2001). We refer to this genotype as M×g '1993-1780' in ref-erence to its accession number in the Kew Living Collection (Hodkinson & Renvoize, 2001), and it was first imported from Yokohama Japan to Denmark in the 1930s (Greef & Deuter, 1993;Linde-Laursen, 1993); in previous studies, we called this genotype 'Illinois', as we obtained our initial stock plant from the Chicago Botanic Garden (Maughan et al., 2011). Natural populations of Miscanthus are found from tropical and subtropical areas of East Asia and Oceania to ~50°N in eastern Russia (Clifton-Brown et al., 2008;Sacks, Juvik, Lin, Stewart, & Yamada, 2013). Such a wide distribution of Miscanthus provides a genetically diverse and valuable gene pool from which to breed new cultivars, including those of M×g.
Overwintering ability is an important selection criterion for perennial bioenergy crops in temperate environments (Burner, Tew, Harvey, & Belesky, 2009;Clifton-Brown et al., 2001). Insufficient overwintering ability is a consistent limitation of M×g '1993-1780' in regions with cold winters, such as those in USDA hardiness zone 4 (average annual minimum air temperature of −34.4°C to −28.9°C; U.S. Department of Agriculture, 2012) and lower, and an intermittent problem in hardiness zone 5 (average annual minimum air temperature of −28.9°C to −23.3°C) and zone 6 (average annual minimum air temperature of −23.3°C to −17.8°C) Dong, Green, et al., 2019;Dong et al., 2018;Dong, Liu, et al., 2019;Heaton et al., 2008;Jørgensen, Mortensen, Kjeldsen, & Schwarz, 2003;Lewandowski, Clifton-Brown, Scurlock, & Huisman, 2000). Clifton-Brown and  observed variation among two M×g, one M. sacchariflorus, and two M. sinensis grown at four field trial locations in Europe for ability to survive the first winter, and this was not associated with plant size or early senescence in the first autumn. Moreover, controlled environment freeze tests on dormant rhizomes (cold-acclimated) removed from the field during the winter indicated that the lethal temperature at which 50% of the rhizomes died (LT 50 ) was −3.4°C for M×g '1993-1780' and M. sacchariflorus, whereas the most winter-hardy M. sinensis genotype had an LT 50 of −6.5°C, which were consistent with prior observations from field trials . Similarly, Friesen, Peixoto, Lee, and Sage (2015) reported that the LT 50 for dormant rhizomes of M×g '1993-1780' was −4°C and it had only 10% survival at −8°C based on controlled environment tests. Heaton et al. (2008) noted that M×g '1993-1780' could survive in the field when air temperatures during the establishment winter were as low as −8°C, with only a 14% loss of plants observed at a northern Illinois site (Shabbona, IL) in the first year, and no losses observed at two more southern locations in year 1 or at any of the locations in the following years.
Miscanthus is typically most sensitive to winter damage during the establishment year (i.e., first winter). Christian and Haase (2001) hypothesized that the first-year plants of M×g '1993-1780' did not go dormant early enough in the across 50 fivefold cross-validations was high (~0.73) after accounting for population structure. Common genomic regions for overwintering ability were detected by GWA analyses and a previous parallel QTL mapping study using three interconnected biparental F 1 populations. One QTL on Miscanthus LG 8 encompassed five GWA hits and a known cold-responsive gene, COR47. The other overwintering ability QTL on Miscanthus LG 11 contained two GWA hits and three known cold stress-related genes, carboxylesterase 13 (CEX13), WRKY2 transcription factor, and cold shock domain (CSDP1). Miscanthus accessions collected from high latitude locations with cold winters had higher rates of overwintering, and more alleles for overwintering, than accessions collected from southern locations with mild winters.

K E Y W O R D S
breeding, genome-wide association analysis, genomic prediction, germplasm, Miscanthus sinensis, Overwintering ability season to avoid damage from cold temperatures. Boersma, Dohleman, Miguez, and Heaton (2015) found that when autumn temperature fell below 10°C in Ontario, Canada, CO 2 assimilation rate and photosystem II efficiency for first-year stands of M×g '1993-1780' were almost four times higher than for the third-year plants, and leaf [N] was about 2.4 times greater, suggesting that the first-year plants were still actively growing before a killing frost. From this photosynthesis data, Boersma et al. (2015) concluded that limited translocation of nutrients to rhizomes in the first-year plants of M×g '1993-1780' contributed to their poorer overwintering ability than the third-year plants in Ontario, though no measurements of rhizome nutrient contents were made. Genetic variation in Miscanthus for overwintering ability in the establishment year could be due to differences in freeze tolerance of dormant rhizomes, differences among genotypes in timing of dormancy and translocation of nutrients belowground, differences in depth of dormant buds in the soil to facilitate avoidance of freezing temperatures, or combinations of these traits.
Prior studies of overwintering ability and freeze tolerance have been constrained by access to only a narrow subset of Miscanthus genetic diversity. In studies of population structure, Clark et al. (2014Clark et al. ( , 2015 identified seven distinct genetic groups of M. sinensis from Asia, including three from Japan, three from China, and one from Korea and northern China; additionally, two M. sinensis groups from the US that were derived from the southern and central Japan groups, including US ornamental cultivars and US naturalized populations, and a group of natural diploid M. sacchariflorus × M. sinensis hybrids (i.e., diploid M×g) found in China were identified. Given that the M. sinensis germplasm that has been available in Europe and North America was derived entirely from the southern and central Japan populations (Clark et al., 2014, studies of hardiness have been limited to just a subset of the southern Japanese germplasm, which itself is a subset of the entire species' diversity. For M. sacchariflorus, the situation is similar, with prior studies of hardiness conducted almost entirely on tetraploids from southern Japan, with the exception of just one diploid genotype (M. sacchariflorus 'Robustus') that originated from near the Amur River around the border of northeastern China and southeastern Russia . Moreover, prior studies of M×g have focused on only one or a few genotypes because different genotypes were not available prior to recent efforts by our research group and others to collect and breed new ones.
Thus, an important objective for Miscanthus breeders who target temperate growing environments, such as the northern US, Canada, and northern Eurasia, is to identify species, populations within species, and genotypes within populations that are best adapted to the severe winters in these locations yet are also able to produce high biomass yields. We hypothesize that Miscanthus accessions collected from high latitude and/or high altitude locations with cold winters will have higher rates of overwintering in hardiness zone 5 environments such as central Illinois, and more alleles for overwintering, than accessions collected from southern locations with mild winters.
Genome-wide association (GWA) analysis is a powerful approach to determine the genetics of complex traits by exploiting linkage disequilibrium between a marker allele and the causative QTL allele (Lipka et al., 2015). Similarly, genomic prediction can be an effective method to facilitate marker-assisted selection but without the need to identify statistically significant associations between genomic regions and phenotypes. Given that Miscanthus is a long-lived perennial, marker-assisted selection is expected to enable plant breeders to efficiently select superior genotypes and greatly improve breeding efficiency. However, few GWA or genomic prediction studies on Miscanthus have been published and these have focused primarily on yield traits of M. sinensis (Nie et al., 2016;Slavov et al., 2014Slavov et al., ,2018 and M. sacchariflorus (Clark et al., 2016); to the best of our knowledge, none GWA studies have been published on overwintering ability, and only one genetic mapping study reported QTLs for overwintering ability (Dong, Liu, et al., 2019). To facilitate the use of a previously characterized M. sinensis germplasm panel (Clark et al., 2014 and identify superior genotypes for breeding, we evaluated the overwintering ability of M. sinensis genotypes in this germplasm panel in North America and Asia. The objectives of this study were to (a) quantify phenotypic and genetic variation for overwintering ability in M. sinensis; (b) identify potential associations between overwintering ability and environmental factors such as latitude and hardiness zone at origin; (c) identify molecular markers associated with overwintering ability, dissect their allelic effects, and explore the potential of marker-assisted selection using these detected markers; and (d) assess the potential of genomic prediction for overwintering ability in M. sinensis using tens of thousands of markers.

| Plant materials and experimental design
We studied the overwintering ability of 565 Miscanthus genotypes, including 561 M. sinensis, 3 closely related M. floridulus, and the M×g '1993-1780' control. Hereafter, we refer all 561 M. sinensis and 3 closely related M. floridulus as M. sinensis since Clark et al. (2014) found that the M. floridulus were part of the southeastern China M. sinensis group. These 564 genotypes were previously assigned to the following eight genetic groups (Clark et al., 2014: US ornamental cultivars (76), US naturalized populations (38), Southern Japan (28), Northern Japan (83) (25), Yangtze-Qinling (73), Southeastern China plus tropical (87). Six of the aforementioned genetic groups were previously determined based on discriminant analysis of principal components (DAPC; Jombart, Devillard, & Balloux, 2010) and by the software Structure (Falush, Stephens, & Pritchard, 2003); the US ornamental cultivars and US naturalized populations were found to be derived from southern Japan but were labeled independently to denote their provenance (Clark et al., 2014. Moreover, about half of the accessions in the US ornamental cultivars group have ≤30% ancestry from diploid M. sacchariflorus, presumably the result of past efforts by ornamental grass breeders in Germany to increase winter hardiness and obtain earlier flowering (Clark et al., 2014). Detailed information about these 564 M. sinensis accessions is listed in Data S1. All accessions were maintained as clonal stock plants in pots in a greenhouse at New Energy Farms in Ontario, Canada and vegetatively propagated. Ramets of each accession were distributed to each of five field trial locations during the winter of 2012.
In the spring of 2012, field trials were planted at two locations in East Asia (Table 1, Figure 1; HU = Hokkaido University, Sapporo, Japan; KNU = Kangwon National University, Chuncheon, South Korea) and three in North America (UI = University of Illinois at Urbana-Champaign, Urbana, IL, USA; CSU = Colorado State University, Fort Collins, CO, USA; NEF = New Energy Farms, Leamington, ON, Canada). Additional field trials were planted in 2013 at three locations (UI, CSU, and KNU; Table 1). Thus, a total of eight field trials were conducted. Nearly a full complement of the entire panel of 564 M. sinensis genotypes was planted in each of the trials at HU, KNU, and NEF. However, due to U.S. quarantine restrictions, only a subset of M. sinensis genotypes was planted in the field trials at UI andCSU (141 in 2012 trial, 163 in 2013 trial). Each field trial was arranged in a randomized complete block design with four replications except for KNU, where only one block included the entire panel and two additional small blocks included a same subset of 141 genotypes due to space constraints. Each plot contained one plant. Spacing between and within rows was 1.5 m. In each trial, irrigation was applied as needed during the first year to ensure good establishment and discontinued in subsequent years. The USDA hardiness zones of the field locations ranged from 5 to 7 (Table 1), representing a range of low temperature stresses during the winter (30-year average annual minimum temperature from −26.1°C to −12.2°C).

| Phenotypic data collection
Overwintering ability was calculated as: 0, plant was alive in previous year's autumn but was dead in current year's spring; 1, plant was alive in previous year's autumn and also regrew in current year's spring. For each trial location, data on air temperature (daily mean, maximum, and minimum) and precipitation were compiled from nearby weather stations ( Figure 1). For the IL location, daily records of soil temperature at 10 cm below the surface were also obtained. Overwintering ability of each plant in each of the eight field trials was evaluated for the first and second winters post-establishment. For the 2012 trials at all five locations, overwintering ability was evaluated after the 2012-2013 and 2013-2014 winters. For the 2013 trials at UI, CSU, and KNU, overwintering ability was evaluated for 2013-2014 winter only at UI and CSU, and for both the 2013-2014 and 2014-2015 winters at KNU. In summary, a total of 14 overwintering ability evaluations were performed for these eight field trials across five locations.

| Statistical analysis of phenotypic data
Analyses of variance (ANOVAs) were conducted in SAS 9.4 procedure MIXED (SAS Institute Inc., Cary, NC, USA) to determine how overwintering ability during the first and T A B L E 1 Field trials of a Miscanthus sinensis diversity panel and controls planted at five locations and evaluated for overwintering ability.
Field trials were established at each location in 2012, and three additional field trials were established at UI, CSU, and KNU in 2013 second winters was each affected by genetic group, genotype, location, year of establishment, number of growing seasons, and their interactions using the mixed linear model: where OWA represents overwintering ability, μ is the grand mean, DAPC represents the genetic groups determined in Clark et al. (2014Clark et al. ( , 2015, G(DAPC) is genotype nested in genetic group, L is location, Y is the year the trial was planted, S is number of growing season, B(LY) is block nested in location by year (i.e., in field trial), DAPC*L represents genetic group by location interaction, DAPC*Y represents genetic group by year interaction, DAPC*S represents genetic group by growing season interaction, G(DAPC)*L represents genotype nested within genetic group by location interaction, G(DAPC)*Y represents genotype nested within genetic group by year interaction, G(DAPC)*S represents genotype nested within genetic group by growing season interaction, and ε is a random error term. All model terms were set as fixed except for block. Traditionally, logistic regression is preferred over linear regression for dichotomous dependent variable modeling (i.e., overwintering ability in this study). However, it has been shown that, in most practical cases, linear regression can be as accurate as or even better than logistic regression, especially in high-dimensional data (Hellevik, 2009). To confirm the reliability of mixed linear model in analyzing binary overwintering ability data, we ran both logistic regression and a mixed linear model for significance tests. Results of the two tests were identical, and the correlation coefficient between the p-values of two testes was 0.99998 (Table S1). Therefore, this conformed to the conclusion of Hellevik (2009) that linear regression could be reliably used in analyzing dichotomous dependent variables. To estimate best linear unbiased predictors (BLUPs) of overwintering ability for each combination of winter, trial location, and year of establishment, ANOVAs were conducted using the following linear mixed model:

| BioClim analysis
In order to identify associations between overwintering ability and environmental factors, we first explored the relationship between overwintering ability and latitude and hardiness zone at origin.

| Marker development and missing data imputation
Restriction site-associated DNA sequencing (RAD-seq) was performed according to a previously described protocol (Clark et al., 2014). In brief, genomic DNA was digested with MspI and either PstI-HF or NsiI-HF (New England Biolabs). Digested DNA was then ligated to a barcoded adapter with a PstI/NsiI overhang and a universal adapter with an MspI overhang. Ninety-five barcoded samples were then pooled into one library, and 200-500 bp PCR products were selected on 2% agarose. Libraries were sequenced on an Illumina HiSeq 2500 at the Roy J. Carver Biotechnology Center at the University of Illinois. Nine PstI libraries from a previous study (Clark et al., 2014; data available at the NCBI Sequence Read Archive, accession SRP026347), as well as eight additional PstI libraries and 13 NsiI libraries (Clark et al., 2016; data available at NCBI, accession SRP063572), were included in the analysis. Because diploid M. sinensis has a large genome of ~5.3 Gb (Chae et al., 2014;Rayburn, Crawford, Rayburn, & Juvik, 2009), multiple sequencing runs can improve read depth and reduce missing data, which are especially important for correctly calling heterozygous loci in obligately outcrossing species such as this. Therefore, every individual in the study was included in at least two PstI libraries and two NsiI libraries. (1) (2) OWA ij = + G i + B j + ij F I G U R E 1 Environmental conditions during multiyear evaluations of a Miscanthus sinensis diversity panel at five field trial locations. The five locations included three in North America (University of Illinois at Urbana-Champaign = UI, Colorado State University = CSU, New Energy Farms in Leamington, Ontario = NEF) and two in Asia (Hokkaido University = HU and Kangwon National University = KNU). Planting date at each location was highlighted using dashed vertical blue lines (UI, CSU, and KNU had two trials planted in consecutive years). Red lines indicate daily average air temperature in °C, black lines represent daily precipitation in mm, gray-shaded areas indicate when plants were dormant, or field trial was not performed, and unshaded areas represent the growing season. For CSU trials, irrigation was applied as indicated by cyan solid bars Single-nucleotide polymorphisms (SNPs) were mined from all M. sinensis entries included in the field trial using the UNEAK pipeline (version 3.0 standalone; Lu et al., 2013) with an initial minimum call rate of 0.04 and a minimum minor allele frequency of 0.002. Three doubled haploid M. sinensis individuals (Głowacka, Kaczmarek, & Jeżowski, 2012;Swaminathan et al., 2012) were also included in the analysis to identify and remove any SNPs that appeared heterozygous and thus likely the result of paralogs. The dataset was then filtered to only include SNPs that had at least a 70% call rate in at least one of the genetic groups previously identified by discriminant analysis of principal components (Clark et al., 2014;Jombart et al., 2010) and a minor allele frequency of at least 0.05 in at least one of the nine groups, resulting in 70,327 SNPs retained, with an overall missing data rate of 39%. Imputation was performed with an estimation-maximization (EM) method based on relatedness (Poland, Brown, Sorrells, & Jannink, 2012) implemented in the R package rrBLUP (Endelman, 2011). To obtain genomic positions of SNPs for Manhattan plots and identify candidate genes, sequence tags were aligned to the Sorghum bicolor 3.0 reference genome (Paterson et al., 2009) using Bowtie2 (Langmead & Salzberg, 2012).

| Genome-wide association analysis
GWA analyses for overwintering ability were performed on the 564 M. sinensis genotypes (including the three M. floridulus) for each combination of growing season (i.e., winter) and field trial separately, for a total of 14 analyses. GWA analyses were performed using the multilocus mixed model (MLMM) implemented in R (Segura et al., 2012). MLMM adopts a forward-backward model selection procedure to potentially improve power (the ability to find true positives) and reduce false-positive rates. We allowed up to 10 forward selection and backward elimination steps for each GWA analysis. The optimal model was determined based on the extended Bayesian information criterion, defined as the largest model in which all cofactors have a p-value below the 5% Bonferroni-corrected threshold. To control for the cofounding effects of cryptic relatedness and population structure, we incorporated into the MLMM model the additive relationship matrix that was calculated using the EM imputation method in rrBLUP (Endelman, 2011;Poland et al., 2012) and the principal components of population structure calculated in R package adegenet (Jombart et al., 2010). Significant markers associated with overwintering ability were dissected for additive and dominance effects. In brief, additive effect was calculated as the phenotypic difference between two homozygous genotypic classes, and dominance effect was calculated as the phenotypic difference between the average of two homozygous genotypes and the heterozygous genotype. Then, for each M. sinensis accession, using all markers that were significantly associated with overwintering ability detected across 14 GWA analyses, the genotypic value of overwintering ability was calculated based on genotypic effects (additive and dominance), and breeding value was calculated based on only additive effects.

| Genomic prediction
Genomic prediction for first-winter overwintering ability was conducted in the R package rrBLUP (Endelman, 2011) based on the 70,327 SNPs. Predictive ability was assessed by a fivefold cross-validation strategy. Briefly, M. sinensis accessions were randomly split into five sets, of which four sets were used as training population to estimate marker effects, and the remaining set was used as validation population to calculate the genomic estimated breeding values (GEBVs). This process was repeated for each of the five sets. Thus, GEBVs were obtained for all M. sinensis accessions. Predictive ability (r) was defined as the Pearson's correlation between GEBVs and trait BLUPs (Daetwyler, Calus, Pong-Wong, de Los Campos, & Hickey, 2013). This fivefold cross-validation process was iterated 50 times. Mean and standard deviation of these 50 correlation coefficients were reported as predictive ability and its standard deviation.
As the presence of closely related individuals in both training and validation populations could artificially inflate prediction accuracies, it is recommended to control for population structure in designing the cross-validation scheme (Ly et al., 2013;Riedelsheimer et al., 2013;Spindel et al., 2015). To do so, three statistical methods were investigated Overview photos of UI, NEF, HU, KNU fields. Note the severe overwintering losses in the UI field trial in genomic prediction for first-winter overwintering ability across the whole panel (N = 564). These three statistical methods differed in controlling the previously determined genetic structure in this M. sinensis diversity panel. The M. sinensis accessions were previously assigned to one of the following eight groups: Southern Japan (28), Northern Japan (83), Korea/N China (154), Yangtze-Qinling (73), Sichuan Basin (25), Southeastern China plus tropical (87), US ornamental cultivars (76), and US naturalized populations (38). For Method 1, as a control comparison, population structure was not accounted for, and best linear unbiased predictor (BLUP) values were calculated using the following equation: where y is the first-winter OWA raw data, μ is grand mean, G is genotype, L is location, Y is year of field trial establishment, B(LY) is block nested within location by year (i.e., in field trial), G*L is genotype by location interaction, G*Y is genotype by year interaction, and ɛ is random error. All model terms were set as random. BLUPs of the G term were used in genomic prediction.
For Method 2, the BLUPs obtained from Equation 3 were used as the independent variables, and then, we fitted a mixed model by including DAPC group as fixed effect: where y is the genotype BLUPs calculated from Equation 3, μ is grand mean, DAPC represents the previously determined where OWA represents overwintering ability, μ is the grand mean, DAPC is genetic group, G(DAPC) is genotype nested in genetic group, L is location, Y is year, S is growing season, B(LY) is block nested in location by year (field trial), DAPC * L represents genetic group by location interaction, DAPC * Y represents genetic group by year interaction, DAPC * S represents genetic group by growing season interaction, G(DAPC) * L represents genotype nested within genetic group by location interaction, G(DAPC) * Y represents genotype nested within genetic group by year interaction, G(DAPC) * S represents genotype nested within genetic group by growing season interaction, and ε is random error. All model terms were set as fixed except for block.  Clark et al., 2014Clark et al., , 2015. The eight genetic groups included US ornamental cultivars (1), US naturalized populations (2), Southern Japan (3), Northern Japan (4), Korea/N China (5), Sichuan Basin (6), Yangtze-Qinling (7), Southeastern China plus tropical (8). On the X-axis, the number of genotypes with data for each DAPC group in each location is indicated in parentheses after the group number. The Y-axis represents overwintering ability ranging from 0 to 1, where 0 indicates plants failed to survive the winter and 1 means plants survived the winter. Overwintering ability of M. ×giganteus '1993-1780' is shown as horizontal dashed line (only shown where <100%). Boxes span from the first to third quartile for each group. Whiskers extend to the minimum and maximum values or to the first and third quartile ±1.5 times the box length, respectively, whichever is shorter. Points indicate genotypes with values outside the range spanned by the whiskers. Internal bar shows the median value genetic groups in this M. sinensis diversity panel (Clark et al., 2014. The residuals from Equation 4 were used as phenotypic data in genomic prediction. For Method 3, we fitted a model by nesting genotype within DAPC group using the following equation: where model terms were defined as in Equations 3 and 4, except that G(DAPC) represents genotype nested within DAPC group. All model terms were set as random effect except for DAPC, which was set as fixed effect. BLUPs of the G(DAPC) term were used in genomic prediction.
Additionally, we explored the effectiveness of genomic prediction for first-winter OWA within each genetic group. Within-genetic-group genomic prediction was performed for each of the following M. sinensis genetic groups: Southern Japan (N = 142), Northern Japan (N = 83), Korea/N China (N = 154), Yangtze-Qinling (N = 73), and Southeastern China plus tropical (N = 87). Within each of these five genetic groups, BLUPs were calculated using Equation 3, and genotype BLUPs were used in genomic prediction. The Sichuan Basin group (N = 25) was not included in the within-genetic-group genomic prediction analysis due to its small size. For the within-genetic-group genomic prediction analysis, the Southern Japan group studied was a composite of its subset groups, US naturalized population (N = 38), US ornamental cultivars (N = 76), and Southern Japan group (N = 28; Clark et al., 2014Clark et al., ,2015; DAPC analysis had included these genotypes in a single group and keeping them together allowed for a larger population size, which is typically advantageous for genomic prediction.

| RESULTS
In the M. sinensis diversity panel, highly significant differences in overwintering ability were observed among genetic groups, genotypes within genetic groups, locations, year of establishment, number of growing seasons, and their interactions except for genetic group by year interaction (Table 2, Figures 2, 3, 4a). The control genotype M×g '1993-1780' survived in all field trials except for the first winter of the CSU 2012 trial (Figure 3). Among the six M. sinensis genetic groups from Asia, first-winter overwintering ability, averaged over all of the field trial locations, was highest for the Korea/N China group (0.98; Figure 4a, Table 3), and relatively high for the Northern Japan group (0.96), the Southern Japan group (0.93), and the Yangtze-Qinling group (0.88). The Korea/N China group showed consistently superior overwintering ability at each of the field trial locations (0.86-0.99). The lowest first-winter overwintering ability was observed for the Southeastern China plus tropical group (0.36), followed by the Sichuan Basin group (0.40), which was consistent with their natural adaptation to low-latitude environments that have mild winters (Figure 4a, Table 3). Such differences in overwintering ability among different genetic groups were further delineated by a strong association between first-winter overwintering ability and latitude of origin (R 2 = 0.49, Figure 4a). Moreover, a similar association was also observed between first-winter overwintering ability and hardiness zones from which M. sinensis accessions were collected (R 2 = 0.34, Figure 4d), indicating M. sinensis accessions collected from high latitude locations with cold winters had higher rates of overwintering than accessions from southern locations with mild winters. In addition, BioClim analysis indicated that BIO1 (Annual mean temperature), BIO6 (Minimum temperature of coldest month), BIO10 (Mean temperature of warmest quarter), and BIO11 (Mean temperature of coldest quarter) were significant at α = 0.05 (Table S2). BioClim PC1 explained 92.4% of total variance among these four significant variables. Intriguingly, a similarly strong association between first-winter overwintering ability and BioClim PC1 was observed (R 2 = 0.50; Table S2). Nevertheless, some exceptional individuals from the Southeastern China plus tropical group and the Sichuan Basin group survived first and second winters at HU, KNU, and/or NEF (Figure 3, Table 3). Among the five field trial locations, first winter overwintering ability was typically lowest at the two US locations, UI (0.40-0.49, two trials) and CSU (0.22-0.76, two trials), though these field trial locations primarily tested entries from the US ornamental and US naturalized groups, and included few of the entries from Asia due to quarantine regulations that delayed their importation beyond the timeframe of this study (Table 3). In contrast to the US locations, first-winter overwintering was typically higher at HU (0.99), KNU (0.69-0.82, two trials), and NEF (0.73), even when the comparison was limited to just the US ornamental and US naturalized groups (Table 3). The lower overwintering rates observed at the US field trial locations relative to the other sites was consistent with a combination of lower minimum air temperatures and/or less snow cover for the US trials ( Figure 1, Table 3). The winter of 2013-2014 was especially cold at the US sites, with minimum air temperatures of −25.3°C at UI and −24.5°C at CSU, whereas the NEF site was moderated by Lake Erie and reached only −19°C, and the two locations in Asia were warmer still at about −9°C (Figure 1, Table 3).
At a given location, overwintering ability (i.e., the proportion of plants that were alive the previous growing season and also survived the subsequent winter) was typically lower after the first winter than after the second winter (Table  3), because the least hardy plants typically died during the first winter, leaving the more adapted genotypes and better established plants to be challenged by the second winter. For example, at NEF in Ontario, overwintering ability was 0.99 after the second winter (Table 3; min. air temp.: −19.0°C) but was only 0.73 after the first winter (min. air temp.: −15.0°C), indicating nearly all the plants that survived the first winter also survived the second colder winter. However, a notable exception to greater overwintering after the second winter was observed for the 2012 trial at HU in Japan, in which nearly all plants from each of the eight genetic groups survived the first winter (Table 3; overwintering ability: 0.99; min. air temp.: −8.4°C), whereas after the second winter, the overall overwintering ability was only 0.80 (min. air temp.: −9.2°C). In particular, two genetic groups, including Sichuan Basin and Southeastern China plus tropical, had substantial losses (0.60-0.61 OWA) during the second winter at HU. This atypical outcome at HU was likely the result of high and early snow cover during the first winter, which insulated the plants from low air temperatures (Figure 1).
Genome-wide association analyses across the five field trial locations detected a total of 73 significant markertrait associations for overwintering ability ( Figure S1), of which 42 could be aligned to physical positions on S. bicolor genome version 3.0 (Table 4). For the first-winter overwintering ability, additive effects among these 73 SNPs ranged from 0.01 to 0.42, and dominance effects ranged from −0.59 to 0.44 (Data S2). Among the 564 M. sinensis accessions, the frequency of desirable alleles (up to two desirable alleles per locus per individual) over all 73 overwintering ability loci ranged from 0.51 (accession PMS-044, Sichuan Basin genetic group) to 0.88 (accession JM0047.002, Northern Japan genetic group) (Data S2, Figure 4b). The predicted genotypic values for overwintering ability based on the 73 SNPs accounted for 56% of the variation observed for first-winter survival among M. sinensis accessions averaged over the five field trial locations (Figure 4c). Similarly, the breeding values for overwintering ability explained 55% of the variation for observed first-winter overwintering ability ( Figure  S2). Moreover, a strong association was observed between frequency of desirable alleles and latitude (R 2 = 0.52, Figure  4b), and a moderate association was observed between frequency of desirable alleles and hardiness zones (R 2 = 0.27, Figure 4e), indicating that M. sinensis accessions collected from high latitude locations with cold winters had more alleles for overwintering than accessions from southern locations with mild winters. Notably, observed first-winter overwintering ability was also significantly correlated with the frequency of desirable alleles (R 2 = 0.35, Figure S3).
Genomic prediction analyses for first-winter overwintering ability resulted in high and stable estimates, as indicated by the small standard error and range of genomic prediction abilities (Table 5). For Method 1, which did not control for population structure, genomic prediction accuracy for the entire panel averaged 0.81. In contrast, for Methods 2 and 3, which controlled population structure (i.e., DAPC groups), average genomic predictive ability was 0.73 and 0.71, respectively. However, genomic prediction within DAPC group varied greatly by group (Table 5). The Southeastern China plus tropical group, and the Southern Japan group had moderate genomic predictive abilities, with average values of 0.34 and 0.48, respectively. For the Yangtze-Qinling group, genomic predictive ability was low, with an average of 0.04. For the Northern Japan and Korea/N China groups, genomic predictions could not be obtained because nearly all accessions in these two groups survived the first winter ( Figure 4a, Table 3), resulting in zero variance among genotypes for this trait.

M. sinensis for overwintering ability: Implications for breeding
Like most prior studies of overwintering ability in Miscanthus  Year 1 OWA indicates the proportion of genotypes that were alive in the autumn of the first year that regrew in the spring of the second year (i.e., survived the first winter). Year 2 OWA indicates the proportion of genotypes that were alive in the autumn of the second year that regrew in the spring of the third year (i.e., survived the second winter). T A B L E 4 (Continued) panel (N = 564) that represented most of the species' natural geographic range, and evaluated the results in the context of previously ascertained population structure (Clark et al., 2014. Consistent with prior studies of overwintering ability in Miscanthus (Clifton-Brown et al., 2001), we found that M. sinensis plants were typically more at risk of dying during the first winter after planting than during the second winter. However, the one exception we observed to this progression was associated with high snowfall early season at HU in year 1 that resulted in nearly all plants surviving the first winter (0.99; Table 3), but a lower overwintering rate was observed in year 2 (0.80), suggesting that avoidance of killing temperatures (in this case by insulating snow cover) can be an important contributor to overwintering in Miscanthus. First-year plants of Miscanthus, which typically produce much less aboveground biomass than second-year plants (Dong, Green, et al., 2019;Gifford, Chae, Swaminathan, Moose, & Juvik, 2015), likely also produce less belowground biomass and more shallow rhizomes than more mature plants, and thus may have few or no buds sufficiently deep in the soil to avoid killing temperatures. Especially as Miscanthus plants mature, those genotypes that are able to produce deep rhizomes will be able to avoid temperatures sufficiently cold to damage or kill dormant buds. Great variation among and within M. sinensis genetic groups was observed for overwintering ability at five field trial locations, with 73 significantly associated SNPs identified via GWA analyses. Genetic groups and genotypes that originated from temperate environments at relatively high latitudes typically had greater overwintering ability than those from subtropical and tropical environments, as expected (Table 3, Figures 2, 3, 4a). We observed positive association between latitude of origin and overwintering ability (R 2 = 0.49, Figure  4a). Similarly, Yan et al. (2012) observed that first-year overwintering ability of 31 M. sinensis accessions collected in China was positively associated with latitude of origin when grown at two northern field trial sites, and the correlation was greatest at the most northern site with the coldest winter (Xilinhot). In addition, BioClim analysis revealed that BIO1 (annual mean temperature), BIO6 (minimum temperature of coldest month), BIO10 (mean temperature of warmest quarter), and BIO11 (mean temperature of coldest quarter) had significant effects on first-winter overwintering ability (Table  S2), indicating that adaptation to low temperatures in winters as well as to a certain amount of heat units in the growing season (as determined by location of origin) both influenced actual overwintering at the trial sites. For example, if a plant is adapted to a long growing season but is then grown in a place with a short growing season, we would expect it to be less likely to survive the winter because it will not be at the right physiological state (e.g., insufficient storage of carbohydrates , which fitted DAPC group as a fixed effect and the derived residuals were used in genomic prediction; Method 3 was based on Equation 5, which controlled population structure by nesting genotype within DAPC group. For within-DAPC-group analysis, the Southern Japan group (N = 142) included genotypes directly collected in southern Japan (N = 28), US ornamentals (N = 76), and US naturalized populations (N = 38), given that the latter two were derived from southern Japan (Clark et al., 2014 and larger sample sizes are advantageous for genomic prediction. b Genomic prediction statistics were calculated based on 50 fivefold cross-validations for each analysis. c Detailed ANOVA results provided in Data S3. d Genomic prediction within Northern Japan and Korea/N China groups were not estimable because nearly all accessions survived (Figure 4a, Table 3), which resulted in zero-value estimates for the percentage of variance explained by genotype for overwintering ability.  (Figure 2a). This event highlighted the importance of selecting for the coldest winter expected in a target environment, especially given that Miscanthus is a long-lived perennial crop and establishment of new commercial plantings is expensive. Thus, one strategy to breed M. sinensis for temperate environments is to select concurrently for overwintering ability and high biomass yield within the genetic groups that naturally have the greatest adaptation to cold winters (e.g., Korea/N China, Northern Japan group, Southern Japan group, and the Yangtze-Qinling group). Additionally, the rare relatively hardy individuals identified within the subtropicaland tropical-adapted Sichuan Basin group (e.g., PMS-005; Data S1) and Southeastern China plus tropical group (e.g., PMS-008; Data S1) represent a potentially valuable breeding opportunity because M. sinensis accessions from lower latitudes typically have higher biomass yield potential than those from higher latitudes . Moreover, the presence of advantageous overwintering alleles even within the subtropical-and tropical-adapted M. sinensis genetic groups, albeit at lower frequencies than in the temperateadapted groups (Figure 4b, Data S2), should facilitate rapid genetic gains via marker-assisted selection. The opportunity to conduct marker-assisted selection for advantageous overwintering alleles at a large number of loci within any of the M. sinensis genetic groups is expected to greatly improve breeding efficiency, especially if the breeder must otherwise rely on unusually cold winters to screen populations.

| Introgression of genes for overwintering ability from northern-adapted M. sacchariflorus into M. sinensis is another potentially useful breeding strategy
Introgression of winter hardiness alleles from cold temperate-adapted M. sacchariflorus may be another viable strategy for improving overwintering ability in M. sinensis. We have shown previously that about half of the ornamental cultivars sold as M. sinensis in the US and Europe are in fact hybrids between M. sinensis and diploid M. sacchariflorus that have been backcrossed to M. sinensis one or more times. Notably, the M. sacchariflorus ancestry in these hybrids is from northern China or eastern Russia (Clark et al., 2014, which are cold temperate environments (hardiness zones 3-5). In this study, we included 76 accessions of predominantly M. sinensis ancestry from the US ornamental cultivars genetic group, and M. sacchariflorus ancestry in these cultivars ranged from 0.00 to 0.36 (Data S1). A significant association was observed between first-winter overwintering ability of the ornamental cultivars and their proportion M. sacchariflorus ancestry (R 2 = 0.17, p = 0.0003, Figure S4). Yan et al. (2012) observed that among 48 accessions of M. sacchariflorus and 31 of M. sinensis from China grown at two northern field trial sites in China, the former was more winter-hardy than the latter. Similarly, in eastern Russia, Clark et al. (2016) did not find wild populations of M. sinensis in areas colder than hardiness zone 5, but did find abundant populations of M. sacchariflorus through hardiness zone 3. Thus, there appears to be a good case for using cold tolerance genes from M. sacchariflorus to improve M. sinensis.
However, we should not assume that all M. sacchariflorus genotypes are equally good sources of winter hardiness. Previously, we have identified three tetraploid and three diploid M. sacchariflorus genetic groups, each associated with different geographical regions in eastern Asia, ranging from the mild Yangtze River region to the cold winters in eastern Russia along the Amur River . For example, Clifton-Brown et al. (2001) observed that first-winter (1997-1998) survival for a tetraploid M. sacchariflorus (EMI no. 5, collected along the Nagara River in Gifu Prefecture, Japan) was only 50% in Sweden and 33% in Denmark, whereas an M. sinensis from Hokkaido (Northern Japan group) had 95%-99% survival at the same locations. Similarly, Clifton-Brown et al. (2001) found that nearly all of the triploid M×g tested in Sweden and Denmark died after the first winter, but survival of four M. sinensis from Honshu Japan had survival rates of 84%-99%.
If the ultimate goal of breeding M. sinensis is to use it as a parent for making improved triploid M×g by crossing it with tetraploid M. sacchariflorus, then it would also be desirable to have an understanding of potential heterotic groups within M. sacchariflorus, and how introgressed genes from a diploid M. sacchariflorus accession into M. sinensis might interact with a tetraploid M. sacchariflorus from a different genetic group. However, there is currently little information on heterotic groups and gene interactions in Miscanthus to guide such breeding strategies.

| Trait association comparisons between this GWA and a parallel study of three interconnected biparental populations
In our parallel study of three interconnected diploid Miscanthus populations (Dong, Liu, et al., 2019), we identified nine QTL for overwintering ability via joint population analysis, whereas in the current GWA analyses, we identified 73 significant SNPs for overwintering ability ( Figure S1, Table 4), an eightfold difference. Greater sampling of genetic diversity in the M. sinensis germplasm panel (564 genotypes from throughout the species' geographic range) in comparison to the four parents of the interconnected populations, likely contributed to the greater number of significant associations identified. Additionally, the five field trial locations in years 1 and 2 for phenotyping in the GWA study was likely advantageous over the one location in year 3 for phenotyping in the interconnected populations study. Moreover, in linkage analysis of biparental populations, there are only a few opportunities for recombination to occur within each population, and QTLs are only detected based on recombination events that occurred during population development, resulting in relatively low mapping resolution (typically 10-20 cM, Doerge, 2002;Holland, 2007). In contrast, historical recombination and greater genetic diversity can be exploited for high resolution mapping in association analysis (Lipka et al., 2015;Yu & Buckler, 2006), and this appears to have been the case for our studies of overwintering ability in Miscanthus.
Two QTLs identified in the interconnected biparental population study included a total of seven SNP-trait associations from the GWA study ( Figure 5). One QTL on

Miscanthus
LG 8 encompassed five marker-trait associations from GWA that aligned to S. bicolor chromosome 4 (Figure 5a), and the other QTL on Miscanthus LG 11 encompassed two marker-trait associations from GWA that aligned to S. bicolor chromosome 6 ( Figure 5b). The identification of significant SNPs via GWA within independently identified QTL from three interconnected biparental populations lends strong support to the conclusion that these regions of the genome are important for overwintering ability in Miscanthus.

| Marker-assisted selection with 73 significant markers from GWA and genomic prediction for overwintering ability
Using the 73 markers detected in GWA analyses, both the estimated genotypic values (considering both additive and dominance effects) and the estimated breeding values (considering only additive effect) explained 55%-56% of the variation for observed first-winter survival in this M. sinensis germplasm panel (Figure 4c, Figure S2). However, the predictions for overwintering ability based on these 73 markers F I G U R E 5 Correspondence between QTLs for overwintering ability detected across three interconnected diploid F 1 Miscanthus populations (Dong, Liu, et al., 2019) and genome-wide associations (GWA) for overwintering ability in an M. sinensis diversity panel. Black dashed bars represent Miscanthus linkage groups (LGs), and blue bars represent S. bicolor (version 3.0) LGs. Orange lines represent corresponding genomic regions between Miscanthus and Sorghum bicolor. (a) An overwintering ability QTL on Miscanthus LG 8 identified via genetic mapping in three interconnected biparental populations (purple bar) encompassed five marker-trait associations from the GWA analysis (black text). Text inside the parenthesis represents gene symbols. The cold tolerance gene COR47 (red text; previously identified in Arabidopsis by Puhakainen et al. (2004) and Bozovic et al. (2013)) was located within this QTL and was 0.02 Mb away from one GWA hit (UIMiscanthus016552). (b) A QTL for overwintering ability on Miscanthus LG 11 identified via genetic mapping in three diploid interconnected biparental populations (purple bar) corresponded to two marker-trait associations from the GWA analysis (black text). Three additional candidate genes (red text) including carboxylesterase 13 (CEX13), WRKY transcription factor (WRKY2), and cold shock domain protein 1 (CSDP1) were also located inside this QTL could be overestimated because the genotypic and breeding values were calculated from the same dataset that was used to identify the significant SNPs. This phenomenon of inflated prediction is known in as 'inside trading' in genomic selection (Arruda et al., 2016). Moreover, given that the dominance effects across these 73 markers ranged from −0.59 to 0.44 (Data S2), the similar prediction accuracy between estimated genotypic values and estimated breeding values suggested that the dominance effects of alleles at multiple loci might have cancelled out each other.
Genomic prediction could be more efficient than markerassisted selection in breeding cold-hardy Miscanthus. GWA analyses cannot detect with statistical significance all causative loci underlying a trait of interest. Moreover, in both linkage mapping and GWA analyses, the separation of marker-trait association detection from marker effect estimation results in biased effect estimation (Beavis, 1994;Jannink, Lorenz, & Iwata, 2010;Lande & Thompson, 1990). In contrast, genomic prediction uses all markers for modeling the performance of individuals. Nevertheless, GWA enables identification of candidate genes that can help elucidate gene networks. Thus, GWA and genomic prediction are complementary strategies. Even though the equal variance assumption (i.e., all markers contribute equally to the genetic variances of trait) in ridge regression is unrealistic in a breeding context, superior predictive ability of genomic selection relative to marker-assisted selection has been demonstrated (Arruda et al., 2016;Heffner, Jannink, & Sorrells, 2011;Jannink et al., 2010;Meuwissen, Hayes, & Goddard, 2001).
Population structure affected genomic prediction of overwintering ability for M. sinensis. For whole panel analysis, Method 1, which did not account for population structure, resulted in prediction accuracies that were biased upward relative to Methods 2 (by 8 points) and 3 (by 10 points) which did account for population structure; the observed bias was consistent with prior studies on other crops (Fiedler et al., 2018;Riedelsheimer et al., 2012;Spindel et al., 2015). To effectively differentiate the overwintering potential of individuals within the Northern Japan and Korea/N China groups, from which nearly all genotypes survived the first winter, higher resolution phenotyping schemes such as a 1-10 ordinal system to capture relative vigor, rather than the current binary system that captured only survival data could be helpful; additionally, challenging these genotypes in colder winter environments would be another promising strategy.

| Candidate genes
GWA analyses provided information for dissecting the genetic mechanism of overwintering ability in Miscanthus. Of the 42 marker-trait associations that aligned to the S. bicolor genome, many were located near known cold-responsive genes (Table 4). Moreover, some of the candidate genes associated with significant SNPs were found within QTL identified in our parallel study of three interconnected biparental populations ( Figure 5). For example, within the overwintering ability QTL on LG 8, one of the significant SNPs from GWA (UIMiscanthus005328, aligned to S. bicolor chromosome 4) was within a MYB gene (MYB43, Table 4, Figure  5a), which is a transcription factor that has been shown to regulate stress response in Arabidopsis (Barah et al., 2016). The known cold tolerance gene, COR47, was also within the QTL on LG 8 and 0.02 Mb downstream from one GWA hit (UIMiscanthus016552; Figure 5a, Table 4). COR47 protects thylakoid membranes during freezing by encoding dehydrin proteins, thereby enhancing cold tolerance in Arabidopsis (Bozovic, Svensson, Schmitt, & Kohn, 2013;Puhakainen et al., 2004). Upregulation of COR genes has also been reported in perennial ryegrass (Lolium perenne L. cv. 'Caddyshack') during cold acclimation (Zhang, Fei, Warnke, Li, & Hannapel, 2009). Within the overwintering ability QTL on LG 11, three candidate cold-tolerance genes were identified, including carboxylesterase 13 (CEX13), WRKY transcription factor (WRKY2), and cold shock domain protein 1 (CSDP1). Two of these candidate cold-tolerance genes (CEX13 and WRKY2) were closely linked to the significant SNPs from GWA (CEX13 was 0.53 Mb from UIMiscanthus014960; WRKY2 was 0.88 Mb from UIMiscanthus118375; Figure 5b, Table  4). Carboxylesterase is an enzyme that regulates biological activities through hydrolysis (Gershater & Edwards, 2007). Although the functional details of carboxylesterase are still unknown, several members of this gene family have been detected following cold treatment in Arabidopsis (Wagstaff et al., 2010) and in grape (Vitis vinifera) (Xin et al., 2013). WRKY is a large transcription factor family that plays a broad spectrum regulatory role in plant defense regulation, response to abiotic stresses, growth, and development (Agarwal, Reddy, & Chikara, 2011;Seki et al., 2002). WRKY2 has been shown to have elevated expression in response to osmotic stress in Arabidopsis (Jiang & Yu, 2009). In plants, cold shock domain proteins play essential roles in acquiring freezing tolerance (Sasaki & Imai, 2012); they are among the most evolutionarily conserved nucleic acid-binding domains, predating the divergence of prokaryotes and eukaryotes (Graumann & Marahiel, 1998;Karlson & Imai, 2003).
One theme to emerge from this analysis of candidate genes was the identification of multiple members of gene families known to be involved with plant response to cold stress in or near newly identified SNPs throughout the Miscanthus genome that were associated with overwintering ability, which is consistent with the importance of gene duplication in plant evolution (Hu et al., 2003;Ming et al., 2002;Paterson, Schertz, Lin, Liu, & Chang, 1995). For example, in addition to the cold shock domain locus on S. bicolor chromosome 6 that was within the QTL on Miscanthus LG 11 and near three significant SNPs for overwintering ability (Figure 5b), another cold shock domain locus on a S. bicolor chromosome 3 was 0.01 Mb from marker-trait association UIMiscanthus118370 (Table  4). Similarly, in addition to the WRKY2 transcription factor on Miscanthus LG 11/ S. bicolor chromosome 6, another WRKY locus (WRKY54) was located downstream of markertrait association UIMiscanthus097991, which aligned to S. bicolor chromosome 8 (Table 4). Similar to COR47 within the QTL on Miscanthus LG 8/S. bicolor chromosome 6, another locus associated with protecting thylakoid membranes was identified, COR314, which was located 0.07 Mb from UIMiscanthus012595 and aligned to S. bicolor chromosome 9 (Grundy, Stoker, & Carré, 2015;Li et al., 2018). Additionally, three marker-trait associations (UIMiscanthus118368 aligned to S. bicolor Chromosome 2, UIMiscanthus017647 aligned to S. bicolor Chromosome 3, and UIMiscanthus005328 aligned to S. bicolor Chromosome 4) were all located within or near MYB transcription factor genes (Table 4). However, some likely candidate genes were only observed once. For example, the marker-trait association, UIMiscanthus118366 aligned to S. bicolor chromosome 2, was only 177 bp downstream from the JAZ8 gene (jasmonate-zim-domain protein 8), which has increased expression under cold stress in wheat (Triticum aestivum L.); moreover, this gene is part of a gene family that regulates responses to biotic and abiotic stresses typically by interacting with MYB transcription factors (Wang et al., 2017), which we have also found to play a role in the overwintering ability of Miscanthus. Thus, the candidate genes identified here represent testable hypotheses about the genes underlying overwintering ability of Miscanthus. Moreover, the functions of the newly identified candidate genes suggest that tolerance to freeze stress is an important component overwintering ability in Miscanthus, in addition to likely avoidance strategies.
In summary, screening a large panel of M. sinensis germplasm enabled us to determine that the Korea/N China genetic group would likely be a valuable gene pool for cold tolerance, and cold-tolerant genotypes are also frequent in the Yangtze-Qinling, Southern Japan, and Northern Japan groups. We detected 73 marker-trait associations for overwintering ability using the large M. sinensis germplasm panel and observed that Miscanthus accessions collected from high latitude locations with cold winters had higher rates of overwintering, and more alleles for overwintering, than accessions collected from low latitude locations with mild winters. Consistency between QTLs and GWA hits suggested that these genomic regions are important for response to cold stress. Many of the candidate genes underling peak marker-trait associations provided interesting hypotheses for further testing. Similarly, genomic prediction abilities for overwintering ability were high. Given these results, we expect both GWA and genomic prediction to substantially improve breeding efficiency for winter hardiness in M. sinensis. Thus, this study represents a significant step toward the development of new Miscanthus cultivars that are optimally adapted to temperate regions with cold winters.