Homoeolog-specific activation of genes for heat acclimation in the allopolyploid grass Brachypodium hybridum

Abstract Background Allopolyploid plants often show wider environmental tolerances than their ancestors; this is expected to be due to the merger of multiple distinct genomes with a fixed heterozygosity. The complex homoeologous gene expression could have been evolutionarily advantageous for the adaptation of allopolyploid plants. Despite multiple previous studies reporting homoeolog-specific gene expression in allopolyploid species, there are no clear examples of homoeolog-specific function in acclimation to a long-term stress condition. Results We found that the allopolyploid grass Brachypodium hybridum and its ancestor Brachypodium stacei show long-term heat stress tolerance, unlike its other ancestor, Brachypodium distachyon. To understand the physiological traits of B. hybridum, we compared the transcriptome of the 3 Brachypodium species grown under normal and heat stress conditions. We found that the expression patterns of approximately 26% and approximately 38% of the homoeolog groups in B. hybridum changed toward nonadditive expression and nonancestral expression, respectively, under normal condition. Moreover, we found that B. distachyon showed similar expression patterns between normal and heat stress conditions, whereas B. hybridum and B. stacei significantly altered their transcriptome in response to heat after 3 days of stress exposure, and homoeologs that were inherited from B. stacei may have contributed to the transcriptional stress response to heat in B. hybridum. After 15 days of heat exposure, B. hybridum and B. stacei maintained transcriptional states similar to those under normal conditions. These results suggest that an earlier response to heat that was specific to homoeologs originating from B. stacei contributed to cellular homeostasis under long-term heat stress in B. hybridum. Conclusions Our results provide insights into different regulatory events of the homoeo-transcriptome that are associated with stress acclimation in allopolyploid plants.

genome duplication events occurred widely in plants at the Cretaceous-Paleogene boundary, which is a major extinction event in the earth's history. Evolutionary views represent a hypothesis that plants with duplicated genomes might have a better chance for survival under global adverse conditions [11,12].
Interspecific hybridization and subsequent genome duplication led to evolutional changes in hybrid species that represented fixed heterozygosity. Allopolyploid plants generally show better growth vigor and stress tolerance than their ancestors [13][14][15]. For example, relative to their ancestors, allopolyploid Arabidopsis (Arabidopsis suecica) shows more vigorous growth, and allopolyploid Spartina (Spartina anglica) shows better tolerance to reducing conditions and sulfite-rich sediments [16][17][18]. Both natural and synthetic wheat (Triticum aestivum) have higher fitness under salt stress than their diploid and tetraploid ancestors [19,20]. It has been suggested that allopolyploidization may have contributed to the adaptation to a wide range of environmental conditions [21,22].
The expression levels of homoeologs in allopolyploid species often show nonadditive expressions that are deviated from parental additivity compared to their ancestors [23]. The nonadditive expression has been widely observed in various allopolyploid species [24][25][26] and can be explained by at least 3 possible scenarios: total gene expression of a homoeolog group in an allopolyploid species is similar to that of 1 of its parental species (expression level dominance [ELD]), total gene expression is lower or higher than in both parents (transgressive expression), and uneven contribution of homoeologs to gene expression (homoeolog expression bias) [23,27]. Through comparative gene expression analysis of total expression levels of homoeolog groups in an allopolyploid species and "mid-parental values" as the average expression levels of genes in its parents, parental additivity and nonadditivity of homoeologs can be examined. By discrimination of homoeologous gene expression, homoeolog expression bias can be identified across tissues, developmental stages, and environmental conditions [28][29][30].
The contribution of expression plasticity of duplicated genes observed in diverse allopolyploid species to adaptive and ecological fitness has been debated. Comparative transcriptome analysis of the natural allopolyploid Coffea arabica and its ancestral species (Coffea canephora and Coffea eugenioides) revealed evidence of a genomic ELD that depends on growth temperature [31]. Comparative transcriptome analysis of hexaploid Triticum aestivum (AABBDD) and its ancestors (Triticum turgidum [AABB] and Aegilops tauschii [DD]) demonstrated differential enrichment of overrepresented gene functions between the ELD genes of both parents, which suggests they may have differentially contributed to particular biological functions in hexaploid wheat [32]. Comparative gene expression analysis of a synthetic tetraploid wheat and its parental species showed transgressive expression of genes that are involved in particular biological functions such as transport, modification, and uroporphyrinogen decarboxylase activity, suggesting that transgressive gene expression may have rapidly occurred following allopolyploidization [33]. Studies on the allotetraploid Arabidopsis kamchatica and its ancestral species demonstrated homoeolog expression bias in genes related to adaptation to cold and heavy metal environments, which were inherited from the respective ancestors that are tolerant against cold and heavy metal stresses, respectively [34,35].
Brachypodium hybridum (2n = 30) is a natural allopolyploid that is derived from a cross between Brachypodium distachyon (2n = 10) and Brachypodium stacei (2n = 20) that occurred approximately 1 million years ago (MYA; Additional file 1, Fig. S1) [36][37][38][39][40][41][42][43]. Although these species inhabit a circum-Mediterranean region, their environmental niches are clearly different. Brachypodium distachyon grows in higher-altitude, cooler, and wetter areas, whereas B. stacei is found in lower-altitude, warmer, and drier areas. The hybrid species B. hybridum grows in areas that are overlapping as well as specific to its ancestors; this growth pattern suggests that its speciation is associated with particular environmental conditions and adaptations to diverse environmental conditions [37]. The broader growth area of B. hybridum compared to its ancestors led us to the hypothesis that its duplicated genes may have contributed to its adaptive and ecological fitness through gene expression changes following allotetraploidization. Specifically, in the B. hybridum transcriptome, each of the homoeologs may show particular expression patterns in response to abiotic stimuli, which are associated with the ancestral growth habitats and gene expression. Therefore, this trio of species has been recently proposed as a model for grass speciation via adaptation and polyploidization [39,40,43].
In the current study, we conducted a comparative global transcriptome analysis of B. hybridum and its ancestors. The analysis revealed nonadditive transcriptome changes in the leaf and root tissues of B. hybridum. Additionally, we performed a homoeolog-specific transcriptome analysis by discriminating RNA-sequencing (RNA-seq) reads of each homoeolog of B. hybridum and determined the nonancestral gene expression patterns. Finally, we assessed homoeolog-specific transcriptome changes in response to heat stress in B. hybridum. Here, we discuss the differential regulation of the homoeo-transcriptome that is associated with heat stress tolerance in B. hybridum inherited from B. stacei.

