Transcriptional Aneuploidy Responses of Brassica rapa-oleracea Monosomic Alien Addition Lines (MAALs) Derived From Natural Allopolyploid B. napus

Establishing the whole set of aneuploids, for one naturally evolved allopolyploid species, provides a unique opportunity to elucidate the transcriptomic response of the constituent subgenomes to serial aneuploidy. Previously, the whole set of monosomic alien addition lines (MAALs, C1-C9) with each of the nine C subgenome chromosomes, added to the extracted A subgenome, was developed in the context of the allotetraploid Brassica napus donor “Oro,” after the restitution of the ancestral B. rapa (RBR Oro) was realized. Herein, transcriptomic analysis using high-throughput technology was conducted to detect gene expression alterations in these MAALs and RBR. Compared to diploid RBR, the genes of all of the MAALs showed various degrees of dysregulated expressions that resulted from cis effects and more prevailing trans effects. In addition, the trans-effect on gene expression in MAALs increased with higher levels of homology between the recipient A subgenome and additional C subgenome chromosomes, instead of gene numbers of extra chromosomes. A total of 10 trans-effect dysregulated genes, among all pairwise comparisons, were mainly involved in the function of transporter activity. Furthermore, highly expressed genes were more prone to downregulation and vice-versa, suggesting a common trend for transcriptional pattern responses to aneuploidy. These results provided a comprehensive insight of the impact of gene expression of individual chromosomes, in one subgenome, on another intact subgenome for one allopolyploid with a long evolutionary history.


