Skip to main content
Download PDF

The role of non-additive gene action on gene expression variation in plant domestication

EvoDevo volume 14, Article number: 3 (2023) Cite this article

Abstract

Background

Plant domestication is a remarkable example of rapid phenotypic transformation of polygenic traits, such as organ size. Evidence from a handful of study cases suggests this transformation is due to gene regulatory changes that result in non-additive phenotypes. Employing data from published genetic crosses, we estimated the role of non-additive gene action in the modulation of transcriptional landscapes in three domesticated plants: maize, sunflower, and chili pepper. Using A. thaliana, we assessed the correlation between gene regulatory network (GRN) connectivity properties, transcript abundance variation, and gene action. Finally, we investigated the propagation of non-additive gene action in GRNs.

Results

We compared crosses between domesticated plants and their wild relatives to a set of control crosses that included a pair of subspecies evolving under natural selection and a set of inbred lines evolving under domestication. We found abundance differences on a higher portion of transcripts in crosses between domesticated-wild plants relative to the control crosses. These transcripts showed non-additive gene action more often in crosses of domesticated-wild plants than in our control crosses. This pattern was strong for genes associated with cell cycle and cell fate determination, which control organ size. We found weak but significant negative correlations between the number of targets of trans-acting genes (Out-degree) and both the magnitude of transcript abundance difference a well as the absolute degree of dominance. Likewise, we found that the number of regulators that control a gene’s expression (In-degree) is weakly but negatively correlated with the magnitude of transcript abundance differences. We observed that dominant-recessive gene action is highly propagable through GRNs. Finally, we found that transgressive gene action is driven by trans-acting regulators showing additive gene action.

Conclusions

Our study highlights the role of non-additive gene action on modulating domestication-related traits, such as organ size via regulatory divergence. We propose that GRNs are shaped by regulatory changes at genes with modest connectivity, which reduces the effects of antagonistic pleiotropy. Finally, we provide empirical evidence of the propagation of non-additive gene action in GRNs, which suggests a transcriptional epistatic model for the control of polygenic traits, such as organ size.

Background

Plant domestication has been widely employed as a model to study phenotypic evolution [1]. One of the most remarkable morphological features that distinguish a domesticated from its wild progenitor is the gigantism of vegetative and reproductive organs [1,2,3]. Studies employing genetic crosses between domesticated and wild plants have revealed this gigantism is a recessive feature, as the domesticated phenotype is often masked in F1 hybrids [4,5,6]. This evidence suggests that genetic factors underlying domestication-related phenotypes might include loss-of-function (LOF) mutations whose function can be complemented by a wild functional allele. Furthermore, it has been reported that these domestication-related phenotypes are often disadvantageous or even lethal in natural populations, suggesting these recessive phenotypes are likely the product deleterious mutations [7]. Empirical work revealed these loss-of-function mutations often affect the regulation of genes rather than their coding products [2, 6, 8]. Which suggests a relevant role of transcript abundance variation in the modulation of domestication-related phenotypes. Indeed, several studies have reported that variation on transcript abundance of genes involved in organ size tends to show a recessive pattern of gene action for the domesticated phenotype [9, 10]. Transcript abundance is controlled by regulatory interactions between cis-acting non-coding DNA sequences and trans-acting factors, such as transcription factors and chromatin remodeling factors [11]. These interactions between regulatory elements act coordinately within the cells forming genetic circuits or gene regulatory networks (GRNs) to determine quantitative phenotypes, such as organ size [12].

Genes coding for cell cycle and cell fate determination are among relevant members of GRNs controlling meristem and, therefore, organ size, as they control the rates of cell division and cell expansion [13]. It is likely that mutations affecting the regulatory activity trans-acting factors with a larger number of transcriptional targets produce greater pleiotropic effects than genes connected with a few genes [14]. Hence, it is possible that the former experience greater selective constraint to avoid antagonistic pleiotropy (negative effects on fitness due to regulatory variation on multiple target genes). Similarly, genes that are regulated by a large number of trans-acting factors would experience selective constraint as they have key roles in development [15].

In diploid organisms, cis-acting mutations often produce additive gene action. A mechanistic explanation for this pattern is that cis-acting mutations act in an allele-specific fashion, so the expression of one allele does not affect the expression of the other and the resulting phenotype is the average of the parental alleles. Conversely, trans-acting factors bind to cis-regulatory elements of the two alleles, modulating the phenotype towards one parent, and thus producing non-additive (dominant-recessive) gene action [16, 17]. As mentioned above, recessive phenotypes are often associated with loss-of-function mutations. It is, therefore, possible that trans-acting factors highly connected in GRNs would display dominant recessive gene action less often to avoid antagonistic pleiotropy. It is also possible thus, that trans-acting factors showing dominant-recessive gene action would produce dominant-recessive gene action on their targets downstream. Some simulation-based studies have revealed that dominant-recessive gene action within GRNs is highly propagable [18]. The role of non-additive gene action on transcript abundance differences between wild and domesticated plants has been assessed in independent studies [9, 10, 19,20,21,22,23]. However, the use of different analytical approaches to both quantify transcript abundance and measure gene action has yielded contrasting results that are difficult to compare among them. Furthermore, besides a few studies on Drosophila and yeast that have assessed the association between GRN connectivity properties and transcript abundance differences, no studies have approached such associations in plant models.

In this study, we evaluated the contribution of gene action to transcript abundance divergence in plants evolving under domestication using species that exhibit gigantism of reproductive organs. These species include maize (Zea mays), chili pepper (Capsicum annuum) and sunflower (Helianthus annuus). Although they exhibit different morphology, human-mediated selection has resulted in the enlargement of their reproductive organs, such as flowers, fruits and seeds [2]. Because they were domesticated in what it is nowadays Mexico and the United States, they were likely exposed to a common set of cultural preferences and agricultural practices [1]. As control scenarios of selective pressures, we included data from crosses of two switchgrass (Panicum hallii) subspecies evolving under natural selection as well as two inbred lines of rice (Oryza sativa) and thale cress (Arabidopsis thaliana).

We hypothesized additive gene action on transcript abundance to be less frequent in the genetic crosses evolving under domestication than in the crosses evolving under natural selection. We particularly predicted this scenario for genes associated with cell cycle and cell fate determination, because they are directly associated with the gigantism of reproductive organs. Using a co-expression based GRN, we also investigated the association between both the effect size of transcript abundance differences and gene action to GRN connectivity properties (In-degree and Out-degree) in Arabidopsis thaliana F1 hybrids, to test for evidence of evolutionary constraint in the expression of genes with high Out-degree and high In-degree. In addition, to assess the role of transcriptional epistasis in propagating gene expression variability in GRNs, we analyzed the similarity of gene action between trans-acting regulators and their target genes. We observed an overall reduction of additive gene action on the domesticated transcript abundance phenotype, with a marked pattern for transcripts associated with cell cycle and cell fate determination. We also observed a negative correlation between GRN connectivity properties and both effect size of transcript abundance differences and gene action. Finally, we observed that dominant-recessive gene action is highly propagable through GRNs and that trans-acting factors showing additive gene action might be driving transgressive gene action on genes downstream.

Results

The genome-wide patterns of inheritance of transcript accumulation unravel an overall reduction of additive gene action in domesticated phenotypes

