Original Research ARTICLE
A Co-Expression Network in Hexaploid Wheat Reveals Mostly Balanced Expression and Lack of Significant Gene Loss of Homeologous Meiotic Genes Upon Polyploidization
- 1John Innes Centre, Norwich Research Park, Norwich, United Kingdom
- 2School of Biosciences, University of Birmingham, Birmingham, United Kingdom
- 3Computational and Analytical Sciences, Rothamsted Research, Harpenden, United Kingdom
Polyploidization has played an important role in plant evolution. However, upon polyploidization, the process of meiosis must adapt to ensure the proper segregation of increased numbers of chromosomes to produce balanced gametes. It has been suggested that meiotic gene (MG) duplicates return to a single copy following whole genome duplication to stabilize the polyploid genome. Therefore, upon the polyploidization of wheat, a hexaploid species with three related (homeologous) genomes, the stabilization process may have involved rapid changes in content and expression of MGs on homeologous chromosomes (homeologs). To examine this hypothesis, sets of candidate MGs were identified in wheat using co-expression network analysis and orthology informed approaches. In total, 130 RNA-Seq samples from a range of tissues including wheat meiotic anthers were used to define co-expressed modules of genes. Three modules were significantly correlated with meiotic tissue samples but not with other tissue types. These modules were enriched for GO terms related to cell cycle, DNA replication, and chromatin modification and contained orthologs of known MGs. Overall, 74.4% of genes within these meiosis-related modules had three homeologous copies which was similar to other tissue-related modules. Amongst wheat MGs identified by orthology, rather than co-expression, the majority (93.7%) were either retained in hexaploid wheat at the same number of copies (78.4%) or increased in copy number (15.3%) compared to ancestral wheat species. Furthermore, genes within meiosis-related modules showed more balanced expression levels between homeologs than genes in non-meiosis-related modules. Taken together, our results do not support extensive gene loss nor changes in homeolog expression of MGs upon wheat polyploidization. The construction of the MG co-expression network allowed identification of hub genes and provided key targets for future studies.
Meiosis is a specialized mode of cell division which generates haploid gametes. Prior to meiosis, chromosomes are replicated. On entry into meiosis, homologous chromosomes (homologs) locate each other and intimately align (synapse) along their length. Within this paired structure, chromosomes recombine and crossover before being accurately segregated (Kleckner, 1996; Mercier et al., 2015; Zickler and Kleckner, 2015). This complex and dynamic process is essential to maintain genome stability and integrity over sexual life cycles and to generate genome variation, which is a major evolutionary driving force (Capilla et al., 2016; Melamed-Bessudo et al., 2016). The genetic variation created by meiotic recombination underpins plant breeding to improve crop species (Wijnker and de Jong, 2008; Choi, 2017; Fernandes et al., 2018). Polyploidization has played an important role in the evolution and speciation of flowering plants (Comai, 2005; Alix et al., 2017), although the resultant multiplicity of related genomes poses a major challenge for the meiotic process. Segregation of the chromosomes to produce balanced gametes requires correct pairing, synapsis, and recombination between only true homologs, rather than any of the other highly related chromosomes (homeologs) (Ramsey and Schemske, 2002; Comai, 2005; Stenberg and Saura, 2013).
In the last two decades, there have been significant advances in our understanding of plant meiosis. Since the isolation of the ﬁrst meiosis-speciﬁc cDNA from lily in the mid-1990s (Kobayashi et al., 1993; Kobayashi et al., 1994), more than 110 plant meiotic genes (MGs) have been identified, mainly from studies of the model diploid plants Arabidopsis and rice (Mercier and Grelon, 2008; Luo et al., 2014; Mercier et al., 2015). Although 25–30% of flowering plants are extant polyploids (Alix et al., 2017), the meiotic mechanisms responsible for their stabilization remain poorly understood. An exception is hexaploid wheat (Triticum aestivum L.), where there is now better understanding of these processes (Blary and Jenczewski, 2019). Despite possessing multiple related genomes, durum wheat, a tetraploid (AABB), and bread wheat, a hexaploid (AABBDD) behave as diploids during meiosis. Thus, most of the meiotic studies conducted in hexaploid wheat have focused on providing better understanding of the meiotic processes required to stabilize this polyploid species (Riley and Chapman, 1958; Riley, 1960; Holm, 1986; Martinez-Perez et al., 2000; Martinez-Perez et al., 2003). An emphasis has been to characterize the role of the Ph1 locus in the suppression of recombination between homeologs (Riley et al., 1968; Holm, 1988; Holm and Wang, 1988; Feldman, 1993; Aragon-Alcaide et al., 1997; Martinez-Perez et al., 2001; Prieto et al., 2004; Colas et al., 2008; Martín et al., 2017). Recent studies have defined this phenotype to a ZIP4 gene which duplicated and diverged on polyploidization (Martín et al., 2017; Rey et al., 2017; Martín et al., 2018; Rey et al., 2018). This event resulted in the suppression of homeologous crossover, and promotion of homologous synapsis.
Although all flowering plants have undergone at least one event of whole genome duplication during their evolutionary history (Jiao et al., 2011), it has been suggested that MG duplicates return to a single copy following whole genome duplication, more rapidly than the genome-wide average (Lloyd et al., 2014). Therefore, it has been assumed that the stabilization process upon the polyploidization of wheat also involved rapid changes in the content and expression of the genes on homeologs. This process would facilitate the correct pairing and synapsis of homeologs during meiosis. The recent development of an expression atlas for hexaploid wheat revealed that 70% of homeologous genes in syntenic triads showed balanced expression (Ramírez-González et al., 2018); however, this study did not include analysis of the genes expressed during meiosis.
Here, we assessed whether the level of expression of all genes in triads was balanced between homeologs during meiosis. Analysis indicated similar balanced expression to that observed in other wheat tissues. However, it could be argued that only meiotic specific genes might show differential expression between homeologs. Sets of candidate MGs were identified using co-expression network analysis and orthology informed approaches, allowing us to evaluate the effect of polyploidization on wheat MG copy number and expression. The combination of co-expression network analysis, in conjunction with orthologue information, will now contribute to the discovery of new MGs and greatly empower reverse genetics approaches (such as wheat TILLING and CRISPR) that can be used to validate the function of the identified candidate genes in wheat.
Results and Discussion
An initial assessment of the homeolog expression pattern in triads during meiosis in hexaploid wheat was undertaken. Relative expression abundance of 19,801 triads (59,403 genes) was calculated for 8 tissues, including meiotic anther tissue, according to published criteria (Ramírez-González et al., 2018). This analysis revealed that the percentage of balanced triads was slightly higher in meiotic anther tissue (77.3%) than in other type of tissues (ranging from 67.3% in floral organs and 76.6% in leaves) (Figure S1). The copy number of genes expressed during meiosis was also investigated. This involved the definition of 19,801 triads (59,403 genes), 7,565 duplets (15,130 genes), 15,109 monads (single-copy genes), and 18,250 genes from the “others” group with various copy numbers, based on the Ensembl Plants database for the high confidence (HC) genes of hexaploid wheat (International Wheat Genome Sequencing Consortium, 2018) (IWGSC v1.1 gene annotation; Table S1). Comparison of copy number of genes expressed in the eight different tissues showed that 70.9% of the genes expressed during meiosis belonged to triads. This percentage ranged between 66.5 and 72.5% for the genes expressed in floral organs and stem tissues, respectively (Figure S2). The overall results were consistent with a previous study, which analyzed a range of tissues, but not meiotic anther tissue, and reported significant balanced expression between homeologous genes in these tissues (Ramírez-González et al., 2018). However, our analysis revealed a slightly higher percentage of balanced triads in meiotic anther tissue than in other types of tissues. These observations did not support the hypothesis that stabilization of polyploidization in wheat involved significant changes in gene content and expression between homeologs (Freeling, 2009; Lloyd et al., 2014; Edger et al., 2017). Considering that not all genes expressed in meiotic anther tissue are directly involved in the meiosis process, it is possible that meiotic specific genes exhibit a different pattern. Therefore, a co-expression gene network was developed to compare the expression pattern of homeologous genes in meiosis-related modules, which potentially represent meiosis specific genes, and other tissue-related modules.
Weighted Co-Expression Network Construction
Network-based approaches have been proved useful in systems biology, to mine gene function from high-throughput gene expression data. Gene co-expression analysis has become a powerful tool to build transcriptional networks of genes involved in common biological events in plants (Usadel et al., 2009; He and Maslov, 2016; Krishnan et al., 2017; Ma et al., 2018; Yu et al., 2018; Liu et al., 2019). The use of co-expression networks has uncovered candidate genes to regulate biological processes in many plants including wheat (International Wheat Genome Sequencing Consortium, 2018), rice (Aya et al., 2011; Tan et al., 2017), and Arabidopsis thaliana (Costa et al., 2015; Silva et al., 2016). The recently published high-quality genome reference sequence (International Wheat Genome Sequencing Consortium, 2018) and a developmental gene expression atlas (Borrill et al., 2016; Ramírez-González et al., 2018), together with the gene expression data collected from meiotic samples, were used to build a co-expression gene network. One hundred and thirty samples from different tissue types were included in this co-expression analysis (Table S2; Figure 1A). A set of 60,379 genes out of the total 107,892 HC genes was considered expressed during meiosis [transcript per million (TPM) > 0.5 in at least one meiosis sample (biological replicate)] and used to run the co-expression analysis. Using the “WGCNA” package in R (Langfelder and Horvath, 2008; Langfelder and Horvath, 2012), genes with similar expression patterns were grouped into modules via the average linkage hierarchical clustering of normalized count expression values (Figure 1E). The power of β = 7 (scale free topology R2 = 0.91) was selected as the soft threshold power to emphasize strong correlations between genes and penalize weak correlations to ensure a scale-free network (Figures 1B, C). Based on this analysis, 50,387 out of 60,379 genes (83.5% of expressed genes) could be assigned to 66 modules. Module size ranged from 52 to 7,541 genes (mean 763 genes; median 429 genes) (Figure 1D). The expression patterns of all genes within a single module were summarized into a module eigengene (ME; representative gene of the module) to minimize data size for subsequent analyses. Expression patterns of modules are shown as a heatmap by plotting ME values in relation to tissue samples (Figure 2).
Figure 1 The weighted gene co-expression networks analysis (WGCNA). (A) Clustering dendrogram of 130 samples from different types of tissues. The sample clustering was based on the expression data of the genes expressed during meiosis. (B) Analysis of the scale-free fit index for various soft-thresholding powers (β). A scale-free network is a network with the property that the number of connections k originating from a given node exhibits a power law distribution P (k) ∼ k−ɤ (where 2 < ɤ < 3). The red line indicates the scale-free topology fit index (of 0.9) for which the network obeys scale free network property. The soft power threshold used in constructing the weighted gene co-expression networks was chosen as the first power to exceed the red line (then β = 7). (C) Analysis of the mean connectivity for various soft-thresholding powers. (D) Number of genes in the modules with their frequency. (E) Dendrogram of the analyzed genes (60,379 genes) clustered based on a dissimilarity measure of topological overlap matrices (1-TOM). Blockwise dendrogram was obtained using average linkage hierarchical clustering with maximum block size of 46,000 genes. Modules were identified using height cut off equal to 0.15 with minimum module size of 30 genes. The different color modules correspond to the branches of the dendrogram.
Figure 2 Heatmap plotting of MEs values in relationship to tissue samples. n indicates number of samples.
Identification of Meiosis-Related Modules
A correlation analysis was conducted using the 66 MEs and the 8 different tissue types. A module was considered as meiosis-related when there was a strong correlation (r) with the 17 meiosis samples, and a weak or negative correlation with other tissue types. Accordingly, three meiosis-related modules were identified: module 2 (containing 4,940 genes), module 28 (544 genes), and module 41 (313 genes). Module 41 showed the strongest correlation with meiotic tissue (r = 0.73, FDR = 2.7 x 10−20), compared to module 2 (r = 0.61, FDR = 9.2 x 10−13) and module 28 (r = 0.52, FDR = 2.1 x 10−8) (Figure 3).
Figure 3 Co-expression network modules in relationship to tissues samples. Each row corresponds to a module; each column corresponds to a tissue type. Each cell contains the correlation value and, in parentheses, its corresponding FDR adjusted P value. n indicates number of samples. Only modules that have correlation value > 0.5 with meiotic anther tissue are shown.
Two other modules (modules 11 and 25) also showed significantly positive correlation with meiotic samples (r = 0.68 and 0.65, respectively); however, they were not considered meiosis-related because they also correlated with samples from floral organs (at stages other than meiosis) and spike tissues, as shown in Figure 3. Therefore, our analysis focused on the three modules (2, 28, and 41), exhibiting a strong correlation with meiotic tissues and not with other floral organs, while modules 11 and 25 were considered as non-meiosis specific modules (referred to in this paper as non-meiotic modules). Other tissue-related modules (the top three correlated modules) were also identified to be used as controls for the meiosis-related modules in the subsequent analysis. These modules were grain-related modules 5, 13, and 32 (r = 0.89, 0.89, and 0.85, respectively); leaves-related modules 1, 45, and 60 (r = 0.72, 0.68, and 0.71, respectively); and roots-related modules 7, 9, and 64 (r = 0.70, 0.76, and 0.85, respectively) (Figure S3).
Biological Significance of Expression Similarity in Modules
Several approaches were undertaken to validate the meiosis-related modules. The three modules (2, 28, and 41), strongly correlated with meiotic tissue expression, were found to be significantly enriched with the gene ontology (GO) slim terms “cell cycle,” “DNA metabolic process,” “nucleobase-containing compound metabolic process,” and “nucleus” (Figure 4). Among the top five enriched GO slim terms in each of the 66 modules, the term “cell cycle” was significant only in the three meiosis-related modules, suggesting this was not a general property of all modules and was instead specific to the meiosis modules. Module 2 in particular was significantly enriched with GO terms related to many biological processes occurring during meiosis such as “DNA replication,” “histone methylation,” “cytokinesis,” “nucleosome assembly,” and “chromatin silencing” (Table 1; Table S3). The term “double-strand break repair via homologous recombination,” an important process during meiosis, was the primary enriched Biological Processes GO term in module 41 (FDR < 0.05). The biological processes mediated by genes in module 28 included “protein deneddylation,” “positive regulation of G2/M transition of mitotic cell cycle,” “COP9 signalosome,” and other terms related to protein deneddylation and cell cycle control (Table 1). GO terms of meiosis-related modules were compared with those of modules highly correlated with other tissues. The GO terms “chloroplast,” “plastids,” “thylakoid,” and “photosystem” were significantly enriched in module 1, the most highly correlated module with leaves. The terms related to protein ubiquitination and protein binding were enriched in module 5 (the most highly correlated module with grain), while the terms “lignin biosynthetic process,” “phenylpropanoid metabolic process,” and “response to wounding” were enriched in module 7, the largest module correlated with roots (Figure S4). This indicated that our co-expression module-tissue correlation was meaningful both from the biological and physiological points of view. Detailed information of the enriched GO and GO slim terms in all modules is listed in the supplementary table Table S3. In summary, GO analysis confirmed that the three modules (2, 28, and 41) were enriched for genes associated with meiotic processes.
Figure 4 Enriched GO slim terms in the meiosis-related modules. Top five enriched GO terms in each module are shown. BP indicates biological processes and CC cellular component. No molecular function (MF) GO terms appear among the top five GO slim terms. Black bars indicate the number of genes in the GO term.
Enrichment of Meiosis-Related Modules for Wheat Orthologs of Known MGs
An assessment was undertaken to confirm that the meiosis-related modules were enriched for wheat orthologs of known MGs. Although the first wheat meiotic cDNA clones were isolated concurrently with the early discoveries of MGs in other plants (Ji and Langridge, 1994), the identification of MGs in wheat has been hampered by the large wheat genome size, its polyploid nature, and the absence of a complete genome sequence. Thus, in comparison to model plants (Arabidopsis and rice), few MGs have been functionally characterized in wheat. Characterized wheat MGs include TaASY1 (Boden et al., 2007; Boden et al., 2009), TaMSH7 (Lloyd et al., 2007), TaRAD51 (Khoo et al., 2008; Devisetty et al., 2010), TaDMC1 (Khoo et al., 2008; Devisetty et al., 2010), TaPSH1 (Khoo et al., 2012), TaZIP4 (Martín et al., 2017; Martín et al., 2018; Rey et al., 2018), and RecQ-7 (Gardiner et al., 2019). Given that, the assessment of whether the three modules contain known MGs was undertaken using orthology informed approaches. A set of 1,063 candidate MGs in wheat was identified and categorized based on the method used to identify the genes:
- The “orthologs” group contained 407 genes (Table S4; Sheet 1) that correspond to wheat orthologs of 103 functionally characterized MGs in model plant species. The majority of these genes (97.5%) were retrieved from the Ensembl Plant ortholog database, while for MGs with no wheat orthologs identified using this method, potential wheat orthologs were identified by searching for amino acid sequence similarity using BLASTP (see Materials and Methods). A list of 10 wheat gene IDs that are potentially orthologs of 4 MGs (AM1, ATM, FANCM, and ZYP1) were identified using the BLASTP methods. These four genes were included in our analysis due to their importance for meiosis in other plant models. However, using homology methods (like BLASTP) to infer orthology has a considerably high false positive error rate (Chen et al., 2007). Thus, the corresponding proteins of those 10 genes together with any conclusions drawn from their identification are tentative and must be treated with caution. There were no apparent wheat orthologs for 14 plant MGs (AtDFO, AtNACK2, AtPRD2, AtRECQ4B, AtSGO1, CDKG1, GIG1, JASON, MMD1, MS5, OsMEL1, OsMOF, PANS1, and XRI).
- The “meiotic GO” group contained 927 wheat genes annotated with one or more meiotic GO terms (Table S4; Sheet 2).
There were 271 genes overlapping between the two groups (Figure 5A), which were considered in the “orthologs” group when undertaking gene enrichment analysis. The presence of each gene in the different modules was determined. A set of 848 genes was assigned to modules in the co-expression network (Figure 5B), including 340 genes in meiosis-related modules. Genes from both groups were significantly over-represented (P < 0.05) in four modules, including two meiosis-related modules (2 and 28). Module 2, in particular, was the most enriched for these genes, possessing more than one third of the total candidate MGs assigned to modules. Module 2 had 142 wheat orthologs of MGs and 155 genes with meiotic GO terms, compared to the expected numbers (based on module size) of 27 and 42, respectively. Module 41 (the third meiosis-related module) was enriched only with genes from the “orthologs” group, having 15 orthologs of known MGs, whereas the expected number was 2 (Figure 5C). Consistent with this, genes from the “orthologs” and/or “meiotic GO” groups were significantly under-represented in modules strongly correlated with other types of tissue (modules 1, 7, and 9), and in modules with negative or no correlation with meiotic tissue (modules 3, 8, 14, and 17) (Table S4; Sheet 3).
Figure 5 Enrichment of meiosis-related genes in the co-expression network modules. (A) Venn diagram of total number of genes in the three groups: meiosis-related modules (genes in modules 2, 28, and 41), orthologs (wheat orthologs of MGs in model plant species), and meiotic GO (genes with meiotic GO terms). N indicates number of genes in each group. Numbers in brackets refer to number of genes not included in the WGCNA analysis because they are not expressed in meiotic anther tissue. (B) Total number of genes assigned to modules from orthologs of MGs (orthologs) and genes with meiotic GO terms (meiotic GO). (C) Gene enrichment in modules. Statistical significance of gene enrichment in a module is color coded (red indicates over-represented, blue under-represented, and gray not significant; P < 0.05). Rhombus shape indicates the expected number of genes in module.
In this study, three co-expression gene modules were identified that are strongly correlated to meiotic anther tissue and highly enriched with GO terms related to many processes occurring during meiosis, orthologs of known MGs, and genes having meiosis-specific GO terms. Although 67 (65%), out of the 103 wheat orthologs, had at least one gene copy assigned to one of the three meiosis-related modules, there were 36 orthologs whose gene homeologs were assigned to other modules (Table S5). Some of those genes have essential meiotic functions, like ASY1, which encodes a protein essential for homologous chromosome synapsis (Caryl et al., 2000; Armstrong et al., 2002; Boden et al., 2007); and DMC1, a gene encoding a recombination protein that acts only in meiosis (Klimyuk and Jones, 1997; Devisetty et al., 2010). Others are known to have both meiotic and mitotic functions, like BRCA2, a DNA repair gene required for double strand breaks repair by homologous recombination (Trapp et al., 2011) and SMC1 and SMC3, chromosome cohesion genes (Lam et al., 2005); thus, they are expressed in both reproductive and non-reproductive tissues. Assessment of these 36 orthologs showed that the expression patterns of their gene copies did not allow them to be clustered in any of the meiosis-related modules (or allocated to module 0 that is composed of genes not forming part of a co-expressed module), either because they were expressed in most samples from all types of tissues, or because they were expressed in a few samples of a specific tissue type (like meiotic anther tissue). The expression values (TPMs) of all gene copies of those 36 orthologs are summarized in Table S5. The number of meiotic anther samples (17) used in the present study might not be enough to identify all MGs being expressed in a specific meiotic stage. Such genes might be identified by WGCNA analysis when a larger number of meiotic samples become available. However, the analysis confirmed that the meiosis-related modules were indeed enriched for orthologs of known MGs, and for GO terms associated with processes involved in meiosis.
Copy Number of MGs
It has previously been suggested that MG duplicates return to a single copy following whole genome duplication more rapidly than the genome-wide average in angiosperms (Lloyd et al., 2014). The analysis of 19 meiotic recombination genes in hexaploid wheat and oilseed rape showed no evidence of gene loss after polyploidization. However, a recent study in tetraploid oilseed rape showed that reducing the copy number of MSH4, a key meiotic recombination gene involved in the ZMM (an acronym stands for the MGs Zip1/Zip2/Zip3/Zip4, Msh3/Msh5, and Mer3 identified initially in yeast) pathway, prevents meiotic crossovers between non-homologous chromosomes (Gonzalo et al., 2019). This led to the suggestion that meiotic adaptation in polyploids could involve “fine-tuning” the progression or the effectiveness of meiotic recombination, which could be achieved through the loss of MG duplicates in the newly formed polyploids (Lloyd et al., 2014; Gonzalo et al., 2019). This hypothesis was evaluated in hexaploid wheat. The gene copy number was assessed for the genes in the three meiosis-related modules and compared with genes in all modules and in other tissue-related modules. Analysis showed that the percentage of genes belonging to triads was 74.4% in the meiosis-related modules, which was similar to this percentage in other tissue-related modules (72.5, 74.0, and 76.1% in leaf-, grain-, and root-related modules, respectively); however, it was significantly higher than those of the non-meiotic modules (57.4%). The highest percentage of genes with three homeologs (83.3%) and lowest percentage of genes with single copy (2.7%) were observed in the group of genes identified as MG orthologs and/or possessing a meiotic GO term (Figure 6A).
Figure 6 Copy number and homeolog expression pattern for genes from meiosis-related and other tissue-related modules. (A) Proportion of genes in each copy number category (triads, duplets, monads, and others) for different sets of expressed genes during meiosis including: “meiotic modules” refers to the three meiosis-related modules 2, 28, and 41. “Non-meiotic modules” refers to the modules 11 and 25 that showed high correlation with meiotic anther but were not considered meiosis-related because they were also correlated with spike and floral organs tissues. The top three correlated modules with each of leaves (modules 1, 45, and 60), grain (modules 5, 13, and 32), and roots (modules 7, 9, and 64) tissues. “All modules” contains all genes assigned to modules in the co-expression network, and “orthologs and meiotic GO” refers to the set of genes that are orthologs of known MGs in other plant species and/or have meiotic GO terms. n number of genes in each set. (B) Proportion of genes from each homeolog expression pattern category (balanced, dominant, and suppressed) calculated for triads in the previously mentioned sets of genes. n number of genes in each set. (C) Ternary plot showing relative expression abundance in meiotic anther tissue of 2,366 triads to which the genes of meiosis-related modules (2, 28, and 41) belong. Each circle represents a gene triad with an A, B, and D coordinate consisting of the relative contribution of each homeolog to the overall triad expression. Triads in vertices correspond to single-subgenome-dominant categories, whereas triads close to edges and between vertices correspond to suppressed categories. Box plots indicate the relative contribution of each subgenome based on triad assignment to the seven categories (balanced, A dominant, B dominant, D dominant, A suppressed, B suppressed, D suppressed). Percentages between brackets indicate the percentage of triad number in each category to the total number of triads.
The high percentage of meiosis-related genes present as triads provides evidence that polyploid wheat did not experience significant gene loss (gene erosion) after polyploidization. However, this assumes that these genes were originally present as single copy genes in each of the A-, B-, and D-genome progenitor species which gave rise to polyploid wheat. Therefore, the copy number of the 103 wheat MG orthologs in wheat progenitor species was investigated. All possible orthologs (high and low confidence predicted orthologs) were retrieved from Ensembl Plants Genes 43 database for Triticum urartu (ASM34745v1; Ling et al., 2013), the diploid progenitor of the wheat A-genome; the D-genome ancestor Aegilops tauschii (Aet_v4.0; Luo et al., 2017), the diploid progenitor of the wheat D-genome; and Triticum dicoccoides (WEWSeq_v.1.0; Avni et al., 2017), the tetraploid progenitor of the hexaploid wheat (genome AABB). There was no change in copy number of 78.4% of genes, while 6.3 and 15.3% of genes had a lower and greater number of copies, respectively (Table 2 and Table S6). Regardless of genome of origin, the percentage of MGs with more copies was always greater than the percentage of genes with fewer copies. Comparing the A-genome MG copy number in hexaploid wheat with the relevant orthologs copy number in the corresponding A-genome ancestor, 86 genes (84.5%) had the same gene copy number in T. dicoccoides as in hexaploid wheat, while only 64 genes (63.1%) had the same gene copy number in T. urartu. This is consistent with the evolutionary history of hexaploid wheat, with T. dicoccoides being a more recent wheat progenitor (∼10,000 years) than T. urartu (> 5 million years) (Marcussen et al., 2014).
Table 2 Changes in copy number of wheat MGs in comparison with their orthologs in wheat progenitors.
An analysis on a subset of wheat genes, which were expected to be involved in meiotic recombination based on the function of their orthologs in model plants (64 genes; Table S6), was conducted. Again, results showed that the majority (94.9%) of those genes had greater or no change in number of copies (Table 2). Given it has been suggested that the reduction in the copy number of ZMM pathway genes could stabilize meiosis in Brassica (Gonzalo et al., 2019), the copy number of the wheat orthologs of seven ZMM genes was evaluated. Five of the seven ZMM genes (MER3, MSH5, ZIP4, PTD, and SHOC1) had equal or greater number of copies. However, TaMSH4 gained one A-genome copy (comparing with T. urartu) and lost one D-genome copy (compared with Ae. tauschii), while TaHEI10 lost A-genome copy and gained a D-genome copy (Table S6). In conclusion, our findings did not support any significant gene loss upon the polyploidization of hexaploid wheat, as suggested for other polyploids (Scannell et al, 2006; Lloyd et al., 2014; Gonzalo et al., 2019).
Homeolog Expression Patterns in Triads of MGs
Initial analysis revealed that most genes expressed during meiosis showed balanced expression between homeologs (Figure S1). The analysis was repeated using gene expression within the validated meiosis modules. Genes from all modules were assigned to three categories (balanced, dominant, and suppressed). Homeolog expression patterns in triads showed that meiosis-related modules 2, 28, and 41 had the highest percentage (87.3%) of genes with balanced expression (belong to balanced triads), compared to the top three tissue-related modules for grain, leaves, and roots (Figure 6B). Surprisingly, the group of candidate MGs selected for being orthologs of known MGs in other plant species and/or having meiotic GO terms had a higher percentage of genes from balanced triads (88.3%), whereas the modules not considered meiosis-specific (having high correlation with meiotic anther tissue and with spike and floral organs) contained only 68.3% genes with balanced expression (Figure 6B). The majority (84.19%) of triads with genes in meiosis-related modules (2,366 triads) showed balanced expression in meiotic anther tissue (Figure 6C).
In wheat, meiotic recombination and gene evolution rates are strongly affected by chromosome position, with relatively low recombination rates in the interstitial and proximal regions (genomic compartments R2a, C, and R2b) but notably higher rates toward the distal ends of the chromosomes (genomic compartments R1 and R3) (Akhunov et al., 2003; Choulet et al., 2014). The lack of significant changes in gene content and more balanced expression between homeologs suggested that these genes might be more prevalent in the proximal genomic compartments (Ramírez-González et al., 2018; International Wheat Genome Sequencing Consortium, 2018). The distribution of MGs was therefore assessed across the genomic compartments compared with the distribution of all HC genes across chromosomes. Analysis showed that genes from the meiotic modules (modules 2, 28, and 41) were significantly over-represented in the genomic compartments R2a, C, and R2b (P = 2.4 x 10−5, 3.1 x 10−6, and 1 x 10−5, respectively), while they were under-represented in the R1 and R3 genomic compartments (P = 1.7 x 10−8 and 1.7 x 10−10, respectively) (Figure 7A). Enrichment in the R2a genomic compartment region was not observed for genes from any of the other top three tissue-related modules (Figure 7B), since 21.7% of genes from the meiosis-related modules were assigned to R2a, while this percentage ranged between 18.2 and 19.5% in other tissue-related modules (Table 3). Interestingly, the set of the genes identified through orthology approaches and MG ontology approaches had also similar high percentage (21%) of genes assigned to R2a genomic compartment (Table 3).
Figure 7 Enrichment of genes from different tissue-related modules in the wheat genomic compartments. (A) Number of genes (actual and expected) from the three meiosis-related modules in each genomic compartment. (B) Comparison of number of genes from different tissue-related modules in genomic compartments. Statistical significance of gene enrichment in modules is color coded (red indicates enriched, blue depleted, and gray not significant; P < 0.05). Black dots indicate the expected number of genes in groups.
Our analysis reveals that homeologous MGs on homeologs mostly show balanced expression and lack a significant change in MG content following polyploidization. The majority of homeologous genes (not only MGs) on homeologs also show over 95% sequence identity to each other (Schreiber et al., 2012; International Wheat Genome Sequencing Consortium, 2018). Given these observations, such homeologs could synapse and recombine during meiosis. However, in allohexaploid wheat, homologs rather than homeolog synapse and recombine during meiosis ensuring the stability and fertility of this species, and the Ph1 locus, in particular to the TaZIP4 gene copy inside this locus, has been identified as the main locus controlling this process. The wheat ZIP4, an ortholog of ZIP4/Spo22 in A. thaliana and rice, is a member of the ZMM genes involved in the synaptonemal complex formation and class I crossover pathway (Chelysheva et al., 2007; Shen et al., 2012). Moreover, wheat lacking Ph1 exhibits extensive genome rearrangements, including translocation, duplications, and deletions (Martín et al., 2018). Thus, the evolution of Ph1 during wheat polyploidization is likely to explain why wheat has largely maintained a similar gene content and balanced expression of its homeologs. How meiosis has adapted to cope with allopolyploidy in other species is still to be resolved; however, it has been suggested that reduction in the copy number of MGs may stabilize the meiotic process after polyploidization (Lloyd et al., 2014; Gonzalo et al., 2019). The present study shows that this is not the case in wheat. It is likely that the presence of Ph1 in wheat enabled the retention of multiple copies of MGs as an alternative mechanism to ensure proper segregation of chromosomes during meiosis. The identification of the TaZIP4 gene within the Ph1 locus as the gene responsible for the Ph1 effect on recombination and the observed effects of Ph1 in wheat suggests that it may have more of a central role in meiosis than originally suspected from studies on model systems (Chelysheva et al., 2007; Shen et al., 2012). It has recently been suggested that ZIP4 might act as a scaffold protein facilitating physical interactions and assembly of different proteins complexes (De Muyt et al., 2018). Therefore, our co-expression network was used to identify the wheat orthologs of known MGs connected with TaZIP4. The analysis indicates that the three TaZIP4 homeologs on group 3 chromosomes (TraesCS3A02G401700, TraesCS3B02G434600, and TraesCS3D02G396500) were clustered in module 2, the largest meiosis-related module, and strongly connected to many orthologs of MGs with various meiotic functions (Figure 8). However, the TaZIP4 copy responsible for Ph1 phenotype (TraesCS5B02G255100) did not cluster in the same module, reflecting its different expression profile from the other homeologs, being expressed in most tissues (Martín et al., 2017; Martín et al., 2018; Rey et al., 2018). The analysis reveals that TaZIP4 is connected to several genes involved in centromere function. Studies on budding yeast have suggested that ZIP4 may affect centromere pairing during meiosis (Kurdzo et al., 2017). Moreover, the Ph1 locus has been shown to affect centromere pairing during meiosis in hexaploid wheat (Martinez-Perez et al., 2001). Therefore, it will be important to assess whether TaZIP4 within the Ph1 locus is responsible for this centromere effect. TaZIP4 was connected to wheat orthologs of genes known to be involved in crossover formation such as MSH2, SHOC1, FANCM, FLIP, EME1B, and MUS81 (Mercier et al., 2015). This suggests that there may be an interplay between TaZIP4 and genes from the anti-crossover pathway. This may be important as knockouts of genes involved in the anticrossover pathway have been shown to increase crossovers in other crops (Mieulet et al., 2018). However, on wheat’s polyploidization, TaZIP4 has been duplicated and diverged in order to improve homolog pairing and prevent homeolog crossover (Martín et al., 2018; Rey et al., 2018). This event may also affect the action of these anticrossover genes and the effects of their knockouts. Thus, TaZIP4 sub-network analysis supports a more central role of ZIP4 in meiosis than originally suspected from studies on model species.
Figure 8 The wheat MG orthologs connected to TaZIP4. The alluvial diagram shows the connected genes to the TaZIP4 homeologs TaZIP4-A1, TaZIP4-B1, and TaZIP4-D1. Edge weight > 1 was used as threshold to visualize connected genes. Black bars indicate the number of homeologs for each connected gene.
Further Characterization of the Wheat Meiotic Co-expression Network
Identification of Hub Genes in the Meiosis-Related Modules
Hub genes were identified within our meiosis-related modules by calculating the correlation between expression patterns of each gene and the ME: the most highly correlated genes to the eigengene being the hub genes. The top 10 hub genes of each module with their functional annotation are shown in Table 4. The top 10 hub genes in module 2 were core histone genes, supporting the strong contribution of histones in this meiosis-related module. For further verification of histone involvement in module 2 and other modules in general, all wheat genes annotated as core histones or having GO terms related to histone modification were retrieved for enrichment analysis. Analysis showed that the five types of histones (H1, H2A, H2B, H3, and H4) were enriched only in module 2 (P = 3.6 x 10−4, 1.2 x 10−22, 1.1 x 10−19, 9.4 x 10−21, and 3.3 x 10−26, respectively), having 433 genes (85% of all core histone genes in all modules), compared to an expected number of genes of 39 (Figure 9). Similar results were obtained for histone modification genes. Module 2 was the most enriched module with this group of genes (P = 9.3 x 10−52), containing 438 genes (30% of all histone modification genes in all modules). The histone modification genes were also enriched in 11 other modules, including the other meiosis-related modules (modules 28 and 41), however, with much lower numbers of enriched genes (Figure 9). Detailed information about genes included in this analysis is provided in Table S7. The strong enrichment of histone modification genes in module 2 (the largest meiosis-related module) supports the important role of histone modifications in meiosis (Maleki and Keeney, 2004; Oliver et al., 2013; Hu et al., 2015; Luense et al., 2016; Xu et al., 2016; Wang et al., 2017).
Figure 9 Histone genes enrichment in the gene co-expression network modules. The analysis included all the genes annotated as core histones (H1, H2A, H2B, H3, and H4) in the wheat genome and the genes with GO terms related to histone modification. Statistical significance of gene enrichment in a module at P < 0.05 is color coded (red indicates enriched, blue depleted, and gray not significant). Rhombus shape indicates the expected number of genes in module.
Hub genes such as “poor homologous synapsis 1” (PHS1) were also identified with module 41, the module most highly correlated to meiotic samples. This gene has been previously reported to have a key role in homologous chromosome pairing, synapsis, DNA recombination, and accurate chromosome segregation during meiosis in maize (Pawlowski et al., 2004), Arabidopsis (Ronceret et al., 2009), and wheat (Khoo et al., 2012). Other hub genes identified in modules 28 and 41 encoded for the high mobility group proteins (Tessari et al., 2003; Di Agostino et al., 2004; Pedersen et al., 2011; Antosch et al., 2015; Alonso-Martin et al., 2016), histone deacetylase (Akiyama et al., 2004; Wang et al., 2006; Magnaghi-Jaulin and Jaulin, 2006; Getun et al., 2017), and F-box proteins (Zheng et al., 2016).
Analysis of Transcription Factors Within the Meiosis-Related Modules
Many transcription factors (TFs) have been reported as key regulators of meiosis from studies on animals (Bolcun-Filas et al., 2011; Yan et al., 2016), yeast (Xu et al., 1995; Horie et al., 1998; Pierce et al., 2003; Beaudoin et al., 2018), and protozoa (Zhang et al., 2018). However, very little is known about the involvement of TFs in plant meiosis. The meiotic co-expression network was therefore exploited to identify potential meiosis-specific TFs. An assessment was undertaken of the enrichment of previously identified TF families in hexaploid wheat in the meiosis-related modules 2, 28, and 41. A total of 4,889 HC genes belonging to 58 TF families were predicted from the annotation of the wheat genome sequence. Of these, 2,439 TFs from 57 families could be assigned to the 66 modules in the gene co-expression network. Modules 2, 28, and 41 (meiosis-related modules) had 225, 25, and 17 TFs belonging to 31, 13, and 9 TF families, respectively (Table S8). Compared to the expected number of TF family genes in each module, only five TF families were significantly enriched in module 2: mitochondrial transcription termination factor (mTERF), growth-regulating factor (GRF), abscisic acid-insensitive protein 3/viviparous1 (ABI3/VP1), forkhead-associated domain (FHA), and E2F/dimerization partner (DP). On the other hand, four TF families were significantly depleted (Figure 10). The TF family NAC was the only TF family significantly enriched in module 41, containing 5 NAC genes (expected number 0.6; P < 0.05). Module 28 was not enriched with any TF family, although E2F/DP TFs were enriched in this module with borderline statistical significance (P = 0.06), with 4 genes in this module (the expected number was 0.2). Except in module 2, E2F/DP and FHA TF families were not enriched in any other modules in the gene co-expression network (Table S8). E2F/DP plays an important role in regulating gene expression necessary for passage through the cell cycle in mammals and plants (Zwicker et al., 1996; Zheng et al., 1999; Sozzani et al., 2006). Members of FHA contain the forkhead-associated domain, a phosphopeptide recognition domain found in many regulatory proteins. Genes belonging to the FHA group are reported to have roles in cell cycle regulation (Hollenhorst et al., 2000; Zhu et al., 2000; Kim et al., 2002), DNA repair (Sun et al., 1998; Bashkirov et al., 2003; Palmbos et al., 2008; Liang et al., 2015), and meiotic recombination and chromosome segregation (Pérez-Hidalgo et al., 2003; Cigliano et al., 2011; Crown et al., 2013). A previous meiotic transcriptome study identified up-regulation of TFs belonging to the MADS-box, bHLH, bZIP, and NAC families in Arabidopsis and maize meiocytes at early meiosis (Dukowic-Schulze et al., 2014). Zinc finger-like TFs have also been suggested to be regulators of maize MG expression (Ma et al., 2008). The present study indicates that TF families known to have roles in cell cycle, and meiosis processes are over-represented in the meiosis-related modules (module 2 in particularly). Those TF families contain about 20 meiosis-specific candidate TF genes whose function can be validated using the available reverse genetics resources in polyploid wheat (Krasileva et al., 2017).
Figure 10 Transcription factor families in the meiosis-related modules. Statistical significance of gene enrichment in the modules is color coded (red indicates over-represented, blue under-represented, and gray not significant; P < 0.05). Rhombus shape indicates the expected number of genes in module.
Visualization of Networks and Identification of Candidate MGs
Having identified meiosis-related modules, the networks within such modules can be visualized, highlighting genes for future studies. Edge files were created with gene annotation for the three meiosis-related modules 2, 28, and 41. Those files can be used to investigate the relation between orthologs of MGs within a module and ranked based on the strength of the connection (weight value). Another application of co-expression networks is the identification of previously uncharacterized genes regulating biological processes (Usadel et al., 2009; Aya et al., 2011; Costa et al., 2015; Silva et al., 2016; Krishnan et al., 2017; Tan et al., 2017; International Wheat Genome Sequency Consortium, 2018; Ma et al., 2018; Yu et al., 2018; Liu et al., 2019). Cytoscape 3.7.1 software (Shannon et al., 2003) was therefore used to visualize our network and to show connections between different orthologs of MGs in meiosis-related modules. Wheat MG orthologs in meiosis-related modules were used as “guide genes” to construct co-expression subnetworks containing only genes with direct connections to the guide genes. One such subnetwork is shown in Figure 11, where the following wheat orthologs of MGs in module 41 were selected and used to construct a meiotic subnetwork: poor homologous synapsis 1 (TaPHS1; Khoo et al., 2012), argonaute (AGO9/AGO104; Durán-Figueroa and Vielle-Calzada, 2010; Singh et al., 2011), replication protein A2c (OsRPA2c; Li et al., 2013), meiotic nuclear division protein 1 (AtMND1; Domenichini et al., 2006; Kerzendorfer et al., 2006), MMS and UV sensitive 81 (AtMUS81; Hartung et al., 2006; Higgins et al., 2008), and parting dancers (AtPTD; Wijeratne et al., 2006; Lu et al., 2014) (guide genes; red circles in Figure 11). The network complexity was reduced using an edge weight > 0.05. The visualized subnetwork contained 53 gene IDs including 9 guide gene copies. The gene TraesCS7A02G233900 (TaPHS1), a hub gene in module 41, was central in the network having the highest number of direct edges (41 direct edges; connected with 77.4% of the genes in the subnetwork). This subnetwork allowed identification of other genes with putative roles in meiosis (pink circles): (a) RNA recognition motif-containing gene (TraesCS5A02G319000) similar to Mei2, a master regulator of meiosis and required for premeiotic DNA synthesis as well as entry into meiosis in Schizosaccharomyces pombe (Watanabe and Yamamoto, 1994; Watanabe et al., 1997); (b) the gene TraesCS4D02G050000 showed similarity to Male meiocyte death 1 (MMD/DUET), a PHD-finger protein plays role in chromatin structure and male meiotic progression in A. thaliana (Reddy et al., 2003); and (c) the gene TraesCS5D02G454900, a possible TF belonging to the FHA family known to have function in cell cycle regulation (Hollenhorst et al., 2000; Zhu et al., 2000; Kim et al., 2002), DNA repair (Sun et al., 1998; Bashkirov et al., 2003; Palmbos et al., 2008; Liang et al., 2015), and meiotic recombination and chromosome segregation (Pérez-Hidalgo et al., 2003; Cigliano et al., 2011; Crown et al., 2013). The meiotic subnetwork contained genes with similarity to cell cycle like F-box family proteins, high mobility family proteins, and chromatin remodeling genes. The subnetwork also contained a group of genes connected to most of our guide genes, which thus might be involved with them in similar biological processes. Examples of such genes are TraesCS3A02G101000, TraesCS1A02G292700, and TraesCS1D02G291100 which encode for zinc finger CCCH domain-containing proteins (Figure 11). Other meiotic subnetworks were also constructed using other guide genes from modules 2 and 28.
Figure 11 A meiotic co-expression subnetwork in hexaploid wheat. This subnetwork was constructed using 9 guide genes in module 41. Guide genes are wheat orthologs of MGs in other plant species (red circles); pink circles represent genes with putative meiosis function. Edge weight 0.05 was used as threshold to visualize genes in the subnetwork using Cytoscape 3.7.1 software.
The Meiotic Co-Expression Network Is Accessible in a Larger Biological Contest
Our WGCNA co-expression network and GO enrichment data have been integrated with the wheat knowledge network (Hassani-Pak et al., 2016) to make it publicly accessible and searchable through the KnetMiner web application (http://knetminer.rothamsted.ac.uk; Hassani-Pak, 2017). KnetMiner can be searched with keywords (incl. module ID and GO terms) and wheat gene identifiers. The gene knowledge graphs generated contain many additional relation types such protein–protein interactions, homology, and links to genome wide association studies and associated literature placing the co-expression networks generated here in a wider context.
In summary, the present study shows that most MGs in wheat are retained as three homeologous genes, which are expressed during meiosis at similar levels, suggesting that they have not undergone extensive gene loss nor sub/neo-functionalization. Meiosis-related modules have been used to create networks and identify hub genes providing targets for future studies. The network containing the ZIP4 gene, recently defined as Ph1 (Martín et al., 2017; Martín et al., 2018; Rey et al., 2018)—for example, highlights potential interacting partners. Finally, the networks highlight genes such as ZIP4 and “poor homologous synapsis 1,” which may play a more central role in meiosis than previously thought. The co-expression network analysis combined with orthologue information will contribute to the discovery of new MGs and greatly empowers reverse genetics approaches to validate the function of candidate genes (Krasileva et al., 2017). Ultimately, this will lead to better understanding of the regulation of meiosis in wheat (and other polyploid plants) and subsequently improve wheat fertility.
Materials and Methods
RNA-Seq Data Collection
For co-expression network analysis, we included 130 samples, containing 113 samples previously described in Ramírez-González et al. (2018) and 17 samples from anthers during meiosis (9 samples from Martín et al. (2018), and 8 samples downloaded from https://urgi.versailles.inra.fr/files/RNASeqWheat/Meiosis/). Samples were selected to represent all main tissue types: grain (n = 37 samples), leaves (n = 21 samples), roots (n = 20 samples), anther at meiosis (n = 17 samples), spike (n = 12 samples), floral organs (anther, pistil, and microspores) at stages other than meiosis (n = 10 samples), stem (n = 7 samples), and shoots (n = 6 samples). All samples were under nonstress conditions and mostly from the reference accession Chinese Spring. Detailed information about the used samples are listed in the Supplementary Materials (Table S2).
Mapping of RNA-Seq Reads to Reference
Kallisto v0.42.3 (Bray et al., 2016) was used to map all RNA-Seq samples to the Chinese Spring transcriptome reference IWGSC RefSeq Annotation v1.1 (International Wheat Genome Sequencing Consortium, 2018), following default parameters previously shown to result in accurate homeolog-specific read mapping in polyploid wheat (Borrill et al., 2016; Ramírez-González et al., 2018). Tximport v1.2.0 was then used to summarize expression levels from transcript to gene level (Text S1; Part 1).
Co-Expression Network Construction
The WGCNA package in R (Langfelder and Horvath, 2008; Langfelder and Horvath, 2012) was used to construct the scale-free co-expression network. Metadata for all samples were assigned with eight tissue types (average 16.25; median 14.5 replicates per factor). Only HC genes (International Wheat Genome Sequencing Consortium, 2018) with expression > 0.5 TPM in at least one meiosis sample were retained for co-expression network construction using the R Package WGCNA (version 1.66). Using the varianceStabilizingTransformation() function from DESeq2 (Love et al., 2014), the count expression level of selected genes was normalized to eliminate differences in sequencing depth between different RNA-Seq studies (Text S1; Part 2). To select a soft power threshold (β) for adjacency calculation (as aij = |sij|β; where sij is the correlation between gene i and gene j), the scale-free topology criterion was used (Zhang and Horvath, 2005). The soft thresholding means suppressing low correlations in a continuous (“soft”) manner by using β value to power the correlation of the genes to that threshold, which reduces the noise of the correlations in the adjacency matrix. Using the pickSoftThreshold() function to calculate β values, the soft power threshold emphasizing strong correlations between genes and penalizing weak correlations was selected as the first power to exceed a scale-free topology fit index of 0.9 (Ramírez-González et al., 2018) (Text S1; Part 3). The correlation type used to calculate adjacency matrices was biweight midcorrelation (bicor). The adjacency matrices were transformed into a topological overlap matrix (TOM), measuring the network connectivity of a gene deﬁned as the sum of its adjacency with all other genes for network generation. The blockwiseModules() function was used to calculate matrices and construct blockwise networks considering the following parameters: network type (NetworkType) = “signed hybrid,” maximum block size (maxBlockSize) = 46,000 genes, soft power threshold (power) = 7, correlation type (corType) = “bicor” (biweight midcorrelation with maxPOutliers set to 0.05 to eliminate effects of outlier samples), topological overlap matrices type (TOMType) = “unsigned” with the mergeCutHeight = 0.15, and the minModuleSize = 30 to classify genes with similar expression proﬁles into gene modules using average linkage hierarchical clustering, according to the TOM-based dissimilarity measure with a minimum module size of 30 genes (Text S1; Part 4). MEs, summarizing the expression patterns of all genes within a given module into a single characteristic expression profile, were calculated as the first principal component in the principal component analysis (PCA) using the moduleEigengenes() function (Text S1; Part 4).
Identifying Meiosis-Related Modules
The MEs were used to test correlations between gene modules and traits (eight tissue types) using the cor() function. To assess the significance of correlations, Student asymptotic P values for correlations were calculated using the function corPvalueStudent() and corrected for multiple testing by calculating FDR (false discovery rate) using a p.adjust() function following the Benjamini and Yekutieli (2001) method. We considered a module meiosis-related when its correlation was strong with meiosis samples (r > 0.5 and FDR < 0.05) and weak (r < 0.3) or negative with other type of tissues (Text S1; Part 5).
Analysis of GO Term Enrichment in Modules
GO term enrichment was calculated using the “goseq” package (Young et al., 2010). Gene ontology (GO) annotations of IWGSC RefSeq v1.0 genes were retrieved from the file “FunctionalAnnotation.rds” in https://opendata.earlham.ac.uk/wheat/under_license/toronto/Ramirez-Gonzalez_etal_2018-06025-Transcriptome-Landscape/data/TablesForExploration/FunctionalAnnotation.rds (Ramírez-González et al., 2018) by filtering for ontology “IWGSC+Stress.” GO data was then converted to IWGSC RefSeq Annotation v1.1 by replacing “01G” by “02G” in the IWGSC v1.0 gene IDs and retaining only genes > 99% similar with > 90% coverage in the v1.0 and v1.1 annotation versions (as determined by BLASTn of the cDNAs) (called “all_go”). P values for GO term enrichment were calculated using the goseq() function (using the following parameters: the pwf object was created using the nullp() function which calculated a probability weighting function for the genes v1.1 based on their length, the gene2cat = all_go, and test.cats = “GO : BP,” to specify the biological process GO term category to test for over representation amongst the inquired genes) and corrected using the FDR method (Benjamini and Hochberg, 1995). A GO term was considered enriched in a module when FDR adjusted P value < 0.05 (Text S1; Part 6). All figures shown for enriched GO terms in the modules were produced using RAWGraphs software (Mauri et al., 2017).
Orthologs of MGs in Wheat
A comprehensive literature search was performed for MGs in model plant species (mainly A. thaliana and rice; Table S9), identifying gene IDs based on the “Os-Nipponbare-Reference-IRGSP-1.0” for rice (Oryza sativa Japonica Group) and “TAIR10” for A. thaliana. Wheat orthologs of MGs were then retrieved from Ensembl Plants Genes 43 database through BioMart (in; http://plants.ensembl.org/biomart) where orthologs calculated according to Vilella et al. (2009) using the following gene datasets: “Triticum aestivum genes (IWGSC),” “Oryza sativa Japonica Group genes (IRGSP-1.0),” and “A. thaliana genes (TAIR10)” for wheat (International Wheat Genome Sequencing Consortium, 2018), rice (Kawahara et al., 2013; Sakai et al., 2013), and A. thaliana (Waese et al., 2017), respectively. For MGs with no wheat orthologs identified using this method, potential wheat orthologs were identified by searching for amino acid sequence similarity using BLASTP (Korf et al., 2003) in Ensembl Plants according to the following criteria: e-value < 1e-10; ID% > 25% with Arabidopsis and > 70% with rice. By blasting the amino acid sequences of rice and Arabidopsis MGs against wheat proteins, a list of genes (that do not have rice and Arabidopsis MG orthologs in the Ensembl Plants database) was identified. For this list of wheat genes, we checked whether they have any other rice or Arabidopsis orthologs. Then, only wheat genes that did not have any rice or Arabidopsis orthologs were considered as orthologs of MGs. Finally, 407 wheat gene IDs were identified as orthologs of 103 plant MGs (listed in Table S4; Sheet 1). This group of genes was referred to in this study as “orthologs.”
Wheat Genes With MG Ontology (GO)
A total number of 46,909 GO terms used by Ramírez-González et al. (2018) to calculate GO term accessions for wheat genes (IWGSC v1.0 gene annotation) were filtered for meiosis-related GO terms, using 15 meiosis-specific keywords (“meiosis,” “meiotic,” “synapsis,” “synaptonemal,” “prophase I,” “metaphase I,” “anaphase I,” “telophase I,” “leptotene,” “zygotene,” “pachytene,” “diplotene,” “chiasma,” “crossover,” and “homologous chromosome segregation”). A total of 284 meiosis GO accessions were identified and used to retrieve 927 wheat genes with potential roles during meiosis (Table S4; Sheet 2). All genes identified by gene orthologs and gene ontology methods were then filtered to retain only genes had expression > 0.5 TPM in at least one meiosis sample. This group of genes was referred to in this paper as “meiotic GO.” Enrichment analysis for the genes from “orthologs” and “meiotic GO” groups in all module was conducted (Text S1; Part 7). The number of genes from each group was assessed in all modules and compared with the expected number based on the module size. There was a set of genes overlapping between “orthologs” and “meiotic GO” groups, which was considered in the “orthologs” group when undertaking gene enrichment analysis. Fisher’s exact test was used to calculate significant enrichment in the modules. Gene group considered over- or under-represented in a module when P < 0.05.
Identifying Highly Connected Hub Genes
Hub genes within each module were identified using the WGCNA R package function signedKME() to calculate the correlation between expression patterns of each gene and the ME. Hub genes were considered those more highly correlated to the eigengene (Text S1; Part 8).
Assessment of TF Families in Modules
A total of 4,889 wheat HC genes (IWGSC RefSeq Annotation v1.1; International Wheat Genome Sequencing Consortium, 2018) belonging to 58 TF families were predicted from the annotation of the wheat genome sequence (https://github.com/Borrill-Lab/WheatFlagLeafSenescence/blob/master/data/TFs_v1.1.csv). The number of TFs from each family was assessed in all modules and compared with the expected number based on the module size. Fisher’s exact test was used to calculate significant enrichment of TFs in the modules. TF family considered over- or under-represented in a module when P < 0.05 (Text S1; Part 9).
Defining Gene Categories Based on Number of Homeologs
A list of homeologs for all HC hexaploid wheat genes (IWGSC v1.1 gene annotation; International Wheat Genome Sequencing Consortium, 2018) was retrieved from Ensembl Plants Genes 43 database through BioMart (in; http://plants.ensembl.org/biomart). Based on number of homeologs from each of the A-, B-, and D-sub-genomes, genes were assigned to four groups: triads that refer to 1:1:1 triads (with a single copy from each of the A-, B-, and D-sub-genomes); duplets referring to 1:1:0, 1:0:1, and 0:1:1 duplets; monads group containing genes with no homeologs (e.g. 0:0:1); and “others” containing genes with more than two homeologs, in conjunction with genes from the homeologous groups 0:1:2, 0:2:1, 1:0:2, 2:0:1, 1:2:0, 2:1:0, 2:0:0, 0:2:0, and 0:0:2. Accordingly, 19,801 triads (59,403 genes), 7,565 duplets (15,130 genes), 15,109 monads (single-copy genes), and 18,250 genes from the “others” group were identified (Table S1).
Defining Gene Categories Based on Homeolog Expression Patterns in Triads
Homeolog expression pattern in triads was determined for each of the eight tissue types (Text S1; Part 10). For triads, it was calculated according to Ramírez-González et al. (2018) where a triad can be described as balanced, A dominant, A suppressed, B dominant, B suppressed, D dominant, or D suppressed, based on the relative expression contribution of its A, B, and D homeologs. Briefly, the relative contribution of each gene in a triad was calculated with the following formula:
Relative expression contribution (A) = TPM(A)/[TPM(A) + TPM(B) + TPM(D)]
where A, B, and D represent the gene corresponding to the A, B, and D homeologs in the triad. Each category is defined by the following ideal relative expression contributions:
The triad was assigned to the category with the shortest Euclidean distance to its relative contribution. Triads were defined as expressed when one of its homeologs was expressed according to the criterion used in our WGCNA analysis (Table S10; Sheet 1). This insured that all triads contain genes from modules were included in the homeolog expression bias analysis (Text S1; Part 11). Genes from a triad might not belong to the same module due to dissimilarity of their expression patterns. Thus, to allow the assessment of the expression pattern of genes in each module, each homeolog (A, B, and D homeologs) in a triad was assigned to one of the three categories “balanced,” “dominant,” and “suppressed” based on the homeolog origin (A, B, and D sub-genome) and the triad description (balanced, A dominant, A suppressed, B dominant, B suppressed, D dominant, or D suppressed) as shown in Table S10 (Sheet 2). The values of the relative contributions of each homeolog per triad were used to plot the ternary diagrams using the R package ggtern (Hamilton, 2016).
Co-Expression Gene Network Visualization
Cytoscape software (version 3.7.1; Shannon et al., 2003) was used to visualize the network described in this study. Firstly, the “exportNetworkToCytoscape” function was used to create edge files which could be used to visualize the network, then depending on network complexity, different weight value thresholds were used to filter genes to be visualized (Text S1; Part 12). The term “weight value” in the input files for Cytoscape refers to the connection strength between two nodes (genes) in terms of correlation value obtained from the topological overlap matrices (TOM). The co-expression network data has also been integrated with the wheat knowledge network (Hassani-Pak et al., 2016) to make it publicly accessible through the KnetMiner web application (http://knetminer.rothamsted.ac.uk; (Hassani-Pak, 2017). The data was semantically modeled as nodes of type gene, co-expression-module, co-expression-study, and GOterm, connected by relations of type part-of and enriched. Each module was given a unique identifier composed of the module number and the prefix “AKA.” KnetMiner can be searched with keywords (incl. module ID and GO terms) and wheat gene identifiers.
Data Availability Statement
All datasets used in this study can be found in in the Earlham Institute repository (https://opendata.earlham.ac.uk/wheat/under_license/toronto/Martin_etal_2018_Alabdullah_etal_2019_wheat_meiosis_transcriptome_and_co-expression_network/). All R scripts are provided as.text files in the Supplementary Material. The gene network data is accessible and searchable through the public KnetMiner website (http://knetminer.com).
AA, PB, AM, CU, and GM contributed conception and design of the study. AA, AM, RR-G, and KH-P organized and curated the data; AA performed the formal analysis; AA, PB, and RR-G wrote the R scripts; GM and PS secured the funding and supervised the work. AA wrote the original draft. All authors contributed to manuscript revision, read and approved the submitted version.
This work was supported by the UK biotechnology and Biological Sciences Research Council (BBSRC) through a grant as part of Designing Future Wheat (DFW) Institute Strategic Programme (BB/P016855/1) and response mode grant (BB/R007233/1).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This manuscript has been released as a pre-print at https://www.biorxiv.org/ (Alabdullah et al., 2019). The authors would like to acknowledge the UK biotechnology and Biological Sciences Research Council (BBSRC; www.bbsrc.ukri.org) for supporting this work.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.01325/full#supplementary-material
Figure S1 | Homeolog expression patterns of expressed triads in hexaploid wheat. Homeolog expression pattern was calculated for 19,801 triads (59,403 genes) across 8 tissue types according to published criteria (Ramírez-González et al., 2018), where triad defined as expressed when the sum of the A, B, and D subgenome homeologs was > 0.5 TPM. (A) Proportion of triads in each homeolog expression pattern across the 8 tissues. n is number of expressed triads. (B) Ternary plot showing relative expression abundance of 14,837 expressed triads (44,511 genes) in the meiotic anther tissue. Each circle represents a gene triad with an A, B, and D coordinate consisting of the relative contribution of each homeolog to the overall triad expression. Triads in vertices correspond to single-subgenome–dominant categories, whereas triads close to edges and between vertices correspond to suppressed categories. Box plots indicate the relative contribution of each subgenome based on triad assignment to the seven categories. Percentages between brackets indicate the percentage of triad number in each category to the total number of triads.
Figure S2 | Proportion of genes in each homeologs number category. Expressed genes across 8 tissues were assigned to four categories (triads, duplets, monads and others). n indicates number of expressed genes.
Figure S3 | Module-tissue relationship. Each row corresponds to a module; each column corresponds to a tissue type; Each cell contains the correlation value (r) and, in brackets, its corresponding FDR adjusted P value. n indicates number of samples. Only modules that have correlation value > 0.5 are shown.
Figure S4 | Enriched GO terms in the meiosis-related and other tissue-related modules. Top 5 GO terms are shown for each module. Black bars indicate the number of genes in the GO term.
Akhunov, E. D., Akhunova, A. R., Linkiewicz, A. M., Dubcovsky, J., Hummel, D., Lazo, G., et al. (2003). Synteny perturbations between wheat homoeologous chromosomes caused by locus duplications and deletions correlate with recombination rates. Proc. Natl. Acad. Sci. U.S.A. 100, 10836–10841. doi: 10.1073/pnas.1934431100
Alabdullah, A. K., Borrill, P., Martin, A. C., Ramirez-Gonzalez, R. H., Hassani-Pak, K., Uauy, C., Shaw, P., Moore, G. (2019). A co-expression network in hexaploid wheat reveals mostly balanced expression and lack of significant gene loss of homeologous meiotic genes upon polyploidization. BioRxiv. doi: 10.1101/695759
Alix, K., Gérard, P. R., Schwarzacher, T., Heslop-Harrison, J. S. (2017). Polyploidy and interspeciﬁc hybridization: partners for adaptation, speciation and evolution in plants. Ann Bot. 120, 183–194. doi: 10.1093/aob/mcx079
Alonso-Martin, S., Rochat, A., Mademtzoglou, D., Morais, J., de Reynies, A., Aurade, F., et al. (2016). Gene expression profiling of muscle stem cells identifies novel regulators of postnatal myogenesis. Front. Cell Dev. Biol. 4, 58–58. doi: 10.3389/fcell.2016.00058
Antosch, M., Schubert, V., Holzinger, P., Houben, A., Grasser, K. D. (2015). Mitotic lifecycle of chromosomal 3xHMG-box proteins and the role of their N-terminal domain in the association with rDNA loci and proteolysis. New Phytol. 208, 1067–1077. doi: 10.1111/nph.13575
Armstrong, S. J., Caryl, A. P., Jones, G. H., Franklin, F. C. H. (2002). Asy1, a protein required for meiotic chromosome synapsis, localizes to axisassociated chromatin in Arabidopsis and Brassica. J. Cell Sci. 115, 3645–3655. doi: 10.1242/jcs.00048
Avni, R., Nave, M., Barad, O., Baruch, K., Twardziok, S. O., Gundlach, H., et al. (2017). Wild emmer genome architecture and diversity elucidate wheat evolution and domestication. Science 357 (6346), 93–97. doi: 10.1126/science.aan0032
Aya, K., Suzuki, G., Suwabe, K., Hobo, T., Takahashi, H., Shiono, K., et al. (2011). Comprehensive network analysis of anther-expressed genes in rice by the combination of 33 laser microdissection and 143 spatiotemporal microarrays. PLoS ONE 6, e26162. doi: 10.1371/journal.pone.0026162
Bashkirov, V. I., Bashkirova, E. V., Haghnazari, E., Heyer, W. D. (2003). Direct kinase-to-kinase signaling mediated by the FHA phosphoprotein recognition domain of the Dun1 DNA damage checkpoint kinase. Mol. Cell Biol. 23, 1441–1452. doi: 10.1128/MCB.23.4.1441-1452.2003
Beaudoin, J., Ioannoni, R., Normant, V., Labbé, S. (2018). A role for the transcription factor Mca1 in activating the meiosis-specific copper transporter Mfc1. PLoS ONE 13, e0201861. doi: 10.1371/journal.pone.0201861
Benjamini, Y., Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B (Methodol.) 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x
Boden, S. A., Langridge, P., Spangenberg, G., Able, J. A. (2009). TaASY1 promotes homologous chromosome interactions and is affected by deletion of Ph1. Plant J. 57, 487–497. doi: 10.1111/j.1365-313X.2008.03701.x
Boden, S. A., Shadiac, N., Tucker, E. J., Langridge, P., Able, J. A. (2007). Expression and functional analysis of TaASY1 during meiosis of bread wheat (Triticum aestivum). BMC Mol. Biol. 8, 65. doi: 10.1186/1471-2199-8-65
Bolcun-Filas, E., Bannister, L. A., Barash, A., Schimenti, K. J., Suzanne, A., Hartford, S. A., et al. (2011). A-MYB (MYBL1) transcription factor is a master regulator of male meiosis. Development 138, 3319–3330. doi: 10.1242/dev.067645
Caryl, A. P., Armstrong, S. J., Jones, G. H., Franklin, F. C. (2000). A homologue of the yeast HOP1 gene is inactivated in the Arabidopsis meiotic mutant asy1. Chromosoma. 109, 62–71. doi: 10.1007/s004120050413
Chelysheva, L., Gendrot, G., Vezon, D., Doutriaux, M. P., Mercier, R., Grelon, M. (2007). Zip4/Spo22 is required for class I CO formation but not for synapsis completion in Arabidopsis thaliana. PLoS Genet. 3 (5), e83. doi: 10.1371/journal.pgen.0030083
Choulet, F., Alberti, A., Theil, S., Glover, N., Barbe, V., Daron, J., et al. (2014). Structural and functional partitioning of bread wheat chromosome 3B. Science 345, 1249721. doi: 10.1126/science.1249721
Cigliano, R. A., Sanseverino, W., Cremona, G., Consiglio, F. M., Conicella, C. (2011). Evolution of parallel spindles like genes in plants and highlight of unique domain architecture. BMC Evol. Biol. 11, 78. doi: 10.1186/1471-2148-11-78
Colas, I., Shaw, P., Prieto, P., Wanous, M., Spielmeyer, W., Mago, R., et al. (2008). Effective chromosome pairing requires chromatin remodeling at the onset of meiosis. Proc. Natl. Acad. Sci. U.S.A. 105, 6075–6080. doi: 10.1073/pnas.0801521105
Costa, M. C., Righetti, K., Nijveen, H., Yazdanpanah, F., Ligterink, W., Buitink, J., et al. (2015). A gene co-expression network predicts functional genes controlling the re-establishment of desiccation tolerance in germinated Arabidopsis thaliana seeds. Planta 242, 435–449. doi: 10.1007/s00425-015-2283-7
Crown, K. N., Savytskyy, O. P., Malik, S., Logsdon, J., Williams, R. S., Tainer, J. A., et al. (2013). A mutation in the FHA domain of Coprinus cinereus Nbs1 Leads to Spo11-independent meiotic recombination and chromosome Segregation. G3 (Bethesda) 3, 1927–1943. doi: 10.1534/g3.113.007906
De Muyt, A., Pyatnitskaya, A., Andréani, J., Ranjha, L., Ramus, C., Laureau, R., et al. (2018). A meiotic XPF-ERCC1-like complex recognizes joint molecule recombination intermediates to promote crossover formation. Genes Dev. 32, 283–296. doi: 10.1101/gad.308510.117
Devisetty, U. K., Mayes, K., Mayes, S. (2010). The RAD51 and DMC1 homoeologous genes of bread wheat: cloning, molecular characterization and expression analysis. BMC Res Notes 3, 245. doi: 10.1186/1756-0500-3-245
Di Agostino, S., Fedele, M., Chieffi, P., Fusco, A., Rossi, P., Geremia, R., et al. (2004). Phosphorylation of high-mobility group protein A2 by Nek2 kinase during the first meiotic division in mouse spermatocytes. Mol. Biol. Cell. 15, 1224–1232. doi: 10.1091/mbc.e03-09-0638
Domenichini, S., Raynaud, C., Ni, D. A., Henry, Y., Bergounioux, C. (2006). Atmnd1-delta1 is sensitive to gamma-irradiation and defective in meiotic DNA repair. DNA Repair (Amst) 5, 455–464. doi: 10.1016/j.dnarep.2005.12.007
Dukowic-Schulze, S., Harris, A., Li, J., Sundararajan, A., Mudge, J., Retzel, E. F., et al. (2014). Comparative transcriptomics of early meiosis in Arabidopsis and maize. J. Genet. Genomics 41, 139–152. doi: 10.1016/j.jgg.2013.11.007
Durán-Figueroa, N., Vielle-Calzada, J. P. (2010). Argonaute9-dependent silencing of transposable elements in pericentromeric regions of Arabidopsis.Plant Signal. Behav. 5 (11), 1476–1479. doi: 10.4161/psb.5.11.13548
Edger, P. P., Smith, R., McKain, M. R., Cooley, A. M., Vallejo-Marin, M., Yuan, Y., et al. (2017). Subgenome dominance in an interspecific hybrid, synthetic allopolyploid, and a 140-year-old naturally established neo-allopolyploid monkeyflower. Plant Cell 29, 2150–2167. doi: 10.1105/tpc.17.00010
Fernandes, J. B., Séguéla-Arnauda, M., Larchevêquea, C., Lloyda, A. H., Mercier, R. (2018). Unleashing meiotic crossovers in hybrid plants. Proc. Natl. Acad. Sci. U.S.A. 115 (10), 2431–2436. doi: 10.1073/pnas.1713078114
Freeling, M. (2009). Bias in plant gene content following different sorts of duplication: tandem, whole-genome, segmental, or by transposition. Annu. Rev. Plant Biol. 60, 433–453. doi: 10.1146/annurev.arplant.043008.092122
Gardiner, L.-J., Wingen, L. U., Bailey, P., Joynson, R., Brabbs, T., Wright, J., et al. (2019). Analysis of the recombination landscape of hexaploid bread wheat reveals genes controlling recombination and gene conversion frequency. Genome Biol. 20, 6. doi: 10.1186/s13059-019-1675-6
Getun, I. V., Wu, Z., Fallahi, M., Ouizem, S., Liu, Q., Li, W., et al. (2017). Functional roles of acetylated histone marks at mouse meiotic recombination hot spots. Mol. Cell Biol. 37, e00942–e00915. doi: 10.1128/MCB.00942-15
Gonzalo, A., Lucas, M.-O., Charpentier, C., Sandmann, G., Lloyd, A., Jenczewski, E. (2019). Reducing MSH4 copy number prevents meiotic crossovers between non-homologous chromosomes in Brassica napus. Nat. Commun. 10, 2354. doi: 10.1038/s41467-019-10010-9
Hamilton, N. (2016). ggtern: An extension to ‘ggplot2’, for the creation of ternary diagrams. R package version 2.2.1 [software]. Available from: https://CRAN.R-project.org/package=ggtern.
Hartung, F., Suer, S., Bergmann, T., Puchta, H. (2006). The role of atmus81 in dna repair and its genetic interaction with the helicase atrecq4a. Nucleic Acids Res. 34, 4438–4448. doi: 10.1093/nar/gkl576
Hassani-Pak, K., Castellote, M., Esch, M., Hindle, M., Lysenko, A., Taubert, J., et al. (2016). Developing integrated crop knowledge networks to advance candidate gene discovery. Appl Transl Genom. 11, 18–26. doi: 10.1016/j.atg.2016.10.003
Higgins, J. D., Buckling, E. F., Franklin, F. C. H., Jones, G. H. (2008). Expression and functional analysis of atmus81 in Arabidopsis meiosis reveals a role in the second pathway of crossing-over. Plant J. 54, 152–162. doi: 10.1111/j.1365-313X.2008.03403.x
Hollenhorst, P. C., Bose, M. E., Mielke, M. R., Müller, U., Fox, C. A. (2000). Forkhead genes in transcriptional silencing, cell morphology and the cell cycle: overlapping and distinct functions for FKH1 and FKH2 in Saccharomyces cerevisiae. Genetics. 154, 1533–1548.
Holm, P. B. (1986). Chromosome pairing and chiasma formation in allohexaploid wheat, Triticum aestivum analyzed by spreading of meiotic nuclei. Carlsberg Res.Commun. 51, 239–294. doi: 10.1007/BF02906837
Horie, S., Watanabe, Y., Tanaka, K., Nishiwaki, S., Fujioka, H., Abe, H., et al. (1998). The Schizosaccharomyces pombe mei4 + gene encodes a meiosis-specific transcription factor containing a forkhead DNA-binding domain. Mol. Cell Biol. 18, 2118–2129. doi: 10.1128/MCB.18.4.2118
Hu, J. L., Donahue, G., Dorsey, J., Govin, J., Yuan, Z. F., Garcia, B. A., et al. (2015). H4K44 acetylation facilitates chromatin accessibility during meiosis. Cell Reports. 13, 1772–1780. doi: 10.1016/j.celrep.2015.10.070
International Wheat Genome Sequencing Consortium. (2018). Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science 361, eaar7191. doi: 10.1126/science.aar7191
Jiao, Y., Wickett, N. J., Ayyampalayam, S., Chanderbali, A. S., Landherr, L., Ralph, E., et al. (2011). Ancestral polyploidy in seed plants and angiosperms. Nature 473, 97–100. doi: 10.1038/nature09916
Kawahara, Y., de la Bastide, M., Hamilton, J. P., Kanamori, H., McCombie, W. R., Ouyang, S., et al. (2013). Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice (N Y). 6, 4. doi: 10.1186/1939-8433-6-4
Kerzendorfer, C., Vignard, J., Pedrosa-Harand, A., Siwiec, T., Akimcheva, S., Jolivet, S., et al. (2006). The arabidopsis thaliana mnd1 homologue plays a key role in meiotic homologous pairing, synapsis and recombination. J. Cell Sci. 119, 2486–2496. doi: 10.1242/jcs.02967
Khoo, K. H. P., Able, A. J., Able, J. A. (2012). Poor homologous synapsis 1 interacts with chromatin but does not colocalise with ASYnapsis 1 during early meiosis in bread wheat. Int. J. Plant Genomics. 23, 514398. doi: 10.1155/2012/514398
Khoo, K. H. P., Jolly, H. R., Able, J. A. (2008). The RAD51 gene family in bread wheat is highly conserved across eukaryotes, with RAD51A1 upregulated during early meiosis. Funct. Plant Biol. 35, 1267–1277. doi: 10.1071/FP08203
Kim, M., Ahn, J. W., Song, K., Paek, K. H., Pai, H. S. (2002). Forkhead-associated domains of the tobacco NtFHA1 transcription activator and the yeast fhl1 forkhead transcription factor are functionally conserved. J. Biol. Chem. 277, 38781–38790. doi: 10.1074/jbc.M201559200
Klimyuk, V. I., Jones, J. D. G. (1997). AtDMC1, the Arabidopsis homologue of the yeast DMC1 gene: characterization, transposon-induced allelic variation and meiosis-associated expression. Plant J. 11, 1–14. doi: 10.1046/j.1365-313X.1997.11010001.x
Kobayashi, T., Hotta, Y., Tabata, S. (1993). Isolation and characterization of a yeast gene that is homologous with a meiosis-specific cDNA from a plant. Mol. Gen. Genet. 237, 225–232. doi: 10.1007/BF00282804
Kobayashi, T., Kobayashi, E., Sato, S., Hotta, Y., Miyajima, N., Tanaka, A., et al. (1994). Characterization of cDNAs induced in meiotic prophase in lily microsporocytes. DNA Res. 1, 15–26. doi: 10.1093/dnares/1.1.15
Krasileva, K. V., Vasquez-Gross, H. A., Howell, T., Bailey, P., Paraiso, F., Clissold, L., et al. (2017). Uncovering hidden variation in polyploid wheat. Proc. Natl. Acad. Sci. U.S.A. 114, 6. E913–E921. doi: 10.1073/pnas.1619268114
Krishnan, A., Gupta, C., Ambavaram, M. M. R., Pereira, A. (2017). RECoN: rice environment coexpression network for systems level analysis of abiotic-stress response. Front. Plant Sci. 8, 1640. doi: 10.3389/fpls.2017.01640
Kurdzo, E. L., Obeso, D., Chuong, H., Dawson, D. S. (2017). Meiotic centromere coupling and pairing function by two separate mechanisms in Saccharomyes cerevisiae. Genetics 205, 657–671. doi: 10.1534/genetics.116.190264
Lam, W. S., Yang, X., Makaroff, C. A. (2005). Characterization of Arabidopsis thaliana SMC1 and SMC3: evidence that AtSMC3 may function beyond chromosome cohesion. J. Cell Sci. 118, 3037–3048. doi: 10.1242/jcs.02443
Li, X., Chang, Y., Xin, X., Zhu, C., Li, X., Higgins, J. D., et al. (2013). Replication protein a2c coupled with replication protein a1c regulates crossover formation during meiosis in rice. Plant Cell 25, 3885–3899. doi: 10.1105/tpc.113.118042
Liang, J., Suhandynata, R. T., Zhou, H. (2015). Phosphorylation of Sae2 mediates forkhead-associated (FHA) domain-specific interaction and regulates its DNA repair function. J. Biol. Chem. 290, 10751–10763. doi: 10.1074/jbc.M114.625293
Liu, W., Lin, L., Zhang, Z., Liu, S., Gao, K., Lv, Y., et al. (2019). Gene co-expression network analysis identifies trait-related modules in Arabidopsis thaliana. Planta 249, 1487. doi: 10.1007/s00425-019-03102-9
Lloyd, A. H., Milligan, S. A., Langridge, P., Able, J. A. (2007). TaMSH7: a cereal mismatch repair gene that affects fertility in transgenic barley (Hordeum vulgare L.). BMC Plant Biol. 7, 67. doi: 10.1186/1471-2229-7-67
Lloyd, A. H., Ranoux, M., Vautrin, S., Glover, N., Fourment, J., Charif, D., et al. (2014). Meiotic gene evolution: can you teach a new dog new tricks? Mol. Biol. Evol. 31, 1724–1727. doi: 10.1093/molbev/msu119
Lu, P., Wijeratne, A. J., Wang, Z., Copenhaver, G. P., Ma, H. (2014). Arabidopsis ptd is required for type i crossover formation and affects recombination frequency in two different chromosomal regions. J. Genet. Genomics 41, 165–175. doi: 10.1016/j.jgg.2014.02.001
Luense, L. J., Wang, X. S., Schon, S. B., Weller, A. H., Shiao, E. L., Bryant, J. M., et al. (2016). Comprehensive analysis of histone post-translational modifications in mouse and human male germ cells. Epigenet. Chromatin. 9, 24. doi: 10.1186/s13072-016-0072-6
Luo, M. C., Gu, Y. Q., Puiu, D., Wang, H., Twardziok, S. O., Deal, K. R., et al. (2017). Genome sequence of the progenitor of the wheat D genome Aegilops tauschii. Nature 551, 498–502. doi: 10.1038/nature24486
Ma, J., Skibbe, D. S., Fernandes, J., Walbot, V. (2008). Male reproductive development: gene expression profiling of maize anther and pollen ontogeny. Genome Biol. 9, R18110. doi: 10.1186/gb-2008-9-12-r181
Ma, X., Zhao, H., Xu, W., You, Q., Yan, H., Gao, Z., et al. (2018). Co-expression gene network analysis and functional module Identiﬁcation in bamboo growth and development. Front. Genet. 9, 574. doi: 10.3389/fgene.2018.00574
Magnaghi-Jaulin, L., Jaulin, C. (2006). Histone deacetylase activity is necessary for chromosome condensation during meiotic maturation in Xenopus laevis. Chromosome Res. 14, 319–332. doi: 10.1007/s10577-006-1049-2
Marcussen, T., Sandve, S. R., Heier, L., Spannagl, M., Pfeifer, M., The International Wheat Genome Sequencing Consortium, et al. (2014). Ancient hybridizations among the ancestral genomes of bread wheat. Science 345 (6194), 1250092. doi: 10.1126/science.1250092
Martín, A. C., Borrill, P., Higgins, J., Alabdullah, A., Ramírez-González, R. H., Swarbreck, D., et al. (2018). Genome-wide transcription during early wheat meiosis is independent of synapsis, ploidy level, and the Ph1 Locus. Front. Plant Sci. 9 (1), 791. doi: 10.3389/fpls.2018.01791
Martinez-Perez, E., Shaw, P., Aragon-Alcaide, L., Moore, G. (2003). Chromosomes form into seven groups in hexaploid and tetraploid wheat as a prelude to meiosis. Plant J. 36, 21–29. doi: 10.1046/j.1365-313X.2003.01853.x
Mauri, M., Elli, T., Caviglia, G., Uboldi, G., Azzi, M. (2017). RAWGraphs: A Visualisation Platform to Create Open Outputs. CHItaly ‘17. 2017; https://doi.org/10.1145/3125571.3125585. [Software] Available at: http://app.rawgraphs.io/. doi: 10.1145/3125571.3125585
Palmbos, P. L., Wu, D., Daley, J. M., Wilson, T. E. (2008). Recruitment of Saccharomyces cerevisiae Dnl4-Lif1 complex to a double-strand break requires interactions with Yku80 and the Xrs2 FHA domain. Genetics 180, 1809–1819. doi: 10.1534/genetics.108.095539
Pawlowski, W. P., Golubovskaya, I. N., Timofejeva, L., Robert, B., Meeley, R. B., William, F., et al. (2004). Coordination of meiotic recombination, pairing, and synapsis by PHS1. Science 303 (5654), 89–92. doi: 10.1126/science.1091110
Pedersen, D. S., Coppens, F., Ma, L., Antosch, M., Marktl, B., Merkle, T., et al. (2011). The plant-specific family of DNA-binding proteins containing three HMG-box domains interacts with mitotic and meiotic chromosomes. New Phytol. 192, 577–589. doi: 10.1111/j.1469-8137.2011.03828.x
Pérez-Hidalgo, L., Moreno, S., San-Segundo, P. A. (2003). Regulation of meiotic progression by the meiosis-specific checkpoint kinase Mek1 in fission yeast. J. Cell Sci. 116, 259–271. doi: 10.1242/jcs.00232
Pierce, M., Benjamin, K. R., Montano, S. P., Georgiadis, M. M., Winter, E., Vershon, A. K. (2003). Sum1 and Ndt80 proteins compete for binding to middle sporulation element sequences that control meiotic gene expression. Mol. Cell Biol. 23, 4814–4825. doi: 10.1128/MCB.23.14.4814-4825.2003
Ramírez-González, R. H., Borrill, P., Lang, D., Harrington, S. A., Brinton, J., Venturini, L., et al. (2018). The transcriptional landscape of polyploid wheat. Science 361, eaar6089. doi: 10.1126/science.aar6089
Reddy, T. V., Kaur, J., Agashe, B., Sundaresan, V., Siddiqi, I. (2003). The duet gene is necessary for chromosome organization and progression during male meiosis in Arabidopsis and encodes a phd finger protein. Development 130, 5975–5987. doi: 10.1242/dev.00827
Rey, M. D., Martín, A. C., Higgins, J., Swarbreck, D., Uauy, C., Shaw, P., et al. (2017). Exploiting the ZIP4 homologue within the wheat Ph1 locus has identified two lines exhibiting homoeologous crossover in wheat-wild relative hybrids. Mol. Breed. 37, 95. doi: 10.1007/s11032-017-0700-2
Rey, M. D., Martín, A. C., Smedley, M., Hayta, S., Harwood, W., Shaw, P., et al. (2018). Magnesium increases homoeologous crossover frequency during meiosis in ZIP4 (Ph1 Gene) mutant wheat-wild relative hybrids. Front. Plant Sci. 9, 509. doi: 10.3389/fpls.2018.00509
Riley, R., Chapman, V., Johnson, R. (1968). The incorporation of alien disease resistance in wheat by genetic interference with the regulation of meiotic chromosomal synapsis. Genet Res. 12, 199–219. doi: 10.1017/S0016672300011800
Ronceret, A., Doutriaux, M. P., Golubovskaya, I. N., Pawlowski, W. P. (2009). PHS1 regulates meiotic recombination and homologous chromosome pairing by controlling the transport of RAD50 to the nucleus. Proc. Natl. Acad. Sci. U.S.A. 106, 20121–20126. doi: 10.1073/pnas.0906273106
Sakai, H., Lee, S. S., Tanaka, T., Numa, H., Kim, J., Kawahara, Y., et al. (2013). Rice Annotation Project Database (RAP-DB): an integrative and interactive database for rice genomics. Plant Cell Physiol. 54, e6. doi: 10.1093/pcp/pcs183
Scannell, D. R., Byrne, K. P., Gordon, J. L., Wong, S., Wolfe, K. H. (2006). Multiple rounds of speciation associated with reciprocal gene loss in polyploid yeasts. Nature 440, 341–345. doi: 10.1038/nature04562
Schreiber, A. W., Hayden, M. J., Kerrie, L., Forrest, K. L., Kong, S. L., Langridge, P., et al. (2012). Transcriptome-scale homoeolog-specific transcript assemblies of bread wheat. BMC Genomics. 13, 492. doi: 10.1186/1471-2164-13-492
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Shen, Y., Tang, D., Wang, K., Wang, M., Huang, J., Luo, W., et al. (2012). ZIP4 in homologous chromosome synapsis and crossover formation in rice meiosis. J. Cell Sci. 125, 2581–2591. doi: 10.1242/jcs.090993
Silva, A. T., Ribone, P. A., Chan, R. L., Ligterink, W., Hilhorst, H. W. (2016). A predictive co-expression network identifies novel genes controlling the seed-to-seedling phase transition in Arabidopsis thaliana. Plant Physiol. 170, 2218–2231. doi: 10.1104/pp.15.01704
Singh, M., Goel, S., Meeley, R. B., Dantec, C., Parrinello, H., Michaud, C., et al. (2011). Production of viable gametes without meiosis in maize deficient for an argonaute protein. Plant Cell 23, 443–458. doi: 10.1105/tpc.110.079020
Sozzani, R., Maggio, C., Varotto, S., Canova, S., Bergounioux, C., Albani, D., et al. (2006). Interplay between arabidopsis activating factors E2Fb and E2Fa in cell cycle progression and development. Plant Physiol. 140, 1355–1366. doi: 10.1104/pp.106.077990
Tan, M., Cheng, D., Yang, Y., Zhang, G., Qin, M., Chen, J., et al. (2017). Co-expression network analysis of the transcriptomes of rice roots exposed to various cadmium stresses reveals universal cadmium-responsive genes. BMC Plant Biol. 17, 194. doi: 10.1186/s12870-017-1143-y
Tessari, M. A., Gostissa, M., Altamura, S., Sgarra, R., Rustighi, A., Salvagno, C., et al. (2003). Transcriptional activation of the cyclin A gene by the architectural transcription factor HMGA2. Mol. Cell Biol. 23, 9104–9116. doi: 10.1128/MCB.23.24.9104-9116.2003
Usadel, B., Obayashi, T., Mutwil, M., Giorgi, F. M., Bassel, G. W., Tanimoto, M., et al. (2009). Coexpression tools for plant biology: opportunities for hypothesis generation and caveats. Plant Cell Environ. 32, 1633–1651. doi: 10.1111/j.1365-3040.2009.02040.x
Vilella, A. J., Severin, J., Ureta-Vidal, A., Heng, L., Durbin, R., Birney, E. (2009). EnsemblComparaGeneTrees: complete, duplication-aware phylogenetic trees in vertebrates. Genome Res. 19, 327–335. doi: 10.1101/gr.073585.107
Waese, J., Fan, J., Pasha, A., Yu, H., Fucile, G., Shi, R., et al. (2017). ePlant: visualizing and exploring multiple levels of data for hypothesis generation in plant biology. Plant Cell 29, 1806–1821. doi: 10.1105/tpc.17.00073
Watanabe, Y., Shinozaki-Yabana, S., Chikashige, Y., Hiraoka, Y., Yamamoto, M. (1997). Phosphorylation of RNA-binding protein controls cell cycle switch from mitotic to meiotic in fission yeast. Nature 386, 187–190. doi: 10.1038/386187a0
Watanabe, Y., Yamamoto, M. (1994). S. pombe mei2+ encodes an RNA-binding protein essential for premeiotic DNA synthesis and meiosis I, which cooperates with a novel RNA species meiRNA. Cell. 78, 487–498. doi: 10.1016/0092-8674(94)90426-X
Wijeratne, A. J., Chen, C., Zhang, W., Timofejeva, L., Ma, H. (2006). The Arabidopsis thaliana parting dancers gene encoding a novel protein is required for normal meiotic homologous recombination. Mol. Biol. Cell. 17, 1331–1343. doi: 10.1091/mbc.e05-09-0902
Xu, L., Ajimura, M., Padmore, R., Klein, C., Kleckner, N. (1995). NDT80, a meiosis-specific gene required for exit from pachytene in Saccharomyces cerevisiae. Mol. Cell Biol. 15, 6572–6581. doi: 10.1128/MCB.15.12.6572
Xu, Z., Song, Z., Li, G., Tu, H., Liu, W., Liu, Y., et al. (2016). H2B ubiquitination regulates meiotic recombination by promoting chromatin relaxation. Nucleic Acids Res. 44, 9681–9697. doi: 10.1093/nar/gkw652
Yan, Z., Fan, D., Meng, Q., Yang, J., Zhao, W., Guo, F., et al. (2016). Transcription factor ZFP38 is essential for meiosis prophase I in male mice. Reproduction 152, 431–437. doi: 10.1530/REP-16-0225
Yu, H., Jiao, B., Lu, L., Wang, P., Chen, S., Liang, C., et al. (2018). NetMiner-an ensemble pipeline for building genome-wide and high-quality gene coexpression network using massive-scale RNA-seq. PLoS ONE 13 (2), e0192613. doi: 10.1371/journal.pone.0192613
Zhang, J., Yan, G., Tian, M., Ma, Y., Xiong, J., Miao, W. (2018). A DP-like transcription factor protein interacts with E2fl1 to regulate meiosis in Tetrahymena thermophila.Cell Cycle. 17, 634–642. doi: 10.1080/15384101.2018.1431595
Zheng, N., Fraenkel, E., Pabo, C. O., Pavletich, N. P. (1999). Structural basis of DNA recognition by the heterodimeric cell cycle transcription factor E2F-DP. Genes Dev. 13, 666–674. doi: 10.1101/gad.13.6.666
Zhu, G., Spellman, P. T., Volpe, T., Brown, P. O., Botstein, D., Davis, T. N., et al. (2000). Two yeast forkhead genes regulate the cell cycle and pseudohyphal growth. Nature 406, 90–94. doi: 10.1038/35017581
Keywords: wheat, meiosis, hexaploid, polyploidization, co-expression, network
Citation: Alabdullah AK, Borrill P, Martin AC, Ramirez-Gonzalez RH, Hassani-Pak K, Uauy C, Shaw P and Moore G (2019) A Co-Expression Network in Hexaploid Wheat Reveals Mostly Balanced Expression and Lack of Significant Gene Loss of Homeologous Meiotic Genes Upon Polyploidization. Front. Plant Sci. 10:1325. doi: 10.3389/fpls.2019.01325
Received: 25 July 2019; Accepted: 24 September 2019;
Published: 18 October 2019.
Edited by:Ingo Ebersberger, Max F. Perutz Laboratories GmbH, Austria
Reviewed by:Jin Koh, University of Florida, United States
Kai Wang, Fujian Agriculture and Forestry University, China
Copyright © 2019 Alabdullah, Borrill, Martin, Ramirez-Gonzalez, Hassani-Pak, Uauy, Shaw and Moore. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Abdul Kader Alabdullah, AbdulKader.Alabdullah@jic.ac.uk