INTRODUCTION
As a deviation from the normal genome, by either gaining or losing entire chromosomes or chromosomal segments, aneuploidy disrupts the genome balance, which generally results in profoundly and severely impaired phenotypes (Antonarakis et al., 2004;Birchler and Veitia, 2007;Huettel et al., 2008;Henry et al., 2010;Siegel and Amon, 2012;Zhu et al., 2015). Some syndromes and diseases, which lead to severe developmental defects and even death in humans, have been associated with variants of entire chromosomes or chromosomal segments (Antonarakis et al., 2004;Morrow, 2010), such as trisomy 21 (Down syndrome) (Letourneau et al., 2014). Moreover, typical chromosomal structural changes, which occur extensively in the majority of solid tumor cells, are believed to induce tumor development (Matzke et al., 2003;Sheltzer and Amon, 2011;Gordon et al., 2012). Since aneuploid studies that involve a single or few chromosomes are limited, genome-wide transcriptional response patterns to aneuploidy remains unclear. Next-generation sequencing, which has been extensively applied to transcriptional analysis, provides a robust means to elucidate gene expression in aneuploid organisms.
Researchers initially believed that phenotypic consequences of aneuploidy are mainly influenced by gene dosage effects of altered chromosomes, which is supported by the correlation between chromosome copy numbers and relative gene expression levels in Down syndrome (Mao et al., 2003), aneuploid yeast (Torres et al., 2008), as well as various types of cancer cells (Gao et al., 2007;Xu et al., 2010;Lu et al., 2011). However, studies on aneuploids in Drosophila and plants showed that changes in the gene expression levels were observed along altered copy numbers (cis-effects) as well as in unaltered disomic chromosomes (trans-effects) (Guo and Birchler, 1994;Birchler and Veitia, 2007;Huettel et al., 2008;Malone et al., 2012;Zhu et al., 2015Zhu et al., , 2018Zhang et al., 2017;Rey et al., 2018). In addition, the phenotypic consequences of aneuploidy are believed to result from either the contribution of altered dosages or the impact of a genome imbalance involving the rest of the genome (Guo and Birchler, 1994;Birchler and Veitia, 2007). Aneuploids exhibit a more severe degree of perturbation in gene expression than changes in ploidy (Birchler and Newton, 1981;Guo et al., 1996), indicating a disruption in the dosage-sensitive gene products (Birchler et al., 2001;Henry et al., 2006).
Plants have better tolerance to chromosomal changes than animals do (Siegel and Amon, 2012). The frequency of aneuploidy increases in polyploids, and the whole set of aneuploids with the loss of each chromosome was established in allohexaploid wheat (Sears, 1954), which may be attributed to a significantly higher capacity to buffer dosage imbalance. Moreover, the aneuploidy with whole-chromosome and structural alterations extensively occurs in nascent allopolyploids (Mestiri et al., 2010;Xiong et al., 2011;Chester et al., 2012;Zhang et al., 2013), implying that aneuploidy is a continuously occurring process until a stable polyploid status is attained. Recent studies have suggested that aneuploidy serves as a rapid mechanism for the adaptation to rare cases of strong selective pressure (Pavelka et al., 2010;Siegel and Amon, 2012).
The whole set of aneuploids in one species is of great value to elucidate the genome structure, relationships and functional interplay from a classical genetic analysis to genome sequencing (Sears, 1954;Mayer et al., 2014;Zhang et al., 2017). However, the establishment of a complete stock of aneuploids is relatively time consuming and has only been achieved for a few species, such as bread wheat (Sears, 1954) and tobacco (Clausen and Cameron, 1944). The allotetraploid Brassica napus L. (2n = 4x = 38, AACC), which originated only 7500 years ago from natural hybridization between the ancestors of two extant diploids, B. rapa L. (2n = 2x = 20, AA) and B. oleracea L. (2n = 2x = 18, CC) (Nagaharu, 1935;Chalhoub et al., 2014), has been used as a model system to investigate the genomic alterations and interactions during allopolyploidization (Albertin et al., 2006;Szadkowski et al., 2010;Zhang et al., 2015). In particular, the ancestral B. rapa (RBR Oro) was restituted from natural B. napus genotype "Oro" by successively inducing the preferential elimination of the C subgenome chromosomes in intertribal crosses with another crucifer (Tu et al., 2010). Then, by crossing and backcrossing RBR Oro with B. napus donor "Oro, " the C subgenome was dissected in situ by adding each of its nine chromosomes to RBR Oro and establishing the whole set of B. rapa-oleracea monosonic alien addition lines (MAALs) (Zhu et al., 2016). These MAALs, with the genetic context of the same natural B. napus, provide a unique avenue to decipher the gene expression interaction of two constituent subgenomes under aneuploidy situations.
In view of this, we performed a RNA-Seq to detect the changes in the global gene expression in these nine MAALs, compared to RBR Oro, to reveal transcriptional aneuploidy response patterns associated with the disomic A subgenome chromosomes. The results of this study provide new insights into the impact of individual chromosomes, from one subgenome on another complete subgenome, for one natural allopolyploid.

Plant Materials
The ancestral B. rapa (RBR Oro) was restituted from a natural B. napus cultivar "Oro, " by inducing the preferential elimination of the C subgenome chromosomes in intertribal crosses between B. napus and another crucifer Isatis indigotica (Tu et al., 2010). Then, a complete set of nine B. rapa-oleracea MAALs (2n = 21, AA+1C 1−9 ) was established to dissect the C subgenome in the A subgenome background through crossing and successive backcrossing between Oro and RBR Oro (Zhu et al., 2016), which was designed, for convenience, as C1-C9. Due to severe self-incompatibility and variable transmission rates of extra chromosomes in MAALs, MAALs plants were pollinated by RBR Oro to generate sufficient progeny seeds. Seedlings were grown in the greenhouse at the Huazhong Agriculture University (Wuhan, China). Fifty progeny plants were established for C1, C2, C3, C5-C8, respectively, and 100 plants for MAALs C4 and C9, due to the much lower transmission rate of these two C subgenome chromosomes via female gametes (Zhu et al., 2016). At least six additional progeny plants for each of the nine MAALs were identified for further analysis.
To screen the additional plants with extra chromosomes of the C subgenome, these progeny plants were first investigated by PCR amplification of the chromosome-specific gene markers that were distributed on both arms of the target chromosomes, and 6 to 21 target plants were identified. Then, a chromosome counting and fluorescence in situ hybridization (FISH) analysis with the C subgenome-specific probe (Zhu et al., 2016) was conducted to determine the karyotype integrity of the target plants (Supplementary Figure S1). The ovaries of young flower buds were collected and treated with 2 mM 8-hydroxyquinoline for 3 h at ∼20 • C, and then fixed in Carnoy's solution (3:1 ethanol: glacial acetic acid, v/v) and stored at −20 • C for cytological analysis. A cytological observation and FISH analysis was conducted according to Li et al. (1995) and Cui et al. (2012), respectively.

