We produced a nascent synthetic allotetraploid wheat (AADD) through hybridization and artificial WGD of two diploid parental species of bread wheat, Triticum urartu (AA), and Aegilops tauschii (DD). To avoid previously reported noise effects caused by dramatic changes in chromosome number and structure in AT2 (Zhang et al., 2013; Gou et al., 2018), all individuals of AADD wheat (accession AT2) were initially karyotyped. Karyotyping was performed by sequential FISH with two DNA probes (pSc119.2, green; pAS1, red) followed by GISH, which enabled the identification of all A- and D-subgenome homoeologous chromosome groups in AT2 individuals (Figure 1A). Individuals without variations in chromosome number and structure (euploids, 2n = 28) were used for subsequent tissue collection and RNA-seq analyses.

We sampled various tissues from AT2 and the diploid AA and DD parents at three developmental stages: 4 days after germination (primary roots and shoots), the trefoil stage (roots and leaves), and the inflorescence stage (inflorescences, 0.5–2.5 cm in length) (Figures 1A,C, Materials and Methods). To obtain a global transcriptomic atlas of all tissues, developmental stages, and genotypes, we performed RNA-seq analysis of 81 samples (three genotypes × nine tissues/stages × three replicates). In total, we generated ~390 Gb of clean data, with an average of 22.9 million reads covering each subgenome (Supplementary Table S1). An in silico parental mix (hereafter denoted Mix) was constructed by mixing the diploid parental RNA-seq data at a ratio of 1:1 to represent for the orthologous gene expression divergence in those two parents. After aligning the sequencing reads and retaining only genes with >1 TPM (transcripts per million reads) in at least one sample, we obtained an expression matrix of 50,194 genes in AT2 and Mix for use in principal component analysis (PCA) to generate an overview of gene expression in the AADD synthetic allotetraploid lineage. The first two principal components (PC1 and PC2) explained 85.76% of the total variance and grouped all samples into four distinct clusters according to tissue type (leaf, shoot, two root tissues, and five inflorescence tissues). Tissue type, therefore, appeared to be the main factor responsible for differences in gene expression between samples. This was consistent with the grouping pattern reported in other studies (Appels et al., 2018) (Figure 1D). To further examine the transcriptome atlas between homoeologous subgenomes, 9,642 homoeologs from the A and D subgenomes expressed in at least one tissue were used to perform hierarchical clustering analysis. Consistent with the PCA results, hierarchical clusters primarily reflected tissue types (Figure 1E). However, within each tissue cluster, subgenomes in the allopolyploid and corresponding parental genomes tended to form sub-clusters (Figure 1E), indicating that interspecific differences and parental legacy were the most important factors that influenced global transcription patterns in the synthetic allopolyploid lineage.

To summarize and compare differences in transcriptomic expression between the parental diploid genomes and their respective allopolyploid subgenomes, we categorized parental diploid orthologs (Ad and Dd) and allopolyploid homoeologs (At and Dt) in terms of their expression conservation across tissues and developmental stages. In brief, genes were categorized based on the ubiquity of their expression: genes expressed in all or most tissues (8–9 tissues/stages) were placed in the “conserved” group, genes with intermediate expression (2–7 tissues/stages) in the “intermediate” group, genes expressed in a single tissue in the “specific” group, and genes with no expression in the “silent” group (Figure 2A). Intriguingly, most orthologs and homoeologs in both the parental diploids and the allotetraploids consistently displayed either a “conserved” or a “silent” expression pattern (34.02–38.83% and 32.55–41.02%, respectively; Figure 2A), giving rise to a bimodal pattern for the number of tissues in which genes were expressed. The “intermediate” and “specific” groups contained 19.43–22.57% and ~6% of all homoeologs (Figure 2A). Major differences in “conserved,” “intermediate,” and “silent” gene numbers between the allotetraploid At and Dt subgenomes were mainly inherited from their diploid parents (Ad vs. Dd; Figure 2A). We found majority of genes (91.21 and 89.63% in A and D subgenomes, respectively) inherited parental types in AT2. Specifically, in A genome, the “intermediate” group showed maximum change (2.86%), whereas the “conserved” group showed minimum change (1.3%); similar to A genome, the “intermediate” group showed maximum change (3.73%), whereas the “silent” group showed minimum change (1.38%) in D genome (Supplementary Figure S1). Such a result suggests that interspecific differences (cis effects) rather than allopolyploidy itself have an important role in regulating transcriptome differences between the two allopolyploid subgenomes. Furthermore, when we compared the expression level group compositions between At and Ad and between Dt and Dd, the At and Ad genomes exhibited almost identical expression patterns (χ²-test, p-value = 0.45; Figure 2A; Supplementary Figure S2); in contrast, the Dt subgenome exhibited greater differences in gene activation and silencing relative to the Dd genome at the whole-tissue scale after allopolyploidization (χ²-test, p-value < 0.001; Figure 2A; Supplementary Figure S2).

