Genetic loci associated with winter survivorship in diverse lowland switchgrass populations

High winter mortality limits biomass yield of lowland switchgrass (Panicum virgatum L.) planted in the northern latitudes of North America. Breeding of cold tolerant switchgrass cultivars requires many years due to its perennial growth habit and the unpredictable winter selection pressure that is required to identify winter‐hardy individuals. Identification of causal genetic variants for winter survivorship would accelerate the improvement of switchgrass biomass production. The objective of this study was to identify allelic variation associated with winter survivorship in lowland switchgrass populations using bulk segregant analysis (BSA). Twenty‐nine lowland switchgrass populations were evaluated for winter survival at two locations in southern Wisconsin and 21 populations with differential winter survivorship were used for BSA. A maximum of 10% of the individuals (8–20) were bulked to create survivor and nonsurvivor DNA pools from each population and location. The DNA pools were evaluated using exome capture sequencing, and allele frequencies were used to conduct statistical tests. The BSA tests revealed nine quatitative trait loci (QTL) from tetraploid populations and seven QTL from octoploid populations. Many QTL were population‐specific, but some were identified in multiple populations that originated across a broad geographic landscape. Four QTL (at positions 88 Mb on chromosome 2N, 115 Mb on chromosome 5K, and 1 and 100 Mb on chromosome 9N) were potentially the most useful QTL. Markers associated with winter survivorship in this study can be used to accelerate breeding cycles of lowland switchgrass populations and should lead to improvements in adaptation within USDA hardiness zones 4 and 5.


