Subgenomic Stability of Progenitor Genomes During Repeated Allotetraploid Origins of the Same Grass Brachypodium hybridum

Abstract Both homeologous exchanges and homeologous expression bias are generally found in most allopolyploid species. Whether homeologous exchanges and homeologous expression bias differ between repeated allopolyploid speciation events from the same progenitor species remains unknown. Here, we detected a third independent and recent allotetraploid origin for the model grass Brachypodium hybridum. Our homeologous exchange with replacement analyses indicated the absence of significant homeologous exchanges in any of the three types of wild allotetraploids, supporting the integrity of their progenitor subgenomes and the immediate creation of the amphidiploids. Further homeologous expression bias tests did not uncover significant subgenomic dominance in different tissues and conditions of the allotetraploids. This suggests a balanced expression of homeologs under similar or dissimilar ecological conditions in their natural habitats. We observed that the density of transposons around genes was not associated with the initial establishment of subgenome dominance; rather, this feature is inherited from the progenitor genome. We found that drought response genes were highly induced in the two subgenomes, likely contributing to the local adaptation of this species to arid habitats in the third allotetraploid event. These findings provide evidence for the consistency of subgenomic stability of parental genomes across multiple allopolyploidization events that led to the same species at different periods. Our study emphasizes the importance of selecting closely related progenitor species genomes to accurately assess homeologous exchange with replacement in allopolyploids, thereby avoiding the detection of false homeologous exchanges when using less related progenitor species genomes.


Figure S3 .
Figure S3.Comparison of gene structure characteristics and gene family in the newly sequenced B. hybridum Bhyb-ECI reference genome from Evolution Canyon I (Israel) to those in other assemblies (reference genomes of B. stacei Bsta-ECI, Bsta-ABR114, B. distachyon Bdis-Bd21, Bdis-Bd1-1, B. hybridum Bhyb-ABR113, Bhyb-26, and Oryza sativa.)(Gordon et al. 2020; Mu et al. 2023, Scarlett et al 2022).A) gene length; CDS length; exon length; intron length; exon number.B) Pan-genomic comparison of the number of gene families in different levels, the numbers in petal indicated amount of unique families of each (sub-)genome, the circle in inner showed core gene families number in different comparison levels.C) Classification of genes into large gene types in Brachypodium and O. sativa.Color codes for genomes and for gene types are indicated in the respective charts.

Figure S4 .
Figure S4.Inter-and intra-(sub-)genome comparisons of ks density using orthologous nuclear gene pairs between compared S and D genomes and subgenomes.The heatmap shows the number of gene pair that had no synonymous substitutions (Ks = 0, excluded in ks density analysis) between the pairwise compared genomes/subgenomes (chisq.test).

Figure S6 .
Figure S6.Phylogenomic tree based on whole plastome sequences of Brachypodium hybridum, B. stacei and B. distachyon accessions and plastome lengths summary of newly assembled plastomes and other published plastomes.A) Maximum-likelihood tree of B. hybridum individuals from ECI plus available plastomes of other eastern, central and western Mediterranean B. hybridum and progenitor species B. distachyon and B. stacei accessions.Boostrap support values are indicated on branches.B) Plastome lengths of the reference genomes and other accessions of the three annual Brachypodium species; different plastotypes and species are marked by color and shape (see Sancho et al. (2018) and Gordon et al. (2020).

Figure S7 .
Figure S7.Bidimensional PCA plots of observed data and of simulated data for each of the four tested evolutionary scenarios on the alternative origins of Brachypodium hybridum allotetraploids (see fig. 2C and supplementary table S5).A) D genomes/subgenomes.B) S genomes/subgenomes.4,000 simulated data sets per scenario were generated with DIYABC-RF.The observed data fell within the clouds of simulated data.

Figure S8 .
Figure S8.Genome structure analysis based on nuclear syntenic SNP data of B. hybridum (BhS subgenome) and its progenitor species B. stacei.A) Bidimensional Principal Component Analysis (PCA) plots of nuclear syntenic SNPs of S-genomes discriminating the B. hybridum S-plastotype plus B. stacei eastern vs western Mediterranean populations along PCA1 (10.1% accumulation of variance) and the B. hybridum D-plastotype along PCA2 (4.32%).B) Population structure analysis plot of S-genomes for the best (K=3) genomic groups, showing differentiating the eastern B. hybridum populations group, the eastern B. stacei populations group, and western B. stacei + B. hybridum populations group.C) The nucleotide diversity (π) and pair-wised population divergence (FST) were examined across the three main populations groups.The value within each circle represents the calculated π for that particular group, while the value along each line indicates the FST between two adjacent groups.The distribution of π and FST values are also shown in subplots.**p-value < 0.01 (Wilcoxon test).