RNA Extraction and Preparation of cDNA Libraries
Since there were minimal phenotypic differences between RBR and MAAL plants at the three-leaf stage, the third newly expanded leaves, without petioles, from target plants were collected and immediately stored in liquid nitrogen for RNA extraction. Three biological replicates, per MAAL and RBR Oro, were prepared to assess gene expression. Briefly, young leaves from two target plants per MAAL were pooled equally and ground in liquid nitrogen. Total RNA was extracted using a commercial RNA kit (Tiangen, China) according to the manufacturer's protocol. The quality and quantity of the isolated RNA was assessed using a Qubit (Invitrogen, Life Technologies), and the RNA integrity number (RIN) was then assessed with an Agilent Technologies 2100 Bioanalyzer (Agilent). Samples with RIN < 8 were excluded from further analysis. For the cDNA library construction, 1.5 µg of total RNA per biological replicate was prepared according to the TruSeq RNA Sample Prep v2 protocol (Illumina). Subsequently, 30 libraries were sequenced on an Illumina HiSeq TM 3000 platform (Illumina), to generate 150bp paired-end reads. The raw sequence data were subjected to data quality control checks based on per base sequence quality scores, per sequence quality scores, per sequence GC content, sequence length distribution, and overrepresented sequences, which are implemented in the FastQC 1 .

Differentially Expressed Genes (DEGs) Analysis
Trimmomatic version 0.33 2 (Bolger et al., 2014) was used to remove adapters and low-quality reads, to generate clean reads with the following parameters: LEADING: 3, TRAILING: 3, SLIDINGWINDOW: 4:15, MINLEN: 36, LEADING: 3, TRAILING: 3, SLIDINGWINDOW: 4:15, MINLEN: 36, and TOPHRED: 33. Clean reads were then aligned to the B. napus reference genome with 101,040 predicted genes 3 (Chalhoub et al., 2014) using HISAT version 0.1.6 4 (Kim et al., 2015). Considering the homologous regions between the A sub-genome and by adding chromosomes of the C subgenome, we performed a strict mismatch tolerance (one-base mismatch) to determine the affiliation of these clean reads according to the following criterion: specifying -L, 0, −0.15, and only using the uniquely mapped reads, in gene expression profiling. The fragments per kilobase per million mapped reads (FPKM) value was calculated using Cufflinks version 2.2.1 5 with default parameters (-p -Gmulti-read-correct) to assess gene expression (Trapnell et al., 2012). To identify differences in gene expression between MAALs and euploid RBR Oro, the package of Deseq2 in the R-project version 3.4.3 6 , with a fold change of at least 2 in the gene expression and false discovery rate of q < 0.05 (Trapnell et al., 2013), was employed to identify DEGs. Because one replicate for C3 and C8 showed a low pairwise correlation coefficient (0.44 and 0.67 in C3, 0.47 and 0.65 in C8) with the remaining two replicates, respectively, only two replicates for C3 and C8 were performed to calculate the DEGs. For GO enrichment, because of the lack of gene annotation information in B. napus, we used Blast2GO 7 to retrieve the B. napus GO annotation file, by searching the protein database of the GenBank (Conesa et al., 2005). The gene ontology (GO) enrichment analysis was then performed by hypergeometric distribution in R, with an adjusted p-value under 0.05 as a cutoff, to determine significantly enriched GO terms (Conesa et al., 2005 8 ). Another tool for GO analysis, AgBase GoAnna version 2.00 9 (McCarthy et al., 2006) was used to annotate the co-regulated trans-effected DEGs (adjusted p < 0.05).