Plant materials
Three Brachypodium species were used in this study: the allotetraploid B. hybridum Bd14-1 and the diploid ancestors B. distachyon Bd21 and B. stacei ABR114. The accessions of these species were provided by the National Plant Germplasm System of the US Department of Agriculture-Agricultural Research Service (David F. Garvin; Plant Science Research Unit, University of Minnesota, USA) and Pilar Catalán (Department of Agriculture and Environment Science, High Polytechnic School of Huesca, University of Zaragoza, Spain). Dry seeds were incubated on wet filter paper in a Petri dish at 4 • C in the dark for 6-7 days to synchronize germination. The germinated seeds were grown in a growth chamber at 25 • C under a 16-h day photoperiod (60 μmol·m -2 ·s -1 ) for 4 days. The plants were transplanted to pots filled with autoclaved Pro-Mix R BX Mycorrhizae TM (Premier Tech, Quebec, Canada). The potted plants were grown in a growth chamber at 22 • C (normal condition) or 32 • C (heat stress condition) under a 20-h day photoperiod (100 μmol·m -2 ·s -1 ) and watered with 5000-fold diluted Professional Hyponex 10-30-20 (Hyponex Japan, Osaka, Japan) every 3 or 4 days.

Whole genome sequence data
Genomic DNA from the leaf tissues of B. hybridum and B. stacei were extracted using the DNeasy Plant Mini Kit (Qiagen, Tokyo, Japan). Libraries for single-end DNA sequencing were obtained using Ion Xpress TM Plus Fragment library kits (Life Technologies Japan, Tokyo, Japan); semiconductor chips that were used for sequencing were prepared using the Ion OneTouch 2 System (Life Technologies Japan) and Ion P1 TM Chip v2 (Life Technologies Japan). The sequencing analyses were performed using an Ion Proton sequencer (Life Technologies Japan). The genome sequence data of B. hybridum and B. stacei were archived at DNA Data Bank of Japan (DDBJ) under the accession number DRA005717 (Additional file 2, Table S1).