To investigate the genetic nature of phenotypic variation at the transcript abundance level in the context of plant domestication, we compiled a large data set that included genetic crosses evolving under natural selection and under domestication (Additional file 1: Table S1). Our data set includes the genome-wide transcriptomes of both F0s and F1 hybrids of two switchgrass (Panicum hallii) subspecies diverging under natural selection; of cultivated maize (Zea mays), chili pepper (Capsicum annuum), and sunflower (Helianthus annuus) genotypes as well as their wild progenitors; of two rice (Oryza sativa) inbred lines and two reciprocal crosses of thale cress (Arabidopsis thaliana). We were particularly interested in transcripts showing abundance differences between the F0s of each genetic cross, because as we see it, they are a proxy of organismal phenotypic variation. We observed a wide range (1–44%) in the portion of the transcriptome showing significant abundance differences between F0s (p-adjusted < 0.05, |Log2(FC)|> 0.5) (Additional file 1: Table S1). This difference was substantially larger in the comparisons between domesticated-wild F0s (average of 35%) compared to the switchgrass F0s (1%), which diverged under natural selection. This result could be partially explained to the nature of the grass itself, or by differences in sequencing coverage. However, a random sampling-based coverage analysis showed a uniform median transcript coverage of ~ 60X for all the libraries used in this study (Additional file 1: Figure S1), discarding this feature as a possible source of bias in our results. Therefore, this result is instead consistent with the polygenic nature of domestication-related phenotypes, such as the overall increase in organ size, which requires the coordinated expression of large constellations of developmentally regulated genes [2, 24]. We employed the degree of dominance (k) as a proxy of gene action to investigate the nature of genetic variation controlling transcript abundance variation across the twelve genetic crosses. The degree of dominance (k) quantitatively describes the extent to which the F1 hybrid phenotype deviates from a midpoint between the two F0 progenitors. Usually, k values near zero indicate purely additive gene action or equal contributions of the two alleles to the F1 hybrid phenotype. k values around |1| reflect dominant-recessive gene action, in which the phenotype produced by one allele masks the phenotypic expression of the other allele. While k values above |1.25| are interpreted as transgressive gene action. In this mode of gene action, the F1 hybrid displays a phenotype above (overdominance) or below (underdominance) the phenotypic values of its F0 progenitors. We hypothesized additive gene action on transcript abundance to be less frequent in the genetic crosses evolving under domestication than in the switchgrass genetic cross. Overall, our results show that the vast majority of the transcripts show k values within −5 and 5 (Fig. 1A), this range is consistent with k estimates on plant morphology and transcriptional phenotypes in maize and chili pepper, respectively [9, 25]. Supporting our hypothesis, the data of switchgrass show a symmetrical shape with a mode centered near to zero, within the range of additive gene action (|k|< 0.25). Whereas for the crosses of domesticated-wild plants and inbred–inbred genotypes, the distributions show non-symmetrical shapes with modes outside the additive gene action range but below the transgressive gene action range (|k|> 0.25, |k|< 1.25), suggesting a relevant role of partially and fully dominant-recessive allelic interactions in shaping transcript abundance under domestication (Fig. 1A).

Fig. 1
figure 1

Inheritance of transcript abundance unveils a relevant role of non-additive gene action in domesticated plants. A Distribution of the degree of dominance (k) across a series of genetic crosses including switchgrass, maize, chili pepper, sunflower, rice and Arabidopsis. |k|> 0.25 indicate non-additive inheritance. The color of circles indicates the category of the genetic cross. B Distribution of k for each category of genetic cross. Natural Selection: wild-wild; domestication: cultivated-wild; Inbreds: cultivated-cultivated. C Empirical cumulative distribution analysis revealed the switchgrass species evolving under natural selection are distinguishable from those evolving directional selection for both domestication (Kolmogorov–Smirnov: D = 0.21; p = 2.8E-14) and Inbreds (Kolmogorov–Smirnov: D = 0.28; p = 2.2E-16) genetic crosses

Full size image

In particular, we noticed that in the domesticated-wild crosses of chili pepper and sunflower, the domesticated phenotype tended to be recessive, although the opposite pattern was shown in maize–teosinte crosses [26]. Finally, transgressive gene action (|k|> 1.25) was a common feature for all the genetic crosses.

Data for k were aggregated into three discrete categories of genetic crosses, namely: natural selection, domestication, and inbreds, and their distributions were compared using a probability density function (Fig. 1B). In general, we observe a pattern of wider distributions in domestication and inbreds compared to natural selection, supporting our prediction of less of additive gene action under domestication.

To formally test this observation, the distributions were compared using a cumulative distribution approach (Fig. 1C). We observed a significantly larger portion of transcripts showing additive gene action (35% of total) in the cross of natural selection compared to the domestication and inbred crosses (25% and 20%; Fisher’s exact test: p = 2.905e-5, p = 5.996e-11, respectively).

By employing a pairwise Kolmogorov–Smirnov test to quantitatively compare the distributions, a substantial and significant difference was found for the three comparisons (Natural selection vs Domestication: p = 2.864e-14; Natural selection vs Inbreds: p = 2.2e-16; Domestication vs Inbreds: p = 2.2e-16). A similar result of less additive gene action in domesticated plants was observed by comparing gene action between cultivated maize and wild teosinte populations [27].

Cell cycle and cell fate determination-related genes tend to display non-additive gene action in genetic crosses involving domesticated plants

A gene ontology (GO) enrichment analysis was performed across the twelve genetic crosses to gain insight into the biological functions associated with transcriptional differences and mode of gene action. Instead of testing for the enrichment of GO terms in the sets of differentially abundant transcripts in each genetic cross, we tested whether a set of specific GOs, likely related to the domestication syndrome were enriched or present in the discrete categories of gene action. Overall, we observed that the majority of the GOs are represented by sets of transcripts showing non-additive gene action This supports the role of non-additive gene action in producing phenotypic divergence at the transcript abundance level, particularly in domesticated-wild genetic crosses. Among the most relevant ontologies showing this pattern, we found cell cycle (GO:0007049), cell fate determination (GO:0001709), flower development (GO:0009909), response to auxin (GO:0009733), and seed dormancy (GO:0009793), which could be related to the domestication syndrome (Additional file 1: Figure S2). We found weak enrichment for some of those categories, but this result was expected, because the number of genes in each category is a subsample of the total DE transcripts and hypergeometric tests are sensitive to small sample sizes [28]. We consider the presence–absence variation of GOs in each category of gene action provides relevant evidence of which biological functions were altered in each divergence type. These results highlight the role of non-additive gene action in modulating cellular and development-related phenotypes through transcriptional divergence under domestication. The genetic cross of switchgrass, evolving under natural selection, shows abundance differences only on transcripts associated with nitrogen metabolism (GO:0006807), which reflects local adaptation to the environmental variation between the two studied ecotypes [28].

The continuous distributions of the degree of dominance (k) were analyzed for genes associated with meristematic activity, such as the cell cycle, cell fate determination, flower development, auxin response, and seed dormancy, which are strongly related to the domestication syndrome. In the switchgrass cross, in addition to have found a small number of genes with differential expression performing these biological functions, these genes show values of k closer to zero, indicating additive gene action (Fig. 2).

Fig. 2
figure 2

Inheritance of meristem-related genes in cultivated plants suggest selection under domestication employs non-additive gene action. Jitter + violin plots show the distribution of the absolute degree of dominance (k) for genes annotated to each biological process. Labels denote specific genes with experimental evidence of their roles in each biological process

