Genome Assembly. The Omni-C chromosomal assembly produced a near complete assembly for B. tectorum. The assembly was 2,482 megabases (Mb) in total length (Supplementary Table S1). The contig N50 was 19.4 Mb, the scaffold N50 was 357.4, and the resulting assembly contained 92.1% of the BUSCO genes with 259 gaps. The final assembly had a L50 of four and a L90 of seven corresponding to the seven chromosomes (x = 7) expected for members of the Pooideae in the the Poeae, Aveneae, Bromeae and Triticeae tribes27.
Based on CDS comparisons using MCScanX28, the seven largest scaffolds of the B. tectorum reference genome were found to have a one-to-one syntenic relationship with the seven chromosomes of the barley genome – not surprising given their phylogenetic proximity. As expected, synteny between B. tectorum and H. vulgare was highest in the chromosome arms where gene density is known to be high and lowest though the centromeric region where gene density is substantially reduced (Fig. 2). A translocation is present between chromosomes two and five where the first half of Bt5 has synteny with chromosome two in barley and the second half of Bt2 has synteny with chromosome 5 in barley (Fig. 2).
Phenotypic Analysis. Broad Sense Heritability (Reliability), the mean, standard deviation, min, max and correlation Best Linear Unbiased Estimates (BLUEs) were calculated to characterize the variation in traits and obtain a measure of the effect of each genotype. The distributions of BLUEs for reproductive phenology traits indicated a bimodal distribution (Fig. 1), where more rapidly flowering and taller plants consisted of accessions from Washington and GRIN collections. In contrast, Montana accessions flowered later and were shorter (Supplementary Table S3). Broad sense heritability was high for all the traits and ranged from 0.94 for tiller number to 0.99 for days to first visible panicle (VPN) (Table 1). Genotypes had a wide range of BLUEs for all traits measured: plant height (PH) ranged from 36.5 to 89 cm, number of tillers ranged from 5.5 to 38.5, and days to first ripe seed (FRS) ranged from 43.2 to 112.5 d (Table 1). Spearman correlations between BLUEs of different reproductive phenology traits were all above 0.95 (Fig. 1). Reproductive phenology and PH traits were moderately negatively correlated, with Spearman correlations ranging from the − 0.44 (PH and FRS) to -0.51 (PH and J1).
Genome-wide Association Mapping for Height, Tiller Number and Phenology Traits. To identify regions of the genome associated with variation in adaptive traits, a GWAS was performed on 121 genotypes to find QTL for PH, VPN, days to first visible joint (J1), FRS, days to 50% ripe seed (AWN50) and tiller number using BLINK with 3 principal components and significance threshold of 0.05 after multiple testing correction. Nineteen QTLs were significantly (q < 0.05) associated with PH and reproductive phenology, with one QTL for PH (1), nine QTL for VPN (9), three QTL for J1 (3), nine QTL for FRS (9), and three QTL for AWN50 (3). A QTL on Bt6 (Bt6:8628087) was significant (q < 0.05) for all the reproductive phenology related traits except J1, where a second significant (q = 2.070E-07) QTL for J1 was at a nearby position (Bt6:8417488) on Bt6 (Table 2). The QTL on Bt1 (Bt1:276755111) was significant (q < 0.05) for J1, FRS, and AWN50 (Table 2). QTLs on Bt2 (Bt2:9403921) and Bt7 (Bt7:347158582) were significantly (q < 0.05) associated with VPN and FRS (Table 2).
Plant Height. We identified a QTL that was significantly (q = 1.360E-05, Table 2) associated with PH on Bt6:301800092. The QTL explained 16.9% (Table 2) of the phenotypic variation. The MAF at the PH QTL on Bt6 was 0.44, indicating both allelic states are common in our panel. Searching the area flanking the QTL associated with PH at Bt6:301800092 on both sides by up to 500 Kb revealed a homolog of Xanthine Dehydrogenase (XDH) 29 Kb from the QTL and a homolog of Indole-3-pyruvate monooxygenase YUCCA6 (YUC6) 242 Kb from the QTL (Table 3). XDH and YUC6 are promising candidate genes for the QTL associated with PH we identified on Bt6 as they both are well document to be associated with changes in PH, senescence and response to drought29,30. In rice (Oryza sativa L.) overexpression of the XDH homolog led to increased PH while under-expression of the XDH homolog resulted in reduced PH29, indicating that homologs of XDH would be fitting candidate genes for PH. Indeed, in rice, a GWAS identified XDH as a candidate gene for coleoptile length in response to flooding31 indicating XDH homologs may be involved with stem elongation in grasses and thus PH. Furthermore, in Arabidopsis knock-out mutations of the XDH gene led to reduced PH32. Knock-out mutations of the YUC6 gene in Arabidopsis was found to increase auxin production in the shoots leading to reduced PH33. Further investigation, including gene expression experiments and/or knockout of the XDH gene homolog in B. tectorum, is needed to validate and further understand the underlying genetics controlling the large effect QTL for PH.
Unexpectedly, PH and reproductive phenology timing were negatively correlated, contradicting previous findings that PH and reproductive phenology were positively correlated (i.e., earlier flowering plants do not have as much time to grow)34,35. Lolium perenne (L.) was also found to have a negative genetic correlation between PH and flowering time36. The negative correlation was thought to be the result of selection imposed by grazing, where plants biological fit plants remained short until they were ready to flower at which point they elongated and flowered quickly37. Bromus tectorum may also be under similar grazing selection pressure38. Competition with crops could also select for taller, fast-growing phenotypes to facilitate competition for space.
Phenology Traits. Flowering time is an important adaptive trait for ensuring the survival of plants in a broad range of climates, often driving local adaptation39,40,41. In predominantly self-fertilizing species, large effect QTLs control the variation in flowering time42,43, in contrast to maize where flowering time is controlled by many QTL of small effect44. Here, we phenotyped four traits (VPN, J1, FRS, AWN50) associated with flowering time to facilitate a GWAS to identify candidate reproductive phenology genes. Days to first visible panicle (VPN) was the first reproductive phenology stage observed and had the highest heritability at 0.99 (Table 1). The largest effect estimated for a significant (q = 3.630E-07) QTL for VPN was found at Bt7:70795764, explaining 16.3% of the variation in VPN with a MAF of 0.26. The eight-remaining significant (q < 0.05) QTL on Bt3 (Bt3:2025814; q = 8.210E-07), Bt4 (Bt4:382636922; q = 1.920E-05), Bt6 (Bt6:8628087; q = 4.830E-03), Bt7 (Bt7:347158582; q = 1.046E-02), Bt2 (Bt2:9403921; q = 3.485E-02), Bt3 (Bt3:383345273; q = 3.797E-02), Bt3 (Bt3:324738958; q = 3.797E-02) and Bt3 (Bt3:324739133; q = 3.797E-02) explained 4.3, 4.3, 2.9, 4.7, 3.3, 3.8, 3.4 and 3.4 percent of the phenotypic variation, respectively, with MAF ranging from 0.2 to 0.49 (Table 2).
Developmentally, days to first joint (J1) was the second reproductive phenology trait to occur. Three significant (q < 0.05) QTL were associated with J1. The significant QTL at Bt6 (Bt6:8417488; q = 2.070E-07) was 12 Kb from a homolog of HDR1 and explained 6.7% of the phenotypic variation with an MAF of 0.35 (Table 2 and Table 3). The QTL at Bt1 (Bt1:276755111; q = 3.870E-05) explained 8.1% of the phenotypic variation with a MAF of 0.2. The other QTL on Bt1 (Bt1:169814345; q = 1.909E-03) explained 4.2% of the phenotypic variation and had an MAF of 0.35.
The next developmental stage was first ripe seed (FRS). Days to first ripe seed was associated with nine significantly (q < 0.05) associated QTL. The leading QTL located on Bt7 (Bt7:347158582; q = 1.360E-05) explained 7.2% of the phenotypic variation of FRS and was common (MAF = 0.29) (Table 2). The remaining QTLs on Bt7 (Bt7:347158582; q = 5.280E-09), Bt1 (Bt1:173453655; q = 1.582E-04), Bt1(Bt1:39144935; q = 1.786E-04), Bt3 (Bt3:48964562; q = 1.786E-04), Bt1 (Bt1:276755111; q = 2.823E-04), Bt2 (Bt2:9403921; q = 4.166E-04), Bt4 (Bt4:382201414; q = 1.026E-03), Bt3 (Bt3:4974562; q = 1.836E-03) and Bt6 (Bt6:8628087; q = 4.110E-02) explained 1.4, 1.8, 1.2, 4.1, 1.9, 0.6, 4.0 and 2.3 percent of the variation, respectively (Table 2). The MAF of QTL significantly associated with FRS were lower than those identified for the other reproductive phenology traits, varying from 0.058 to 0.29.
The final developmental reproductive phenology trait to occur was AWN50, resulting in three significant associations (q < 0.05). The leading QTL on Bt1 (Bt1:276755111; q = 6.760E-05) explained 6.5% of the phenotypic variation and had a MAF of 0.2 (Table 2). The other two QTL were located on Bt5 (Bt5:348443767; q = 1.637E-04) and Bt6 (Bt6:8628087; q = 5.973E-03) explaining 5.6% and 4.8% of the phenotypic variation with MAFs of 0.37 and 0.2, respectively.
Survey of Candidate Genes Associated with Phenology. Seventeen genes were identified as punitive candidates for reproductive phenology traits, based on their proximity to QTL (within a 500Kb interval on either side of the significant SNP) and previous functional characterization related to specific phenologies. The most promising of these genes is a homolog of Heading Date Repression 1 (HDR1) The proximity of a homolog of HDR1 with QTL from all the maturity related trait analyzed suggests that HDR1 significantly influences maturation in B. tectorum. In O. sativa, HDR1 promotes Heading date 1 (Hd1) and represses Early heading date 1 (Ehd1) delaying flowering45. Knockout mutations or RNA interference of HDR1 resulted in rice plants that flowered 30 days earlier in long day light conditions, making HDR1 a promising candidate gene for maturity traits45.
Multiple candidate genes were identified for several reproductive phenology associated QTL, indicating multiple genes may underlie reproductive phenology in B. tectorum. The QTL on Bt2 (9403921) associated with VPN and FRS was 13Kb, 29Kb and 133Kb from homologs of ABC transporter B family member 19 (ABCB19), Cullin-3A (CUL3A), and Gibberellin 20 oxidase 2 (GA20OX2), respectively. Loss of function mutations in ABCB19, CUL3A, and GA20OX2 led to longer flowering times, indicating all three promote advancement of reproductive phenology46,47,48. FRS5 is a punitive transcription factor regulating far-red light control of development49, associated with adaptation to photoperiod50. Dehydration-responsive element-binding protein 1F (DREB1F) is a putative transcription factor, when upregulated, represses GA biosynthesis catalyzed by GA20OX genes, causing shorter plants with delayed flowering51. A likelihood ratio test of the interaction between the FRS QTLs on Bt2 (9403921) and Bt1 (173453655) revealed epistasis is likely (X2 = 3.71 df = 1, p = 0.054), indicating DREB1F and GA20OX are interacting to control FRS in an epistatic manner.
The QTL associated with VPN on Bt3 (2025814) was 120Kb, 177Kb and 240Kb from homologues of BTB/POZ and MATH domain-containing protein 1 (BPM1), BTB/POZ and MATH domain-containing protein 2 (BPM2) and PHYTOCHROME-DEPENDENT LATE-FLOWERING (PHL), respectively. PHL triggers flowering under long-day conditions by repressing the phytochrome b (PHYB) and the constans (CO) genes52. BPM1 and BPM2 are transcription factors of the BTB/POZ and MATH domain-containing protein (BPM) gene family and are involved in regulating flowering time by making proteins that are part of the Cullin E3 ubiquitin ligase complexes that include CUL3A53. A likelihood ratio test did not detect epistasis (X2 = 0.02, df = 1, p = 0.901) between the QTL on Bt2 (9403921) near CUL3A and the QTL on Bt3 (2025814) near BPM1 and BPM2.
A QTL on Bt3 (324738958) associated with VPN was located 64Kb, 113Kb and 189Kb from homologs of UPSTREAM of FLC (UFC), Ultraviolet-B receptor 8 (UVR8) and Early Flowering 7 (VIP2), respectively. UVR8, VIP2 and UFC are involved with the regulation of Flowering Locus C (FLC) in A. thaliana54,55,56, suggesting one mechanism maintaining genetic variation in B. tectorum flowering time is regulation of FLC expression. Homologs of four genes in the FHY3/FAR1 gene family were near QTL found associated with maturity traits indicating the FHY3/FAR1 gene family is part of a mechanism maintaining variation in maturity traits in B. tectorum. Homologues of FHY3 and FRS6 were found 147 and 149 Kb from a QTL on Bt7 (Table 3), respectively. A homologue of FRS7 was near the QTL on Bt3 associated with VPN (Table 3) FRS5 discussed earlier is also member of the FHY3/FAR1 gene family. The FHY3/FAR1 gene family is comprised of 14 homologous genes, regulating transcription as a response to far red light (Lin and Wang, 2004). FHY3, FRS6, FRS7 have been demonstrated to regulate flowering time in A. thaliana48,57,58 making members of the FHY3/FAR1 good candidate genes for maturity traits.
Additionally, a Cytokinin dehydrogenase 2 (CKX2) homolog was identified as a candidate gene for the QTL at Bt3:48964562 associated with FRS. CKX2 catalyzes the oxidation of cytokinins59. Cytokinins have been shown to regulate flowering time via transcriptional activation of Twin Sister of FT (TST)60 and overexpression in members of the Cytokinin dehydrogenase (CKX) gene family have shown to delay flowering time in long day conditions59 indicating that the CKX2 homolog could be controlling maturity traits through the regulation of cytokinins in B. tectorum.