RNA-seq data
Shoots and roots from each species grown in a growth chamber for 4 days after synchronized germination were sampled to elucidate the global homoeolog expression patterns in B. hybridum. Shoots and most young leaf blades from each species grown in a growth chamber under different temperature conditions (22 • , normal condition, and 32 • C, heat stress condition) for 3 and 15 days, respectively, after being transplanted were sampled to elucidate the homoeolog-specific transcriptional response to heat stress in B. hybridum. Total RNA was extracted from each sample using the RNeasy plant mini kit (Qiagen). Poly(A) RNAs were purified using the NEBNext R Poly(A) mRNA Magnetic Isolation Module (New England Biolabs, MA, USA). The libraries for single-end strand-specific RNA-seq were obtained using the Ion Total RNA-seq kit v.2 (Life Technologies Japan). Size-selected libraries were purified using the Agencourt AMPure XP (Beckman Coulter, CA, USA). Semiconductor chips that were used for sequencing were prepared using the Ion PI TM Hi-Q TM OT2 200 kit (Life Technologies Japan), Ion PI TM Hi-Q TM Sequencing 200 kit (Life Technologies Japan), and Ion P1 TM Chip v3 (Life Technologies Japan). The sequencing analyses were performed using an Ion Proton sequencer (Life Technologies Japan) with 3 biological replicates. The RNA-seq data were archived at DDBJ under the accession number DRA005699 (Additional file 2, Table S1).

B. hybridum and B. stacei show significant tolerance to long-term heat stress
Unlike B. distachyon, B. hybridum and B. stacei showed significant tolerance to long-term heat stress. The allopolyploid B. hybridum and its ancestor B. stacei grow in warmer regions, in contrast to the other ancestor B. distachyon, which suggests that B. hybridum and B. stacei, but not B. distachyon, might have adapted to hightemperature conditions. To test this hypothesis, we compared the plant biomasses of B. hybridum, B. stacei, and B. distachyon grown under normal (22 • C) and heat stress (32 • C) conditions at 2 time points (after 3 and 15 days of heat stress exposure). After 3 days of exposure to heat stress, the 3 species showed no differences in growth (Fig. 1A). Conversely, after 15 days of heat exposure, B. distachyon showed a significant decrease in fresh weight (P < 0.01, t test), whereas B. hybridum and B. stacei maintained their growth (Fig. 1B)  and B. hybridum Bd14-1 and mapped the reads to the reference genome sequence of B. distachyon Bd21. We found that 85% and 91% of the genomic reads of B. stacei and B. hybridum mapped to the reference genome and covered 89% and 98% of the genic region of the reference genome, respectively, which suggests high similarity among the homoeologous genomes ( Fig. 2 and Additional file 2, (1) identified in the genic regions (Additional file 1, Fig. S4A). Specifically, 45.4% and 17.6% of SNPs are distributed in intronic and exonic regions, respectively, and 12.5% and 1.6% of SNPs are distributed in untranslated regions and splice sites, respectively. We assessed the types of polymorphism found in genic regions and found that 26 986 genes (79% of the annotated genes in the Bd21 genome) contain nonsynonymous variations. Of these, 4558 genes (13% of the annotated genes in the Bd21 genome) contain several variations in the B. stacei genome; 683, 2980, and 895 genes show start codon loss, stop codon gain, and stop codon loss, respectively (Additional file 1, Fig. S4B and Additional file 2, Table S4). Based on the synonymous substitution rate (Ks) between the homoeologs, we estimated a divergence time of approximately 6.4 MYA between B. distachyon and B. stacei, which is in accordance with the divergence time of B. stacei of 5.8-16.4 MYA previously estimated by Catalan et al. [38]. Moreover, we computed the nonsynonymous substitution rate (Ka) and Ka/Ks ratio and identified 1714 homoeologs whose Ka/Ks values were greater than 1, supporting their accelerated evolution. Interestingly, we found that genes involved in DNA mismatch repair with gene ontology (GO) terms such as proliferating cell nuclear antigen complex and DNA polymerase processivity factor activity (GO:00 43626 and GO:00 30337, respectively, P < 0.01) are significantly enriched in those genes, which suggested that DNA damage tolerance may reflect selection and adaptive divergence between B. distachyon and B. stacei.   3A and B). Then, we classified these expressed homoeolog groups according to additive or nonadditive expression and found that 3994 and 4681 homoeolog groups (23% and 26% of the analyzed homoeolog groups in leaf and root tissues, respectively) were nonadditively expressed in the B. hybridum leaf and root transcriptomes, respectively (the remaining 13 262 and 13 340 homoeolog groups in leaf and root tissues, respectively, were additively expressed; Fig. 3C and D). Moreover, we assessed ELD in the nonadditively expressed homoeolog groups. We found that 611 and 668 homoeolog groups in leaf and root tissues, respectively (3.5% and 3.7% of the analyzed homoeolog groups), in B. hybridum showed ELD, with expression levels similar to those in B. distachyon (ELD-Bd) (Fig. 3C and D). On the other hand, 1263 and 1283 homoeolog groups in leaf and root tissues, respectively (7.3% and 7.1% of analyzed homoeolog groups), in B. hybridum showed ELD of B. stacei (ELD-Bs) ( Fig. 3C and D). These findings indicate that the number of homoeolog groups of ELD-Bs is nearly twice that of ELD-Bd, suggesting that gene expression patterns of B. stacei are preferentially inherited to the B. hybridum transcriptome in leaf and root tissues. We also identified 627 and 850 homoeolog groups in B. hybridum leaf and root tissues, respectively (3.6% and 4.7% of the analyzed homoeolog groups), that showed transgressive expression, i.e., higher expression than in both ancestors. Various GO terms, such as primary metabolic processes, primary cellular processes, stress response, protein transport, localization, and translation, were enriched in the set of homoeolog groups that showed transgressive expression (Additional file 2, Tables S5 and S6), suggesting that the overproduction of transcripts may affect various biological functions through the tetraploidization.

