Homoeolog expression in polyploid wheat mutants shows limited transcriptional compensation

Whole-genome duplications (WGDs), hybridisations, and small-scale duplications have resulted in an abundance of gene duplicates in plant genomes – on average, 65% of genes in plant genomes have a paralog (Panchy et al., 2016). These duplications can provide opportunities for genetic innovation when mutations cause one of the copies to adopt a novel function (neo-functionalisation; Birchler & Yang, 2022). Alternatively, the presence of gene duplicates can result in functional redundancy, whereby the effect of a loss-of-function mutation in one gene copy is masked by the remaining functional copy. It is debated whether this phenotypic robustness is evolutionarily advantageous or whether selection for lower expression noise results in retention of duplicated genes (Pires & Conant, 2016; Iohannes & Jackson, 2023). Nevertheless, this buffering of deleterious mutations has hindered the functional characterisation of genes (e.g. through gene knockouts) and reduced the variation available for crop breeding (Uauy et al., 2017). Hybridisation and polyploidisation events are widespread in the evolutionary history of the angiosperms and have given rise to major cereal crops such as bread wheat (Triticum aestivum) and pasta wheat (Triticum turgidum ssp. durum; Matsuoka, 2011). Bread wheat is a hexaploid, formed by the merger of three diploid progenitor species. The hexaploid wheat genome therefore consists of three subgenomes (A, B, and D; El Baidouri et al., 2017; IWGSC et al., 2018). Pasta wheat is a tetraploid, and its genome consists of the A and B subgenomes (El Baidouri et al., 2017). As the wheat progenitor species were closely related, many genes are present in highly similar copies across the subgenomes – these duplicate genes formed by polyploidisation are known as homoeologs and share an average nucleotide sequence identity of 97.2% (Schreiber et al., 2012; IWGSC et al., 2018). Similar to paralogs, functional redundancy can exist between homoeologs, and loss-of-function mutations in multiple homoeologs may be required before a phenotype is observed. For example, in hexaploid wheat, loss-of-function mutations in all three homoeologs of the Ms26 gene are required to confer male sterility, while single mutants have no reduction in fertility (Singh et al., 2017). Similarly, only triple mutants in Qsd1 have extended seed dormancy (Abe et al., 2019). This redundancy may result in hidden variation in single homoeologs that cannot be observed until mutations in multiple homoeologs are combined (Borrill et al., 2015). However, it remains unclear how phenotypic compensation between gene duplicates occurs. Active transcriptional compensation, the transcriptional upregulation of genes with a high degree of sequence similarity (such as paralogs or homoeologs) in response to a loss-of-function mutation in a gene, has been proposed as an explanation for functional redundancy (Sztal & Stainier, 2020). Such transcriptional compensation has been observed in multiple species, including Caenorhabditis elegans (Serobyan et al., 2020) and Danio rerio (zebrafish; El-Brolosy et al., 2019; Ma et al., 2019). The degree of transcriptional compensation between paralogs is also hypothesised to underlie the penetrance of mutations – for example, in zebrafish, a greater degree of craniofacial distortion was observed due to mutation of MEF2CA in the absence of upregulation of the MEF2C paralogs (Bailon-Zambrano et al., 2022). Transcriptional compensation between paralogs has also been documented in numerous plant species, for example in the CLE gene family. Transcriptional upregulation of CLE9 (a CLV3 paralog) is observed in clv3 mutants in Solanum lycopersicum (tomato), Petunia hybrida (petunia), and Physalis grisea (groundcherry; Rodriguez-Leal et al., 2019; Kwon et al., 2022). As such, clv3 mutants in these species have weaker phenotypes than the clv3 mutant in Nicotiana benthamiana (tobacco), in which CLE9 has become pseudogenised (Kwon et al., 2022). However, active compensation is not always observed between paralogs in plants – for example, in Arabidopsis thaliana, there was no transcriptional upregulation of RPL23aA in response to the knock-out of its paralog RPL23aB, or vice versa (Degenhardt & Bonham-Smith, 2008; W. Xiong et al., 2020). It is unknown whether active transcriptional compensation occurs in wheat or other polyploid plants, buffering the effects of mutations in individual homoeologs. Here, we assessed whether transcriptional compensation is prevalent between homoeologs in hexaploid and tetraploid wheat using ethyl methanesulfonate (EMS) mutagenised Targeting Induced Local Lesions in Genomes (TILLING) lines (Krasileva et al., 2017). Each wheat TILLING line has a large number of mutations (> 5000 in the hexaploid cultivar Cadenza; Krasileva et al., 2017), allowing us to simultaneously screen the effect of many mutations. We performed RNA-sequencing and differential gene expression analysis to determine whether genes are frequently upregulated in response to a loss-of-function mutation in one of their homoeologs. We found no evidence for widespread active transcriptional compensation between homoeologs, with a rate of c. 3% in hexaploid and tetraploid wheat, indicating that such a mechanism is unlikely to be the primary cause of functional redundancy between homoeologs in polyploid wheat. To test whether active transcriptional compensation occurs between homoeologs in hexaploid wheat, we screened four representative EMS-mutagenised TILLING lines (cv Cadenza) to identify homoeolog groups with a premature termination codon (PTC) mutation in one of the homoeologs. PTC mutations are expected to cause truncation of the protein encoded by the gene, thus resulting in loss of function that could be associated with upregulation of the gene's homoeologs to provide a buffering effect. Each line underwent two generations of single seed descent (SSD) to increase homozygosity. Differential expression of the gene affected by the PTC mutation and its homoeologs in the SSD TILLING line relative to wild-type (WT) was determined by analysis of RNA-sequencing data (Fig. 1a). Across the four mutagenised SSD lines (C0604, C0895, C1015, and C1704), we identified 158 unique homoeolog groups with a PTC mutation in one homoeolog. In 20.6% (C0604) to 41.7% (C0895) of cases, the homoeolog with the PTC mutation was downregulated in the mutagenised line relative to WT (false discovery rate (FDR) adjusted P-value < 0.05; Fig. 1b). However, despite our relaxed threshold for detection, we did not observe widespread upregulation of the homoeologous genes – across the mutagenised lines, we identified only four homoeolog groups (2.5%) in which at least one of the homoeologs was upregulated (FDR-adjusted P-value < 0.05; Fig. 1b; Supporting Information Table S1). The PTC mutations associated with homoeologous upregulation affected groups with various homoeolog expression patterns in WT, suggesting that this effect is not associated with the loss of dominantly expressed homoeologs (Table S1). Furthermore, the proportion of homoeolog groups with upregulated homoeologs was not significantly different for homoeolog groups affected by PTC mutations or synonymous mutations (2.5% of homoeolog groups with a PTC mutation, 3.8% of homoeolog groups with a synonymous mutation; chi-squared test P-value = 0.401; Fig. 1c). Active transcriptional compensation is not expected to occur in response to synonymous mutations as they are not expected to affect protein function. This suggests that active transcriptional compensation is not the default mechanism for buffering single homoeolog mutations. We found a similar lack of active transcriptional compensation in an independent RNA-sequencing dataset of two EMS-mutagenised lines produced by H. C. Xiong et al. (2020). Thirty homoeolog groups with a PTC mutation were identified, of which only one (3.3%) showed upregulation of the homoeologs (Fig. S1a; Table S1). Again, there was no significant difference between the proportion of upregulated homoeologs in groups affected by a PTC mutation compared with those affected by a synonymous mutation (3.5% of homoeolog groups with a synonymous mutation, chi squared test P-value = 0.951; Fig. S1b). Next, we investigated whether the location of the PTC mutation within the transcript affects whether active transcriptional compensation occurs. We hypothesised that PTC mutations occurring earlier within the coding sequence would have a greater impact on protein function and therefore promote more compensatory upregulation. We identified genes with multiple PTC mutations across the four Cadenza mutagenised lines, differing in their position within the gene. Accordingly, we found TraesCS5A02G119300 (with an early PTC in C1015 and a late PTC in C1704) and TraesCS6A02G155200 (with an early PTC in C0895 and a late PTC in C0604; Fig. 1d,e; Table 1). The relative expression of TraesCS5A02G119300, TraesCS6A02G155200, and their homoeologs in the mutagenised lines compared with WT was determined by reverse transcription polymerase chain reaction (RT-qPCR). We found no significant effect of genotype on expression of any of the homoeologs (ANOVA, P-value > 0.05; Fig. 1d,e). This was further confirmed by the RNA-sequencing results (FDR-adjusted P-value > 0.05; Fig. S2a,b). Moreover, in those cases where homoeologous upregulation was observed, the PTC mutations were distributed throughout the coding sequence (Table S1). Overall, we found no evidence for widespread active transcriptional compensation between homoeologs in hexaploid wheat. Next, we investigated whether active transcriptional compensation takes place in tetraploid wheat. As each homoeolog group consists of only two homoeologs (rather than three), a loss-of-function mutation in one homoeolog would result in a 50% reduction in functional transcript, compared with 33% in hexaploid wheat. Thus, we hypothesised that active transcriptional compensation would be more likely to occur in tetraploid wheat. From six EMS-mutagenised TILLING lines (cv Kronos; K2619, K2864, K3239, K0427, K4533, and K0774), a total of 101 unique homoeolog groups containing a PTC mutation were identified. In line with the results from hexaploid wheat, between 8.3% (in K0427) and 32.4% (in K2864) of the homoeologs with a PTC mutation were downregulated in the mutagenised lines compared with WT (FDR-adjusted P-value < 0.05; Fig. 2a). We found three homoeolog groups (3.0%) affected by a PTC mutation in which the nonmutated homoeolog was upregulated (FDR-adjusted P-value < 0.05; Fig. 2a; Table S1). The proportion of homoeolog groups with upregulation of the nonmutated homoeolog was not significantly different between those affected by a PTC mutation compared with those affected by a synonymous mutation (1.5% of homoeolog groups with a synonymous mutation, chi-squared test P-value = 0.234; Fig. 2b). Despite the lack of widespread active transcriptional compensation between homoeologs in hexaploid and tetraploid wheat, we hypothesised that active transcriptional compensation may be responsible for the functional redundancy observed within particular homoeolog groups. The PHS1 homoeolog group, encoding the plastidial α-glucan phosphorylase, has been characterised as functionally redundant by Kamble et al. (2023) – the phs1-1 double mutant in tetraploid wheat has altered starch granule morphology (it has a decreased proportion of small starch granules and larger B-type granules), while neither the phs1-A1 single mutant nor the phs1-B1 single mutant has altered granule morphology. By contrast, the CSP41a homoeolog group, encoding a chloroplast RNA-binding protein, is nonredundant in tetraploid wheat – mutation of the A homoeolog alone (in the csp41-A single mutant) is sufficient to increase resistance to yellow rust (Corredor-Moreno et al., 2022). We used mutants in each of the homoeologs of PHS1 and CSP41a to investigate whether the level of redundancy between homoeologs affected the degree of active transcriptional compensation observed (Table 2). For both homoeologs of PHS1, the presence of a PTC mutation was associated with downregulation of the affected homoeolog (for PHS1-A1, PTC mutation in K4533, fold change = 0.18, FDR-adjusted P-value = 7.8 × 10−19; for PHS1-B1, PTC mutation in K2864, fold change = 0.15, FDR-adjusted P-value = 4.3 × 10−7; Fig. 2c). However, there was no upregulation of the nonmutated homoeolog in either case (FDR-adjusted P-value > 0.05; Fig. 2c). Similarly, for CSP41a, the presence of a PTC mutation resulted in downregulation of the affected homoeolog (for CSP41a-A, PTC mutation in K3239, fold change = 0.045, FDR-adjusted P-value = 1.9 × 10−61; for CSP41a-B, PTC mutation in K2619, fold change = 0.093, FDR-adjusted P-value = 7.6 × 10−15; Fig. 2d). No upregulation of the nonmutated homoeolog was observed in either case (FDR-adjusted P-value > 0.05; Fig. 2d). This indicates that active transcriptional compensation is not necessary for functional redundancy between homoeologs. We have shown, for 188 homoeolog groups in hexaploid across two independent experiments and 101 homoeolog groups in tetraploid wheat, that active transcriptional compensation occurs at a low frequency – only 2.5–3.3% of PTC mutations were associated with upregulation of the homoeologous genes across three independent datasets. These proportions were not significantly different from the proportion of synonymous mutations associated with homoeolog upregulation, suggesting that apparent compensatory upregulation of gene expression is not necessarily linked to the functional impact of the mutation (Figs 1, 2, S1). Mechanisms to compensate for paralog loss-of-function have previously been shown to be nonuniversal – a similar study in yeast (Saccharomyces cerevisiae) found increased levels of c. 11% of proteins when their paralogs were deleted (DeLuna et al., 2010). DeLuna et al. (2010) also showed that compensatory upregulation occurred almost exclusively for proteins with essential functions, suggesting that upregulation is induced in response to a physiological deficit, rather than as a direct response to the mutation. Such ‘needs-based’ responses could explain the upregulation observed in the homoeolog group encoding a subunit of Respiratory Complex I (Table S1; NADH dehydrogenase). However, this did not extend to any transcription factors, which are generally considered to be dosage-sensitive (Birchler & Veitia, 2010), and might therefore be expected to be transcriptionally compensated to maintain gene–dosage balance. The low level of active transcriptional compensation observed in wheat is not in keeping with the functional redundancy that is often observed between homoeologs. Although almost all of the mutagenised genes in our study are functionally uncharacterised, the PHS1 homoeologs in tetraploid wheat, previously shown to be functionally redundant, also did not show active transcriptional compensation (Fig. 2c; Kamble et al., 2023). Future experiments to characterise the few genes that showed transcriptional compensation may provide additional insight into the connection to functional redundancy. It is possible that a different active compensatory mechanism facilitates functional redundancy between homoeologs – for example, Diss et al. (2014) propose translational upregulation or changes in protein localisation as modes of compensation. Our study is also unlikely to capture cell-type or tissue-specific transcriptional compensation, as highlighted by Iohannes & Jackson (2023). However, current models of transcriptional compensation in mouse, zebrafish and C. elegans indicate that expression of the PTC-containing transcript can cell-autonomously trigger compensation and paralog compensation is shared between different cell states, suggesting that tissue-specificity may not be a major limitation (Sztal et al., 2018; Serobyan et al., 2020; Mellis et al., 2024). Alternatively, there may be significant passive compensation between homoeologs, with two-thirds (in hexaploid wheat) or half (in tetraploid wheat) of the WT levels of functional transcript being sufficient to maintain the WT phenotype. Passive compensation between homoeologs in wheat is probable given that these gene duplications are relatively young and have had less time to diverge in expression levels – the hybridisation events occurred between 500 000 (forming A and B subgenomes) and 10 000 yr ago (addition of the D subgenome; El Baidouri et al., 2017). Increased time since gene duplication is likely to result in hypofunctionalisation, as the paralogs acquire mutations potentially decreasing their expression, and compensatory drift, in which the expression of one paralog decreases over time and function is maintained by its highly expressed copy (Iohannes & Jackson, 2023). When hypofunctionalisation or compensatory drift have occurred, expression from one paralog is insufficient to maintain gene function and an active compensatory mechanism would be required for functional redundancy (Iohannes & Jackson, 2023). Accordingly, Cusack et al. (2021) found that the age of duplication is associated with functional redundancy in A. thaliana, with duplicates formed during the α-WGD more likely to be functionally redundant than those that arose during the more ancient β- or γ-WGD events. Furthermore, the HBEGF paralogs exhibiting active transcriptional compensation in zebrafish are much older than homoeologs in wheat, having arisen during the fish-specific genome duplication event c. 320 million years ago (Ma; Vandepoele et al., 2004; Laisney et al., 2010; El-Brolosy et al., 2019), and are thus more likely to have experienced hypofunctionalisation/compensatory drift. Similarly, the CLV3 paralog in the Solanaceae, CLE9, arose at least 30 Ma and shows active transcriptional compensation (Rodriguez-Leal et al., 2019). Overall, this underscores active transcriptional compensation as an evolutionary innovation that is only favourable when passive compensation is insufficient (whether to maintain functional redundancy or to adequately control gene expression noise, per Pires & Conant, 2016), and is otherwise unfavourable due to the cost of excess mRNA synthesis (with excess mRNA synthesis being selected against, as shown by Hausser et al., 2019). Passive compensation between homoeologs suggests that targeting active transcriptional compensation mechanisms characterised in other species is unlikely to break functional redundancy in wheat and would not be a viable strategy for breeding. Rather, mechanisms to simultaneously introduce mutations in multiple/all homoeologs, such as gene editing, or to reduce the expression of multiple/all homoeologs, such as RNAi, are advantageous to tackle functional redundancy. Target genes were selected based on the availability of multiple EMS-mutagenised hexaploid wheat (Triticum aestivum L. cv Cadenza) lines (Krasileva et al., 2017) with PTC mutations either early or late within the coding sequence of the same gene. Further screening of publicly available gene expression data (Borrill et al., 2016; Ramírez-González et al., 2018) was used to select genes with high expression in seedling leaf tissues. Based on these conditions, the genes TraesCS5A02G119300 and TraesCS6A02G155200 were chosen. Selected TILLING lines with PTC mutations in these genes are in Table SSD lines were produced from the Cadenza TILLING seed for the selected lines (C0604, C0895, C1015, and Krasileva et al., 2017). For each a single seed was and in The resulting was in the in to The were for the RNA-sequencing For tetraploid wheat, genes were selected from previously characterised redundant and nonredundant were based on the availability of EMS-mutagenised tetraploid wheat turgidum ssp. cv with PTC mutations in the coding sequence (Krasileva et al., 2017) and high expression in seedling leaf tissues. Selected TILLING lines with PTC mutations in these genes are in Table In to and the TILLING lines and were to for the lower mutation in the TILLING for the selected lines were from the were on at for and to to promote for a further were into cell containing were in at the with and to maintain a of with and No phenotypic or were observed in the mutant lines compared with WT at the seedling For the TILLING lines were to that they were for the PTC mutation of was from seedling leaf from plants the from The used for are shown in Table was using to the were using the leaf were from plants in and at was using the was using the synthesis was using reverse to the was in a using the I with the of for of 10 at at 30 at and for at are shown in Table expression was determined by the for each transcript relative to the The effect of genotype on relative expression was using a and in were from the leaf of in and at For the Cadenza was using the from three for each TILLING For the was using a and with the & and were by were to the IWGSC et al., 2018) and using et al., 2016). For the were to the with from the D average of per of the total were for the Cadenza TILLING lines and an average of per of the total were for the TILLING from were into using et al., 2015). expressed genes between the TILLING lines and WT were identified using et al., a threshold of adjusted P-value < genes were as those with fold change > 1, and downregulated genes as those with fold change < To identify present in the TILLING lines, the were with et al., and to the IWGSC sequence et al., 2018) using et al., and et al., was used for 10 & and effect was using et al., 2016). were to only those to cause a PTC or synonymous mutation and were Further was to identify in homoeologous genes homoeolog groups in et al., and those present in one of the TILLING lines, not WT to for cultivar The of expressed genes were used to identify changes in expression in the genes affected by a mutation and their homoeologs. expression were from per million for each of the homoeolog groups in the WT the in Ramírez-González et al. For the lines, the was to for the tetraploid genome (with expression of for and for The same analysis was on RNA-sequencing data from H. C. Xiong et al. two independent EMS-mutagenised lines and to the WT We for on RNA-sequencing and for during the This was by the through This was by the and through and the and the in was through a the and the and the with from and the Cadenza SSD and performed the the analysis of RNA-sequencing data and the have and the The Cadenza SSD lines and are available from the through the with and The TILLING lines can be from the and from RNA-sequencing have been in the used for the analysis of the RNA-sequencing data can be found on at Fig. of widespread active transcriptional compensation in an independent dataset of EMS-mutagenised hexaploid wheat from H. C. Xiong et al. (2020). Fig. PTC location not affect whether transcriptional compensation occurs between homoeologs. Table groups in which at least one nonmutated homoeolog is upregulated relative to the WT Table of all used in this is not responsible for the or of any Supporting Information by the than be to the The is not responsible for the or of any by the than be to the for the The remains with to in and in any