Data Statistics and Visualization
To check the statistical significance of the dysregulated expression between each MAAL vs. RBR Oro, a Chi-square test was used with a 0.05 q-value as the cutoff. To determine whether the trans-acting effects on gene expression were associated with the gene number of extra chromosomes, Pearson correlation coefficient (PCC) with an FDR-adjusted p-value < 0.05 as the cutoff was used. To check whether the different expression levels of the genes were similarly impacted by these extra chromosomes of the C subgenome, a Wilcoxon signed-rank test was used with 0.05 q-value as the cutoff. Statistical analyses were conducted using the R-project. All of the RNA-Seq data, including raw data and process documentation used in this study is available at the Gene Expression Omnibus (GEO), accession number GSE111510.
Fluorescence in situ hybridization images were captured using a computer-assisted fluorescence microscope with a CCD camera (Axio Scope A1, Zeiss, Oberkochen, Germany). Images of the gene expression analysis were generated by plotting functions of the R-project and then composed by Adobe Illustrator version CC.

Pairwise Comparisons of Cis-and Trans-Effects on Gene Expression
In the present study, 14.76−60.78% of the DEGs were directly attributed to the extra chromosomes of the C subgenome (ciseffects), and the majority consisted of upregulated genes, primarily due to the genome dosage. However, several downregulated genes on additional chromosomes of the C subgenome were observed in eight MAALs, except for C7, probably as a result of the failure to distinguish noise caused by the high expression of their homeologs in RBR.

Trans-Effect Dysregulation Genes Not Associated With the Gene Number of Extra Chromosomes, and Different Transcriptional Responses of Disomic Chromosomes to Aneuploidy
Analysis of the Pearson correlation coefficient (PCC) did not detect any correlations (R 2 = 0.11, p = 0.39) between the number of genes along the additional C chromosomes and corresponding dysregulated genes on the A subgenome. Furthermore, no strong relationship between the expressed genes (FPKM > 0.1) of the extra chromosomes of the C subgenome and changes in gene expression profiles in the A subgenome, were observed (R 2 = 0.27, p = 0.15). Inversely, the gene expression along each chromosome of the A subgenome was distinctly impacted by different additional chromosomes in MAALs, as the extent of the trans-effects on each of its 10 chromosomes, were uneven in all comparisons except for RBR vs. C3, and the genes on some chromosomes were more prone to transcriptional perturbation, indicating that the identity of additional chromosomes could affect gene expression throughout the A subgenome in a specific manner (Supplementary Table S2). However, MAAL C3 exhibited a similar change in the expression of genes in the A subgenome chromosomes (1.10−2.03% of the genes were .08 * * The group of up-/down-DEGs is significantly higher than oppositely regulated genes (chi-square, p < 0.01). A,B,C,D Different groups were calculated by chi-square (p < 0.001). A,B,C,D,E,F Different groups were calculated by chi-square (p < 0.001). * * The group of up-/down-DEGs is significantly higher than oppositely regulated genes (χ 2 , * p < 0.05, * * p < 0.01).
dysregulated). To further test the relationship between transeffects on the A subgenome and extra chromosomes of the C subgenome, we then calculated the PCC values of all of the genes in the A subgenome, between any two of the nine MAALs and RBR Oro, and screened for pairwise transcriptome differences (Supplementary Table S3). Compared with other MAALs, C1 always showed significantly lower pairwise correlation coefficients, followed by C2 (Supplementary Table S3). As the homoeologous pairs (A1/C1 and A2/C2) showed the most extensive homoeology, which could frequently be replaced and compensated (Xiong et al., 2011;Chalhoub et al., 2014), we concluded that a higher similarity between the additional chromosome and the recipient genome (chromosome) induced more severe changes in gene expression of the recipient genome. The differential impacts of different additional chromosomes on gene expression in the A subgenome were also reflected by a hierarchical analysis of the similarities in terms of overall expression patterns ( Figure 1A). The differences of a global average of all gene expressions showed that MAALs C1 and C2 were most divergent in the overall expression patterns from the rest, which was in line with the result of the pairwise PCC values. Furthermore, the two MAALs were clustered together into an independent group (Figure 1A). Box plots of all expressed gene profiles (FPKM > 0.1) on the recipient A subgenome, in all of the pairwise comparisons, showed that medians of the gene expressions in MAALs C1 and C2 were far from that in RBR Oro and other MAALs, confirming the observed divergence in gene expression ( Figure 1B). In addition, the remaining seven MAALs and RBR Oro were divided into another group (Figure 1A). It was determined that MAAL C9 was closely related to RBR Oro and a subgroup was formed. It did however, exhibited a divergent phenotype from RBR Oro (Zhu et al., 2016).
Another interesting result was that among these transeffected DEGs in the pairwise comparisons of each MAAL vs. RBR Oro (Table 2), only the DEGs in RBR vs. C6 showed significantly enhanced expression (χ 2 test, p < 0.01). Considering that C6 was the only MAAL that flowered earlier than RBR Oro (Zhu et al., 2016), we questioned whether flowering genes constituted a large fraction of genes that were found to be highly expressed. However, only three genes related to the flowering time in B. napus, from the public Brassica Database 10 were dysregulated, including one downregulated (BnaA03g39820D) and two upregulated DEGs (BnaA01g15350D and BnaA07g13990D).