Homoeolog expression bias in B. hybridum
To identify homoeolog expression bias in B. hybridum, we sorted the B. hybridum RNA-seq reads to the 2 homoeologous genomes of B. distachyon (Bd-subgenome) and B. stacei (Bs-subgenome), based on their sequence identities after alignments with both genomes, and quantified the expression levels of each of the Bdand Bs-homoeologs in the B. hybridum transcriptome. Through homoeolog-specific RNA-seq read sorting, we discriminated 75-79% of the RNA-seq reads between the Bd-subgenome origin reads and Bs-subgenome origin reads (Additional file 1, Fig. S5 and S6). We found that 6512 and 6296 homoeolog groups show nonancestral expression patterns in the B. hybridum leaf and root transcriptomes, respectively (38% and 35% of the analyzed homoeolog groups; orange, light blue, and purple boxes in Fig. 3E and F), and the remaining 10 744 and 11 725 retained the ancestral expression patterns (gray and red boxes in Fig. 3E and F). In the homoeolog groups showing ancestral expression patterns, we found that similar proportions of the homoeolog groups retained the expression bias of each of the ancestral transcriptomes (27% of the analyzed homoeolog groups in leaf and root tissues; red boxes in Fig. 3E and F). In the homoeolog groups showing nonancestral expression patterns, we observed bias loss in most cases in the B. hybridum transcriptomes (25% and 23% of the analyzed homoeolog groups in leaf and root tissues, respectively; light blue boxes in Fig. 3E and F). These findings suggest that more than 30% of the homoeologous gene pairs in B. hybridum had changed their expression pattern toward nonancestral expression during speciation after the genomescale gene duplication via allopolyploidization.