Full size image

Congruently, for the genes associated with the domestication syndrome, we observed that crosses involving domesticated and inbred plants showed higher k values compared to the cross evolving under natural selection. Furthermore, we found that the median k values for transcripts annotated with those biological functions showed gene action that is distinguishable from additive gene action (|k > 0.25|) (Fig. 2, Additional file 1: Table S2). This result suggests that genes associated with meristematic activities tend to display dominant-recessive or transgressive gene action. CDKD1, CDC22, EL2, and TSS were amongst the most relevant genes annotated as cell cycle regulators. They are part of the machinery that ensures the transitions between different cell cycle stages. For example, CDKD1 encodes for a cyclin-dependent kinase that controls the transition between the G2 and M phases of the cell cycle [29]. CDC22, a cofactor of the anaphase-promoting complex (APC) controls meristem size by regulating mitosis [30]. While EL2 inhibits the progression of the cell cycle by directly binding to CDKD1, thus having a role in the endoreduplication cell cycle and cell size [31]. Finally, TSS indirectly controls the cell cycle and organ size by modulating sugar availability, therefore, having a relevant role in meristem and organ size [32].

The pattern of non-additive gene action was also notable for transcripts associated with cell fate determination, flower development, response to auxin, and seed dormancy (Fig. 2). These biological functions are well known for having a relevant role in the domestication syndrome [8]. For instance, ZWIP2 is a zinc finger transcription factor that controls organ size by transcriptionally regulating the activity of fruit development-related genes, such as FUL, BP, and RPL [33]. ELF3 is implicated in the regulation of flowering time, a known feature of the domestication syndrome in some crops [34].

Auxin-mediated developmental transitions such as germination and growth are strongly associated with environmental cues. Amongst other genes associated with such functions, there are remarkable examples, such as PKL and TOR. PKL is a chromatin remodeling factor whose loss-of-function mutant shows delayed vegetative-to-reproductive transition and overall smaller vegetative and reproductive structures [35].

By repressing the embryonic state and modulating the seed cellular program into a seedling one, PKL also determines the transition between the embryonic and vegetative life stage. This process involves cell division, cell expansion, and the activation of photosynthetic machinery [36].

TOR is implicated in cell growth and cell division and possibly in seed dormancy [37]. These results highlight the role of non-additive gene action in modulating meristematic activity and ultimately producing the domestication syndrome characterized by gigantism in organs and delayed flowering times.

Trans-acting factors with high connectivity avoid antagonistic pleiotropy

Natural selection acts on phenotypic variability. Genetic variation in individual genes can produce phenotypic variability. However, cumulative evidence supports an omnigenic model that posits that it is the coordinated action of constellations of genes which ultimately produces phenotypic variation on continuous traits, such as organ size [38]. In such constellations or gene regulatory networks (GRNs), the trans-acting elements would produce the majority of the phenotypic variability via expression variation on genes downstream. Therefore, genetic variation at transcription factors and chromatin remodeling factors could be enough to modulate phenotypic variation at the whole GRN level. There is empirical evidence that suggests that to avoid antagonistic pleiotropy, trans-acting factors regulating a large number of genes (large Out-degree) tend to experience small expression changes [39]. Similarly, we reasoned that the number of regulators by which a gene is regulated (In-degree) could produce evolutionary constraint on its expression level, as mutations at some regulators could compensate for the phenotypic effects of other regulators. We hypothesized a negative correlation between the magnitude of expression variation and both, the Out-degree and the In-degree. To formally test this hypothesis, an A. thaliana co-expression-based GRN was employed. This GRN represents 330,775 regulatory interactions among 1,883 trans-acting factors (regulators) and 22,478 non-trans-acting genes (targets). From those interactions, 1,828 are self-regulatory. We then obtained the sub-network for a set of 3584 transcripts that display expression differences between the Col and the C24 A. thaliana accessions and that had detectable expression in the F1 Col-C24 hybrids. These accessions differ substantially in vegetative growth and reproductive yields [40]. This sub-GRN represents 16,226 regulatory associations between 443 regulators and 3,141 targets, including 416 self-regulatory interactions.

Overall, we found a weak but significant negative correlation between the magnitude of transcript abundance differences (|Log2(Col/C24)|) and both the Out-degree and the In-degree (ρ = −0.12, p = 0.02; ρ = −0.033, p = 0.04, respectively) (Fig. 3A left). These correlation values are similar to those found in a study employing different Drosophila species [39].

Fig. 3
figure 3

Analysis of the A. thaliana GRN showing connectivity properties of regulators and their targets. A Association between the Out-degree and the In-degree to the magnitude of divergence on transcript accumulation between A. thaliana Col-0 and C24 as well as their absolute degree of dominance (k) are shown as scatter plots. The line depicts the loess regression. B Violin plots showing the absolute Log2 foldchange in the degree of dominance of TFs and their target genes for each category of inheritance. Boxplot depict the median and quartile values. Values of |Log2(KREG/KTAR)| near to 0 indicate the TF and its target have highly similar values of k. C Gene regulatory network for genes with transcript accumulation divergence between A. thaliana Col-0 and C24. Edges are color coded according to the (|Log2(KREG/KTAR)|). Darker edges denote REG–TAR interaction with high similarity on their degrees of dominance. D Co-expression network for selected regulatory genes. Nodes are color-coded accordingly with the degree of dominance in a discrete fashion

Full size image

Recessive phenotypes are often the result of deleterious mutations [41]. Therefore, if we consider the non-additive mode of gene action as a proxy of dominant-recessive variation, it is feasible to predict a negative correlation between the Out-degree and the degree of dominance (k). For the case of the In-degree, genes with a large number of regulators could display stronger non-additive effects than their counterparts due to the cumulative effects of their regulators. We tested such hypotheses employing the aforementioned GRN and the degree of dominance estimates in the Col x C24 F1 hybrid. We found a weak but significant negative correlation between the degree of dominance (k) and the Out-degree (ρ = −0.08, p = 0.08). Conversely, the degree of dominance showed a weak but significant positive correlation with the In-degree (ρ = 0.073, p = 1.19e-5) (Fig. 3A right). Despite its relevance, to our knowledge, no previous studies have addressed the association between these two connectivity GRN properties and the degree of dominance.

Non-additive gene action propagates through GRNs

Under the scenario of GRNs being modulated by the activity of key trans-acting regulators, it is natural to think that different modes of gene action of trans-acting factors would have different phenotypic effects on the genes they regulate downstream. For example, if a trans-acting factor shows non-additive (dominant-recessive) expression in an F1 hybrid background, it would be expected for its target genes to show dominant-recessive gene action. This is because trans-acting factors bind to the two homologous chromosomes in diploid organisms, thus completely recapitulating the phenotype of one parent, unless the trans-acting factor is haploinsufficient. Conversely, for a trans-acting factor showing additive gene action (both functional alleles), it would be expected for its targets to show transgressive gene action, because they experience the cumulative trans-acting regulation of the two functional alleles. We tested these hypotheses by establishing a similarity index between the degree of dominance of trans-acting regulators and their targets (|Log2(KREG/KTAR)|) in our GRN.

