From reciprocal crosses between inbred lines of B. rapa L. (AA, 2n = 20, genotype 3H120) and B. oleracea L. (CC, 2n = 18, genotype Chijielan), B. napus F₁ hybrid (AC, 2n = 19), allotetraploid (CCAA, 2n = 38) and triploid hybrid (CCA, 2n = 28) were produced (Cui et al., 2012) and used for this study, together with another triploid (AAC, 2n = 29) synthesized here. The allotetraploid CCAA originated directly from the cultured embryo-plantlet obtained from the B. oleracea × B. rapa cross without colchicine treatment, probably the spontaneous chromosome doubling occurred in vitro (Cui et al., 2012). The triploid hybrid CCA also from B. oleracea × B. rapa cross likely resulted from the fusion of unreduced gamete (CC) by the female parent B. oleracea and reduced gamete (A) by the male parent (Cui et al., 2012). The allotriploid (AAC) was produced by pollinating the clonal plants of the allotetraploid (AACC) with B. rapa, with the aid of immature embryo culture on MS agar medium without hormones (Murashige and Skoog, 1962; Figure 1A). All these materials were maintained and propagated by subculturing the young buds on MS medium with 1.5 mg/l⁻¹ 6-benzyl aminopurine (6-BA) and 0.25 mg/l⁻¹ α-naphthalenacetic acid (NAA), to produce enough plants for study, following the previous procedure (Cui et al., 2012). Plantlets grew on MS medium in the growth chamber at 25°C and a 14/10 h (day/night) photoperiod, and the newly emerged and expanded young leaves were collected and immediately frozen in liquid nitrogen for RNA extraction. We adopted a mixed sampling strategy with three plantlets of each material and two biological replicates for each sample.

Young ovaries were collected and treated with 2 mM 8-hydroxyquinoline for 3–4 h at 22°C, and subsequently fixed in Carnoy's solution I (3:1 ethanol: glacial acetic acid, v/v) for 24 h, and stored at −20°C. The C-genome specific probe (BAC BoB014O06) was used to identify the C-genome chromosomes. The procedures of FISH analyses followed the protocol of Cui et al. (2012) and Zhu et al. (2015).

Total RNA was isolated with TRIzol reagent (Invitrogen, Life Technologies) following standard protocol from two biological replicates. RNA quality and purity were assessed with the Agilent Technologies 2100 Bioanalyzer (Agilent) according to the RNA Integrity Number (RIN) value. RNA-Seq library construction was processed following TruSeq RNA Sample Prep v2 protocol. Subsequently, the 100 bp paired-end reads were generated via Illumina HiSeq 2000.

We used NGSQCToolkit (v2.3.3) (Patel and Jain, 2012) to check and visualize the quality of the raw data, in order to trim and filter the pair-end reads containing Ns, reads containing adapters, and low quality reads (Reads of low quality base were greater than 20%). Then the clean reads were aligned to the reference genome of Brassica napus (Brassica _napus.annotation_v5.gff3.) using HISAT (HISAT version 0.1.6-beta) (Kim et al., 2015) with the default parameters, except for setting the minimum alignment score of L, 0, −0.18. To provide sensitive and accurate results, only unique mapped reads were used in further study. FPKM (Fragments per Kilobase of transcript per Million mapped reads) method was used to predict the gene expression levels (Trapnell et al., 2010).

To study the effects of the genome-dosage regulatory between A- and C-subgenome genes, we used the 31,526 homoeologous gene pairs between A and C subgenomes according to the reference genome sequence data of B. napus (Chalhoub et al., 2014). About 54% of these homoeologous gene pairs were removed due to no expression or low expression (FPKM < 1) in diploid parents, finally we selected 14380 homoeologous gene pairs with FPKM values greater than 1.0 in both diploid progenitors. All further analyses were performed on these 14,380 homoeologous gene pairs.

Pearson correlation tests between homoeologous gene expression and genotype dosage (1: 2/3: 1/2:1/2: 1/3) for each A- and C-subgenome gene were performed, respectively. Pearson test and multiple test correction were calculated using adjustment method in R. The P values in the analysis were adjusted for the multiple test correction by the Benjamini–Hochberg method (Benjamini and Hochberg, 1995). Genes with significant expression and dosage correlation (p < 0.05) were defined as dosage dependent, whereas genes with no significant expression and dosage correlation were defined as dosage independent.