Homoeolog-specific transcriptional response to heat stress in B. hybridum
Brachypodium distachyon transcriptomes were not significantly different between normal (22 • C) and heat stress ( Fig. 4A), suggesting that B. stacei and B. hybridum respond more strongly to heat stress than B. distachyon. We also found that after 3 days of stress exposure, 5649 genes and 3725 homoeolog groups showed significantly higher expression in B. stacei and B. hybridum, respectively, than in B. distachyon (false discovery rate [FDR], ≤0.001). We then dissected the expression patterns of the 3725 homoeolog groups from B. hybridum into Bd-and Bshomoeologs and found that 2088 homoeolog groups (56% of the 3725 homoeolog groups) were preferentially expressed by the Bs-homoeologs (Fig. 4B). The distribution of fold-change expression values between Bd-and Bs-homoeologs in all expressed homoeolog groups in B. hybridum under heat stress for 3 days indicated that most of the homoeologs were evenly expressed by both subgenomes (Additional file 1, Fig. S7A), whereas the 3725 homoeolog groups showed abundant expression in the Bshomoeologs (Additional file 1, Fig. S7B). Additionally, of the 2088 Bs-homoeologs, 1791 genes were shared with the 5649 genes that showed higher expression in B. stacei than in B. distachyon; these genes included those specifically involved in metabolic processes as well as cellular response to stress and damage stimulus ( Fig. 4C and D). These results suggest that the functions of B. hybridum genes that were inherited from B. stacei may have contributed to the transcriptional stress response and associated metabolic changes in the B. hybridum transcriptome during the earlier response to heat stress.

Early transcriptional responses of Bs-homoeologs contribute to the maintenance of cellular functions in B. hybridum under long-term heat stress
After 15 days of exposure to heat stress, B. hybridum and B. stacei maintained their transcriptional states similar to those under normal condition, in contrast to the severely damaged cellular system of B. distachyon. When we compared the expression patterns between B. distachyon, B. stacei, and B. hybridum after 15 days of exposure to heat stress, we found that B. stacei and B. hybridum showed similar expression patterns between stress and normal conditions (PCC, 0.94-0.96 and 0.92-0.96, respectively), whereas B. distachyon showed remarkable changes in its transcriptome (PCC, 0.84-0.89), reflecting the severe decrease in its biomass ( Fig. 1B and 5A). Specifically, genes involved in primary metabolism, such as photosynthesis and metabolite and energy generation, were significantly less represented in B. distachyon than in B. stacei and B. hybridum under long-term heat stress ( Fig. 5B and C and Additional file 2, Tables S7 and S8), which likely indicates the physiological sensitivity of its cellular system to heat stress. The number of Bs-homoeologs that were expressed more abundantly than Bd-homoeologs was reduced in homoeolog groups that were highly expressed in B. hybridum than in B. distachyon at 15 days after exposure to heat stress compared to that found after 3 days of heat stress exposure ( Fig. 4B and Additional file 1, Fig. S8), which suggests that Bs-homoeologs are significantly activated at the earlier phase of transcriptional response to heat stress. We also found significantly higher expression of Bs-homoeologs and genes in B. stacei encoding A2-type heat shock transcription factor (HsfA2) as well as putative HsfA2-targeted genes such as heat shock protein 101 (Hsp101) and ascorbate peroxidase 2-like (APX2-like; homologs of Arabidopsis APX1 and rice APX2), which are known as key factors in response to heat [44][45][46] compared to expression of their counterparts in B. distachyon at 3 days after exposure to 32 • C (Additional file 1, Fig. S9). To validate the homoeologspecific gene expression in B. hybridum observed in our RNAseq data, we quantified the polymorphic nucleotides in the cDNA amplicons of B. hybridum using a TaqMan SNP genotyping assay in combination with a digital polymerase chain reaction (PCR) system. The assay discriminates between 2 homoeoalleles of a specific SNP that are labeled with 2 fluorescent dyes (FAM TM and VIC TM (Thermo Fisher Scientific, MA, USA)) and quantifies the expression of each allele by detecting the fluorescence. Specifically, we assessed the quantitative distribution of the decisive SNPs in the 2 homoeolog groups corresponding to the genes annotated in the Bd21 genome, Bradi1g16510 (APX2-like) and Bradi2g49660 (Hsp101). We found that the proportions of SNPs in the cDNA amplicons from the B. hybridum RNAs sampled after 3 days of exposure to 32 • C were comparable to their expression bias observed in the RNA-seq analysis (Additional file 1, Fig. S10 and S11). These results suggest that the earlier response to heat stress that was specific to the Bshomoeologs likely contributed to the maintenance of their cellular homeostasis under the long heat exposure.