In general, we observed that from the total 443 trans-acting regulators contained in the GRN, one third (33.8%) showed transgressive gene action, about one quarter showed dominant-recessive (24.8%) or partially dominant-recessive gene action (24.8% and 26.6%, respectively), and a smaller portion (14.6%) showed additive gene action. These data recapitulate the genome-wide pattern described above for this genetic cross (Fig. 1A), suggesting that trans-acting variation could be a relevant source of transcriptional variability between the Col and C24 inbred lines. We tested for differences in the similarity index among the four discrete categories of gene action of transcription factors. We found higher similarity in gene action among trans-acting factors and their targets when the trans-acting factor showed dominant-recessive gene action, followed by instances, where the trans-acting factor was partially dominant-recessive (median |Log2(KREG/KTAR)|= 0.63 and 1.03, respectively). These results suggest that trans-acting factors displaying these two modes of gene action tend to produce the same gene action on their targets. Whereas a lower similarity was observed when trans-acting factors showed transgressive gene action (median |Log2(KREG/KTAR)|= 1.28), and even lower for the trans-acting factors showing additive gene action (median |Log2(KTF/KTAR)|= 2.93) (Fig. 3B, C). These data indicate, on one hand, that dominant-recessive allelic interactions are more propagable through GRNs than transgressive gene action. In addition, on the other hand, trans-acting factors showing additive gene action tend to produce transgressive gene action on their targets, supporting our hypothesis of the role of trans-acting variation in heterotic phenotypes. A Wilcoxon rank-sum test in a pairwise fashion found significant differences in all comparisons. Comparisons that included transcription factors with additive gene action had larger effects sizes (Additional file 1: Table S3).

The propagation of gene action through GRNs was analyzed in the subnetworks of GRNs of PKL, TSS, ZWIP2, CDC22, and ELF3 (Fig. 3D). Notably, for PKL, which showed partially dominant-recessive gene action, only three out of its 24 target genes showed additive gene action (binomial test; p = 0.1). Unfortunately, the size of the other sub-networks is too small to be formally analyzed. These results suggest that genetic variation affecting the activity of trans-acting factors propagates to some degree.

Discussion

By comparing transcript abundance variation and its patterns of gene action between a single genetic cross of plants evolving under natural selection to those under domestication, we found that the latter exhibit such variation in a substantially larger portion of their transcriptomes and that this variation appears to be driven by LOF genetic mutations that produce non-additive gene action. This result provides evidence of the role of human mediated domestication in shaping gene expression patterns despite an overall small genome-wide genetic divergence compared to species evolving under natural selection. For instance, the switchgrass subspecies exhibit 1.2% of genome-wide nucleotide differences [42], while these values are 0.3%, 0.2% and 0.1% for the chili pepper, maize and sunflower accessions, respectively [43,44,45].

Despite transcriptional variability due to sampling at varied developmental timepoints or tissues, our normalized measure of gene action (k) allowed us to investigate the genetic nature of allelic variation influencing transcript abundance across different species. Indeed, our study revealed that a substantial portion of the transcripts that showed abundance differences between wild and domesticated plants display partially non-additive (dominant-recessive) gene action. This mode of gene action has been ignored in the majority of experiments assessing the inheritance of transcript accumulation in plants [10, 19,20,21, 23, 28, 46, 47, 48]. Perhaps, due to the use of discrete categories of modes of gene action based on differential expression between F0s and F1 hybrids. Our findings suggest a relevant role of non-additive gene action on modulating a polygenic architecture of domestication-related phenotypes.

The recessive patterns of transcript accumulation observed in the domesticated sunflower and chili pepper suggest that divergence due to selection under domestication is associated with loss-of-function genetic variation. The opposite pattern was shown in maize crosses, which indicates a possible role of gain-of-function genetic changes in the domestication of this crop. However, it is likely that maize-like alleles initially arose in a teosinte genetic background and behaved as loss-of-function. Therefore, the dominance of modern maize phenotypes that we observe in an F1 hybrid could differ from the initial conditions of maize domestication [49]. This finding supports a growing collection of studies that propose human-mediated selection for gigantism of organs, unconsciously selected for genetic variation that would not be selected in naturally evolving populations due to its deleterious effects on fitness [41].

The gigantism of vegetative and reproductive organs is a common feature of domesticated plants. Several studies have proposed a relevant role of loss-of-function variation in cell cycle and cell fate determination in producing such gigantism by altering meristem size [2, 8]. Congruently, our study revealed that transcript abundance of genes associated with cell cycle and cell fate determination showed dominant-recessive gene action. Empirical evidence indicates that genes involved in such functions belong to homologous GRNs, suggesting crop domestication is an outstanding example of parallel evolution. For instance, the modification of plant architecture and inflorescence via the modulation of the cell cycle GRN by the action of the TCP transcriptional repressor tb1 is common for maize, pearl millet, and barley [3]. In our study, both tb1 and transcripts annotated as cell cycle regulators showed strong non-additive effects in the crosses of maize, chili pepper, and sunflower. We could hypothesize that if GRNs are at least partially conserved between these species, the non-additive trans-regulatory activity of their tb1 homologs could contribute to alterations in the cell cycle machinery and meristem size. It would be interesting to experimentally test the aforementioned hypothesis by contrasting the meristem size of both gain and loss-of-function mutants of tb1 in different species.

Using an A. thaliana GRN and differential expression data of two accessions of this species that differ in meristem size, we found a weak negative correlation between the magnitude of expression divergence and both, the In-degree and the Out-degree. These findings suggest that trans-acting factors with high Out-degree avoid antagonistic pleiotropy by constraining their expression levels. Whereas the constraint of expression differences on genes with high In-degree could be due to the compensatory regulation of their multiple trans-acting regulators [39]. Our observation of a weak negative correlation between the degree of dominance (k) and the Out-degree has no previous record in the literature. We speculate, however, that given that mutations producing non-additive phenotypes often result in negative effects on fitness [41], natural selection would purge such deleterious mutations on highly connected trans-acting regulators to avoid antagonistic pleiotropy. It remains to be seen if our findings are supported by GRNs of other domesticated species, yet, they add to our knowledge about the interplay between GRNs and the nature of genetic changes controlling domestication-related phenotypes, such as organ size. Finally, we acknowledge the use of Arabidopsis microarray-based expression profiles to build the co-expression based GRN could be a relevant source of technical bias. Yet, this data set is a valuable resource for the plant biology community.

The nature of genetic changes that produce non-additive gene action within GRNs at two levels (dominant-recessive and transgressive), led us to consider two classical hypotheses that propose mechanistic explanations [17, 50]: (i) the dominance hypothesis posits that the deleterious allelic effects of one parent are complemented in the F1 hybrid by the functional allele of the other parent. (ii) The beneficial allele hypothesis predicts that two functional alleles would produce a transgressive phenotype in the F1 hybrid due to their complementary beneficial phenotypic effects. Our study provides relevant empirical evidence for these two hypotheses. First, in support of the dominance hypothesis, we found that trans-acting factors showing dominant-recessive gene action produce dominant-recessive effects on most of their target genes. In agreement with this finding, computational models propose that dominant-recessive gene action propagates through GRNs via trans-acting factors [18]. One possible explanation for this finding is that, in absence of cis-acting variation in the target genes, the functional trans-acting allele complements the non-functional trans-acting allele of the other parent. Second, in support of the beneficial allele hypothesis, our data revealed that targets genes of trans-acting factors that showed additive gene action tended to display transgressive phenotypes. This indicates that parental allelic trans-acting factors differ in their cis-acting regulatory elements, but the two alleles are fully functional in their respective genetic background, thus reinforcing their regulatory activity on their transcriptional targets. The beneficial allele hypothesis is an example of transcriptional epistasis. A classic example of transcriptional epistasis is the anthocyanin pathway in maize, where the additive interaction of allelic variants of the trans-acting regulator B, produces transgressive transcript abundance on the anthocyanin synthesis genes A1, A2, and Bz1 [51]. The catalog of 65 trans-acting regulators that showed additive inheritance in the Col-C24 A. thaliana cross is a relevant source of candidates to experimentally test the role of functional cis-acting allelic variants in producing transgressive gene action on meristem size. It remains to be seen if our findings are supported by GRNs of other species, yet, they add to our knowledge about the interplay between GRNs and the nature of genetic changes controlling domestication-related phenotypes, such as organ size.