INTRODUCTION
Switchgrass (Panicum virgatum L.) is a native North American perennial grass with potential as a biofuel feedstock due to its high biomass potential on marginal lands with limited inputs. There are two ecotypes of switchgrass, lowland and upland, both of which are relevant to the development of biofuel feedstocks, largely dependent on geographic region (Porter, 1966;Zalapa et al., 2011). The lowland ecotypes are adapted to more southerly climes, whereas upland ecotypes are found in more northern latitudes (Casler, 2005(Casler, , 2012Casler et al., 2007;Casler et al., 2004;Lovell et al., 2021). Lowland and upland ecotypes are phenotypically distinct as exhibited by differences in overall height, leaf blade dimensions, tiller number, and heading/flowering dates (Casler, 2005).
The lowland ecotype of switchgrass has generated considerable interest, largely due to its extremely late flowering characteristics (Grabowski et al., 2017;Tornqvist et al., 2018). In USDA hardiness zones 3-5, lowland switchgrass flowers 4-6 wk later than upland switchgrass, increasing the duration of the biomass accumulation phase and increasing biomass yields by up to 50% (Casler et al., 2004). Multiple-location field evaluations of highly selected late-flowering populations confirmed a 50% increase in biomass yield, with the added benefit of broadening the adaptation range of lowland switchgrass from hardiness zones 6-9 to include hardiness zones 4 and 5 . However, lowland germplasm derived from prairie and savanna remnants in the southern United States (USDA hardiness zones 6-9) have low winter survivorship at the northern latitudes (Casler, 2012;Casler & Vogel, 2014;Casler et al., 2004). Lowland cultivars planted in northern Wisconsin (USDA hardiness zone 3b) had up to fourfold yield and winter survival loss compared with Oklahoma (USDA hardiness zone 7a) (Casler et al., 2004). Lovell et al. (2021) reported 42% winterkill in lowland plants compared with only 2% in upland plants when planted in northernmost sites in the United States (SD, NE, IL, and MI). These authors also reported a strong correlation between extreme minimum temperature of the site of origin with the biomass yield of lowland populations. As yield and winter survival are highly correlated (Bhandari et al., 2011;Lovell et al., 2021), the sustainable and economic use of switchgrass as an energy grass in northern latitudes can be improved by increasing winter survivorship of lowland ecotypes.
Previous experiments have documented that most mortality of switchgrass occurs when the plants are in their winter dormant phase (Sarath et al., 2014). Cold tolerance, generally measured as survivorship following stressful winter conditions (cold temperatures, high winds, frequent freezethaw cycles, and open winters with little or inconsistent snow cover), is a heritable trait (Hopkins et al., 1995) and can be improved by simple phenotypic selection methods Poudel et al., 2020). However, this process requires many years to ensure adequate levels of stress and to accumulate sufficient changes in allele frequency to warrant the release of a new cultivar. Identification of DNA markers associated with winter hardiness in lowland switchgrass would be extremely valuable to support marker-based selection programs and accelerate the assessment and selection of winter-hardy lowland genotypes.
Many traits of agronomic importance in crop plants are quantitatively inherited, meaning that they are controlled by many quantitative trait loci (QTL) with small individual effects (Falconer & Mackay, 1996). Identification of these QTL would provide the basis for marker-based selection approaches that reduce selection cycle time through its effectiveness on early-stage selection and allow selection in the absence of selective forces. Mapping QTL in experimental populations (Morton, 1955) and genome-wide association studies (GWAS) across diverse populations (Klein et al., 2005;Price et al., 2010;Zhu et al., 2008) are the most common approaches used for the identification of trait-associated variants (TAV). Conventional QTL mapping is based on newly generated recombination events which are often populationspecific (Feltus et al., 2006). On the other hand, GWAS exploit historical recombination within a large diversity panel and provide higher resolution to the mapping of the quantitative traits (Price et al., 2010;Zhu et al., 2008). Both of these approaches require genotyping of a large number of individuals, which is labor-intensive, time-consuming, and costly. Researchers working on cold tolerance in switchgrass are still in the initial stages of identification of TAV. To date, six QTL have been identified for freezing tolerance of switchgrass; however, these are likely specific to one particular lowland × upland hybrid population (Poudel et al., 2019b).
Bulk segregant analysis (BSA) is an alternative approach to identify TAV, which involves genotyping pooled DNA from individuals with shared phenotype that significantly reduce genotyping cost (Michelmore et al., 1991). This method relies on the differences of allelic frequencies between the phenotypically distinct pools and can be implemented with any type of marker data (Magwene et al., 2011). Traditionally, BSA studies were performed in biparental mapping studies. However, next-generation sequencing-based BSA approaches are also equally applicable in detecting QTL with small effects in diversity panels, given that effective population size and selection intensity are both adequate (Ehrenreich et al., 2010;Magwene et al., 2011). Hence, with the advent of highthroughput genotyping technologies and the availability of a reference genome, BSA is becoming an increasingly useful method of identifying TAV in both biparental mapping populations and diverse populations.
In this study, we performed BSA to identify genetic loci related to winter survivorship in lowland switchgrass populations. The lowland populations were subjected to natural winter pressure at two locations in southern Wisconsin. We hypothesized that alleles favorable to cold hardiness would be enriched in individuals that survive winters and that these alleles would be underrepresented in individuals that fail to survive winters. Therefore, using the differential rate of survival following winter stress, DNA pools were constructed for each population by bulking samples into survivor and nonsurvivor pools. These DNA pools were genotyped via exome capture sequencing. Using allele frequencies, we identified markers associated with winter survivorship in lowland populations.

Plant materials and experimental design
Twenty-nine diverse lowland switchgrass populations from the United States and Mexico (Table 1), representing a wide range of geographic areas, were used to test winter survival in this study. Eighteen populations were collected as prairie remnants in autumn 2014 (natural populations) and the other 11 populations were derived from breeding populations (BP) previously described Grabowski et al., 2017). These BP were all derived from lowland populations collected across the southern United States and previously subjected to selection for winter survivorship in USDA hardiness zones 4 and 5. These 11 populations were planted in Normal, IL, in 2012, allowed to open pollinate, and their seeds harvested in 2014.
In February 2015, seeds of each population were germinated in plastic seedling tubes containing a 1:1 mixture of topsoil and peat moss in the glasshouse at the U.S. Dairy Forage Research Center, Madison, WI. The seeds were cold stratified for 15 d prior to seeding. Seedlings that were 16 wk of age were transplanted to one of two locations: Arlington, WI (Plano silt loam; fine-silty, mixed, superactive mesic Typic Argiudoll), or Madison, WI (Kegonsa silt loam; finesilty over sandy or sandy-skeletal, mixed, superactive, mesic Mollic Hapludalf), in June 2015. The experimental design was a randomized complete block design with eight blocks at each location. Plots consist of one row of up to 25 plants on a 0.3m spacing with 0.9 m between rows. Although Arlington and Madison testing sites are located in the same USDA hardiness zone 4b, these two locations were chosen because of differences in soil type and 30-yr mean temperature. Weeds were controlled with pre-emergence herbicides and hand weeding as described by Casler (2008). Plants were allowed to grow throughout the 2015 season and their biomass was removed after a killing frost.
There was little winter injury during the first winter, 2015-2016, with >95% survivorship for all populations at both loca-