Figure S9 .
Figure S9.Comparisons of genomic features between the Brachypodium hybridum Bhyb-ECID and Bhyb-ECIS subgenomes and those of other allotetraploids (Bhyb-ABR113, Bhyb-26) and the genomes of the progenitor species B. distachyon (Bdis-Bd1-1, BdisBd21) and B. stacei (Bsta-ECI, Bsta-ABR114).A) GC content.B) Repetitive content coverage nearby gene regions; down indicates gene downstream 0-10k sequence, up gene upstream 0-10k sequence, and gene the encoding region containing CDSs and introns.The dashed lines link homologous chromosome pairs between the subgenomes and the respective progenitor genomes.C) Composition of repetitive sequences per main repeat categories.The coordinates represent the ratio of each category to the total length of the corresponding genome or subgenome.D) Estimated insertion times of Copia and Gypsy retrotransposons in the B. distachyon and B. stacei genomes and BhD and BhS subgenomes of B. hybridum.Color and symbol codes for B. distachyon (BdisBd1-1, BdisBd21) and B. stacei (Bsta-ECI, Bsta-ABR114) genomes, and B. hybridum (Bhyb-ECI, Bhyb-ABR113, Bhyb-26) BhD and BhS subgenomes are indicated in the corresponding charts of each subfigure.

Figure S13 .
Figure S13.The relationship between gene expression level in transcripts per million (TPM) and transposable element (TE) ratio nearby genes of the Bhyb-ECI genome from Evolution Canyon I. TE ratio of nearby genes was calculated in 2kb upstream to 2kb downstream windows of the gene.A) Gene expression is negatively correlated with TE density near the genes in Bhyb-ECI, 1000 genes were randomly selected to be visualized.B) Histogram plot of TE ratio in non-expressed genes (TPM = 0).Comparisons of TE density in genes show stability in all case study groups.

Figure S15 .
Figure S15.Heatmap of differential gene expression levels of stress-related NAC genes in the Brachypodium hybridum Bhyb-ECI BhD and BhS subgenomes.L, leaf tissue, R, root tissue, C control (well-watered) conditions, T treatment (drought) conditions.The chart legend is Log10(TPM).
B)Genome assembly parameters of B. hybridum Bhyb-ECI local reference genome.

Table S3 .
Gene predictions and statistics in the assembled B. hybridum Bhyb-ECI genome from Evolution Canyon I, functional annotations and TF predictions.Gene predictions comparisons of the B. hybridum Bhyb-ECI genome to those of the reference genomes of B. stacei Bsta-ABR114, B. distachyon Bdis-Bd21, Bdis- A)

Table S5 .
Inputs and results of tested evolutionary scenarios for alternative origins of B. hybridum using filtered SCOG SNP data and coalescence and supervised machine learning methods implemented in DIYABC Random Forest v1.1.1.Scenarios for the separate Brachypodium D and S (sub)-genomes datasets are described in the main text and shown in fig 2C.(A)Values of demographic and historical parameters used for each of the four scenarios of each D and S (sub)-genomic datasets.t1-t4, time; ra, admixture rate; N, effective population size.All parameters were set to default priors; each prior value was drawn from a uniform distribution.
(B) Scenario choice inferred by DIYABC-RF.Mean number of votes and posterior probability of each scenario using out-of-bag estimators for each scenario with DIYABC-RF.Values of the best scenario (scenario 1) are highlighted in bold.(sub)genome

Table S6 .
Summary of repetitive elements found in the assembled B. hybridum Bhyb-ECIS and Bhyb-ECID subgenomes and those of other allotetraploidswith respect to those of progenitor species genomes B. stacei (Bsta-ECI, Bsta-ABR114) and B. distachyon (Bdis-Bd1-1, Bdis-Bd21).The first column of each genome shows the length(bp) of genome, and the second shows the percentage (%) of genome covered by repeat type.

Table S8 .
Summary of sequence comparisons between the B. hybridum BhS and BhD subgenomes and diploids' genomes at gene level.Ks=0 and identical gene number indicated that two homeologous genes have no synonymous substitutions and have the same base sequence, respectively.

Table S9 .
Summary of reciprocal homoeologous exchange (HE) with replacement regions detected between the two homeologous subgenomes (B.stacei-type BhS; B. distachyon-type BhD) of A) the newly assembled Brachypodium hybridum Bhyb-ECI genome and in other 21 resequenced eastern Mediterranean B. hybridum allotetraploid genomes, B) the B. hybridum ABR113 reference genome and in other 5 western Mediterranean B. hybridum allotetraploid genomes, C) the B. hybridum Bhyb-26 reference ancestral genome, D) in a total of 30 B. hybridum eastern and western Mediterranean and ancestral reference and resequenced genomes.

Table 11 .
Summary of Homeologous Expression Bias (HEB) of Brachypodium hybridum BhD and BhS genes.A) HEB of Bhyb-ECI genes in leaf and root tissues under control (well-watered) and drought conditions.B) HEB of constitutionally expressed genes (control conditions) in Bhyb-26.C) HEB of constitutionally expressed genes (control conditions) in Bhyb-ABR113.Median calculated from all Log2(TPMBhS-gene/TPMBhD-gene) values of BhS-gene and BhD-gene homeologous gene pairs.TPM: transcripts per kilobase per million mapped reads.

Table S12 .
Summary of the number of biased expressions of dominant Brachypodium hybridum genes from homeologous BhD and BhS gene pairs in

Table S13 .
Number of biased expressions of dominant Brachypodium hybridumBhyb-ECI genes from homeologous BhD and BhS gene pairs in different tissues (leaf, root) and different conditions (well-watered, drought).