Evolutionary nonadditive gene expression in the hybrid grass B. hybridum likely underlies its enhanced invasive behavior
Our transcriptome analysis provides comprehensive evidence for nonadditive gene expression in the hybrid species B. hybridum, suggesting the global transcriptional changes in its leaves and roots evolved through allopolyploidization. The GO analysis of nonadditively expressed homoeolog groups in leaf and root tissues of B. hybridum relative to both ancestral species showed an overrepresentation of genes that were involved in the response to stimulus and abiotic stimulus (Additional file 2, Tables S9 and S10). Such functions were also enriched in the homoeolog groups in leaf tissues of B. hybridum that showed nonadditive expression under heat stress conditions (Additional file 2, Tables S11 and S12), which suggests that increased expression divergence of genes related to such functions enhanced this species' ability to respond to environmental change and to adapt to ecological niches. Increased expression divergence of these classes of genes has also been reported in allotetraploid Arabidopsis and wheat [24][25][26]. Although such examples have been reported in few species, future progress in transcriptome datasets of hybrid species and their ancestors might enable the application of universal rules for determining transcriptional changes when new hybrid species are generated. Allopolyploid species have long been hypothesized to possess greater environmental adaptability to wider niches than their ancestors. Specifically, enhanced heterozygosity and genetic diversity that result from the hybridization of multiple diverged genomes have been thought to upgrade stress tolerances and contribute to the expansion of niches in hybrid species [47][48][49]. Although B. hybridum shows the largest niche overlap compared with its diploid ancestral species, it shows a niche breadth that is smaller than that of B. distachyon and slightly greater than that of B. stacei [37,38]. However, B. hybridum also successfully colonized nonnative regions of the world [37], which suggests it has greater ecological tolerance than other diploids. With the nonadditive transcriptional changes in B. hybridum, the expanded diversity of gene expression might contribute to the colonization of nonnative areas while avoiding inbreeding depression and might boost its diversifying selection [50,51].

Homoeolog-specific gene expression underlies acclimation to heat stress
Although no differences in visible traits were found between the 3 Brachypodium species after 3 days of exposure to heat stress, significant differences in physiological traits in the early stage of heat exposure were revealed from their transcriptomes. The comprehensive list of Bs-homoeolog and B. stacei genes with significantly higher expression than the B. distachyon genes after 3-day exposure to 32 • C included genes involved in the regulation of acclimation to heat in plants, such as heat shock transcription factor, heat shock protein, and DNAJ [52][53][54][55] (1791 genes in Fig. 4C, Additional file 2, Table S13). Thus, the functions of B. hybridum in acclimating to heat are likely specifically inherited from B. stacei via allopolyploidization. Based on our transcriptome analysis using the Brachypodium trio, we hypothesize that the transcriptional response of Bs-homoeologs in B. hybridum and their counterparts in B. stacei during the early phase of heat exposure contributes to the maintenance of their physiological equilibrium under long-term heat exposure and might presumably be associated with their heat stress tolerance (Fig. 6). Previous evolutional and ecological studies suggested that B. distachyon could have adapted to different environments by diverging from B. stacei [37,38]. The distribution areas of B. distachyon and B. stacei suggest that the 2 Brachypodium species have adapted to cooler and wetter areas and to warmer and drier areas, respectively, via diversification [37]. Therefore, the heat acclimation function in B. distachyon might have been lost during its adaptation process after branching from the common ancestor of these diploids. The heat-adaptive trait in the B. stacei genome could influence the survival of

Fresh weight measurement
The aboveground parts of plants of each species grown under normal and heat stress conditions for 3 or 15 days were used to measure fresh weight. Twelve individuals were used for the measurement. A t test was used for statistical comparisons between the plants grown under the different conditions. The significance threshold was set at P < 0.01.

Coverage calculation and SNP calling
To obtain high-quality reads, the genome sequence reads were trimmed by cutting bases off from the start and end of reads if quality ≤20 and by removing the reads with final read length <50 using Trimmomatic v.0.32 (Trimmomatic, RRID:SCR 011848) [56]