Finally, our finding in A. thaliana of the propagation of phenotypic effects of trans-acting regulators to their targets adds empirical evidence to the omnigenic model. According to this model, variability in complex phenotypes such as organ size is the result of variations on the trans-acting activity of “peripheral” TFs on the expression of “core” genes in GRNs. Therefore, the domestication-related phenotypes are likely underlaid by a small number of mutations that work in an omnigenic fashion [38]. As large transcriptomic data sets are available to construct high-quality GRNs it will be possible to test our hypothesis in other plant species.

Conclusions

Our study provides evidence of the role of non-additive gene action in the modulation of transcript abundance of genes associated with domestication-related traits in plants, such as the gigantism of reproductive organs. The non-additive transcript abundance variation on genes associated with cell cycle and cell fate determination machinery are likely driven by loss-of-function mutations, in line with the classic observation of aberrant organ sizes in domesticated plants made by Darwin. One major limitation of this study is the use of a single genetic cross as natural selection control. It remains to be seen if the use of additional genetic crosses as tests and as controls, would shed light on whether non-additive patterns of gene expression play a major role in plant domestication. Using a co-expression-based GRN of the model plant A. thaliana, our data revealed a weak but significant negative correlation between the Out-degree and both the magnitude of transcript abundance variation a well the absolute degree of dominance. Which suggested highly connected genes in GRNs are less prone to be influenced by LOF genetic changes that result in non-additive gene action. If the Out-degree is a good measure of pleiotropy, LOF mutations would be expected to be preferentially eliminated by natural selection, providing a mechanism to avoid antagonistic pleiotropy. Finally, our study provides empirical evidence of the propagation of non-additive phenotypes in GRNs. This is congruent with the omnigenic model and highlights the role of trans-acting regulatory changes in modulating polygenic phenotypes such as fruit size in the evolution of plants under domestication.

Methods

Analysis design, data collection and preprocessing

We reasoned that F1 hybrids derived from crosses between cultivated plants and their wild progenitors would display a large portion of non-additive variation in their transcriptional landscapes compared to genetic crosses between lineages evolving under natural selection. Therefore, our data set includes one genetic cross of two switchgrass (Panicum hallii hallii x P. hallii filipes) ecotypes evolving under natural selection [28]; three genetic crosses of cultivated maize and its wild progenitors (Zea mays x Zea mays parviglumis) [22]; one genetic cross of cultivated and wild chili pepper (Capsicum annuum annuum x C. annuum glabriusculum). [9]; one genetic cross of cultivated and wild sunflower (Heliantus annuus x H. petiolaris) [19]; one reciprocal cross between two inbred lines of cultivated rice (Oriza sativa indica x O. sativa japonica) [52]; and two reciprocal crosses of the Arabidopsis thaliana Columbia to the A. thaliana C24 and A. thaliana Landsberg erecta ecotypes [53], which are domesticated to the laboratory conditions. Due to the availability of raw next generation sequencing NGS-based transcriptome data, reference transcriptomes and the experimental conditions from which RNA was obtained, data for specific genetic crosses were chosen and processed into a unified analytical custom pipeline (Additional file 1: Table S4). Raw NGS sequences were retrieved and converted to fastq format employing the fasterq-dump program of the SRA toolkit 2.9.1 [54] using default parameters. Reads were then pre-processed to remove adapters, and to trim and filter low quality bases using the fastp software [55] with default parameters.

Read mapping, transcript accumulation quantification and differential expression analyses

Genome-wide transcript abundance quantification for each triplet (Both parents + F1 hybrid) of sequencing libraries was estimated by pseudo-aligning the clean reads to their respective reference transcriptome using Kallisto [56] with parameters -b 100, -t 16 and -bias. Data for each species were then consolidated into a single data base that included both parents (F0) and their (F1) hybrid. Transcripts with a per-species median of read counts smaller than five were excluded for further analyses. Transcript abundance estimates were then employed to test for differential transcript abundance between F0 genotypes using a generalized linear model via the DESeq2 R library [57]. Transcripts were considered as differentially expressed using an FDR ≤ 0.05 and a |Log2 fold change|> 0.5.

Mode of gene action assessment

We estimated the degree of dominance (k) as a proxy of mode of gene action for all the transcripts showing evidence of both differential expression between F0 and identifiable expression in the F1 hybrids. Specifically, the degree of dominance (k = d/a) estimates the association between the dominant and additive effects. As we see it, the degree of dominance is a normalized continuous measurement that estimates the deviation of the null hypothesis of both parents contributing equally to the observed phenotype in the F1 hybrids. The dominance effect (d = F1—[(X1 + X2)/2]) measures the deviation of the F1 hybrid phenotype from a midpoint of the F0s. The additive effect (a = (X1−X2)/2) provides a quantitative estimate of the effect of substituting one allele, for example a wild allele for a mutant allele in a cultivated phenotype. Using a custom R function, we employed the normalized transcript abundance values and evaluated k for all the differentially expressed transcripts across 12 genetic crosses. We then sorted each transcript according to four discrete categories of gene action as follows: |k|< 0.25 = additive; |k|≥ 0.25 and |k|≤ 0.75 = partially dominant-recessive; |k|≥ 0.75 and |k|≤ 1.25 = dominant-recessive; |k|≥ 1.25 = transgressive. Beyond discrete categories of transcript abundance gene action, we analyzed the transcriptional landscapes employing probability distributions. This allowed us to compare the breath of allelic interactions in different contexts of selective pressures.

Genetic crosses were then split into Natural Selection, Domestication and Inbreds according to the publications from where data were retrieved (Additional file 1: Table S1). We then performed a Kolmogorov–Smirnov test in a pairwise fashion to asses if plants evolving under natural selection accumulate more additive variation than plants evolving under domestication. Data were then visualized using an empirical cumulative distribution function.

Gene ontology and Co-expression network construction

To obtain biological insight into the portion of the transcriptome showing differential expression and showing different modes of gene action, we employed GO enrichment and co-expression network analyses. To avoid possible annotation artefacts, for instance ambiguity in the protein domains and biological functions, we decided to generate a unified annotation for the set whole set of differentially expressed transcripts. We first blasted [58] the transcripts against the Swiss-Prot database [59] and then selected the best hit employing home-made R scripts. GO terms for each transcript were then retrieved from the UniProt database using the ID mapping function [60]. The rationale of our GO enrichment analysis was to test if genes in different categories of gene action in each cross are enriched for different biological functions. The topGO R package [61] was then used to perform the enrichment analysis employing the weight01 algorithm.