GO enrichment analysis was performed using Cytoscape plug-in BiNGO (Maere et al., 2005). GO terms with corrected p < 0.05 (Benjamini and Hochberg FDR-adjusted P value) were considered to be significantly enriched. The orthologous genes in A. thaliana were used to predict the most probable function of the gene pairs as the high homology between B. napus and A. thaliana. We used the whole Arabidopsis genome gene list as background.

We measured the cis/trans effects on gene transcription by comparing the gene transcription difference between A- and C-subgenome genes as defined in previous studies (Tirosh et al., 2009; Shi et al., 2012), briefly as A = log₂(PA/PC) (both cis and trans effects), B = log₂(F₁A/F₁C) (cis effects) and A − B (trans effects) (P = parents, F₁ = hybrid / polypoid). If A = B and B ≠ 0, genes were classified as “only cis effects,” whereas genes were classified as “only trans effects” if A ≠ B and B = 0. And if A = B and B = 0, genes were “no cis-trans effects,” while genes were “cis-trans effects” if A ≠ B and B ≠ 0. Statistically significant differences were identified using a Fisher's exact test and multiple testing correction (p < 0.05).

The interaction data (Release 3.4.134 compiled on February 25th, 2016) set from BIOGRID (Stark et al., 2006; http://thebiogrid.org/) was used to detect protein–protein interactors for dosage-dependent/-independent genes. The orthologous genes in A. thaliana were used to predict the corresponding genes of Brassica. Only genes with ≥ 1 interactors were analyzed. The Wilcoxon's rank-sum test was used to test whether the numbers of dosage-dependent/-independent genes were significantly different.

The chromosome complements of synthesized B. napus hybrid (AC, 2n = 19, A: C = 1: 1) and allotetraploid (CCAA, 2n = 38, A: C = 1: 1), and two derived allotriploids (AAC, 2n = 29, A: C = 2: 1; CCA, 2n = 28, A: C = 1: 2) were confirmed by fluorescence in situ hybridization (FISH) analyses with C-genome specific probe (Figure 1C), before they were used to study the gene expressions. While the hybrid and allotetraploid had an intermediate phenotype between two parental diploids, the allotriploid (AAC) was more biased to B. rapa and another one (CCA) to B. oleracea (Figure 1B), which suggested the dosage effects of the component genomes on the morphological expression.

To analyze genome-wide gene expression levels of leaves from these synthetics using RNA-seq, 19–45 million sequencing clean reads were obtained from each of two biological replicates (Table S1). The gene expression levels between two replicates were correlated very well (Average R = 0.93, Figure S1). FPKM (Fragments per Kilobase of transcript per Million mapped reads) values were used to represent the gene expression levels (Trapnell et al., 2010). The proportions of the number of expressed genes (FPKM > 0) for A- and C-subgenome were about 48.3 and 51.5% in AC, CCAA, and MPV (mid-parent expression values, representing the progenitors), respectively (Figure 2A, Table S2). The proportions were consistent with those of the reference genome sequence data of B. napus (Chalhoub et al., 2014). However, the proportions of A- and C-subgenome genes were significantly different in AAC (Chi-square test, p < 2.2e-16; 53.4 and 46.4%, respectively) and CCA (Chi-square test, p < 2.2e-16; 45.1 and 54.7%, respectively), in comparison with MPV (Figure 2A, Table S2). In addition, the changes of the proportion between A- and C-subgenome genes were in accordance with dosage alterations, suggesting that dosage variables exert an influence on gene expression.

Then the genome-wide gene expression levels (FPKM > 0) of different materials were investigated, and the results indicated that gene/genome dosage imbalance influenced gene expression levels and global gene expression levels were positively related with genome dosages. As shown in Figure 2B, the global expression level of A-subgenome genes was significantly higher than that of C-subgenome genes in AAC (Wilcoxon signed-rank test, p < 2.2e-16), whereas the expression level of A-subgenome genes was significantly lower than that of C-subgenome genes in CCA (Wilcoxon signed-rank test, p < 2.2e-16). But there were no significant expression differences between A- and C-subgenome genes in dosage-balance materials (AC and CCAA, A: C = 1: 1) and the parents (MPV, A: C = 1: 1) (Wilcoxon signed-rank test, p > 0.05). The results confirmed that gene expression levels could be impacted by genome-dosage balance.

To assess what extents of expression levels have changed among materials with different genome dosages, 14,380 homoeologous gene pairs (see Section Materials and Methods) were used to compare expression levels with MPV (mid-parent expression values) calculated by averaging expression values observed in the diploid progenitors. Fold change (FC) was used to quantify the change of expression levels, compared with progenitors. We found that significant differential expressions occurred between A- and C-subgenome genes when the genome dosages varied. The distributions of FC values were almost the same in AC and CCAA with the balanced genome-dosage, and ranged 0.75–1.25 (interquartile range, 0.75–1.25) and the median was about 1.0 (Figure 2C), indicating the genes of both subgenomes expressed highly consistently with those in two parents. Differently in AAC and CCA with the imbalanced genome dosages, the magnitude of the FC values was significantly greater for A-subgenome genes (median = 1.26-fold) than that of C-subgenome genes (median = 0.75-fold) in AAC, while the opposite situation happened in CCA (Wilcoxon's rank-sum test, p < 2.2e-16; Figure 2C). Furthermore, the magnitude of FC increased with the proportion of subgenome dosage. In other words, the position of the interquartile range and the median of A-subgenome genes in AAC (approximately 1.0–1.6 and 1.25, respectively) and C-subgenome genes in CCA (approximately 1.0–1.5 and 1.25, respectively) were obviously higher than 1.0 for the median of AC and CCAA. Whereas the interquartile range and the median of C-subgenome genes in AAC (approximately 0.5–1.0 and 0.75, respectively) and A-subgenome genes in CCA (approximately 0.6–1.0 and 0.75, respectively) were definitely less than 1.0 (Figure 2C), suggesting single dosage effects responding to genome-dosage changes. Besides, the genome-wide gene expression levels of AC and CCAA were almost the same, for no significant difference was detected (Wilcoxon's rank-sum test, p > 0.05; Figure 2B). The observation supported the view that the effect of dosage was not simply due to the absolute dosage of the copy number, but rather resulted from a change in the relative dosage balance between A and C subgenomes (Birchler et al., 2005; Bekaert et al., 2011).

There was also dosage compensation at the mRNA level. We found that it was not the expected 2-fold change in gene expression levels along with the change of subgenome ratio (A/C = 2) in AAC. Data showed that the positions of the interquartile range and the median of A-subgenome genes were about 1.0–1.6 and 1.25, respectively, but were about 0.5–1.0 and 0.75 for C-subgenome genes (Figure 2C). It was similar in CCA. In other words, genes expressed did not show the expected 2-fold change in gene expression levels if mRNA levels correlated perfectly to gene/genome dosage, which suggested that there were buffering effects or dosage compensation at the mRNA level. Some clues could also be revealed from the global gene expression levels (Figure 2B). The results collectively indicated a complex relationship between gene dosage and expression.

To study correlations between gene expression and dosage, we calculated the correlation coefficient (R-values) between homoeologous gene expression levels (FPKM) and relative dosages of the 14,380 homoeologous gene pairs between A and C subgenomes, respectively. For example, for gene BnaA01g05230D with R = 0.97, an A-subgenome gene, expression levels of AA (all A genome), AAC (2/3 A genome), AC (1/2 A genome), CCAA (1/2 A genome), and CCA (1/3 A genome) were 22.22, 11.37, 7.41, 8.98, and 6.34, respectively. Pearson correlation between expression levels (22.22, 11.37, 7.41, 8.98, and 6.34) and the relative dosage (1, 2/3, 1/2, 1/2, and 1/3) was 0.97 (Table S3). Results showed that about 95% R-values were > 0, meaning positive correlations, and the remaining 5% R-values were < 0 for negative correlation (Figure 3A, Table S3). Moreover, among those 5% genes with R < 0, only 38 A-subgenome genes and 21 C-subgenome genes (0.26 and 0.15%, respectively) had statistical significance. While most genes (average 63%) were statistically significant among those genes with R > 0 (Figure 3A, Table S3). We concluded that the majority of genes showed positive correlations between gene expression and genome dosage in Brassica polyploids and hybrids.

The small group of negatively correlated genes which were called “inverse dosage effect” genes (Birchler and Veitia, 2012; Veitia et al., 2013) could be potentially interesting, due to their negative effect of the regulator to cancel the positive effect of the change in gene dosage on the total expression levels. In the following, we investigated the molecular function and biological pathway of the potentially small group of significantly inverse dosage effect genes with the DAVID Functional Annotation Tool (https://david.ncifcrf.gov/tools.jsp). After the orthologous genes in A. thaliana of those negatively correlated genes were submitted to the DAVID, by using the Arabidopsis genome gene list as background, the genes involved in proteolysis and/or macromolecule catabolic process (GO: 0006508, GO: 0009057) were highestly over-represented, followed by those for translation and/or ribosome (GO: 0006412, GO: 0005840; Table S4).

In order to study the effects of the relative dosage balance, coefficients of determination (R²) between dosage and expression level were used to quantify the strength of dosage effects on expression. High R²-values indicated that the gene expression levels and dosage changes were strongly correlated, and the corresponding genes were defined as dosage dependent. On the contrary, genes with low R²-values were called dosage-independent, suggesting that the expression level was weakly correlated with the dosage changes. For example, for dosage dependent BnaA10g21130D, an A-subgenome gene with R² = 0.9996, expression levels were 6.32: 3.96: 2.87: 2.85: 1.69, close to 1: 2/3: 1/2: 1/2: 1/3. In contrast, for dosage independent gene BnaC05g12440D (R² = 2.39E-9), the gene expression levels were 96.35: 71.47: 94.49: 71.31: 99.86, far from 1: 2/3: 1/2: 1/2: 1/3 (Table S3).

All genes were clustered into two groups based on R²-values using the Pearson correlation test with multiple testing correlation (Benjamini and Hochberg, 1995), as defined in Arabidopsis (Shi et al., 2015). Empirically p < 0.05 was used as the cutoff value for the statistically significant level. Maintaining the significant level at 0.05 in both A- and C- subgenome genes, we got R² > 0.59. So those genes with R² > 0.59 were called dosage dependent A (Ad) and C (Cd) genes, while others with R² < 0.59 were called dosage independent A (Ai) and C (Ci) genes (Figure 3B, Table S3). However, the distribution was continuous, there was no obvious boundaries between dosage dependent and independent genes. Under this cutoff value, about 60% genes (8930 and 8582 genes in A and C subgenomes, respectively) were dosage dependent, and the rest of 40% genes (5450 and 5798, respectively) were dosage independent (Figure 3C, Table S3).

Dosage dependent and independent genes definitely influenced the strength of dosage effects on gene expression. In CCAA, the FC values of most Ad and Cd genes (91 and 92%, respectively) distributed in “0.5–1.5” compared with MPV, significantly greater than those of Ai and Ci genes (71 and 69%, respectively; Chi-square test, p < 2.2e-16). On the contrary, less than 10% of genes appeared in the regions with FC values “> 1.5 or < 0.5” of Ad and Cd, whereas 30% of Ai and Ci genes were significantly higher (Chi-square test, p < 2.2e-16; Figure 3D). Dosage dependent and independent genes contributed differentially to variable gene expressions, likely the former to additive expression and the latter to non-additive expression.

A total of 14,380 homeologous gene pairs were classified into four groups: AdCd (40%), AdCi (22%), AiCd (20%), and AiCi (18%) (Figure 4A, Table S3), which was not significantly different from the expected percentages of AdCd (37%), AdCi (25%), AiCd (23%), and AiCi (15%) (Chi-square test, p > 0.05). The similar dosage dependency of both A- and C-subgenome genes suggested possible generalities for factors controlling which genes were dosage dependent or independent. So genes in groups of AdCd and AiCi definitely represented two opposite types divided by whether the gene expression was dependent or independent of the dosage changes.

As the gene expression was either dosage dependent or independent, we wanted to know if the absolute expression levels of genes worked. We divided the absolute expression levels ranking from low to high into “bins,” each containing 1000 genes and counted the number of dosage dependent (Ad, Cd) and dosage independent (Ai, Ci) genes in each bin. The distribution of dosage-dependent and -independent genes in bins with different absolute expression levels had no significant differences than the expected (about 60 vs. 40%) in CCAA (Chi-square test, p > 0.05; Figure 4B), indicating that the absolute expression levels of dosage-dependent genes did not differ from those of dosage-independent genes. Therefore, the absolute expression levels were not the determining factor of dosage dependence and/or independence.

Gene ontology (GO) enrichment analysis was used to classify genes according to their molecular function and the pathway in which they were involved. In this study, the orthologous A. thaliana Gene ID was used to predict the corresponding genes in B. napus. GO analysis revealed different functional enrichment between genes with dosage dependent and independent. First, the top 3000 of AdCd genes were submitted to BiNGO (Maere et al., 2005), and they were significantly enriched (p < 0.05) in lots of categories, including the most basically process, like “cellular process,” “developmental process,” “biological regulation,” “localization,” “transport,” etc. (Figure 4C, Table S5). However, AiCi genes were mainly clustered into the GO terms of “response to stimulus” in biological process, which consisted of a mass of daughter categories, such as “response to stress,” “response to abiotic stimulus,” “response to chemical stimulus” and so on (Figure 4C, Figure S2, and Table S5). So the dosage independent genes probably played a key role in stress responses. The other main modules contained the GO terms including “structural constituent of ribosome,” “DNA binding” and “transcription factor activity” in molecular function (Figure 4C, Table S5).

Gene expression changes in gene transcription could result from cis and/or trans effects (Rockman and Kruglyak, 2006; Williams et al., 2007; Yoo et al., 2014). In this study, we measured the cis/trans effects on gene transcription by comparing the transcription difference between A- and C-subgenome genes (see Section Materials and Methods), as defined in other studies (Tirosh et al., 2009; Shi et al., 2012). Results showed almost the same gene distribution of “only cis effects,” “only trans effects,” “cis-trans effects,” and “no cis-trans effects” in dosage-balanced AC and CCAA. However, the gene percentage of trans effects differed significantly in AAC (22.5%) and CCA (17.8%), compared with AC (10.7%) or CCAA (11.3%) (Chi-square test, p = 4.395e-11), whereas cis effects remained the same (Chi-square test, p = 0.2727), approximately 17.5% (Figures 5A,B, Figure S3, and Table 1). It was shown that trans effects rather than cis effects were mainly responsible for gene/genome dosage variations.

Then we wanted to know if dosage-balance regulatory influenced gene expression associated with cis/trans effects. We explored the distribution of dosage-dependent (AdCd) and -independent (AiCi) genes in groups of “only cis effects” and “only trans effects.” For cis effects, the distribution of AiCi genes was significantly greater than expected (Chi-square test, p < 0.01), whereas that of AdCd genes was significantly less than expected (Chi-square test, p < 0.01; Figure 5C), indicating that dosage-independent genes more likely influenced expression by cis effects. Inversely, dosage-dependent genes were more inclined to impact expression by trans effects, for the data showed that the distribution of AdCd genes was significantly greater than expected (Chi-square test, p < 0.01), while that of AiCi genes was significantly less than expected (Chi-square test, p < 0.01) for trans effects in dosage-imbalance group (Figure 5D).

The dosage balance hypothesis stated that dosage imbalance of subunits of macromolecular complexes could affect the eventual amount of complexes formed (Birchler et al., 2005; Veitia et al., 2008). Liang et al. (2008) predicted that the duplicates of a highly under-wrapped protein should be more sensitive to dosage imbalance. This meant that dosage-dependent genes might be expected to have fewer protein–protein interactions than dosage-independent ones. From the investigations of the number of interactors of dosage-dependent/-independent genes, the dosage-dependent genes (AdCd) had lower number of protein–protein interactions than that of dosage-independent (AiCi), but the difference was insignificant (Wilcoxon's rank-sum test, p = 0.2302; Figure 6A), possibly due to the lack of reference data and the use of the orthologous genes in A. thaliana to predict the corresponding gene in Brassica. By using previously published data in Arabidopsis (Shi et al., 2015) and reanalysis, it was found that the dosage-independent genes (TiAi) had a significant greater number of protein–protein interactions than that of dosage-dependent (TdAd) (Wilcoxon's rank-sum test, p = 0.04692; Figure 6B).

Gene copy number changes generally translate into changes in gene expression, which was supported by lots of cases in yeast, mouse, human and Arabidopsis (Tang and Amon, 2013). We also found gene expression changes caused by gene dosage or copy number alterations (Figure 2), together with phenotypes alteration driven by the cumulative effects of dosage changes of a large number of genes (Figure 1). Recent works have shown that additive expression was the prevailing gene expression pattern when two parental genomes were present in an allopolyploid nucleus, although non-additively expressed genes represented a small portion (Wang et al., 2006; Chagué et al., 2010; Jiang et al., 2013; Yoo et al., 2013; Zhao et al., 2013). In our results, many genes generally increased their expressions along with the increase of gene/genome dosage (Figure 2), and were defined as dosage dependent, for their expressions showed well correlations with dosage changes. However, some genes defined as dosage independent did not exhibit the simple additivity that was predicted to be buffered, because their expressions had poor correlation with dosage changes. Gene dosage balance was one of the important factors for non-additive gene expression in polyploids, although our understanding remained limited (Yoo et al., 2014).

In our study, both A and C subgenomes genes were subjected to dosage dependent and independent expression regulations in synthetic B. napus, which represented approximately 60 and 40%, respectively. About 58% of A and C-subgenome genes were in the same direction to either dosage-dependent (40% AdCd) or dosage-independent expression (18% AiCi), while 42% of genes (22% AdCi and 20% AiCd) were in different directions (Figure 4A). In comparison with the results of TdAd (54%), TdAi (15%), TiAd (13%), and TiAi (17%) from Arabidopsis (Shi et al., 2015), our data indicated a little different percentage distribution. As to the reasons for the difference, the first one might be the different cutoff values used to distinguish dosage-dependent and -independent genes, as there was no clear distinction between them. The second was perhaps the different genetic backgrounds between Arabidopsis and Brassica species, as diploid Brassica genomes were triplicated compared with A. thaliana, as confirmed by recent whole genome sequencing of B. rapa (Wang et al., 2011), B. oleracea (Liu et al., 2014; Parkin et al., 2014), and B. napus (Chalhoub et al., 2014). Thereafter along with some losses of duplicated copies from the genomes during evolution, it was more complex for gene-dosage alterations in Brassica.

Gene expression changes in transcription could result from cis and/or trans effects (Rockman and Kruglyak, 2006; Williams et al., 2007; Yoo et al., 2014). It was previously suggested that trans-acting factors played a much larger role than cis-factors in causing gene expression variation between different species (Dong et al., 2011; McManus et al., 2014; Wang et al., 2015). Our results also showed that trans effects played a key role in response to gene/genome dosage alterations (Table 1). Besides, we observed that dosage-dependent genes more likely impacted expressions by trans effects (Figure 5D), whereas dosage-independent genes were more inclined to influence expression by cis effects (Figure 5C). The cis and trans regulatory factors differed in influencing the evolution of gene regulation, cis-effect regulatory affected the expression of nearby genes (e.g., changes in promoters and enhancers) on the same chromosome (Wittkopp et al., 2008; Tirosh et al., 2009; Dong et al., 2011; Wittkopp and Kalay, 2012). Thus, the regulation via many complex subunits possibly explained why cis effects overrepresented in dosage-independent genes. But trans- effect regulatory impacted both alleles of the diploid progenitors, dosage variations of alleles belonging to dosage-dependent genes most likely regulated expression via trans effects.

Dosage compensation was a common biological phenomenon that was clearly supported by many cases in a wide range of organisms, including Drosophila, mouse, human, maize, etc. (Veitia and Potier, 2015). The most typical example was X chromosome inactivation in mammals, which was a mechanism that equalized the number of active X chromosomes in eutherian females (XX) and males (XY) (Pessia et al., 2012; Veitia and Potier, 2015). In our study, we also found that there were not only clear dosage effects of gene expression responding to gene/genome dosage changes, but also obviously dosage compensation effects. As showed in Figure 2, the genes expressed did not show the expected 2-fold change in their expression levels if mRNA levels correlated perfectly to gene/genome dosage in AAC and CCA, suggesting dosage compensation at the mRNA level.

Veitia et al. (2013) suggested that the inverse dosage effect genes or negatively operated dosage effectors could lead to dosage compensation, as the positive effect of the dosage changes on expression was canceled by the negative effect of the regulator. In our results, some inverse dosage effect genes were detected (Figure 3A). Such genes with statistical significance involved in proteolysis (GO: 0006508) and ribosome (GO: 0005840) were highly over-represented (Table S4), and might contribute to dosage compensation. The result was consistent with previous study that proteolysis was predicted to play a crucial role in buffering the effects of gene dosage alterations (Veitia et al., 2008; Lundberg et al., 2012; Veitia and Potier, 2015), as induction proteolysis appeared to be a general response to the genomic imbalance due to aneuploidy (Lundberg et al., 2012). The cell might increase the level or activity of proteolysis and/or chaperones to cope with the overexpression of hundreds of proteins, many of which belonged to complexes, caused by the presence of supernumerary chromosomes (Veitia and Potier, 2015). In addition, the result that ribosome (GO: 0005840) was highly represented in the inverse dosage effect genes was consistent with a frequent biological phenomenon of nucleolar dominance, a non-additive or uniparental expression of rRNA genes, observed in many allopolyploids including B. napus (Chen and Pikaard, 1997).

Furthermore, we found that dosage-independent genes were more likely buffered for response to dosage changes than dosage-dependent genes (Figure 2C), which coincided with the result that the greater number of protein–protein interactions was investigated in dosage-independent genes (Figure 6B), suggesting that the existence of dosage-balanced multi-subunit complexes or networks contributed to dosage compensation. Indeed, such findings were reported in Drosophila (Zhou, J. et al., 2011; Malone et al., 2012). Changes in gene copy number resulted in changes of protein levels in the majority of cases in aneuploid budding yeast and human cells by quantitative proteomic analyses (Pavelka et al., 2010; Torres et al., 2010; Stingele et al., 2012), while the proteins that did not show the coordinated increase with gene copy number were found to be predominantly components of large protein complexes (Torres et al., 2010; Stingele et al., 2012). This was possibly explained by the dosage balance hypothesis (Veitia et al., 2008, 2013; Birchler and Veitia, 2010) that the co-variation of a target gene along with a linked controlling gene could lead to dosage compensation at the transcriptional level.

Dosage imbalance could influence protein complexes that resulted in novel phenotypes, leading to great variability for adaptations in evolution (Birchler et al., 2005; Veitia et al., 2008; Veitia and Potier, 2015). In our GO enrichment analysis (Figure 4C, Table S5), dosage-dependent genes were mainly associated with the basic biological processes, although Lloyd et al. (2014) suggested an exceptive example that meiosis genes were often maintained in single copy per genome, revealing that these genes were important to maintain growth and developmental stability, and might provide genetic stability against null mutations and selective advantage by dosage-dependent gene regulation (Ha et al., 2007). In other side, dosage-independent genes were more likely enriched in stress response related processes, suggesting that these genes might accelerate fitness and adaptation. Furthermore, the fact that dosage independent genes with Gene Ontology terms associated with response to important environmental factors, transcription factors and ribosomal protein genes was in accordance well with the fact that genes were most likely retained in WGD duplicates in B. rapa (Wang et al., 2011). This result was often explained by the gene dosage hypothesis in which genes encoding products that interacted with one another should be over retained and genes with products that did not interact with other gene products should be lost (Birchler et al., 2005; Freeling and Thomas, 2006; Birchler and Veitia, 2010; Wang et al., 2011; McGrath et al., 2014). It was similar to previous report of Bekaert et al. (2011) that constraints on the relative dosages of central network genes represented an important force for maintaining duplicates. Constraints on dosage-balance over time could result in preferential rewiring of certain biological pathways to execute novel functionality (De Smet and Van de Peer, 2012), and likely helped them to cope with new ecological opportunities and/or challenges (Schranz et al., 2012; Fawcett et al., 2013).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.