Accordingly, we focused on characterizing expression differences in gene pairs from each homoeologous subgenome vs. its parental diploid genome (At vs. Ad and Dt vs. Dd, abbreviated as AtAd and DtDd hereafter). We identified 3,474 (13.44%) and 5,305 (19.41%) differentially expressed genes (DEGs) (>2-fold changes) in at least one tissue between two A and D genomes, respectively. Specifically, in the pair of At vs. Ad (AtAd), there were 135–455 (0.67–2.2%) upregulated genes and 115–1,133 (0.55–5.47%) downregulated DEGs in At; in the pair of Dt vs. Dd (DtDd), there were 138–624 (0.65–3.31%) upregulated genes and 229–1,548 (1.03–6.77%) downregulated genes in Dt (Figure 2B). In particular, the proportion of downregulated DEGs in spike tissues (from 0.5 to 2.5 cm) was much higher in the DtDd comparison than in the AtAd comparison (Fisher's exact test, p-value < 0.001), and this proportion increased gradually with spike development (from 2.65% in 0.5-cm spikes to 6.49% in 2.5-cm spikes, Figure 2B). We next focused on genes that showed the greatest response to allopolyploidization (>4-fold changes). The expression trends observed across all DEGs in paired comparisons were also observed in this group of allopolyploidy-responsive genes (Figure 2B). GO analysis of DEGs identified in DtDd comparison across tissues showed that the DEGs downregulated in spike tissues were enriched in biological process terms related to photosynthesis, the cellulose catabolic process, and the activities of various enzymes (Supplementary Figure S3). The larger percentage of D subgenome genes downregulated in spike tissues may explain, at least in part, why the spike phenotype of AT2 is more similar to that of T. urartu (Figure 1C).

The DEGs identified in paired comparisons (AtAd and DtDd) were also categorized into groups based on the number of tissues in which they were expressed, as described above. We found that 71.48% (1,967) of all DEGs in the AtAd comparison were tissue specific (expressed in a single tissue), and the remaining 28.53% (785) were shared among at least two tissues; these included 26.6% (732) intermediate DEGs (2–7 tissues) and 1.93% (53) conserved DEGs (8–9 tissues) (Figure 2C). The DEG categorization differed markedly in the DtDd comparison: The proportion of tissue-specific DEGs decreased to 50.52% (2,319), and the proportions of intermediate and conserved DEGs increased to 45.4% (2,084) and 4.07% (187), respectively (Figure 2C). To investigate whether the same gene has different expression response (up- or downregulated) to allopolyploidy among different tissues, the gene, which is upregulated (down-regulated) in all checked tissues, was defined as “up-conserved” (“down-conserved”), whereas that with upregulated in some tissues but downregulated in other tissues was defined as “divergent.” For the DEGs expressed in at least two tissues, consistent up- or downregulation in the AtAd (79.95%, 682) and DtDd (88.14%, 2,119) comparisons implied that these genes showed a convergent response to allopolyploidization (in all tissues where they were expressed). The proportion of such DEGs increased gradually as the number of tissues increased in the AtAd comparison but showed little change with tissue number in the DtDd comparison (Figure 2C). The specific up- and downregulation pattern differed between the AtAd and DtDd comparisons. Specifically, there were relatively similar proportions of up- and downregulated DEGs in the AtAd comparison (40.68 vs. 40.56%, χ²-test, p-value = 0.9697), but there were far more downregulated DEGs in the DtDd comparison (16.1% up vs. 72.05% down, χ²-test, p-value < 0.001). A sizable fraction of the DEGs, 18.76% (AA) and 11.86% (DD), showed divergent responses to allopolyploidization; in other words, these DEGs showed a different direction of regulation in at least two tissues, unlike the DEGs that showed a convergent response (Figure 2C). Together, these results show that allopolyploidy can lead to expression differences in given tissues for both subgenomes, but the Dt genome showed greater expression differences than the At genome in response to allotetraploidy.

One of the most general phenomena induced by allopolyploidization is subgenome dominance, which is established by homoeolog expression bias (Grover et al., 2012). Accordingly, we identified 9,642 homoeologous and orthologous gene pairs in AADD and its diploid parents and explored transcriptional subgenome dominance in all allopolyploid samples. To dissect and trace the underlying HEB variation in response to allopolyploidy, we classified homoeolog pairs into four major groups according to their HEB status in the parental mix and AT2 (Materials and Methods; Supplementary Table S3). We found that the majority of homoeolog pairs inherited the parental HEB pattern across tissues (Group i, 90.87–94%), and the proportion of A-biased homoeologs was significantly lower than that of D-biased homoeologs (6.82–8.84% in A subgenome vs. 13.66–14.92% in D subgenome; Fisher's exact test, p-value < 0.001), suggesting that parental legacy is a major determining factor in the transmission of gene expression patterns. Homoeolog pairs that had lost parental expression bias (Group ii, 2.13–5.73%) and pairs with de novo expression bias (Group iii, 2.71–4.02%) were moderately common (Supplementary Table S3). In Group iii, there was a more pronounced de novo expression bias toward the A subgenome (1.39–2.49% for A-HEB, 1.22–1.68% for D-HEB; Supplementary Table S3). Homoeologs with reversed HEB (Group iv) were relatively uncommon and were scattered among different tissues and stages.

Based onthe HEB patterns in the parents and AT2, we calculated the net HEB characterization by comparing the expression of each homoeolog to that of its parental orthologs (Materials and Methods). As shown in Figure 3A, more genes exhibited D-subgenome dominance than A-subgenome dominance in leaf tissues (3.72% D-HEB vs. 1.16% A-HEB). However, other tissues, especially the five spike tissues, showed more A-HEB: 0.99–4.17% A-HEB vs. 0.53–1.17% D-HEB (Figure 3A; χ²-test, p-value < 0.001).

To explore the contribution of allopolyploidy-induced differential gene expression to HEB variation, we divided each of the HEB groups above (except Group iv) into three classes based on whether both, neither, or one of the homoeologs from the A and D subgenomes was DEG (Figure 3B). The latter group of subgenome-specific DEGs was subdivided into A-subgenome DEGs and D-subgenome DEGs. The proportion of subgenome-specific DEGs was lower in Group i (from 1.53% in root tissue to 6.56% in shoot tissue) than in Group ii (from 16.1% in primary root tissue to 40% in 2.-cm spike tissue) and Group iii (from 6.76% in root tissue to 22.93% in 1.5-cm spike tissue), which probably contributed to the significant difference in HEB (Figure 3B; χ²-test, p-value < 0.001). In addition, there were more downregulated homoeologs than upregulated homoeologs in all three groups (Student's T-test, p-value = 0.02 for Group i, p-value = 0.04 for Group ii, and p-value = 0.01 for Group iii). In particular, there were more downregulated D-subgenome homoeologs than A-subgenome homoeologs at all five spike developmental stages (29.01–95.88% vs. 4.12–70.99% in Group i, 11.43–98.17% vs. 1.83–88.57% in Group ii, and 33.33–94.74% vs. 5.26–66.67% in Group iii), suggesting that the overall A-subgenome dominance in AT2 is correlated with dramatic downregulation of gene expression in the Dsubgenome (Figure 3B).

To measure the expression correlation between homoeologs at the network level, we also constructed a homoeolog co-expression network. As shown in Supplementary Figure S4, the homoeologs were divided into three distinct groups (A–C) based on their different expression connectivity scores (abbreviated HEC hereafter) among all tested tissues. The higher the connectivity, the more the homoeologs are assumed to share the same expression profile, whereas lower connectivity indicated greater expression divergence (Materials and Methods). In line with the observed increase in expression divergence, the levels of non-synonymous (Ka) and synonymous substitutions (Ks) and nucleotide divergence in homoeolog promoters gradually increased in the three groups as expression connectivity decreased (leveraged expression divergence) in both Mix and AT2 (Supplementary Figure S4A, ANOVA and Tukey's HSD test, p-value < 0.05). This result indicates that variation in cis-elements and/or amino acid sequences (i.e., cis effects) was likely to underlie the expression divergence between homoeologs.

By characterizing changes in HEC induced by allopolyploidization, we found that most homoeologs inherited the parental expression pattern (55.85%, 28.29%, and 3.94% for the three groups) and maintained relative high HEC scores (84.14% with scores >0.3, Group A + B), implying that the function and expression of most homoeologs remained similar before and after allopolyploidy. For homoeologs that switched between groups, allopolyploidization tended to increase rather than reduce expression relatedness between homoeologs (8.39% from low HEC to high HEC vs. 3.54% from high HEC to low HEC, χ²-test, p-value < 0.001), which implies potential cross-talk between subgenomes via trans factors in the allopolyploid (Supplementary Figure S4B).

To explore cis-trans regulatory mechanisms at the initial stage of allopolyploidy, we characterized homoeolog-specific expression (abbreviated as HSE hereafter) for 9,642 expressed homoeologs in all tested tissues. Based on a published classification method (Bao et al., 2019), we identified seven categories: conserved, cis-only, trans-only, cis + trans, cis × trans, compensatory, and ambiguous (Materials and Methods). An example summary of gene homoeologs in these categories from root tissue is shown in Figure 4A. The overall numbers of genes in a given category were relatively similar in different tissues and stages (Supplementary Figure S5). Notably, with the exception of the conserved category, cis-only homoeologs were most common in all tested tissues (average, 36.59%, Figure 4B), whereas cis × trans homoeologs were least common (average, 0.36%, Figure 4B). When we focused on the frequency of each category, cis-only homoeologs were present in the highest proportion of tissues, but the median proportion for cis-only homoeologs was still <50% (0.44). All these results imply that HSE regulation among different tissues and stages of synthetic AT2 wheat was extremely variable and complex (Figure 4C).

To further explore allopolyploid-induced variations in HEC, we tested whether or not homoeologs had undergone a transition between different modules based on the homoeolog co-expression network described above. For each quadruplet (Ad–Dd ortholog pairs in diploid parents and At–Dt homoeolog pairs in AT2), we defined different rewiring groups based on whether the four members belonged to the same or different co-expression modules. We classified the 9,642 quadruplets into five distinct groups: (i) conserved same: Ad, Dd, At, and Dt belonged to the same module; (ii) conserved different: Ad and Dd werein different modules, but Ad and At were in the same module, as were Dd and Dt; (iii) convergence: Ad and Dd were in different modules, but At and Dt were in the same module; (iv) divergence: Ad and Dd were in the same module, whereas At and Dt were in different modules, and (v) others: any pattern that did not match one of modules i–iv (Figure 4D).

As expected, the majority of quadruplets (4,186, 55.08%) showed a conserved expression pattern in which all four members belonged to the same co-expression module (Group i). Group ii contained 1,409 (18.54%) quadruplets, and we inferred that such a parent-mimic pattern (in which the subgenomes differed and matched their respective parents) may reflect divergence of cis elements in the parents after species differentiation. Indeed, when HSE patterns were considered, Group ii possessed the highest proportion of cis-only homoeologs (Figure 4E; Games-Howell post-hoc test, p-value < 0.001). Furthermore, Group ii had the highest Ka and Ks values and nucleotide diversity between homoeologous promoters, providing further evidence for the role of cis regulation (Figure 4F; Games-Howell post-hoc test, p-value < 0.05). Notably, the number of quadruplets for which only the parent orthologs were in different modules (Group iii; 1,112, 14.63%) was significantly higher than the number for which only the tetraploid homoeologs were in different modules (Group iv; 656, 8.63%) (Figure 4D; χ²-test, p-value < 0.001). The former result may be related to the allopolyploid environment with a shared trans regulator, whereas the latter appears to be associated with subgenome dominance (asymmetry).

A notable phenotypic feature of AT2 is its general morphological mimicry of T. urartu rather than Ae. tauschii during spike development, especially for the spikelet shape and awn length: Relative to T. urartu and AT2 having thin and long spikelet, the spikelet of Ae. tauschii is wider and shorter (Figure 1C). We hypothesized that spike development-related genes (SDRGs) that were downregulated in the D subgenome and/or dominant in the A subgenome may have participated, independently or together, in the regulation of spike development after AADD allopolyploidy.

To explore this hypothesis, we constructed a new gene co-expression network of five spike developmental stages based on the expression of 2,561 SDRGs (1,135 and 1,426 from the A and D subgenomes, respectively; see Materials and Methods). We assigned 2,396 (~94%) SDRGs to nine network modules, and modules 3, 5, and 7 exhibited overall lower expression levels in AT2 (Figure 5A). Notably, there were more D-subgenome genes than A-subgenome genes in these three modules (Figure 5B; 259 vs. 100 in Module 3, 128 vs. 3 in Module 5, and 59 vs. 38 in Module 7; binomial test, p-value < 0.05 for all comparisons). Furthermore, the proportion of stably downregulated DEGs from the D subgenome was further enriched in Modules 3 (from 72.14 to 93.07%), 5 (from 97.71 to 100%), and 7 (from 60.82 to 70.59%) (Figures 5B,C; Materials and Methods). The expression patterns of 108 stably downregulated SDRGs in the D subgenome are shown in Figure 5D. They include a number of genes associated with spike development that have been reported previously. For example, FLOWERING LOCUS T, a major gene that regulates wheat flowering and yield-related traits (Chen et al., 2014; Liu et al., 2020; Isham et al., 2021), showed increased expression throughout early wheat spike development (Figure 5D). In addition, transcription factor (TF) genes, such as MYB and bHLH TF genes, exhibited increased expression from early to late development, suggesting that these TFs play important roles in wheat spike development. Twenty-four of 108 SDRGs in the D subgenome that had counterparts in the A subgenome showed bias toward the D genome in Mix in at least one stage but had similar expression levels in AT2 (Supplementary Figure S6). These results imply that downregulation of D-subgenome homoeologs may contribute to the similar spike developmental phenotypes of T. urartu and AT2.

On the other hand, the overall subgenome dominance of important A-subgenome SDRGs during spike development may also contribute to the phenotypic similarity of AT2 and T. urartu spikes. We further investigated 51 stable A-subgenome HEBs in AT2 (Materials and Methods) and found that the majority of them were intrinsic A-subgenome HEBs (Figure 5E). Intriguingly, spike development-associated transcription factors, including B3, ARF, and MYB factors (Wang et al., 2017; Li et al., 2018), were included in these A-subgenome HEBs, indicating that these TFs tended to be asymmetrically expressed during specific spike developmental stages in the allopolyploid environment. Ubiquitin, F-box protein, and E3 ubiquitin-protein ligase also occurred frequently among the A-subgenome HEB genes, indicating that ubiquitin-related protein processes may play important roles in spike morphology. Moreover, a gene encoding an auxin response factor exhibited A-subgenome HEB (Figure 5E), and this gene has specific roles in the reproductive development of rice, maize, and wheat (Benjamins and Scheres, 2008; Li et al., 2018). These results suggest that some key candidate SDRGs with A-subgenome dominance may also affect spike development in artificially synthesized AT2 wheat.

The plants used in this study included the newly synthesized allotetraploid wheat (AADD genome, accession AT2) along with its diploid parents, Triticum urartu (AA genome, accession TMU38), and Aegilops tauschii (DD genome, accession TQ27). Fluorescent in situ hybridization (FISH) and genomic in situ hybridization (GISH) were conjoined for karyotyping investigation, and only euploid plants were used for further experiments. All the plants were grown under the same controlled growing conditions: 16 h 25°C/8 h 15°C, day/night. We collected samples of different tissues and developmental stages from both AT2 and its diploid progenitors for RNA-seq, which include: shoots and primary roots were collected 4 days after germination; second leaves and mature roots were gathered from trefoil-stage seedlings (Figure 1B); spikes of five successive periods from 0.5 to 2.5 cm in length during stem extension were also obtained (Figure 1C). For each tissue, each of the three biological replicates was pooled from samples collected from three individuals. All collected tissues were immediately frozen in liquid nitrogen. Total RNA was isolated using Trizol (Invitrogen) according to the standard manufacturer's protocol.

Totally, 81 samples (three genotypes × nine tissues/stages × three replicates) were used for RNA-seq profiling. Library construction and sequencing procedure were performed by standard illumine protocols, and c.2.6 billion clean reads (2 × 150 bp) were generated from the Illumina HiSeq2500 (Supplementary Table S1). To estimate and distinguish the accumulated divergence of gene expression (i) between two diploid parental species and (ii) between parents and AT2, an in silico parental mix (hereafter denoted Mix) was constructed by combining the reads of T. urartu and Ae. tauschii at a ratio of 1:1 for all tissues. For each sample, cleaned RNA-seq reads were mapped against the reference genome constructed by integrated or combined T.urartu (http://www.mbkbase.org/Tu/) (Ling et al., 2018) and A.tauschii (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA341983) (Luo et al., 2017) genome sequences using HISAT2 (v.2.0.1- beta) (Kim et al., 2015) with default parameters (Supplementary Table S1). Only uniquely mapped paired-end reads were retained for read counting by feature Counts function of the subread package (v1.6.1) (Liao et al., 2014) to generate the count and normalized expression in TPM (transcripts per million reads) values. Only the genes with TPM > 1 among three biological replicates were used for subsequent analysis. To assess a cross mapping rate (a mis-mapping rate) between the two subgenomes for all test tissues, we completed cross mapping validation by aligning both T. urartu and Ae. tauschii sequencing data to the combined AADD reference genome. We found limited cross mapping rates, which included ~4.9% (4.26–7.51%) of A genome reads mapped to D genome reference and ~0.72% (0.59–1.02%) of D genome reads mapped to A genome reference (Supplementary Table S2). PCA analysis and hierarchical clustering were performed using the prcompfunction and factoextrapackage (http://www.sthda.com/english/rpkgs/factoextra) (Kassambara and Mundt, 2017) in R software (http://www.R-project.org/) (R Core Team, 2013) with default settings.

We examined DEGs (differentially expressed genes) between AT2 and an in silico mix using the DESeq2 (Love et al., 2014) package in R software (R Core Team, 2013). Pairwise comparisons of gene expression were made between expression of Mix and AT2 homologs (At vs. Ad; Dt vs. Dd) at the same developmental stage. Genes with fold change of expression value >2 and FDR adjusted p-value < 0.05 were defined as DEGs in each comparison. The DEGs genes in respective genome pair were categorized according to whether the DEG pattern is observed in only one tissue (tissue specific) in 2–7 tissues (intermediate frequency), or more than 8 tissues (constitutive). To identify spike development-related genes (SDRGs), we initially performed pairwise comparisons of gene expression in spike tissues and non-spike tissues to identify genes specifically upregulated in spike tissues in Mix. Specifically, the genes with at least 2-fold upregulation in at least one spike developmental stage were defined as SDRGs. In addition, the genes downregulated in at least three spike stages compared with parents were defined as stably downregulated SDRGs. GO (gene ontology) enrichment analysis of DEGs was performed using the R package clusterProfiler (v3.10.0) (Yu et al., 2012) with the following parameters settings: pvalueCutoff = 0.05, pAdjustMethod = BH, minGSSize = 15, maxGSSize = 500, qvalueCutoff = 0.05.

To study the features of homoeolog expression patterns in the allotetraploid AT2, the expression levels of 9,642 T.urartu-A.tauschii orthologous gene pairs identified by a python script were monitored. First, CDSs of the T.urartu and A.tauschii genome were compared against each other by BLASTN (E-value < 1e⁻⁵). Then, CDS pairs that met following two criteria were considered to be orthologs, which included: (i) aligned regions cover more than 60% of the CDS length for both paired sequences and (ii) aligned regions harbor minimum 90% sequence similarity. The genes with low RNA-seq reads count (both sequences in a homoeologous pair with relative abundance TPM < 1 in all tissues) were discarded.

For each pair of A- and D-subgenome homoeologs, differential homoeolog expression was calculated using DESeq2 to infer homoeolog expression bias (relative contribution of homoeologs to the transcriptome). The homoeologous gene pairs were initially defined as A-biased homoeologous (A-HEB) and D-biased homoeologous (D-HEB) (fold change > 2, FDR adjusted p-value < 0.05) and non-biased homoeologous gene pairs (nHEB). Homoeologous gene pairs were further divided into four groups: (i) the parental legacy group, in which orthologs and respective homoeologs in Mix and AT2 showed identical HEB in direction; (ii) the group displaying allopolyploidy-induced convergent expression, in which HEB occurred in Mix but diminished in AT2; (iii) the group of allopolyploidy-induced divergent expression, in which de novo HEB was established in AT2, and (iv) the group displaying allopolyploidy-induced HEB reversion. In addition, homoeologs exhibiting A-HEB in at least three spike stages were defined as stably A-biased homoeologs. As for net HEB, the approach described by Flagel and Wendel (2010) was used. Briefly, for the 9,621 expressed ortholog/homoeolog pairs, the TPM values of each pair were converted to the ratio of A to D ortholog/homoeolog. These ratios of the allotetraploid were compared to those of parental mix using the Student's T-test method (a threshold of p-value < 0.05 and fold change >2 was used to assess significance).

The WGCNA (Langfelder and Horvath, 2008) package in R (R Core Team, 2013) was used to build individual weighted undirected co-expression networks for respective datasets with the blockwiseModules function. We input three datasets into the network construction pipeline, including:(i) for constructing homoeolog co-expression networks, the TPM RNA-seq expression values of 9,642 homoeologs in Mix and AT2 of 9 tissues/stages were utilized, respectively; (ii) for constructing quadruplets co-expression network, expression values of 9,642 homoeologs of Mix and AT2 were used; (iii) to build spike-specific network, the TPM values of 2,561 SDRGs in five inflorescence developmental stages of Mix and AT2 were employed. Within the pipeline, the soft thresholding power (β) was set to 18, 18, and 9 for network i, ii, and iii based on the scale-free topology criterion (Zhang and Horvath, 2005). Groups of closely connected genes, namely modules, were identified by clustering genes based on the topological overlap matrix and cutting the clustering tree into branches by the cutreeDynamic method with parameters: deepSplit = 2, pamRespectsDendro = FALSE, minModuleSize = 30, mergeCutHeight = 0.25. The gene expression pattern within a single module was summarized into a module eigengene (ME), corresponding to the first principal component, which was considered the most representative gene expression in such a module (Fuller et al., 2007).

As for the network i, we used the homoeolog expression connectivity (HEC) to assess the degree to which homoeolog pairs diverged in expression. Briefly, the HEC between a pair of homoeologs assesses the degree to which two homoeologs share neighbors in the co-expression network. The higher the connectivity score, the more the homoeologs are assumed to share the same expression profile. According to criteria originally proposed in lotus (Shi et al., 2020), we calculated the hypergeometric-value and the connectivity score, and classified the genes into three groups (Supplementary Figure S4A): Group A with connectivity >0.5 and P-value < 0.01, Group B with connectivity 0.5 > × >0.3 and P-value <0.01, Group C with connectivity <0.3 and P-value >0.99.

Following the method of analyzing classical allele-specific expression (ASE) (McManus et al., 2010; Bao et al., 2019), we completed homoeolog-specific expression (HSE) analysis to explore the potential regulatory insulation and/or interactions among homoeolog genes in allopolyploid. We identified cis- and trans-regulatory divergence by using the procedures as described previously (Bao et al., 2019). Specifically, any significant difference between orthologs in Mix was considered as evidence of cis and trans effects co-regulation [represented by A; A = log2 (Ad/Dd)], whereas any significant difference between homoeologs in the AT2 was considered evidence of cis-regulatory divergence [represented by B; B = log2 (At/Dt)]; accordingly, trans effects were derived by subtracting the expression divergences of gene pairs in AT2 from those of Mix (A – B). Student's T-test was performed to test above three comparisons (comparison A, B, and A – B). All 9,642 expressed homoeologous pairs were specifically classified into seven regulatory categories: (i) cis only: A = B, A ≠ 0, B ≠ 0; (ii) trans only: A ≠ B, A ≠ 0, B = 0; (iii) cis + trans: cis and trans effects of one gene were in the same directions, A ≠ B, A > 0, B > 0 or A ≠ B, A < 0, B < 0; (iv) cis × trans: cis and trans effects of one gene were in the opposite directions A ≠ B, A > 0, B < 0 or A ≠ B, A < 0, B > 0; (v) compensatory: A ≠ B, A = 0, B ≠ 0; (vi) conserved: A = B, A = 0, B = 0; (vii) ambiguous: all other patterns of significant tests, with no clear biological interpretation.