Gene co-expression networks allow for inferring groups of genes that act together to regulate cellular processes. While some widely used approaches such as the WGCNA [62] employ correlation as a gold standard measurement of gene co-expression, other algorithms such as AracNe-AP [63] have been shown to perform better to infer regulatory relationships such as those between transcription factor genes and non-transcription factor genes. Using the AracNe-AP algorithm and 79 public microarray-based transcriptomes across different organs and developmental stages [64], we thus inferred a co-expression-based gene regulatory network (GRN) for the model organism Arabidopsis thaliana. Regulatory associations between TFs and non-TF genes were computed using the mutual information measurement. The AracNe-AP algorithm takes a list of genes defined as TFs and the rest of the genes as possible targets and infers mutual information. TFs are connected by an edge in a network if the mutual information meets a threshold value estimated from random sampling of mutual information among TFs and non-TF genes. We then obtained the specific co-expression or regulatory network for the genes showing differential gene expression between the A. thaliana Col-C24 F0s.

Association between gene regulatory network properties, transcript abundance differences and gene action

Some TFs have thousands of transcriptional targets downstream, it has been reported that there is a negative correlation between the effect size of transcriptional variation (ES) at TFs and the number of genes they regulate (Out-degree). In a similar way, the number of TFs that regulate a single gene (In-degree) has been negatively be correlated with the ES. We thus employed a Pearson correlation analysis and a Loess non-linear regression to investigate the possible association between ES and both the Outdegree and the Indegree.

We were also interested in assessing the association between both Indegree and Outdegree to the degree of dominance (k). We predicted that TFs with high Outdegree would tend to show additive gene action (|k|~ 0). We have no clear prediction for the inheritance of genes with high Indegree. We employed correlation and non-linear regression analyses to investigate the association between those network metrics and the degree of dominance (k).

Gene action propagation in gene regulatory networks

We reasoned that TFs and their targets would display similar degrees of dominance (k), particularly for TFs showing partial or complete dominant-recessive gene action. We thus employed the absolute values of the Log2 fold change between the k values of TFs and their targets (|Log2(KREG/KTAR)|), to generate a normalized measurement of degree of dominance similarity. In addition to visualize the distributions of k similarity among the four categories of gene action of TFs, we explicitly compared the median similarity of TFs showing additive gene action to those showing non-additive gene action using a Wilcoxon rank sum test using the wilcox_test of the Rstatix R package [65]. We then mapped the k similarity values as a color-coded edge attribute and visualized the co-expression networks for each instances of gene action of TFs using the ggraph R package [66]. To visually analyze the similarity of degree of dominance between TFs and their targets, networks for a set of TFs annotated as cell cycle, cell fate and response to auxin regulators were plotted. All data visualizations and data wrangling, unless specified, were done using the Tidyverse R package [67] in the Rstudio programming environment [68, 69].

Availability of data and materials

The data sets generated and analyzed during the current study are available from the figshare repository: https://doi.org/10.6084/m9.figshare.19358249.v1

Abbreviations

APC:

Anaphase-promoting complex

CDC22:

Cell division cycle 20.2

CDKD1:

Cyclin-dependent kinase D1

DE:

Differentially expressed

ELF3:

Early flowering protein

ES:

Effect size

FC:

Fold change

FDR:

False discovery rate

GO:

Gene ontology

GRN:

Gene regulatory network

NGS:

Next generation sequencing

PKL:

Pickle remodeling factor

TOR:

Target of rapamycin

ZWIP2:

Zink finger protein WIP2