Putative Functions of Differently Regulated Gene Groups
To determine whether some core genes or biological pathways of the recipient A subgenome in MAALs responded to additional chromosomes, Venn diagrams were constructed to determine the specific DEGs and co-regulated DEGs. Due to the ambiguous results in one diagram and because of too many pairwise comparisons, the Venn diagrams were divided into four sections, including three groups based on the results of divergent gene expressions and one group comprising of the co-dysregulated genes from the former three groups (Supplementary Figure S2). Each pairwise comparison uniquely detected 22.80−46.29% of trans-effected DEGs (Figure 2). Among these pairwise comparisons, RBR vs. C1 had a significantly higher number of specific DEGs (1,342 DEGs, occupying 46.29% of transeffected DEGs) (χ 2 test, p = 2.60e-10), followed by RBR vs. C8 (225 DEGs, 32.80%), RBR vs. C6 (206 DEGs 31.54%), and RBR vs. C4 (504 DEGs, 28.43%), whereas the remaining five MAALs were within the lower ranges of 22.80−26.75%. However, only 10 trans-effected dysregulated genes, including one downregulated gene and nine upregulated genes, were always observed in all pairwise comparisons. These results indicate that, despite the uniformity of the global transcriptional aneuploidy response patterns, there was no significant overlap in individual dysregulated genes between different pairwise comparisons of each MAAL vs. RBR Oro. To elucidate whether these 10 codysregulated genes were involved in one or specific biological pathway(s), in response to the consequences of aneuploidy on gene expression, we annotated these co-dysregulated genes into the public Gene Ontology (GO) database. However, these genes failed to be clustered into any GO terms or KEGG pathways,  probably due to the low number of genes used in the analysis. Therefore, we had to deduce the putative functions of the codysregulated genes ( Table 3). According to their ortholog gene functions in A. thaliana, the upregulated DEGs were mainly involved in "copper ion binding, " "transporter activity, " as well as "protein transporter activity, " implying that MAALs responded to stress by increasing their levels of gene expression. In addition, the downregulated DEG was involved in "rRNA binding" and "structural constituent of ribosome, " which agreed with the reported decrease in protein synthesis in aneuploids (Sheltzer et al., 2012;Stingele et al., 2012).
We then performed gene annotation to investigate the functions of specific and co-regulated DEGs that were detected in at least two pairwise comparisons in all pairwise comparisons, to reveal whether these different gene groups were involved in distinct or similar biological pathways. We found that 17 GO terms were enriched with these co-regulated DEG groups in all comparisons (p < 0.05, Figure 3A) and were mainly involved in responses to stress (all p < 0.0007), abiotic stimulus (all p < 5.91e −5 ) and oxidoreductase activity (all p < 0.003), particularly in the chloroplast (all p < 4.12e −8 ). As expected, the GO terms enriched by the specific DEGs for each pairwise comparison did not significantly overlap with each other. No GO terms that were enriched in all of the comparisons were observed, and very few terms were detected in most of the comparisons. Since DEGs that were often perturbed in many pairwise comparisons probably represented common features of aneuploidy, 650 trans-affected DEGs were identified in more than four pairwise comparisons. The AgBase GoAnna website tool was then employed to analyze the gene ontology annotation of these DEGs (McCarthy et al., 2006). The majority of GO terms ( Figure 3B) enriched with these DEGs were associated with a photosynthetic process, such as "photosynthesis" (p < 3.81e −19 ), "photosynthesis, light harvesting" (p < 0.00015), "chloroplast" (p < 1.33e −47 ), and "chlorophyll binding" (p < 1.38e −5 ). The terms of "rRNA binding" (p < 0.05) and "ribosome" (p < 0.0027) were also detected, suggesting that additional C chromosomes generally had a consequence on protein synthesis. Altogether, these results indicated that dysregulated genes, in all comparisons belonging to distinct groups, were involved in diverse cellular pathways.