Core Ideas
• Bulk segregant analysis (BSA) is a low-cost alternative to association studies for identifying underlying QTL. • Lowland switchgrass planted in southern Wisconsin exhibited differential winter hardiness. • BSA test revealed 9 QTL in tetraploid and 7 QTL in octoploid populations. • QTL at 88 Mb on chromosome 2N, 115 Mb on 5K, 1 and 100 Mb on 9N are deemed potential. • Fixation of allele for winter survivorship in lowland population was not observed.
tions. Both experiments were kept weed-free in 2016 using herbicides and hand weeding as described above. Prior to senescence, leaf samples (approximately 5 cm) were collected from each plant in August 2016. Identical sized tissue samples were freeze-dried, sealed into airtight plastic bags, and stored at room temperature until the phenotypic scoring of survivorship was recorded following the second winter, 2016-2017. A paper hole-puncher was used to sample an equal amount of leaf tissue of all genotypes. Biomass of all plants was again removed after killing frost. Individual plants were scored twice for winter survival at approximately 10-30 d after the initial spring regrowth by visual assessment of the percentage of living shoots on the scale of 0-20, where 0 = dead, 1 = one green shoot, . . . , and 20 = no winter damage (Poudel et al., 2019a). The winter survivorship scoring system was intended to approximate a linear estimate of crown recovery, where a score of 20 was intended to translate to 100% of crown shoot survival. Based on the phenotypic scores for each population within a location, an equal number of plants was selected and bulked to create survivor and nonsurvivor tissue pools from each location, such that each pool contained no more than 10% (n = 8-20 individuals) of the original population size. The survivor pool consisted of samples with high survival scores and the nonsurvivor pool consisted of randomly selected dead plants (score = 0) while ensuring they represented all blocks within a location. Finally, 40 pairs of DNA pools (19 from Arlington and 21 from Madison) were constructed and used for DNA sequencing and further analysis based on the evidence of differential winterkill. For example, WALS_106 and WALS_109 did not show differential winterkill in the Arlington location and thus were not used in subsequent bulk segregant analysis ( Table 1) c NA, not available because the population were not selected for bulk segregant analysis.