References

  1. Purugganan MD, Fuller DQ. The nature of selection during plant domestication. Nature. 2009;457:843–8. https://doi.org/10.1038/nature07895.

    Article  CAS  PubMed  Google Scholar 

  2. Doebley JF, Gaut BS, Smith BD. The molecular genetics of crop domestication. Cell. 2006;127:1309–21.

    Article  CAS  PubMed  Google Scholar 

  3. Purugganan MD. Evolutionary insights into the nature of plant domestication. Curr Biol. 2019;29:R705–14.

    Article  CAS  PubMed  Google Scholar 

  4. Lester RN. Evolution under domestication involving disturbance of genic balance. Euphytica. 1989;44:125–32.

    Article  Google Scholar 

  5. Alpert KB, Grandillo S, Tanksley SD. fw 2.2:a major QTL controlling fruit weight is common to both red- and green-fruited tomato species. Theor Appl Genet. 1995;91:994–1000.

    Article  CAS  PubMed  Google Scholar 

  6. Goodman RM. Encyclopedia of plant and crop science. London: Routledge; 2004. p. 1069–74.

    Book  Google Scholar 

  7. Pickersgill B. Domestication of plants in the Americas: insights from Mendelian and molecular genetics. Ann Bot. 2007;100:925–40.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Meyer RS, Purugganan MD. Evolution of crop species: genetics of domestication and diversification. Nat Rev Genet. 2013;14:840–52. https://doi.org/10.1038/nrg3605.

    Article  CAS  PubMed  Google Scholar 

  9. Diaz-Valenzuela E, Sawers RH, Cibrian-Jaramillo A. Cis- and trans-regulatory variations in the domestication of the chili pepper fruit. Mol Biol Evol. 2020;37:1593–603. https://doi.org/10.1093/molbev/msaa027.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Escoto-Sandoval C, Ochoa-Alejo N, Martínez O. Inheritance of gene expression throughout fruit development in chili pepper. Sci Rep. 2021;11:1–16. https://doi.org/10.1038/s41598-021-02151-z.

    Article  CAS  Google Scholar 

  11. Wittkopp PJ, Kalay G. Cis-regulatory elements: molecular mechanisms and evolutionary processes underlying divergence. Nat Rev Genet. 2012;13:59–69. https://doi.org/10.1038/nrg3095.

    Article  CAS  Google Scholar 

  12. Hill MS, Vande Zande P, Wittkopp PJ. Molecular and evolutionary processes generating variation in gene expression. Nat Rev Genet. 2020;22(4):203–15. https://doi.org/10.1038/s41576-020-00304-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Cong B, Tanksley SD. FW22 and cell cycle control in developing tomato fruit: a possible example of gene co-option in the evolution of a novel organ. Plant Mol Biol. 2006;62:867–80.

    Article  CAS  PubMed  Google Scholar 

  14. Siegal ML, Promislow DEL, Bergman A. Functional and evolutionary inference in gene networks: does topology matter? Genetica. 2007;129:83–103.

    Article  CAS  PubMed  Google Scholar 

  15. Batada NN, Hurst LD. Evolution of chromosome organization driven by selection for reduced gene expression noise. Nat Genet. 2007;39:945–9. https://doi.org/10.1038/ng2071.

    Article  CAS  PubMed  Google Scholar 

  16. Lemos B, Araripe LO, Fontanillas P, Hartl DL. Dominance and the evolutionary accumulation of cis- and trans-effects on gene expression. Proc Natl Acad Sci. 2008;105:14471–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Springer NM, Stupar RM. Allelic variation and heterosis in maize: How do two halves make more than a whole? Genome Res. 2007;17:264–75. https://doi.org/10.1101/GR.5347007.

    Article  CAS  PubMed  Google Scholar 

  18. Porter AH, Johnson NA, Tulchinsky AY. A new mechanism for mendelian dominance in regulatory genetic pathways: competitive binding by transcription factors. Genetics. 2017;205:101–12. https://doi.org/10.1534/GENETICS.116.195255/-/DC1.

    Article  PubMed  Google Scholar 

  19. Rowe HC, Rieseberg LH. Genome-scale transcriptional analyses of first-generation interspecific sunflower hybrids reveals broad regulatory compatibility. BMC Genomics. 2013;14:342. https://doi.org/10.1186/1471-2164-14-342.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Bell GDM, Kane NC, Rieseberg LH, Adams KL. RNA-seq analysis of allele-specific expression, hybrid effects, and regulatory divergence in hybrids compared with their parents from natural populations. Genome Biol Evol. 2013;5:1309–23.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Combes MC, Hueber Y, Dereeper A, Rialle S, Herrera JC, Lashermes P. Regulatory divergence between parental alleles determines gene expression patterns in hybrids. Genome Biol Evol. 2015;7:1110–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Lemmon ZH, Bukowski R, Sun Q, Doebley JF. The role of cis regulatory evolution in maize domestication. PLOS Genet. 2014;10:e1004745. https://doi.org/10.1371/JOURNAL.PGEN.1004745.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Haas M, Himmelbach A, Mascher M. The contribution of cis-and trans-acting variants to gene regulation in wild and domesticated barley under cold stress and control conditions. J Exp Bot. 2020;71:2573–84. https://doi.org/10.1093/jxb/eraa036.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Studer AJ, Wang H, Doebley JF. Selection during maize domestication targeted a gene network controlling plant and inflorescence architecture. Genetics. 2017;207:755–65. https://doi.org/10.1534/genetics.117.300071.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Yang J, Mezmouk S, Baumgarten A, Buckler ES, Guill KE, McMullen MD, et al. Incomplete dominance of deleterious alleles contributes substantially to trait variation and heterosis in maize. PLoS Genet. 2017;13:1–21.

    Article  Google Scholar 

  26. Lemmon ZH, Bukowski R, Sun Q, Doebley JF, Tanabe M. The role of cis regulatory evolution in maize domestication. PLoS Genet. 2014;10:e1004745. https://doi.org/10.1371/journal.pgen.1004745.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Chen Q, Samayoa LF, Yang CJ, Bradbury PJ, Olukolu BA, Neumeyer MA, et al. The genetic architecture of the maize progenitor, teosinte, and how it was altered during maize domestication. PLoS Genet. 2020;16:1–21. https://doi.org/10.1371/journal.pgen.1008791.

    Article  CAS  Google Scholar 

  28. Lovell JT, Schwartz S, Lowry DB, Shakirov EV, Bonnette JE, Weng X, et al. Drought responsive gene expression regulatory divergence between upland and lowland ecotypes of a perennial C4 grass. Genome Res. 2016;26:510–8. https://doi.org/10.1101/gr.198135.115.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Takatsuka H, Umeda-Hara C, Umeda M. Cyclin-dependent kinase-activating kinases CDKD;1 and CDKD;3 are essential for preserving mitotic activity in Arabidopsis thaliana. Plant J. 2015;82:1004–17.

    Article  CAS  PubMed  Google Scholar 

  30. Kevei Z, Baloban M, Da Ines O, Tiricz H, Kroll A, Regulski K, et al. Conserved CDC20 cell cycle functions are carried out by two of the five isoforms in arabidopsis thaliana. PLoS ONE. 2011. https://doi.org/10.1371/journal.pone.0020618.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Peres A, Churchman ML, Hariharan S, Himanen K, Verkest A, Vandepoele K, et al. Novel plant-specific cyclin-dependent kinase inhibitors induced by biotic and abiotic stresses *. J Biol Chem. 2007;282:25588–96. https://doi.org/10.1074/JBC.M703326200.

    Article  CAS  PubMed  Google Scholar 

  32. Skylar A, Sung F, Hong F, Chory J, Wu X. Metabolic sugar signal promotes Arabidopsis meristematic proliferation via G2. Dev Biol. 2011;351:82–9. https://doi.org/10.1016/j.ydbio.2010.12.019.

    Article  CAS  PubMed  Google Scholar 

  33. Marsch-Martínez N, Zúñiga-Mayo VM, Herrera-Ubaldo H, Ouwerkerk PBF, Pablo-Villa J, Lozano-Sotomayor P, et al. The NTT transcription factor promotes replum development in Arabidopsis fruits. Plant J. 2014;80:69–81.

    Article  PubMed  Google Scholar 

  34. Herrero E, Kolmos E, Bujdoso N, Yuan Y, Wang M, Berns MC, et al. EARLY FLOWERING4 recruitment of EARLY FLOWERING3 in the nucleus sustains the Arabidopsis circadian clock. Plant Cell. 2012;24:428–43. https://doi.org/10.1105/TPC.111.093807.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Li HC, Chuang K, Henderson JT, Rider SD, Bai Y, Zhang H, et al. PICKLE acts during germination to repress expression of embryonic traits. Plant J. 2005;44:1010–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Henderson JT, Li HC, Rider SD, Mordhorst AP, Romero-Severson J, Cheng JC, et al. PICKLE acts throughout the plant to repress expression of embryonic traits and may play a role in gibberellin-dependent responses. Plant Physiol. 2004;134:995–1005. https://doi.org/10.1104/PP.103.030148.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Deprost D, Yao L, Sormani R, Moreau M, Leterreux G, Bedu M, et al. The Arabidopsis TOR kinase links plant growth, yield, stress resistance and mRNA translation. EMBO Rep. 2007;8:864–70. https://doi.org/10.1038/SJ.EMBOR.7401043.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Liu X, Li YI, Pritchard JK. Trans effects on gene expression can drive omnigenic inheritance. Cell. 2019;177:1022-1034.e6. https://doi.org/10.1016/j.cell.2019.04.014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Yang B, Wittkopp PJ. Structure of the transcriptional regulatory network correlates with regulatory divergence in drosophila. Mol Biol Evol. 2017;34:1352–62. https://doi.org/10.1093/MOLBEV/MSX068.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Groszmann M, Gonzalez-Bayon R, Greaves IK, Wang L, Huen AK, James Peacock W, et al. Intraspecific arabidopsis hybrids show different patterns of heterosis despite the close relatedness of the parental genomes. Plant Physiol. 2014;166:265–80. https://doi.org/10.1104/pp.114.243998.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Huber CD, Durvasula A, Hancock AM, Lohmueller KE. Gene expression drives the evolution of dominance. Nat Commun. 2018;9:1–11. https://doi.org/10.1038/s41467-018-05281-7.

    Article  CAS  Google Scholar 

  42. Gould BA, Palacio-Mejia JD, Jenkins J, Mamidi S, Barry K, Schmutz J, et al. Population genomics and climate adaptation of a C4 perennial grass, Panicum hallii (Poaceae). BMC Genomics. 2018;19:1–11. https://doi.org/10.1186/S12864-018-5179-7/FIGURES/4.

    Article  Google Scholar 

  43. Qin C, Yu C, Shen Y, Fang X, Chen L, Min J, et al. Whole-genome sequencing of cultivated and wild peppers provides insights into Capsicum domestication and specialization. Proc Natl Acad Sci U S A. 2014;111:5135–40. https://doi.org/10.1073/pnas.1400975111.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Wang J, Li X, Do Kim K, Scanlon MJ, Jackson SA, Springer NM, et al. Genome-wide nucleotide patterns and potential mechanisms of genome divergence following domestication in maize and soybean. Genome Biol. 2019;20:1–16. https://doi.org/10.1186/S13059-019-1683-6/FIGURES/6.

    Article  Google Scholar 

  45. Hübner S, Bercovich N, Todesco M, Mandel JR, Odenheimer J, Ziegler E, et al. Sunflower pan-genome analysis shows that hybridization altered gene content and disease resistance. Nat Plants. 2018;5:54–62. https://doi.org/10.1038/s41477-018-0329-0.

    Article  CAS  PubMed  Google Scholar 

  46. Zhang C, Lin C, Fu F, Zhong X, Peng B, Yan H, et al. Comparative transcriptome analysis of flower heterosis in two soybean F1 hybrids by RNA-seq. PLoS ONE. 2017. https://doi.org/10.1371/journal.pone.0181061.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Corrinne E. Grover, Mi-Jeong Yoo, Meng Lin, Matthew D. Murphy, David B. Harker, Robert L. Byers, et al. Genetic Analysis of the Transition from Wild to Domesticated Cotton (Gossypium hirsutum L.). Genetics. 2020;10:731–54. https://doi.org/10.1534/g3.119.400909.

  48. Zheng Y, Wang P, Chen X, Sun Y, Yue C, Ye N. Transcriptome and metabolite profiling reveal novel insights into volatile heterosis in the tea plant (Camellia Sinensis). Molecules. 2019;24:3380. https://doi.org/10.3390/MOLECULES24183380.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Stitzer MC, Ross-Ibarra J. Maize domestication and gene interaction. New Phytol. 2018. https://doi.org/10.1111/nph.15350.

    Article  PubMed  Google Scholar 

  50. Mackay IJ, Cockram J, Howell P, Powell W. Understanding the classics: the unifying concepts of transgressive segregation, inbreeding depression and heterosis and their central relevance for crop breeding. Plant Biotechnol J. 2021;19:26–34.

    Article  PubMed  Google Scholar 

  51. Dooner HK, Robbins TP, Jorgensen RA. Genetic and developmental control of anthocy an in biosynthesis. 1991. www.annualreviews.org.

  52. Chodavarapu RK, Feng S, Ding B, Simon SA, Lopez D, Jia Y, et al. Transcriptome and methylome interactions in rice hybrids. Proc Natl Acad Sci U S A. 2012;109:12040–5. https://doi.org/10.1073/PNAS.1209297109/-/DCSUPPLEMENTAL.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Groszmann M, Gonzalez-Bayon R, Lyons RL, Greaves IK, Kazan K, Peacock WJ, et al. Hormone-regulated defense and stress response networks contribute to heterosis in Arabidopsis F1 hybrids. Proc Natl Acad Sci U S A. 2015;112:E6397–406. https://doi.org/10.1073/PNAS.1519926112/-/DCSUPPLEMENTAL.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. NCBI. GitHub - ncbi/sra-tools: SRA Tools. 2021. https://github.com/ncbi/sra-tools. Accessed 24 Feb 2022.

  55. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90. https://doi.org/10.1093/BIOINFORMATICS/BTY560.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34:525–7. https://doi.org/10.1038/nbt.3519.

    Article  CAS  PubMed  Google Scholar 

  57. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1–21.

    Article  Google Scholar 

  58. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    Article  CAS  PubMed  Google Scholar 

  59. Bateman A. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 2019;47:D506–15.

    Article  Google Scholar 

  60. Retrieve/ID mapping. https://www.uniprot.org/uploadlists/. Accessed 24 Feb 2022.

  61. Alexa A, Rahnenführer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006;22:1600–7. https://doi.org/10.1093/BIOINFORMATICS/BTL140.

    Article  CAS  PubMed  Google Scholar 

  62. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:1–13. https://doi.org/10.1186/1471-2105-9-559/FIGURES/4.

    Article  Google Scholar 

  63. Lachmann A, Giorgi FM, Lopez G, Califano A. ARACNe-AP: gene network reverse engineering through adaptive partitioning inference of mutual information. Bioinformatics. 2016;32:2233–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Schmid M, Davison TS, Henz SR, Pape UJ, Demar M, Vingron M, et al. A gene expression map of Arabidopsis thaliana development. Nat Genet. 2005. https://doi.org/10.1038/ng1543.

    Article  PubMed  Google Scholar 

  65. Kassambara A. Package “rstatix.” 2020. https://cran.r-project.org/web/packages/rstatix/index.html. Accessed 24 Feb 2022.

  66. Pedersen T. Package ‘ggraph.’ 2021. https://cran.r-project.org/web/packages/ggraph/index.html. Accessed 24 Feb 2022.

  67. Wickham H, Averick M, Bryan J, Chang W, DMcgowan AL, et al. Welcome to the Tidyverse. J Open Source Softw. 2019;4:1686. https://doi.org/10.21105/JOSS.01686.

    Article  Google Scholar 

  68. RStudio Team. RStudio | Open source & professional software for data science teams - RStudio. RStudio Inc. 2020. https://www.rstudio.com/. Accessed 15 Mar 2022.

  69. Wilson A, Norden N. The R project for statistical computing. 2015. URL: http://www.r-project.org/254. Accessed 15 Mar 2022.