Highly Expressed Genes Are Prone to Downregulation and Lowly Expressed Genes to Upregulation
To analyze the transcriptional aneuploidy response patterns in MAALs, we detected whether the different expression levels of the genes were similarly impacted by these extra chromosomes of the C subgenome. Based on the gene expression level (FPKM > 0.1) in RBR Oro, we classified the genes of the A subgenome into three groups according to low (0.1 < FPKM < 10), medium (10 < FPKM < 100), and high (FPKM > 100) expression levels. About two thirds of the expressed genes (16,714 genes, representing 62.53% of the genes with RPKM > 0.1) belonged to the first group, nearly one third of these (8,699 genes, 32.55%) were in the medium expression level group, and relatively few genes (1,315 genes, 4.92%) were in the high expression level group. Compared to RBR Oro, the lowly expressed genes group gave rise to the lowest proportion of trans-effected DEGs (consisting of 1.58−7.04% of the total number of low expression genes with the average of 2.88%), indicating that the lowly expressed genes group was less susceptible to the impact of extra chromosomes of the C subgenome. Both the highly and moderately expressed genes exhibited escalating contributions to the trans-effected DEGs (comprising 3.03−13.39% of the medium-level expression genes with an average of 5.29% and occupying 5.78−36.2% of high expression genes with an average of 14.04%), suggesting that these genes were more prone to transcriptional perturbation in aneuploidy (Supplementary Table S4 and Figure 4A). Intriguingly, it was observed that upregulated genes predominantly accumulated in lowly expressed genes for all pairwise comparisons (χ 2 test, p < 0.05), whereas the downregulated DEGs were more prevalent in the medium expression and the high expression level (χ 2 test, p < 0.05) groups (Supplementary Table S4). These findings probably indicated some transcriptional patterns in response to extra chromosomes of the C subgenome. To substantiate this hypothesis, we compared the gene expression fold-change distributions of all of the pairwise comparisons for each group, with the expected distribution exhibited by RBR Oro FIGURE 4 | Impact of additional C subgenome chromosomes on changes in gene expression changes in the A subgenome. (A) Proportions of trans-effect DEGs in different gene expression levels. Based on the gene expression level, genes in the A subgenome are divided into three categories: low-expression level genes (0.1 < FPKM < 10), medium-expression level genes (10 < FPKM < 100), and high-expression level genes (FPKM > 100). For all pairwise comparisons, dysregulated genes consistently show the highest frequency in the high-expression level group, second in the medium-expression level group, and last in the low-expression level group. The y axis represents the frequency of the gene groups. (B) Frequency distribution of log 2 (fold change) in four MAALs indicate differences in transcriptional response patterns to aneuploidy in various gene expression groups. High-expression level genes are more prone to be downregulated while the opposite is true for low-expression level genes (Wilcoxon test, q < 0.05). The x axis represents the log 2 (fold-change) in gene expression levels between the MAALs and euploid RBR in three gene groups and the y axis represents frequency distribution of log 2 (fold-change). Figure S3 and Figure 4B). As expected, the highly and moderately expressed genes (except in C7, Wilcoxon signed-rank test, p = 0.88) were generally biased toward the reduced expression in MAALs (all p < 0.05). However, the low expression level genes showed increased expression in all MAALs ( Figure 4A). This result agreed with the observation that highly expressed genes contributed substantially to the downregulated and lowly expressed genes, mainly associated with upregulation in a pair of monozygotic twins, discordant for trisomy 21 (Letourneau et al., 2014), suggesting that there might be a common feature in transcriptional pattern responses to aneuploidy.