Statistical analysis
An association test based on an individual population could result in p-value inflation due to unobserved population structure. Population structure is a concern because of prior reports that showed switchgrass can maintain mixtures of distant ancestral groups even in natural populations Morris et al., 2011). To minimize false-positive associations, the analysis targeted markers that were significant in multiple populations and in both evaluation environments.
To cluster groups of populations with similar genetic backgrounds, the populations were clustered based on genetic similarity using k-means cluster analysis in the R package stats. For each marker, a generalized linear mixed model was cre-ated predicting allele frequencies of all populations within a population group. Each model predicted allele frequencies as a binomial random variable with the read depth (the sum of reference and alternate allele reads) as the number of trials (i.e., logistic regression) (Yang et al., 2015). The binomial success probability was a linear function with the following fixed effects: survival status (survivor vs. nonsurvivor), population, and the interaction between population and survival status. The environment (Madison vs. Arlington) was included as a random variable. The test statistic for TAV significance was the likelihood ratio between models with and without survival status included. The ratio between models with and without the interactive effect between survival status and the population was also collected to estimate the reliability of the marker effect across multiple populations.
As discussed by Yang et al. (2015), the above statistic is prone to inflation and therefore two strategies were implemented to reduce inflation. First, a smoothing tri-cubic kernel with a sliding window was used on the test statistic to reduce the likelihood of spurious relationships (Magwene et al., 2011). The appropriate smoothing window size was estimated through visual observations of test statistic distortion near strong peaks in the dataset. The smoothed test statistic was further corrected for inflation by using a regression approach genomic control strategy (λ [Devlin & Roeder, 1999] and R package gap). Briefly, the tri-cubic smoothed test statistics was divided by the calculated inflation factor before calculating one-tailed chi-square distribution with one degree of freedom. Markers with false discovery rate (FDR) <0.05 were considered significant TAVs (Benjamini & Hochberg, 1995).
The post-hoc analyses assessed the change in allele frequency of TAV within each population. The change in allele frequency between survival and nonsurvival pool (Δ) was calculated as where RC = read count of either the favorable or unfavorable allele, TRC = total read count, S = survivor pool (desirable allele), and NS = nonsurvivor pool (undesirable allele). For each population and TAV, the two-sided p-values were calculated from Δ assuming the t-distribution. The p-values were then adjusted for FDR <0.01. In addition, the change in beneficial allele frequency among populations was assessed using linear regression on the minimum temperature of the coldest month (USDA hardiness zone) of the origin of the population. This test was conducted only in tetraploids. This test provided additional confidence to the identified QTL based on the knowledge that switchgrass exhibits latitudinal and longitudinal adaptability (Casler et al., 2007;Casler et al., 2004;Lovell et al., 2021;Poudel et al., 2019a). The annotation of the TAVs were accessed through Phytozome-13 (https://phytozome-next.jgi.doe.gov/ Goodstein et al., 2012).

Field-based winter mortality
Following the second winter (2016-2017), there was differential mortality among populations at both locations (Table 1). The percentage of winter mortality ranged from 0 to 94% among populations. The lower mortality rate in 2015-2016 compared with 2016-2017 may be attributable to the higher minimum temperature during the early winter months (November and December), providing a longer period of coldacclimation (Guy, 1990). Even though the locations occurred in the same hardiness zone and both locations received similar low-temperature stress (Supplemental Figure S1a), the survival in Madison was much higher than in Arlington. This difference could not be explained by either differential winter temperature or snow cover (Supplemental Figure S1b). Following the second winter, the mean percentage of living shoots was 9% at Arlington and 34% at Madison. Mean survivorship was 26% at Arlington and 61% at Madison (Table 1). Natural populations tended to be lower in mean survivorship compared with the BP (mean difference = 19%, p < .01), perhaps because of their breeding and selection history. The highest winter survival for both locations was observed with population SW788, a breeding population originating from New York, whereas the lowest winter survival was observed in a population from Mexico, largely as expected based on geographic origins.

Sequencing read depth and population differentiation
The average sequencing read count across all samples was 64 (N = 329,033 SNPs), after removal of positions with extreme read counts (15 < read depth for all bulks < 600) and minor allele frequency (< 0.05). A dendrogram of the bulked DNA samples constructed using reference allele frequencies clearly differentiated the populations into two distinct clades and aligned with k-means clustering (Figure 1). Four populations defined as octoploids (Poudel et al., 2019a) were monophyletic and separated from tetraploid populations, with the exception of two tetraploid populations: WALS_106 and WALS_109 from Louisiana and Mississippi, respectively. These two populations formed a distinct clade with a close resemblance to octoploids compared with tetraploids ( Figure 1). The exact reason for these two populations having a higher resemblance to the octoploids is unknown; however, these two populations behaved differently during phenotypic evaluation, with distinctly highest winter mortality together with the population from Mexico in both years of evaluation. As such, those two populations were excluded from the association analysis. The average read depth in octoploids was greater than for tetraploids (Supplemental Figure S2a) in both the survivor and nonsurvivor pools (paired t-test p < .001, mean difference = 4.5) similar to the results from Ramstein et al. (2018) in diverse switchgrass populations. Similarly, the standard error of sequencing depth was greater (paired t-test p < .001) in octoploids in both survivor and nonsurvivor pools (Supplemental Figure S2b).

QTL discovery
Parallel association analyses were conducted using the bulked DNA samples from tetraploid populations (n = 60, with populations 106 and 109 excluded) and octoploid populations (n = 16) to eliminate any impacts of ploidy on the association analysis and provide a constant genetic background for TAVs. Visual evaluation of distortion near significance peaks indicated that test statistics revert to background levels well below 500 bp. The smoothing tri-cubic kernel with a sliding window was therefore set to 250 bp. The genomic inflation factor (λ) of the smoothed function was 2.66 for the tetraploid populations and 1.89 for the octoploid populations (Figure 2a-b) prior to correction. Upon correction, the λ was approximately 1 for both tetraploids and octoploids (Figure 2a-b).
Using the tetraploid populations, 11 significant TAV comprising nine unique QTL were identified ( Table 2). Note that the QTL in this paper are coded as ploidy (4x for tetraploid and 8x for octoploid), followed by chromosome number, and physical base position in Mb (e.g., 4x.9N.1 is a QTL identified using 4x populations at position 1 Mb on chromosome 9N). The post-hoc t-test analyses of these TAV indicated that the change in allele frequency between survival and nonsurvival pools (Δ) was different (p < .01) in at least three and at most seven out of 15 populations irrespective of the site of origin or type of population (Figure 3). Changes in allele frequencies between survival and nonsurvival pools of DNA for three QTL (4x.5K.115, 4x.9N.1, and 4x.9N.100) were different in up to seven populations. However, there was no evidence for fixation of alleles or selection signatures for winter survivorship between survivor and nonsurvivor pools in any of these QTL as indicated by the small differences in allele frequency (Figure 3, Supplemental Table S1). The greatest Δ, equivalent to 0.09, was observed in QTL 4x.9N.1and 4x.9N.9 ( Table 2). Six of these 11 TAV corresponded to four unique QTL (4x.9N.1, 4x.9N.9, 4x.9N.26, and 4x.9N.100) located on chromosome 9N.
The beneficial allele frequency decreased with the increase in minimum temperature of the coldest month of the site of origin of the population in only three out of nine QTL  (Figure 4). This indicates that temperature was not a sole indicator of the change in allele frequencies. However, these three TAV provide additional support for the identified QTL. The association was greatest in QTL 4x.9N.9 (p < .01 and R 2 = 0.45) followed by QTL 4x.2N.88 (p = .02 and R 2 = 0.29) (Figure 4).
For the octoploid populations, nine significant TAV comprising seven unique QTL were identified ( Table 2). The posthoc analysis of these TAV indicated Δ was different (p < .01) in at least two out of four populations for six TAV, with the greatest Δ of 0.16 observed for QTL 8x.4K.62. The population from Georgia (WALS_103) contained only one QTL, but this region was associated with differential mortality in both test locations. As in tetraploids, there was no evidence of a selection signature for winter survivorship, as indicated by very small differences in allele frequency between the survival and nonsurvival pool ( Figure 5, Table 2, Supplemental  Table S1).
A single common QTL region at chromosome 5K has been identified from both tetraploid (115 MB) and octoploid (114 MB) population models. This particular QTL has been found to be associated with seven out of 15 tetraploid populations. The exact SNP of this QTL identified from tetraploid populations (Chr05K_115382984) has not been annotated. This QTL could be considered a promising QTL for improving both tetraploid and octoploid populations simultaneously. F I G U R E 2 Q-Q plot of the chi-square estimates based on (a) 4x populations and (b) 8x populations. The gray dots and black dots represented the -log10(p) before and after control of genomic inflation factor respectively. The green dots are the single nucleotide polymorphisms (SNPs) significant at false discovery rate <0.05

DISCUSSION
Since the introduction of BSA by Michelmore et al. (1991), it has been successfully implemented in several crops for both qualitative and quantitative traits using dominant or codominant markers, mostly in mapping populations. Examples include cold stress resistance in maize (Zea mays L.) (Farooqi et al., 2016), drought resistance in maize (Quarrie et al., 1999), drought resistance in rice (Oryza sativa L.) (Salunkhe et al., 2011;Vikram et al., 2012), salt tolerance in rice (Tiwari et al., 2016), Fusarium head blight in wheat (Triticum aestivum L.) (Shen et al., 2003), xylose utilization from Saccharomyces cerevisiae (Wenger et al., 2010), scab resistance in apple (Malus domestica Borkh.) (Yang et al., 1997), and cold tolerance in rice (Sun et al., 2018;Takagi et al., 2013). All these studies used biparental segregating populations for detecting associations. Among the few studies using a diversity panel for quantitative trait association in BSA was a study of kernel numbers in maize (Yang et al., 2015). Bulk segregant analysis had also been successfully used in diverse populations pertaining to the differentiation of adaptive and reproductive mechanisms among ecotypes within species (Gould et al., 2017;Jones et al., 2012;Twyford & Friedman, 2015).
In this study, we aimed to find winter survivorship QTL in natural and BP of lowland switchgrass. These types of populations had undergone numerous historical recombination events and thus had a different genetic structure than segregating populations (Falconer & Mackay, 1996). The power to detect a QTL in natural and BP using BSA is limited mainly due to rapid linkage disequilibrium (LD) decay, accompa-nied by higher segregation distortion around the QTL as a result of random mating for several generations compared with mapping populations (Quarrie et al., 1999). To account for these spurious associations, smoothing of the G values using kernels is recommended (Magwene et al., 2011). In a biparental or multiparent mapping population, the LD decay rate is slower so the window width of 15-30 cM (Magwene et al., 2011) could be used to smoothen the G values. However, in diverse populations, averaging of the read count by constructing a sliding window with bin size equivalent to LD decay rate is commonly practiced (Gould et al., 2017). In this study, the smoothing of the G-value with a 250 Kb window size was done based on the results from Ramstein et al. (2016) and Poudel et al. (2019a) who suggested that LD decay is less than 5 cM in a similar diverse panel of lowland switchgrass populations.
The power of the BSA test can be increased by two mechanisms: (a) increasing the number of individuals in the DNA pools and (b) using sufficient sequencing depth (Magwene et al., 2011). In whole-genome sequencing, a larger pool size increases the accuracy of 'Pool-seq', and, in most comparisons, a pool size of 50 individuals produces sufficiently accurate allele frequency estimates (Schlötterer et al., 2014). Magwene et al. (2011) recommended to use at least 50 individuals per pool from a diversity panel, but in practice BSA analyses were successful with DNA pools of consisting 10-20 individuals in a segregating population from a biparental mating with recombinant inbred lines (Quarrie et al., 1999;Sun et al., 2018). In this study, we bulked only 10% of the individuals (n = 8-20) for the survivor and nonsurvivor pools. The detection of significant QTL regions could be affected by Type I error (Navabi et al., 2009), but this error is expected to decrease because of high marker density (Sun et al., 2010). The average read depth on the selected marker used for analysis in this study was only 60. However, this sequencing depth is sufficiently large enough to detect even a weak QTL effect given the bulk size of 20 or more DNA samples (Magwene et al., 2011;Yang et al., 2015). The analysis was performed separately on tetraploid and octoploid population and accounted for population and evaluation location effects within the logistic regression model F I G U R E 4 Linear regressions of favorable allele frequencies for three single nucleotide polymorphism (SNP) markers as a function of the mean temperature of the coldest month (30-year mean) for 15 lowland switchgrass populations evaluated at two Wisconsin locations with the goal of isolating QTL that will be of practical value for breeding. Despite these controls, we reported different QTL compared with the prior freezing tolerance study on lowland × upland hybrid biparental populations (Poudel et al., 2019b). In addition to the different populations used in these two studies, the discrepancy in the detected QTL could be attributed to the differences in the type of environment used to screen the genotype. Plants react differently to the varying light spectrum, wind, and available soil resources under controlled environment conditions compared with field conditions, and it is not unusual to detect different QTL in different environments (Dhanaraj et al., 2007). Overall, the QTL from a biparental population would be of interest to the lowland-upland hybrids. However, the QTL identified in the F I G U R E 5 Change in allele frequency for three single nucleotide polymorphism (SNP) markers associated with selection for survivorship (difference between survivor and nonsurvivor groups) for four lowland switchgrass populations evaluated at two Wisconsin locations and classified according to state of origin (** indicates a significant selection effect at p < .01) current study can be useful to improve multiple populations simultaneously.

QTL across ploidy levels
Octoploids have tetrasomic inheritance and accurate dosage call for intermediate heterozygous classes in octoploids (simplex, duplex, or triplex) require genotyping at high resolution (Ramstein et al., 2018;Uitdewilligen et al., 2013). This highlights an advantage of BSA over the GWAS studies for associating the marker effect in octoploids because BSA relies on sequencing depth for analysis which eliminates the need to classify intermediate heterozygous classes. However, the effects of ploidy on marker effect estimates will need to be considered if the individuals from different ploidy levels are bulked together. In such cases, genome doubling and dominance in intermediate heterozygous classes will impact QTL detection (Endelman et al., 2018) and might affect the ability to detect QTL. In this study, we assessed the QTL from the tetraploid and octoploid populations separately. The QTL regions associated with winter survivorship did not overlap among tetraploid and octoploid populations, with the notable exception of one QTL on chromosome 5K around 114-115 Mb. These results attest to the wide genetic divergence between tetraploid and octoploid lowland switchgrass populations at least in terms of winter survival.
In tall fescue (Festuca arundinacea Schreb.), the expression of NDPK in plants treated with 200-mM sodium chloride increased by up to 37-fold compared with untreated plants (Shahabzadeh et al., 2019).The QTL at chromosome 5K at position ∼114-115 Mb were identified in both octoploid and tetraploid populations with maximum Δ of 0.31 in Kanlow, OK (Supplemental Table S1). The exact peak SNP from the octoploid populations is located within the coding region for the protein CULLIN (Phytozome 13.1/P. virgatum 4.1). CULLIN is an important component of the ubiquitin proteasome pathway, which in conjugation with RING E3 ligase degrades ICE1 and is thus regarded as a negative regulator to cold stress (Choi et al., 2014;Dong et al., 2006). Similarly, probable membrane protein DUF221 protein is collocated with these QTL at 115 Mb. DUF221 proteins are a family of osmosensitive calcium-permeable cation channels that regulate cold, drought, and salt tolerance (Ganie et al., 2017;Hou et al., 2014;Sanders et al., 1999). In plants, cold, drought, and salt stresses all stimulate the accumulation of compatible osmolytes and antioxidants (Hasegawa et al., 2000;Xiong et al., 2002). Freezing damage is caused by osmotic dehydration triggered by extracellular ice formation, which leads to cell lysis (Sakai & Larcher, 2012;Sandve et al., 2011;Zhang et al., 2010). As such, the role of DUF221 proteins in response to freezing can be anticipated, but specific research of this protein under cold or freezing stress has not been reported.
Four out of nine QTL identified in tetraploid populations located on chromosome 9N at positions 1, 9, 26, and 100 Mb. An additional QTL identified from the octoploid populations is also located on chromosome 9N at position 41 Mb. The N subgenome is believed to have received genetic introgression of adaptive loci for key biomass-related traits from upland to lowland populations (Lovell et al., 2021). Interestingly, phytochrome B (PHYB) genes are located at positions 9 and 98 Mb (phytozome 13.1/P. virgatum 4.1; Bahri et al., 2018), suggesting the relevance of our identified QTL as an adaptation signature to winter survival within the lowland tetraploid switchgrass population. Previous work by Bahri et al. (2018) reported PHYB, a gene involved in photoperiod response, was responsible for ecotypic differentiation and population adaption. PHYB has been reported to negatively impact cold and/or freezing tolerance in several crops (Franklin & Whitelam, 2007;He et al., 2016;Jiang et al., 2020;Lee & Thomashow, 2012;Wang et al., 2016). The transcript within the QTL 4x.9N.9 encodes for histone H2A. Histone modification through acylation is considered to be involved in several stress responses (Kim et al., 2015) and H2A has been suggested to play an important role in thermal stress response (Boden et al., 2013;Kumar & Wigge, 2010). Similarly, the SNP peak located in QTL 4x.9N.100 has been described as HEAT SHOCK 70 kDa protein 1/8 (Table 2), a regulator of heat stress (Sung et al., 2001) that may have an important role as plants recover from freezing injury. A low temperature and salt responsive protein LTI6 (Pavir.9NG688000) is co-located with this QTL (Phytozome 13.1/P. virgatum 4.1). The candidate genes and proteins colocated with our identified QTL did not show strong selection signature for winter survivorship across multiple lowland populations (data not shown), which could partly be due to the neutralizing allelic effect as a result of bulking populations with different genetic background. These suggestive QTL need further validation for use in marker-based selection schemes but can possibly be incorporated into genomic selection protocols to improve selection accuracy.

Future prognosis
Overall, this study has established a better understanding of the genetic architecture of winter survivorship in lowland switchgrass and we provide a starting point for markerbased selection towards improving winter survival in lowland switchgrass. We did not identify any conserved regions for winter survivorship among the populations and thus infer that winter survival in switchgrass is regulated by many QTL with small additive effects which may often be population specific. We propose that the QTL detected in populations across the broad region of focus in this study, including the Gulf Coast and Southern Great Plains, can be used to improve multiple populations simultaneously. Conversely, QTL detected only in specific populations will be limited to the improvement of those populations only. Understanding and exploiting these QTL should focus on the development of lowland switchgrass for northern latitudes with reduced winter mortality. These results suggest that it should be possible to develop additional and broader-based lowland switchgrass populations adapted to USDA hardiness zones 4 and 5, using germplasm from hardiness zones 7 through 9. Simultaneous stacking of the all 17 putative QTL identified in this study would be logistically challenging, except through genome-wide approaches (Collard et al., 2005). For most marker-based selection approaches, QTL are prioritized based on their proportion of phenotypic variation (Ribaut & Betrán, 1999), but a distinct disadvantage of BSA is the absence of QTL effect estimates. Therefore, we can instead prioritize significant QTL based on biological annotations, the reliability of associations across multiple populations, and relation between allele frequencies and long-term minimum temperature of the coldest month of the site of origin of population under study. Based on source of information, four QTL are considered most promising: 4x.2N.88, 4x.9N.1, 4x.5K.115 or 8x.5K.114, and 4x.9N.100. The first two QTL showed significant associations to minimum temperature and reliable changes in allele frequency across multiple populations. The QTL 4x.5K.115 or 8x.5K.114 that encodes for the CULLIN protein were simultaneously identified from the same genomic region in both tetraploid and octoploid populations. The QTL 4x.9N.100 is considered promising as it encodes for HEAT SHOCK 70 kDa protein 1/8. Designing PCR-based flanking markers to easily identify these QTL could allow rapid identification of superior genotypes across diverse southern lowland populations for adaptation to the northern United States.

D A T A AVA I L A B I L I T Y S T A T E M E N T
Raw reads are available in the National Center for Biotechnology Sequence Read Archive under BioProject PRJNA529229. Due to size, the SNP data matrix containing read counts of 80 bulked samples extracted for 5,596,351 bi-allelic loci is available as WALS_BSA_ReadCount_SNP_Matrix.txt.gz on the Dryad Digital Repository at https://doi.org/10.5061/ dryad.9cnp5hqhv.

A C K N O W L E D G M E N T S
The authors thank the University of Wisconsin Biotechnology Center DNA Sequencing Facility for DNA sequencing services. We thank Thomas Juenger, University of Texas, Austin; DK Lee, University of Illinois, Urbana-Champaign; and Yanqi Wu, Oklahoma State University, Stillwater for providing some of the germplasm used in this experiment. This work was supported by funds from the U. S. Department of Agriculture National Institute of Food and Agriculture (2014-67009-22310) to CRB and MDC.