Download references

Acknowledgements

E.D-V. is grateful to Dulce Valdivia for her help to improve some R scripts and for her time to discuss the statistical approaches.

Funding

E.D-V. is grateful to CONACYT for his PhD scholarship. Endogenomiks Mexico support this research.

Author information

Authors and Affiliations

  1. National Laboratory for Genomics of Biodiversity (UGA-Langebio), CINVESTAV, Irapuato, Guanajuato, Mexico

    Erik Díaz-Valenzuela, Daniel Hernández-Ríos & Angélica Cibrián-Jaramillo

  2. Monterrey Institute of Technology and Higher Education, Monterrey, Nuevo León, Mexico

    Daniel Hernández-Ríos

Authors
  1. Erik Díaz-Valenzuela
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Daniel Hernández-Ríos
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Angélica Cibrián-Jaramillo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

EDV and ACJ conceived and designed the project. EDV and DHR performed the computational analyses. EDV wrote the article. ACJ critically revised the article. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Erik Díaz-Valenzuela or Angélica Cibrián-Jaramillo.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

Not applicable.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1

. GO enrichment analysis showing biological functions altered in the studied genetic crosses. Dot plots show the enrichment, or presence of each GO for each mode of gene action (Additive, Partially-Dominant, Dominant, Transgressive) across the array of genetic crosses. Color of dots indicate the significance of the enrichment; the size depicts the number of annotated genes to each GO. Figure S2. A sampling-based analysis revealed a uniform ~60X sequencing coverage for all the libraries used in the study. The per-transcript sequencing coverage (Log2 normalized read counts) is shown as boxplots for the three genotypes of each species. The line within the boxes depicts the median, whereas the dots outside the box represent outlier values. Table S1. Summary of the per-species analysis of differential expression between parental and F1 hybrids. Table S2. Transcripts associated to domestication-related funcitons showing non-additive gene action. Table S3. Summary of the pair-wise differences between discrete categories of gene action of transcription factors and their target genes. Table S4. Details of the per-species RNA-seq libraries used in the study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Díaz-Valenzuela, E., Hernández-Ríos, D. & Cibrián-Jaramillo, A. The role of non-additive gene action on gene expression variation in plant domestication. EvoDevo 14, 3 (2023). https://doi.org/10.1186/s13227-022-00206-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13227-022-00206-4

Keywords

EvoDevo

ISSN: 2041-9139