No Obvious Dysregulation Domains in MAALs
Gene expression dysregulation domains (GEDDs), wherein dysregulated genes were clustered as domains on chromosomes, were reported in human trisomy 21 (Letourneau et al., 2014) and a nullisomics in B. napus (Zhu et al., 2015). In the present study, we smoothed the distribution of gene expression fold changes (log2 ratios) along each normal chromosome for each MAAL, using the Lowess function of R to identify potential GEDDs. Unexpectedly, only three dysregulated regions on different chromosomes were observed in different comparisons (Figure 5), indicating that these GEDDs might not be a common feature for aneuploidy. Furthermore, we couldn't exclude the possibility that these three dysregulation domains resulted from macroscopic structural variants at the corresponding chromosomes.

DISCUSSION
The advent of high-throughput sequencing technologies has facilitated the analysis of genome-wide gene expressions involving aneuploids. Recently, numerous studies demonstrated that dysregulated genes in aneuploidy are not only restricted to altered chromosome but also disomic chromosomes (Huettel et al., 2008;Makarevitch and Harris, 2010;Letourneau et al., 2014;Zhu et al., 2015Zhu et al., , 2018Zhang et al., 2017), suggesting that this dysregulation in gene expression may be a collective transcriptional response to aneuploidy (Sheltzer et al., 2012). Phenotypic consequences in aneuploidy might result from cis-effects or trans-effects or global alterations of the entire regulatory system (Prestel et al., 2010). This theory of genome imbalance involves disruptions in the dosage-sensitive gene products (Birchler et al., 2001;Henry et al., 2006). It is therefore essential to detect whether changes in transcriptional expression are in direct proportion to alterations in the DNA copy number, or whether dosage compensation minimizes the effects of aneuploidy. Our results clearly showed that when compared to parental RBR, all of the MAALs exhibited prevalent trans-effect gene expression. However, no strong correlations were observed between trans-effects and levels of the chromosome genes number or the expressed gene number. Intriguingly, the MAALs harboring chromosome C1 or C2, showed a significantly higher number of transeffect dysregulated genes ( Table 2) and also provided the closest transcriptional responses based on gene expression divergence ( Figure 1A). From the extensive homoeology FIGURE 5 | Dysregulated domains of gene expression in different chromosomes in some pairwise comparisons. Log 2 (fold-change) value in the gene expression between one of the MAALs and RBR in the chromosomes is depicted in gray. Red lines denote the smoothed distribution for the differentially expressed genes in these selected chromosomes, as representatives, using the Lowess function of R. Dysregulated domains that are exhibited in dotted boxes are clearly observed in only three chromosomes in three pairwise comparisons (chromosome A1 in C5/RBR, chromosome A10 in C4/RBR and chromosome A02 in C6/RBR). The y axis represents the log 2 (fold change) value of FPKM between the MAALs and euploid RBR. The x axis represents the sorted positions of genes on these chromosomes. between chromosome pairs (A1/C1 and A2/C2) (Chalhoub et al., 2014), we proposed that the levels of transcriptional responses to aneuploidy were associated with the similarity of the gene content between the remainder of the genome and the altered chromosome(s), rather than with the extent (gene number) of gene alterations. Another intriguing observation was that significantly different levels of gene expressions were observed in the recipient A subgenome, in response to different extra chromosomes of the C subgenome, wherein a higher number of dysregulated genes were observed in more similar chromosomes, such as A1 and A2, than those of other chromosomes.
A landmark study on gene expression from aneuploidy cells in diverse organisms demonstrated that aneuploids of various chromosomes and in different organisms, impacted similar cellular pathways and resulted in stereotypical anti-proliferative responses (Sheltzer et al., 2012). Although different pairwise comparisons between MAALs revealed no significant overlap of individual dysregulated genes in the A subgenome and only 10 co-regulated genes were detected, their functions were clearly involved in either negative proliferation or enhanced stress. GO analysis of co-regulated DEGs, which were identified in more than one comparison, also suggested that similar cellular pathways were related to changes in transcription profiles due to trans-effects.
In our study it was observed that highly expressed genes, on disomic chromosomes of the A subgenome, were more prone to transcriptional perturbation in all pairwise comparisons. Furthermore, the set of highly expressed genes seemed to be downregulated in response to aneuploidy. Inversely the set of naturally lowly expressed genes was likely upregulated. A similar phenomenon has been observed in one B. napus nullisomics, that lost the chromosome C2 (Zhu et al., 2015) as well as in monozygotic twins discordant for trisomy 21 (Letourneau et al., 2014). Interestingly, in these monozygotic twins, either upregulated or downregulated genes, between the twins, were clustered in the domains along all of the chromosomes (GEDDs), in which the highly expressed genes were mainly distributed in downregulated domains while lowly expressed genes were distributed in upregulated domains (Letourneau et al., 2014).
It was speculated that these dysregulation domains were associated with some phenotypic features of trisomy 21 (Letourneau et al., 2014;Pope and Gilbert, 2014), or that they could be the result of extra chromosomal material (Letourneau et al., 2014). Similar dysregulation domains observed in nullisomic B. napus likely supported the latter hypothesis (Zhu et al., 2015), although a transcriptome analysis of trisomy 5 in Arabidopsis thaliana (Huettel et al., 2008) and a series aneuploidy in hexaploid wheat (Zhang et al., 2017) did not detect these domains. These results, however, were derived from aneuploid organisms involving only one or a few chromosomes. Among our nine MAALs, only three dysregulated regions on different chromosomes were observed in three comparisons, indicating that the dysregulation domains were likely specific to particular aneuploid organisms.
Altogether, transcriptomic analysis of MAALs provided a comprehensive insight into transcriptional aneuploidy response patterns and allowed a deeper understanding of the functional interplay between the intact A subgenome and the extra C subgenome chromosomes in the natural allopolyploid B. napus.

AUTHOR CONTRIBUTIONS
ZB designed this study and wrote the manuscript. PQ analyzed the RNA-Seq data. ZP, CB, and HD phenotyped and identified parental RBR "Oro" and MAALs. GX provided the analysis of FISH for all materials. LZ edited the manuscript. All authors have commented and approved the final manuscript.    High-expression level genes are more prone to be downregulated while the opposite is true for low-expression level genes in all pairwise comparisons (Wilcoxon test, q < 0.05). The x axis represents the log 2 (fold change) value in gene expression between the MAALs and euploid RBR in the three gene groups, and the y axis represents the frequency distribution of log 2 (fold-change).