Read count and reads per million calculations
To obtain high-quality reads, the RNA-seq reads were trimmed by cutting bases off from the start and end of reads if quality ≤20 and by removing the reads with final read length <50 using Trimmomatic with -thread 4 LEADING: 20 TRAILING: 20 MINLEN: 50 commands. The trimmed reads were mapped to Bd21 and the virtual B. stacei genome using TMAP with mapall -n 4 -v -Y -u -o 2 stage1 map4 commands; B. distachyon reads were mapped to the Bd21 genome; B. stacei reads were mapped to the virtual B. stacei genome; and B. hybridum reads were mapped to both of these genomes. The expression levels of the homoeologs in B. hybridum were quantified by classifying the RNA-seq reads into the following 3 groups based on the number of mismatches between read and both genomes using original Perl script (Additional file 5): Bd-subgenome origin reads, Bs-genome origin reads, and unclassified reads (Additional file 1, Fig. S5B). The mapping data from the B. distachyon reads, B. stacei reads, Bdsubgenome origin reads, and Bs-genome origin reads were subjected to read count using featureCounts (v.1.4.6) [66] with the gene structural annotation of Bd21. For the MAV, which is the average of the expression values in both ancestors, the total read count data from the B. distachyon and B. stacei reads were used; for the entire expression value of the homoeolog group in B. hybridum, the total read count data for the Bd-and the Bssubgenome origin reads were used (Additional file 1, Fig. S5A).
The RPM values were calculated for all annotated genes based on the read count data. Genes, MAVs, homoeologs, and homoeolog groups with RPM ≥1 in all 3 biological replicates were defined as expressed.

Validation of homoeolog-specific gene expression in B. hybridum
The ratios of expression between homoeologs were verified using a digital PCR system in combination with the TaqMan R SNP genotyping assay method (Applied Biosystems, CA, USA) [67]. PCR primers were designed at conserved regions between homoeologs, and a targeted homoelogous SNP was included in the amplified region. Two dye-labeled (FAM and VIC) TaqMan probes were designed for quantifying homoeologous SNPs using the TaqMan MGB SNP kit (Applied Biosystems). cDNA was prepared from B. hybridum RNA sampled after 3-day exposure to heat stress (with 3 biological replicates) using SuperScript R IV reverse transcriptase (Invitrogen, CA, USA) and was used for quantitative PCR. The thermocycling conditions were initiated at 96 • C for 10 min, followed by 39 cycles of annealing and extension at 56 • C for 2 min and denaturation at 98 • C for 30 s, followed by a final extension at 56 • C for 2 min. Homoeologous SNPs in the cDNA amplicons were detected with QuantStudio TM 3D (Life Technologies Japan). The PCR primers and TaqMan probes used in this assay are listed in Additional file 2, Table S14.

Differentially expressed genes analysis
Differentially expressed genes between the 3 Brachypodium species, homoeolog groups, and growth conditions were identified using the DESeq2 package (v.1.10.1) [68] in R (v.3.2.4) with the Wald test based on the read count data. The FDR for each comparison was calculated by adjusting the P value using the Benjamini-Hochberg procedure. Genes with an FDR ≤0.001 were defined as differentially expressed.

GO enrichment analysis
GO annotations of the B. distachyon genes were prepared using the method used in Koda et al. [69]. GO terms for the B. distachyon genes were used from the annotation information retrieved from Phytozome (Bdistachyon 314 v3.1.annotation info.txt) [57]. Additional GO terms were associated with the B. distachyon genes based on GO terms of transcripts for Arabidopsis thaliana and rice in the "Best-hit-arabi-name" and "Besthit-rice-name" row in the annotation file. The GO terms for genes annotated in A. thaliana and rice were used from the annotation information downloaded from Phytozome (Athaliana 167 TAIR10.annotation info.txt [70], and Osativa 204 v7.0.annotation info.txt [71]). To reduce bias, GO terms assigned to more than 3000 genes in B. distachyon were excluded. Enriched GO terms for selected genes were identified using Basic Local Alignment Search Tool 2GO (v.3.3.5) [72] with the Fisher exact test. All annotated genes in the Bd21 genome were used as a reference set. The enriched GO terms were summarized using the web-based tool REVIGO [73] if more than 10 enriched GO terms were found. GO annotations of the B. distachyon genes used in this study are provided in Additional file 6.

Availability of supporting data
All sequencing data were archived at DDBJ under accession numbers DRA005717 and DRA005699 (BioProject accessions: PRJDB5654 and PRJDB5657). Supporting data are available in the GigaScience repository GigaDB [74] and via additional files. Additional files 1 and 2 provide Supplementary Figures and Supplementary Tables, respectively. Additional file 3 provides the