ORIGINAL RESEARCH article
Sec. Plant Development and EvoDevo
Volume 11 - 2020 | https://doi.org/10.3389/fpls.2020.619404
Barley Anther and Meiocyte Transcriptome Dynamics in Meiotic Prophase I
- 1Cell and Molecular Sciences, The James Hutton Institute, Dundee, United Kingdom
- 2School of Life Sciences, University of Dundee, Dundee, United Kingdom
- 3School of Agriculture and Wine, University of Adelaide, Adelaide, SA, Australia
In flowering plants, successful germinal cell development and meiotic recombination depend upon a combination of environmental and genetic factors. To gain insights into this specialized reproductive development program we used short- and long-read RNA-sequencing (RNA-seq) to study the temporal dynamics of transcript abundance in immuno-cytologically staged barley (Hordeum vulgare) anthers and meiocytes. We show that the most significant transcriptional changes in anthers occur at the transition from pre-meiosis to leptotene–zygotene, which is followed by increasingly stable transcript abundance throughout prophase I into metaphase I–tetrad. Our analysis reveals that the pre-meiotic anthers are enriched in long non-coding RNAs (lncRNAs) and that entry to meiosis is characterized by their robust and significant down regulation. Intriguingly, only 24% of a collection of putative meiotic gene orthologs showed differential transcript abundance in at least one stage or tissue comparison. Argonautes, E3 ubiquitin ligases, and lys48 specific de-ubiquitinating enzymes were enriched in prophase I meiocyte samples. These developmental, time-resolved transcriptomes demonstrate remarkable stability in transcript abundance in meiocytes throughout prophase I after the initial and substantial reprogramming at meiosis entry and the complexity of the regulatory networks involved in early meiotic processes.
Due to factors such as global population growth and climate change, food security is a major challenge which can be partly addressed by breeding new crop traits. Breeding relies almost entirely upon meiosis, a specialized form of cell division which generates haploid gametes from diploid parental meiocytes. During the first of two cell divisions in meiosis, parental homologous chromosomes condense, pair, synapse, and undergo recombination within a meiosis-specific protein complex called the synaptonemal complex (SC) (Osman et al., 2011; Naranjo, 2012; Zickler and Kleckner, 2015). The meiotic recombination machinery creates and segregates new allelic combinations which may produce novel traits in offspring. In the life cycle of flowering plants, this pivotal event occurs within specialized reproductive male (stamens) and female (carpel) organs. These reproductive organs are almost entirely composed of vegetative tissues embedding a small number of archesporial cells that differentiate to form meiocytes (Goldberg et al., 1993; Hord and Ma, 2007; Kelliher and Walbot, 2011). Understanding the development of these organs and the biology underpinning meiotic recombination could help to improve the efficiency of trait mapping; potentially leading to improvements in crop breeding.
Phylogenomic, mutagenic, and immuno-cytological tools have been extensively used to clone and characterize plant orthologs of meiotic genes (Lambing et al., 2017). Despite these advances, we still know little about the regulatory mechanisms of gene expression during meiosis. The transcriptomic profile of anthers before, during, and after meiosis has been extensively profiled in Arabidopsis (Honys and Twell, 2004; Chen et al., 2010; Yang et al., 2011; Dukowic-Schulze et al., 2014a), rice (Hobo et al., 2008; Fujita et al., 2010; Tang et al., 2010) and maize (Dukowic-Schulze et al., 2014b; Kelliher and Walbot, 2014; Zhang et al., 2014; Yuan et al., 2018; Nelms and Walbot, 2019). These studies provide substantial insights into meiotic transcriptome dynamics in these species with some observations in common. Several of these studies report that many meiotic genes, including ASY1 and DMC1, are both translated and transcribed in pre-meiotic cells—as early as in early archesporial cells in maize (Tang et al., 2010; Kelliher and Walbot, 2014; Zhang et al., 2014; Yuan et al., 2018). This appears to be a highly conserved phenomenon as accumulation of meiotic genes in pre-meiotic cells has also been shown using laser-capture and microdissection of fixed germ cells from human testis (Jan et al., 2017). Early accumulation of meiotic genes appears to be coupled with robust transcriptional changes which may occur at the mitosis-meiosis transition (Yuan et al., 2018), or in a two-step manner during leptotene and early zygotene in maize (Nelms and Walbot, 2019). Female meiocytes are harder to access and isolate in large numbers resulting in fewer studies of gene expression during female meiosis. Instead, Kubo et al. (2013) isolated entire rice ovules at different developmental stages from pre-meiosis to embryo sac and determined their gene expression profiles. This also revealed a major change of the transcriptome in the rice ovule, although later than in the anthers, occurring between zygotene/pachytene and later meiotic stages (Kubo et al., 2013). Other common observations of anther transcriptome dynamics are the significant enrichment of E3 ubiquitin ligases (Fujita et al., 2010; Tang et al., 2010; Yang et al., 2011; Yuan et al., 2018), and mitochondrial processes (Dukowic-Schulze et al., 2014a; Yuan et al., 2018; Nelms and Walbot, 2019) in developing germinal cells.
Although extensively conserved, the meiotic program and gene expression throughout meiosis vary between plant species. Dukowic-Schulze et al. (2014a) compared early meiotic transcripts in Arabidopsis and maize anthers finding low correlation in the expression profile of orthologous genes and a number of transcription factors whose expression was unique to Arabidopsis or to maize. However, higher correlation in expression profiles was observed when looking at well-defined meiotic genes (Dukowic-Schulze et al., 2014a). Recently, Bélanger et al. (2020) examined the transcriptome of anthers before, during, and after meiosis in both barley and wheat. This uncovered apparent divergence in pre-meiotic expression of a group of 24nt phased small interfering RNAs (phasiRNAs), which are expressed in barley and wheat anthers but not in rice and maize (Bélanger et al., 2020). The anther transcriptome has also been probed by RNA-seq in Brassica (Braynen et al., 2017) and sunflower (Flórez-Zapata et al., 2014, 2016) showing a high number of differentially expressed genes when meiocytes were compared to somatic tissue of the same plant or between the cultivated varieties and their wild relative. Meiotic DSB-associated genes were down-regulated in a synthetic autotetraploid Brassica compared to its diploid progenitor resulting in abnormal chromosomal segregation (Braynen et al., 2017). On the other hand, compared to its wild relative, meiotic recombination was shown to be higher in cultivated sunflower and many long non-coding RNAs (lncRNAs) were exclusively up-regulated in its meiocytes (Flórez-Zapata et al., 2016). Similarly, low fertility in autotetraploid rice was attributed to differential expression of lncRNAs during anther and ovule meiosis (Li et al., 2020). LncRNAs are poorly characterized compared to protein-coding genes but these studies highlight their potential importance during development (Ponting et al., 2009; Ariel et al., 2015; Shafiq et al., 2016) including their regulatory role during sexual reproduction (Golicz et al., 2018). The processing of a subset of lncRNAs has been shown to generate small RNAs in animals (Jalali et al., 2012) and plants (Ma et al., 2014). Although the regulation of gene expression by small RNA has been reported in cotton anthers (Chen et al., 2018), Arabidopsis meiocytes (Huang et al., 2019) and rice ovules (Li et al., 2020), the full integration of small RNAs and lncRNAs biogenesis and function remains to be determined.
Meiocytes are represented by a small number of cells (Kelliher and Walbot, 2011; Zhang et al., 2014) embedded within complex tissues of anthers and ovules making their isolation difficult and time consuming. Several groups have isolated these rare cells and studied their transcriptome. However, most of these studies were restricted to pre-meiosis to zygotene, excluding later meiotic stages (Tang et al., 2010; Dukowic-Schulze et al., 2014a; Kelliher and Walbot, 2014; Zhang et al., 2014; Yuan et al., 2018; Nelms and Walbot, 2019). Unlike the studies where male pre-meiotic cells were recovered by means of microdissection (Tang et al., 2010; Kelliher and Walbot, 2014; Zhang et al., 2014) or protoplasts (Nelms and Walbot, 2019), generally careful staging was not performed for prophase I samples resulting in pooled meiotic stages instead (Chen et al., 2010; Yang et al., 2011; Flórez-Zapata et al., 2014, 2016). Although barley and wheat whole anther transcriptomic profiling has been recently carried out, this did not include specific profiling of expression in meiocytes (Bélanger et al., 2020).
To gain further insight into transcriptional and post-transcriptional regulation in barley reproductive tissues we combined precise immuno-cytological staging with RNA sequencing to analyze transcript profiles of barley anther and meiocyte tissues before, during, and after prophase I. Our data reveal that in barley, as in other plant species, there is extensive transcriptional reorganization between pre-meiosis and early prophase I in anthers. In barley, this switch appears to consist predominantly of down-regulation of many lncRNAs. Sizeable down-regulation of lncRNAs was also observed in the comparison of meiocytes to whole anther tissues in early prophase I. This adds further context to the observed enrichment of lncRNAs in sunflower meiocytes compared to somatic tissues (Flórez-Zapata et al., 2016). As prophase I progressed into metaphase I–tetrad, significant changes to transcript abundance were substantially reduced at each successive stage. In addition, meiocytes showed remarkably stable transcript abundance between leptotene–zygotene and pachytene–diplotene. Further, differential and co-expression analyses indicate the enrichment of several regulatory pathway components in meiocytes in prophase I. These included pentatricopeptide repeat (PPR) RNA binding proteins, E3 ubiquitin ligases, lys48 and NEDD8 specific deubiquitinating enzymes, and small RNA processing components, matching previous observations of enrichment of mitochondrial activity and the ubiquitin proteasome system during meiosis in rice and maize (Tang et al., 2010; Yuan et al., 2018). The details captured in our transcript abundance datasets provide significant insights into this highly organized, stepwise biological process and will lead to both greater understanding and possible applied outcomes.
Materials and Methods
Barley cv. Golden Promise plants were grown in cereal compost in a controlled growth room at 18°C for 16h light and 14°C for 8 h dark and 70% humidity. When plants reached the desired stage (pre-meiosis/meiosis) at 5–7 weeks post-germination, they were processed individually by collecting their spikes at the last minute before anther sampling to minimize stress effects.
Male meiotic prophase I in barley progresses over an approximate 36 h time period inside the anthers of developing spikelets (Higgins et al., 2012; Barakate et al., 2014), that are themselves buried within the extending barley culm. To collect appropriately staged meiotic inflorescences, or spikes, we first developed a reliable prediction system based on external morphological characteristics of the main tiller as described for the cultivar Optic (Gómez and Wilson, 2012). In our environmentally conditioned growth room, the barley cultivar Golden Promise plants reached meiosis at 6 weeks post-germination, when the flag leaf was emerging, and the spikes were 0.5–2 cm in length. As we had previously determined a strong correlation between anther length and meiotic progression (Arrieta et al., 2020), anthers were dissected under a stereomicroscope and measured before being stored as length classes in RNAlaterTM (Thermo Fisher Scientific). For each length class, a subset of 10 anthers was stored in PBS buffer and used for chromosomal acetocarmine staining for meiotic staging. This allowed us to determine the precise correlation between anther size and meiotic stage in this material.
Surfaces and tools were cleaned with 70% Ethanol and RNase AWAYTM (Thermo Fisher Scientific) prior to material collection and all steps were done wearing gloves. The collected spikes were placed in a Petri dish with wet filter paper (using RNase-free milli-Q water) on a cool pack or ice during the sampling, including under the stereomicroscope with LED light to avoid RNA degradation. Fine tweezers were used to remove awns and collect anthers of all florets between the fifth and fifteenth position from the bottom under a stereomicroscope resulting in 60 anthers on average, a process that takes about 10 min. All materials were collected in the afternoon between 1:00 p.m. and 3:00 p.m. only to improve sample homogeneity. Four stages were collected: pre-meiosis–G2, leptotene–zygotene, pachytene–diplotene and all stages from metaphase I to tetrad.
Each barley floret contains three synchronized anthers. One anther was staged using acetocarmine staining and the two others were collected in 1.5 ml Eppendorf tube containing 100 μl RNAlaterTM. Anthers of the same length and meiotic stage were pooled together. From this tube, 10 random anthers were transferred into a new tube with 1x PBS buffer for staging validation with immunostaining using meiosis specific antibodies (Supplementary Method 1). The process was repeated until at least three replicates of 60–200 anthers were collected for each meiotic stage.
To prepare meiocyte samples, 30–50 anthers that were staged as leptotene–zygotene and pachytene–diplotene by immuno-cytology were used to release meiocyte clusters onto microscope slides (Supplementary Method 2). A sample of these were stained with DAPI to evaluate the purity of the preparation (Figures 1O,P and Supplementary Videos 1, 2) and the remaining material was transferred into Eppendorf tubes containing TRIzol® for subsequent total RNA isolation. Similar to anther samples, meiocyte clusters at close stages of leptotene and zygotene or pachytene and diplotene were mixed in equal proportions. The duration of entire meiocyte isolation process for each replicate depends on anther stage and took on average 1.5 h.
Figure 1. Anther and meiocyte collection and staging. Anthers were collected according to their size and staged with both acetocarmine (A,C,F,I,K,M) and immuno-cytology (B,D,E,G,H,J,L,N). (A,B) Pre-meiotic stage; (C) leptotene/zygotene stage; (D) leptotene stage; (E) zygotene stage; (F) pachytene stage; (G) pachytene stage; (H) diplotene stage; (I,J) metaphase I; (K,L) anaphase I; (M,N) tetrads; (O,P) isolated fresh meiocyte clusters at leptotene/zygotene and pachytene/diplotene, respectively. Acetocarmine (Gray), DAPI (Blue and White for images O,P), ASY1 (Green), ZYP1 (Magenta). Scale bar 5 μm. Note images (O,P) have associated videos.
For comparison we included four samples of germinating embryo (EMB) as those have been shown to represent a wide range of transcripts [International Barley Genome Sequencing Consortium (IBGSC), 2012]. Barley cultivar (cv.) Golden Promise seeds were germinated in a Petri dish on wet paper for 5 days in dark. Three EMB were collected per sample by removing the residual seed material. In total four samples were collected and used for total RNA extraction.
Cytology and Imaging
Ten anthers of each pool were fixed in 4% formaldehyde (1x PBS/0.5% TritonTM X-100/0.5% Tween 20) for 20 min, rinsed twice in 1x PBS and ground gently in the tube with a disposable pestle. Tubes were centrifuged at 4,500 rpm for 1 min and 20 μl meiocytes suspension was transferred onto a PolysineTM 2 (Poly-L-Lysine coated slides) and left to air dry. The immuno-cytology (Supplementary Method 1) was done according to Colas et al. (2017).
Barley Anther Transcriptome (BAnTr) Construction
Total RNA was extracted from the 12 anther samples, six meiocyte samples and 4 EMB and sequenced using both Illumina NextSeq 550 and PacBio Sequel v1 platforms (Supplementary Method 3).
Illumina Read Mapping
The resulting Illumina sequencing reads were mapped individually against the barley cv. Golden Promise assembly using STAR 2.7.1a (Dobin et al., 2013) in a two-step process based on Veeneman et al. (2016). The first genome index was generated using default settings. The exact parameters for the two-step mapping process can be found in Supplementary Method 4.
PacBio Read Mapping
All three samples of long-read transcriptome data were processed individually using the SMRT Analysis software package v3.1.01 for building of the consensus circular reads (CCS) and demultiplexing which included barcode and primer removal. Afterward all transcripts were merged for the clustering and polishing. Polishing on the merged transcripts was done using proovread v.2.14.1 (Hackl et al., 2014) with the Illumina reads from the 12 anther samples and the 6 meiocyte samples. The polished transcripts were mapped to the Golden Promise reference assembly using GMAP version 2018-07-04 (Watanabe and Wu, 2005).
We used the Mikado pipeline (Mikado v1.2.4; Venturini et al., 2018) to join the different strands of transcriptome evidence, generated above. As additional transcript evidence from a more complete transcriptome we included the BaRTv1 transcripts (Rapazote-Flores et al., 2019). The BaRTv1 transcripts were mapped to the Golden Promise reference assembly (GMAP version 2018-07-04; -n 0 –min-trimmed-coverage = 0.80 –min-identity = 0.90). Input files for the Mikado file included a splice junction file generated by Portcullis (Mapleson et al., 2018; default parameters), open reading frame identification of the transcripts by TransDecoder2 (default parameters; Haas et al., 2013) and blast results using DIAMOND BLASTx (Buchfink et al., 2015; –evalue 1e-5) against the NCBI non-redundant protein database as evidence. Together with the stringtie.gtf, scallop.gtf, pacbio.gtf and bartv1.gtf those were combined and scored, generating the Barley Anther Transcriptome (BAnTr). This dataset comprises 65,795 genes and 119,778 transcripts (Supplementary File 1: AntherTranscriptomeBAnTr.fasta at http://doi.org/10.6084/m9.figshare.12136773).
A padded version of the transcriptome was generated for the transcript quantification as previous experiments in Arabidopsis showed this to improve quantification and provide more similar results to high-resolution real-time PCRs (Zhang R. et al., 2017; Supplementary File 2: AntherTranscriptomeBAnTrPadded.fasta at http://doi.org/10.6084/m9.figshare.12136773).
To distinguish between protein-coding and long non-coding transcripts, all transcripts were run through the NCBI ORF finder. Transcripts and corresponding genes were assigned as coding, if one hit was identified in a DIAMOND BLASTp (v0.9.24; Buchfink et al., 2015) search [–query-cover 60 –evalue 1e-30] against the NCBI non-redundant protein database (Supplementary File 3: AntherProteomeBAnTr.fasta at http://doi.org/10.6084/m9.figshare.12136773). Additionally, transcripts were checked using CPC2 (coding potential calculator 2, run on both strands) and InterProScan (v5.39; default parameters; Jones et al., 2014). If both supported a coding sequence, those were also included in the coding transcripts; if only one of those provided evidence for coding, the sequences were allocated to the unclassified group. All remaining transcripts and their corresponding genes were assigned as long non-coding RNA. For further classification we used the approach from Zhang Z. et al. (2020) by comparing the location of lncRNA to the location of protein-coding RNA using Cuffcompare (Trapnell et al., 2012). In addition, we also compared the lncRNA location with the transposable element location. Based on this lncRNA were then classified into TE-lncRNA, intergenic lncRNA (cuffcompare class code u), intronic lncRNA (cuffcompare class code i), sense lncRNA (cuffcompare class code o), antisense lncRNA (cuffcompare class code x) and other lncRNA (remaining cuffcompare class codes). Conservation of BAnTr lncRNA transcripts was assessed by alignment using BLAST (Altschul et al., 1990) to: RiceLncPedia transcripts (Zhang Z. et al., 2020); Arabidopsis non-code v5 transcripts (Zhao et al., 2016); maize, wheat, and Brachypodium GreeNC transcripts (Paytuví Gallart et al., 2016); and barley CANTATAdb v2.0 transcripts (Szcześniak et al., 2019).
To study differential gene expression the 12 anther and 6 meiocyte Illumina sequencing samples were mapped against our above described padded BAnTr transcriptome using Salmon v.0.14.1 (Patro et al., 2017; -validateMappings –useVBOpt –seqBias –gcBias –posBias). The subsequent analysis was conducted using the 3D RNA-seq pipeline (Guo et al., 2019; accessed 28/11/2019). This pipeline is a combination of different R packages. Transcript abundance is imported using the lengthScaledTPM method from tximport (Soneson et al., 2016) and transformed into gene level counts or kept as transcript level counts. Low expressed transcripts with a count per million reads (CPM) below 1 in fewer than three samples were removed. This reduced the dataset to 31,918 genes and 50,861 transcripts. Gene and transcript read counts were normalized using the trimmed mean of M values (TMM) method, implemented in the edgeR package (Robinson et al., 2010; Robinson and Oshlack, 2010; McCarthy et al., 2012). Based on the experimental design, six contrast groups were built: anthers at leptotene–zygotene versus anthers at pre-meiosis; anthers at pachytene–diplotene versus anthers at leptotene–zygotene; anthers at metaphase I–tetrad versus anthers at pachytene–diplotene; meiocytes at leptotene–zygotene versus anthers at leptotene–zygotene; meiocytes at pachytene–diplotene versus anthers at pachytene–diplotene; and meiocytes at pachytene–diplotene versus meiocytes at leptotene–zygotene. Differential expression analysis was done using the voom function in the limma R package (Law et al., 2014; Ritchie et al., 2015). Genes and transcripts were rated as differentially expressed with a log2-fold change (log2FC) above 1 or below −1 and an adjusted p-value (by the Benjamini-Hochberg method; Benjamini and Hochberg, 1995) of below 0.01. To analyse the relationship between replicates, a multidimensional scaling plot was generated using the transcript expression data.
Differential Alternative Splicing Analysis
Differential alternative splicing (DAS) analysis is part of the 3D RNA-seq pipeline. It integrates the diffSplice function from the limma R package and compares the log2FC of each individual transcript with the gene level log2FC. Transcripts with a Δ percent spliced (ΔPS) ratio of above 0.1 were assigned as significant alternative spliced transcripts. To determine the alternative splicing events, we used SUPPA version 2.3 (Alamancos et al., 2015; Trincado et al., 2018).
To identify co-expressed networks, we used the weighted gene co-expression network analysis (WGCNA) R package (v1.63; Langfelder and Horvath, 2008). The same reduced dataset of 31,918 expressed genes as introduced above was used as input. For the network construction an approximate scale-free topology of above 0.85 was achieved with a soft power of 10. Network construction and module assignment was done with the following settings: a signed hybrid network, mergeCutHeight of 0.35, a minModuleSize of 30 and a deepSplit of 4.
Gene Ontology and Functional Annotation
Multiple strands of evidence were combined for the functional annotation of the proteins. Mercator4 v.2 (Schwacke et al., 2019) was used to classify the proteins into functional bins. eggNOG (Huerta-Cepas et al., 2017; Huerta-Cepas et al., 2019) was used to add COG (Clusters of Orthologous Groups of proteins) annotation and a protein description. PANNZER2 (Törönen et al., 2018) was used for the GO term identification and a protein description. The GO enrichment analysis was done using the topGO R package using the fisher weight01 algorithm and a P-value cut-off below 0.001 (Alexa et al., 2006). TopGO only outputs an exact P-value until 1e-30. Everything below will be given as P < 1e-30. For visual presentation in the GOplot R package (Walter et al., 2015) those were set to 1e-30.
Accurate Staging of Anthers
As meiosis is not completely synchronized along the length of a spike (Lindgren et al., 1969), only anthers from the middle spikelets were collected. In total, three biological replicates with at least 120 anthers in each were collected for pre-meiosis (0.3–0.4 mm anthers), leptotene–zygotene (0.5–0.9 mm anthers), pachytene–diplotene (1.0–1.2 mm anthers) and metaphase I–tetrad (1.3–1.4 mm anthers) (Figure 1).
Acetocarmine stained pre-meiotic spreads showed small nuclei (Figure 1A) without any labeling by antibodies against the SC proteins HvZYP1 and TaASY1, although a signal can be detected with anti-TaASY1 mainly in the cytoplasm (Figure 1B). Contrary to the pre-meiotic stages, anther size and acetocarmine staining were not entirely reliable for leptotene–zygotene staging (Figure 1C) but immunostaining with antibodies against TaASY1 and HvZYP1 allowed us to easily distinguish between these two close developmental stages. At leptotene, linear ASY1 axes are formed while there is no or little labeling of HvZYP1 at the initiation of synapsis (Figure 1D). Zygotene on the other hand was determined by the presence of HvZYP1 labeling during synapsis progression as previously described (Colas et al., 2017; Figure 1E). Pachytene and diplotene were somewhat easier to find based on anther size and acetocarmine spreads (Figure 1F) and were distinguished using antibodies against TaASY1 and HvZYP1 proteins (Figures 1G,H). Pachytene (Figure 1G) was characterized by the polymerization of HvZYP1 between the ASY1 axes along the entire chromosomes and diplotene (Figure 1H) was identifiable by the appearance of tinsel-like chromosomal structures as previously described (Colas et al., 2017). However, this immuno-cytological work revealed that meiocytes within the same anther were not fully synchronized and we opted for combining developmentally close samples to make mixtures of equal proportions of leptotene and zygotene or pachytene and diplotene. Later stages from metaphase I to tetrads were combined into a single sample (Figure 1). Metaphase I is characterized by the presence of seven ring bivalents (Figures 1I,J) while anaphase I (Figures 1K,L) and anaphase II samples show chromosomes or chromatin segregation, respectively (Colas et al., 2017). Anther samples of these later stages were not as homogenous as those at prophase I (leptotene–diplotene), containing both metaphase I and anaphase I or both anaphase II and telophase II. Very few anthers contained tetrads (Figures 1M,N) and none reached pollen stage.
To prepare meiocyte samples, a subset of anther samples that were staged as leptotene–zygotene and pachytene–diplotene by immuno-cytology were used to release meiocyte clusters onto microscope slides. A sample of these were stained with DAPI to evaluate the purity of the preparation (Figures 1O,P and Supplementary Videos 1, 2) and the remaining material was transferred into Eppendorf tubes containing TRIzolTM for subsequent total RNA isolation. The RNA integrity numbers (RIN) were in the range of 5.8–7.3 and 7.3–9.7 for meiocyte and anther samples, respectively. Similar to anther samples, meiocyte clusters at close stages of leptotene and zygotene or pachytene and diplotene were mixed in equal proportions.
The combination of anther size and immuno-cytology allowed us to delineate different developmental stages. We extracted total RNA and sequenced 18 samples (3 biological replicates of 2 meiocyte and 4 anther stages) generating a total of 2.1 billion reads (Supplementary Table 1). For comparative purposes we extracted and sequenced RNA from 4 replicate germinating embryo (EMB) samples generating over 245 million reads.
Creation of a reproduction-inclusive barley transcriptome reference assembly For RNA-seq and differential gene expression analysis we first developed a reference Barley Anther Transcriptome (BAnTr). We did this for two important reasons. First, we anticipated that anther/meiocyte tissues potentially contain a set of very specific transcripts. So, to complement the Illumina short read RNA-seq data, we generated and included three 5′ cap- and 3′ polyA-captured PacBio IsoSeq datasets from a mixture of the RNA samples mentioned above. Second, we are working with the transformation reference barley cultivar Golden Promise and had recently generated a Golden Promise genome reference assembly (Schreiber et al., 2020). Building a reference-guided transcriptome based on the same cultivar will be more complete and more accurate as cultivar specific transcripts will be included. To include as many barley transcripts as possible, the Illumina RNA-seq and IsoSeq reads of anthers and meiocytes were combined with the short reads from the germinated embryo samples and ultimately with the BaRTv1-Reference Transcript Dataset reported recently by Rapazote-Flores et al. (2019). The final transcriptome, BAnTr, comprised 65,795 genes and 119,778 transcripts.
For each sample, transcript quantification files were generated using Salmon (Patro et al., 2017) in conjunction with BAnTr. The output quantification files were then read into the 3D RNA-seq pipeline (Guo et al., 2019) to generate read counts which were converted to transcript per million (TPM) using the tximport R package (Soneson et al., 2016). The data used as input for the 3D RNA-seq pipeline has been deposited at figshare: https://doi.org/10.6084/m9.figshare.11974182. Read counts and TPMs were pre-processed and filtered to reduce noise and technical variance by excluding low abundance transcripts (keeping transcripts with a count per million reads ≥ 1 in at least 3 samples, and genes where at least one transcript passed the filtering). After filtering, 50,861 transcripts encoded by 31,918 genes remained. A multidimensional scaling plot was generated showing clear separation between samples (Supplementary Figure 1). Using differentially expressed genes (DEGs) criteria described in section “Materials and Methods,” a total of 10,713 DEGs (Supplementary Figures 2A,B) were identified and for further analysis divided into three categories based on protein-coding potential: protein-coding, long non-coding RNA (lncRNA) (Supplementary Figure 2C), and unclassified (see section “Materials and Methods”). phasiRNAs were detected in isolated meiocytes and anthers of maize (Zhai et al., 2015; Dukowic-Schulze et al., 2016). Using several male sterile mutants defective in specific anther cell layers, Zhai et al. (2015) revealed the spatiotemporal dynamics of 21 and 24 nt phasiRNAs highly accumulated at pre-meiotic and meiotic stages, respectively. PhasiRNAs have been shown to function in DNA methylation during anther development playing a critical role in fertility (Liu C. et al., 2020; Zhang M. et al., 2020). More recently, Teng et al. (2020) demonstrated that Dicer-like 5 generates the 24 nt phasiRNAs in meiotic anthers of maize. Production of phasiRNAs is triggered by cleavage of non-coding precursor (PHAS) transcripts directed by micro RNAs (miRNAs) miR2118 and miR2275, to produce 21 and 24 nt phasiRNAs, respectively (Johnson et al., 2009; Song et al., 2012a). To identify putative PHAS transcripts, all lncRNA transcripts were assessed for miR2118 and miR2275 recognition sites using psRNATarget (Dai et al., 2018). This identified 568 lncRNAs with predicted miR2118 recognition sites, 553 lncRNAs with predicted miR2275 recognition sites, and 53 predicted to contain recognition sites for both miRNAs. Matching the overall pattern of lncRNA expression, putative targets of both miR2118 and mi2275 were predominantly expressed in anthers and not in germinating embryos (93 and 86%, respectively). Expression of putative miR2275 targets (Supplementary Figure 2E) was more stable than expression of putative miR2118 targets (Supplementary Figure 2D) with 74 and 46% of each showing differential expression in any comparison, respectively. However, putative targets of these miRNAs which were differentially expressed displayed a broadly similar pattern of expression across sample stages (Supplementary Figures 2F,G), predominantly characterized by down-regulation from pre-meiosis to pachytene–diplotene. A small proportion—1% of each—of these lncRNAs were enriched in meiocytes compared to whole anthers.
Several analyses were used to further probe the functional pathway and developmental expression of this high number of DEGs. However, in our analysis we focused on meiotic recombination and gene regulation pathways.
Co-expression Analysis Detected Four Modules Enriched in Meiocyte Samples
We used the weighted gene co-expression network analysis (WGCNA) tool (Langfelder and Horvath, 2008) to identify co-expression networks within all 31,918 expressed genes in the anther and meiocyte dataset. The analysis identified 17 gene clusters (shown as color-coded modules in Supplementary Figure 3) with an additional un-clustered gray module (Supplementary Table 2). The number of genes in each module ranged from 38 to 6209 (Supplementary Table 2). Variation was also evident in differential eigengene network analysis showing levels of correlation between clusters (Supplementary Figure 4). WGCNA navy and yellow modules were enriched in the meiocyte samples in comparison to anthers. In total, four modules (navy, yellow, green, and olive) showed higher expression in isolated meiocytes in comparison to anthers (Figure 2). Gene ontology (GO) enrichment analysis highlighted mRNA and RNA binding GO terms as significantly enriched in the navy module (P = 6E-16 and P = 3.1E-13); the yellow module was significantly enriched for DNA binding (P < 1E-30), RNA modification, and nucleosome GO terms (Figure 3 and Supplementary Table 3) (P < 1E-30 and P < 1E-30). From the remaining modules, genes clustered in the red module showed a pattern of expression which was approximately inverse to genes clustered in the navy module (Supplementary Figure 3). The red module is significantly enriched in genes located in chloroplast components (Supplementary Figure 5 and Supplementary Table 3) (chloroplast thylakoid membrane, P = 2.2E-13 and chloroplast envelope, P = 3.7E-6). The blue and the light-green modules show a similar downward trend of expression with meiotic progression (Supplementary Figure 3).
Figure 2. WGCNA analysis of co-expressed genes. A total of 17 modules were found and the selected four show an interesting pattern for meiocyte enriched genes. The samples (3 replicates each) are A.PRE, anther pre-meiosis; A.LEP-ZYG, anther leptotene–zygotene; M.LEP-ZYG, meiocyte leptotene–zygotene; A.PAC-DIP, anther pachytene–diplotene; M.PAC-DIP, meiocyte pachytene–diplotene; A.MET-TET, anther metaphase I–tetrad. The prefixes A. and M. in the sample names depict anther and meiocyte samples, respectively.
Figure 3. Gene ontology enrichment analysis. Results from the navy and the yellow modules of the WGCNA. Size of bubbles correspond to total number of proteins associated with the GO term.
Large Scale Transcriptional Re-programming in Anthers at Meiosis Entry Is Followed by Progressive Reduction of Significant Changes Into Metaphase–Tetrad Stage
3D RNA-seq analysis revealed a dramatic alteration in gene expression in anthers between pre-meiotic and leptotene–zygotene stages followed by progressively fewer significant changes in pachytene–diplotene and metaphase–tetrad stages (Figures 4A–C and Supplementary Table 4). The number of DEGs was highest (n = 6,119) in comparison of anthers at the leptotene–zygotene stage versus pre-meiosis. Strikingly, 48% of the total number of DEGs and 65% of the down-regulated DEGs are lncRNAs in this comparison (Figures 4A,B). Only 8.6% (n = 2,024) of these lncRNA transcripts returned an alignment with rice, maize, Arabidopsis, wheat, or Brachypodium lncRNA database sequences. When BAnTr lncRNA transcripts were aligned to barley CANTATAdb lncRNA data transcripts the number returning alignment increased to 30.45%. No clear relationship between lncRNA classification and successful alignment was observed (Supplementary Table 5). Compared with germinating embryos, anthers contain a much larger number of lncRNAs represented mainly in bleu module (Figures 5A,B). As expected, the majority of these lncRNAs are relatively short (Figure 5C) and comprise a single exon (Figure 5D). In anther tissues, there was a further pronounced, but much smaller, change in gene expression at pachytene–diplotene versus leptotene–zygotene (Figures 4A,B). However, gene expression in meiocytes appears stable between leptotene–zygotene and pachytene–diplotene stages, with only 4 DEGs (2 up-regulated; 2 down-regulated) observed in this comparison (Figure 4C). Only 164 genes (163 up-regulated; 1 down-regulated) are significantly differentially expressed in anthers at the metaphase I–tetrad stage compared to the pachytene–diplotene stage (Figures 4A,B). When lncRNAs are excluded from the analysis, the largest number of DEGs occurs in meiocytes at the leptotene–zygotene stage when compared to anthers at the same stage (n = 3,517; Figure 4C); the number of protein-coding DEGs in comparison of anthers at the leptotene–zygotene stage versus pre-meiosis is 2710 (Figure 4B).
Figure 4. Comparisons of differential gene expression. (A) Alluvial plot showing up- or down-regulation of genes significantly differentially expressed in all anther (A) meiotic stage comparisons. The comparisons are: A.LepZyg-A.Pre, leptotene–zygotene versus pre-meiosis; A.PacDip-A.LepZyg, pachytene–diplotene versus leptotene–zygotene; and A.MetTet-A.PacDip, metaphase I–tetrad versus pachytene–diplotene. The bar at each stage represents the total complement of genes which are differentially expressed in any anther comparison, the bar is colored by the proportion of genes in each successive meiotic stage comparison on the x axis which is up-regulated (yellow), down-regulated (pink), or not significant (blue). The lines connecting the bars between x axis comparisons represent movement of genes from one expression category to another between successive meiotic stage comparisons. The vast majority of all DEGs (pink and yellow) at any stage are present in the comparison of leptotene–zygotene anthers to pre-meiotic anthers. A large proportion of DEGs in the first comparison do not show a significant change in expression between pachytene–diplotene and leptotene–zygotene, as represented by the pink and yellow lines connecting to the blue shaded area in the second x-axis comparison. (B) UpSet plot of the same anther meiotic stage comparisons and (C) UpSet plot of meiocyte (M) stage and meiocyte-anther tissue comparisons. UpSet plots are similar in principal to the more commonly used Venn and Euler diagrams in that they represent the intersection of distinct groups (sets) but are better at representing large numbers of groups (Lex et al., 2014). The total number of genes in each group (the set size) is represented by a histogram on the left-hand side of the x-axis. The black dots underneath the x-axis indicate which groups are represented in the histogram above showing the number of genes represented by this group. Where only one dot is present the histogram above represents the number of genes unique to this group. Where two or more dots are connected by a line the histogram above represents the number of genes which occur in the connected groups. The comparisons are: M.PacDip-M.LepZyg, meiocytes at pachytene–diplotene versus meiocytes at leptotene–zygotene; M.LepZyg-A.LepZyg, meiocytes at leptotene–zygotene versus anthers at the same stage; M.PacDip-A.PacDip, meiocytes at pachytene–diplotene versus anthers at the same stage. The prefixes A. and M. in the sample names depict anther and meiocyte samples, respectively.
Figure 5. Differential expression and alternative splicing of different gene categories. (A) Compared with germinating embryos, anthers are enriched in lncRNAs. (B) Distribution of coding, non-coding and undefined genes in different modules. The modules Violet, Turquoise, Navy, Green, Yellow and Olive are enriched in meiocytes. (C) Length and (D) exon number distributions of different transcript categories. (E) Distribution of alternative splicing events for the entire dataset and the DAS (differential alternative spliced) genes. A3, alternative 3′ splice-site; A5, alternative 5′ splice-site; AF, alternative first exon; AL, alternative last exon; MX, mutually exclusive exons; RI, retained intron; SE, skipping exon. (F) Differential alternative spliced genes. The comparisons are: A.LepZyg-A.Pre, anther leptotene–zygotene versus anther pre-meiosis; A.PacDip-A.LepZyg, anther pachytene–diplotene versus anther leptotene–zygotene; A.MetTet-A.PacDip, anther metaphase I–tetrad versus anther pachytene–diplotene; M.LepZyg-A.LepZyg, meiocyte leptotene–zygotene versus anther leptotene–zygotene; M.LepZyg-M.LepZyg, meiocyte leptotene–zygotene versus meiocyte leptotene–zygotene; M.PacDip-A.PacDip, meiocyte pachytene–diplotene versus anther pachytene–diplotene. The prefixes A. and M. in the sample names depict anther and meiocyte samples, respectively.
Differential Alternative Splicing
A striking feature of many of the genes that are putatively involved in meiosis and recombination in barley is the number of introns and exons they contain. The average of our protein-coding transcripts in the BAnTr transcriptome is 7.5 exons per gene. Based on the literature and phylogenetic analysis, we compiled a list of 121 barley genes with confirmed role in meiosis in yeast, animal and plant species. Our selection of 121 potential meiotic genes has an average of 25 exons per gene—more than three times the average for all genes. One possible reason may be that regulation of these genes occurs at the post-transcriptional level through delayed pre-mRNA processing or DAS. Indeed, 84 out of the 121 potential meiotic genes have more than one isoform. However, only two of those genes are significantly differentially alternatively spliced (HvHSP90.7 and HvPTB1a). HvHSP90.7 has two main expressed isoforms throughout meiosis, BAnTr.GP.7HG019194.1 and BAnTr.GP.7HG019194.16 with the first showing no significant expression changes and the second being significantly upregulated in the step from pre-meiosis to leptotene-zygotene. BAnTr.GP.7HG019194.16 is characterized by an alternative 3′ splice site of 11 nucleotides which causes a premature stop codon. HvPTB1a has multiple isoforms expressed throughout meiosis, but only two transcripts (BAnTr.GP.4HG005942.17 and BAnTr.GP.4HG005942.18) are differentially expressed and one BAnTr.GP.4HG005942.1 shows differential transcript usage. BAnTr.GP.4HG005942.18 is characterized by an additional exon and BAnTr.GP.4HG005942.17 retains an intron. Both cases result in premature stop codons.
The analysis of DAS genes across the whole transcriptome revealed that most contain isoforms with different 3′ or 5′ splice-sites, representing 48% of the total identified events (Figure 5E). Only 0.8% of the genes contained isoforms with mutually exclusive exons. The DAS genes were similarly distributed with no splicing event being strongly enriched in comparison to the whole dataset. GO enrichment analysis showed for both groups of DAS genes, i.e., with alternative 3′ or 5′ splice-sites, an enrichment in ribonucleoprotein complex (P = 2.8E-4 and P = 2.6E-5) which contains the small nuclear ribonucleoprotein particles (snRNPs) forming the spliceosome (Supplementary Table 3). We also observed enrichment of the GO terms U2-type prespliceosome (P = 7.3E-5) and U1 snRNP (P = 1.3E-4) and could identify both U2A-factor (BAnTr.GP.5HG019656) and U1-70k (BAnTr.GP.1HG004326, BAnTr.GP.1HG006488) components plus additional genes involved in forming the small nuclear ribonucleoprotein particle (snRNP). The most differentially alternatively spliced genes can be observed in the comparison of meiocytes to anthers in leptotene–zygotene stage (Figure 5F). GO enrichment analysis of anther depleted genes highlighted the same GO terms as above (U2-type prespliceosome, P = 7E-5 and U1 snRNP, P = 1.1E-4). The GO term “regulation of DNA damage checkpoint” (P = 2.5E-4) was enriched in meiocytes (Supplementary Table 3). This GO term was assigned to replication Protein A 2A (RPA2A) and two other WD40 repeat containing proteins, one of which was annotated as DNA damage-binding protein 2 (DDB2). Surprisingly, HvRPA2A was not enriched in meiocytes at either stage while DDB2 and the other WD40 repeat containing protein showed significant (P = 1.5E-4 and P = 1.3 E-3) up-regulation in meiocytes compared to anthers at leptotene, although did not exceed our threshold log fold change value of + 1.5.
Genes Annotated With Roles in Post-transcriptional and Post-translational Modification Are a Major Component of Genes Enriched in Prophase I Meiocytes
GO enrichment analysis of DEGs points to a sizeable role for RNA and post-translational modifications in prophase I (Supplementary Table 3). This was most acutely highlighted in the comparison of meiocytes to anthers at leptotene–zygotene—representing the largest contrast in expression of protein-coding genes. Lys48-specific deubiquitinase activity (GO:1990380; P = 0.00083), histone lysine demethylation (GO:0070076; P = 0.00068), histone H3-K9 demethylation (GO:0033169; P = 0.00074), ligase activity (GO:0016874; P = 0.0002), and RNA modification (GO:0009451; P = 3.60E-10) were all significantly enriched in meiocytes. Enrichment of the RNA modification GO term reflected up-regulation of genes annotated as pentatricopeptide repeat (PPR) proteins, which regulate gene expression at the RNA level (Manna, 2015). Up-regulation of these genes also drives significant enrichment of endonuclease activity (GO:0004519; P = 0.00088) in this comparison. Of the two genes enriched in meiocytes at pachytene–diplotene versus leptotene–zygotene one was annotated as a PPR protein, the other was unannotated, and both depleted genes in this comparison were lncRNAs. Continued enrichment of these genes at pachytene–diplotene in meiocytes versus anthers at the same stage is also reflected in significant enrichment of mRNA binding (GO:0003729; P = 0.00033). Meiocyte enriched genes in this contrast group were enriched for GO terms including helicase activity (GO:0004386; P = 0.00088) and synapsis (GO:0007129; P = 0.00023), reflecting the formation of the SC in meiocytes at this stage. GO enrichment of ligase activity in DEGs in meiocytes at leptotene–zygotene compared to anthers at the same stage is driven by meiocyte enrichment in E3 ubiquitin ligases. In total, 890 ubiquitin or ubiquitin-like (SUMO, NEDD8) E3 ligases were annotated in the BAnTr transcriptomic dataset; in our dataset 133 of these were specific to anther and meiocyte tissues, 63 were unique to germinating embryos, and 589 were expressed in both anther and germinating embryo tissues. 166 of these were differentially expressed in a least one contrast group and 71 were enriched in meiocytes compared to anthers at the same stage (Supplementary Table 6). Notably, those enriched in meiocytes included many E3 ligases—or multi-subunit E3 ligase components—with confirmed roles in meiosis in other organisms. These included: a SKP1 ortholog (BAnTr.GP.5HG012508) and three F-box proteins (BAnTr.GP.2HG002834, BAnTr.GP.3HG014344, BAnTr.GP.3HG014348) which form part of the multi-subunit SKP Cullin F-box (SCF) E3 ligase complex; one RING-H2 component of the anaphase promoting complex (BAnTr.GP.1HG006974); and three seven-in absentia (SINA) E3 ligases (BAnTr.GP.2HG018292, BAnTr.GP.3HG000774, BAnTr.GP.3HG000850). HEI10, a highly conserved E3 ubiquitin/SUMO ligase which is known to be involved in DSB repair (Ziolkowski et al., 2017), was stably expressed in all tissues and stages, including in germinating embryo tissues. Although no NEDD8 E3 ligases are differentially expressed in any stage or tissue comparison, NEDD8-specific protease activity (GO:0019784; P = 0.00029) is enriched in pachytene–diplotene compared to leptotene–zygotene in anthers. Significant up-regulation of Lys48 specific deubiquitinase activity, in parallel to E3 ligase activity, might reflect delicate regulation of the ubiquitination cascade at this stage. However, of the four significantly up-regulated protein-coding genes assigned this GO by PANNZER2, only one (BAnTr.GP.7HG015226) is explicitly annotated as a deubiquitinating enzyme (DUB); two are annotated as B-box transcription factors and the other as an O-fucosyltransferase. Lys48 deubiquitination (GO:0071108; P = 3.20E-5) is also enriched in pachytene–diplotene compared to leptotene–zygotene in anthers reflecting the up-regulation of four protein-coding genes annotated as otubain-like DUBs, active on Lys48-linked polyubiquitin chains and NEDD8 (Edelmann et al., 2009). Histone lysine demethylation and histone H3-K9 demethylation enrichment in meiocytes compared to anthers at leptotene–zygotene is driven by up-regulated protein-coding genes orthologous to argonaute (AGO) proteins and demethylases KDM3 (BAnTr.GP.1HG003366), PKDM9 (BAnTr.GP.UnG004036), JMJ705 (BAnTr.GP.UnG004032 and BAnTr.GP.UnG004334), and JMJ706 (BAnTr.GP.1HG007084). Indeed, of 21 AGO orthologs, 10 are differentially expressed in at least one stage or tissue comparison (Supplementary Table 7 and Supplementary Figure 7). When the log2FC cut-offs for differential expression are disregarded 14 AGO orthologs are significant in at least one contrast group. HvAGO18 is significantly enriched in four out of six contrast groups—most significantly (P = 0.0000108; log2FC + 4.17) in anthers at leptotene–zygotene compared to pre-meiosis. HvAGO18 transcript increases significantly in abundance in anthers again at pachytene–diplotene (P = 0.005; log2FC + 1.5) and declines slightly, although not significantly at metaphase–tetrad. HvDCL5 was significantly up regulated in anthers at leptotene–zygotene stage compared to pre-meiosis (P = 0.0015). HvRDR6 was significantly enriched in meiocytes when compared to whole anther tissue at leptotene–zygotene (P = 0.0047).
Protein-coding genes which are up-regulated in anthers at successive stages and down-regulated in meiocyte versus anther contrast groups of the same stage are largely reflective of water transport, photosynthesis, cell signaling, amino acid transport, and other metabolic processes (Supplementary Table 3). However, down-regulated DEGs in anthers at leptotene–zygotene versus pre-meiosis are more informative. The switch to prophase I in anthers, in addition to the large-scale down-regulation of lncRNAs, coincides with massive down-regulation of protein-coding genes annotated with RNA-directed DNA polymerase activity (GO:0003964; P < 1E-30) all of which are expressed only in anther tissues and not in germinating embryos. Several of these genes were annotated as reverse transcriptase-like genes from copia-type centromeric retrotransposons RE1 and TNT 1 in Arabidopsis and N. tabacum respectively (Grandbastien et al., 1989; Yamada et al., 2014). These also formed a sizeable component of the down-regulated DEGs annotated as nucleic acid binding (GO:0003676; P < 1E-30) and zinc binding (GO:0008270; P < 1E-30) proteins; 89% of both groups were also anther specific. DNA binding transcription factors (GO:0003700; P = 3.20E-5) and transcriptional regulation related GO terms GO:0017053 and GO:0006355 (P = 5.50E-5 and 2.20E-6, respectively) were also enriched in down-regulated DEGs in this contrast group.
Many Meiotic Genes Are Expressed Before Meiosis Begins
We probed pre-meiotic and meiotic cells with antibodies against two SC (ASY1 and ZYP1) proteins and one recombination (DMC1) protein. The protein products of HvASY1 and HvZYP1 were detected during prophase I (Figures 6B,C) but not in pre-meiotic and anaphase I cells (Figures 6A,D). HvDMC1 detection on the other hand shows a diffuse, mostly cytoplasmic, signal in pre-meiotic cells (Figure 6E). During meiosis, HvDMC1 appears as discrete foci in leptotene–diplotene nuclei (Figures 6F,G) before becoming diffuse again at metaphase I–tetrad (Figure 6H). All three genes are transcribed at pre-meiosis (Figure 6I) and belong to the navy WGCNA module showing higher transcript abundance in prophase I meiocytes (Figure 2). We then analyzed the gene-level expression of 121 known meiotic genes which were distributed among 11 different modules. While all were expressed, only 29 were differentially expressed across the samples using our experimental thresholds (Figure 6I and Supplementary Figure 6). It is noticeable that the overall level of expression of some genes, like HvPCNA and members of the HSP90 family, remains high in all stages and tissues while others, like HvMET1b, HvMRE11B, HvXRCC3, and HvDDM1C, remained relatively low.
Figure 6. Expression of selected meiotic genes. Immuno-staining of meiotic nuclei at four developmental stages for HvZYP1 (magenta, A–D), and HvDMC1 (red, E–H) proteins. All samples were stained with anti-ASY1 antibody (green) and counterstained with DAPI (blue). (A,E) Pre-meiosis, (B,F), leptotene—zygotene, (C,G), pachytene—diplotene, (D,H), anaphase I—metaphase I. Scale bar 10 μm. (I) Heatmap expression profile of meiotic genes with a statistically significant log fold change in at least one tissue or stage comparison. The genes are ordered and grouped by WGCNA module on the vertical axis. Genes were extracted from the total dataset, transcript counts log transformed, and plotted using ggplot2 (Wickham, 2016) in R (script available at https://github.com/BioJNO/BAnTr). The samples (3 replicates each) are A.Pre, anther pre-meiosis; A.LepZyg, anther leptotene–zygotene; A.PacDip, anther pachytene–diplotene; A.MetTet, anther metaphase I–tetrad; M.LepZyg, meiocyte leptotene–zygotene; M.PacDip, meiocyte pachytene–diplotene. The prefixes A. and M. in the sample names depict anther and meiocyte samples, respectively.
Of the 28 meiotic genes, HvASY1, HvDMC1, HvZYP1, HvBUB3b, HvSDS, HvSWI, HvRPA2c, HvRPA1c.1, and HvHSP90.7 were significantly up-regulated from pre-meiosis to leptotene–zygotene. Further gene-level expression comparisons showed 19 putative meiotic genes (HvASY1, HvATM, HvATR, HvBRCA2b, HvBUB3b, HvMER3, HvPCH2, HvPHS1a, HvPRD3a, HvPRD3b, HvPS1, HvREC8.1, HvREC8.2, HvRPA2C, HvSCC2, HvSWI, HvTOP3alpha, HvZYP1, and HvZIP4) were significantly enriched in meiotic cells compared to anthers at the same stage (Supplementary Table 8). In comparison to our control dataset of germinating embryo samples, expression of only 6 of these meiotic genes are specific to anthers. However, while the rest are expressed in the EMB samples, 16 show lower gene expression in this tissue. The other nine genes show EMB gene expression similar to at least one anther developmental stage.
A Large Number of Transcription Factors Are Differentially Expressed During Anther and Meiocyte Development
Anther and meiocyte transcriptomes have been determined in several plant species but the mechanisms of gene regulation in these tissues are still poorly understood. The varied gene-level expression patterns of modules generated by WGCNA analysis in this dataset could be explained by steady-state gene expression regulation, including at the transcriptional level. Using the Plant Transcription Factor (TF) Database (Pérez-Rodríguez et al., 2010), we found a total of 1,353 annotated TFs of which 382 (28.2%)—belonging to several different TF families (Supplementary Table 8)—were differentially expressed in at least one stage or tissue comparison. 140 of these were expressed in anther tissues and not in germinating embryos. The largest families in the significant log2FC subset are bHLH (n = 29) and MYB (n = 29) (Figure 7A). The distribution of TFs families in different expression modules was determined with the red module containing the highest number (n = 154; Supplementary Table 9), a module showing down-regulation in meiocytes when compared to anthers at both stages of leptotene–zygotene and pachytene–diplotene (Supplementary Figure 3). This is reflected in the significant GO enrichment of transcriptional regulation in anthers at leptotene–zygotene versus pre-meiosis. 254 (66.5% of transcription factors significant in any contrast group) showed differential expression in comparisons between meiocytes and anthers at either pachytene–diplotene or leptotene–zygotene, with 170 and 84 TFs significantly different in one or both stages, respectively (Figure 7B and Supplementary Table 10). This strong transcriptional activation, accompanied by meiocyte enrichment in E3 ubiquitin ligases for likely protein degradation, might underly the vast transcriptome reprogramming at the meiosis onset.
Figure 7. Expression of transcription factor (TF) families in anthers and meiocytes. (A) Total number of TFs per family that are expressed in anthers and meiocytes. (B) Number of differentially expressed TFs determined by comparing their transcript levels in anthers at different stages. The comparisons are: A.LepZyg-A.Pre, anther leptotene–zygotene versus anther pre-meiosis; A.PacDip-A.LepZyg, anther pachytene–diplotene versus anther leptotene–zygotene; A.MetTet-A.PacDip, anther metaphase I–tetrad versus anther pachytene–diplotene. The prefixes A. and M. in the sample names depict anther and meiocyte samples, respectively.
Despite technical advances, and the importance of the reproductive phase in crop breeding, the regulation of the meiotic transcriptome in plants is still poorly understood (Zhou and Pawlowski, 2014). In this study, we sought to establish the most comprehensive transcriptome of barley anthers and meiocytes spanning pre-meiosis–tetrad stages. The barley cultivar Golden Promise (GP) was chosen as it is the reference cultivar for barley transformation which is relevant for future functional studies. In addition, we recently sequenced its genome (Schreiber et al., 2020), which will assist with interpretation of the transcriptome data, and developed a TILLING population using ethyl methanesulfonate (Schreiber et al., 2019). We used immuno-cytology to carefully stage all anther and meiocyte samples and analyzed a combination of Illumina and PacBio reads using the most up to date bioinformatic pipeline. This allowed us to identify novel transcripts and study transcriptome dynamics throughout meiosis.
lncRNAs Are Enriched in Pre-meiotic Anthers but Not in Meiocytes
The number of DEGs in staged anther tissues declines from its highest level between leptotene–zygotene and pre-meiosis with the largest down-regulation of genes of any comparison. This supports previous studies in Arabidopsis (Chen et al., 2010; Yang et al., 2011), maize (Dukowic-Schulze et al., 2014b; Nelms and Walbot, 2019) and sunflower (Flórez-Zapata et al., 2014, 2016) showing large scale transcriptional re-organization at meiosis entry. A strikingly large component of differentially expressed genes at this transition were lncRNAs—98% of which were expressed in anthers but not in germinating embryos. High expression of lncRNAs in plant reproductive tissues has been described in both sunflower (Flórez-Zapata et al., 2016) and in maize (Li et al., 2014). Flórez-Zapata et al. (2016) reported large scale differential expression of lncRNAs in sunflower meiocytes when contrasted with somatic tissues; and between meiocytes in comparison of genotypes exhibiting significantly different recombination rates. Our results indicate that anther specific lncRNAs may play an important role in prophase I entry in barley. However, our results also indicate that expression of many of these lncRNAs is down-regulated during the transition to prophase I and in meiocytes when compared to whole anther tissues—indicating that their most important functional roles may be prior to meiosis or as pre-cursors to meiotic small RNAs (sRNAs).
Rapid evolution of lncRNAs is well documented; however, in most cases function is conserved across large evolutionary distances (Ulitsky, 2016). This is reflected in the low proportion of BAnTr lncRNAs which align to published lncRNA database sequences. LncRNAs are emerging as key plant development regulators (Ariel et al., 2015; Shafiq et al., 2016) and their role in animal and plant sexual reproduction is now well established (reviewed in Golicz et al., 2018), including transcriptional gene silencing (Wierzbicki et al., 2008), nucleosome positioning (Garcia et al., 2010), chromosome looping (Ariel et al., 2014), and pairing of homologous chromosomes (Ding et al., 2012). LncRNAs are also precursors of phasiRNAs—cleaved double-stranded RNAs either 21 or 24 nt in length. Both 21 and 24 nt phasiRNAs are associated with male reproductive development in angiosperms, the former predominantly and the latter exclusively. We have identified 568 and 553 PHAS transcripts with predicted miR2118 and miR2275 recognition sites, respectively, and only 53 lncRNAs contain recognition sites for both miRNAs. These numbers are larger than previously reported for maize anthers using small RNA sequencing in maize (Zhai et al., 2015). RDR6 is involved in processing both classes, converting single-stranded precursor RNA into double-stranded RNA (dsRNA) following AGO-catalyzed cleavage (Song et al., 2012b). DCL4 and DCL5 then cuts this dsRNA in regular intervals producing 21 and 24 nt phasiRNAs, respectively (Song et al., 2012a; Zhang M. et al., 2020). Both RDR6 and DCL5 are present in co-expression networks enriched in meiocytes in the BAnTr dataset. RDR6 showed significant enrichment at leptotene–zygotene in meiocytes compared to anthers and DCL5 was up-regulated in anthers at leptotene–zygotene compared to pre-meiosis (Supplementary Figure 7 and Supplementary Table 7). DCL5 has been shown to generate 24-nt phasiRNAs, which in turn silence the 24-PHAS loci in cis by CHH-type DNA methylation during meiotic anther development in maize (Zhang M. et al., 2020). This may reflect the importance of 24-nt phasiRNAs in prophase I. Several AGO genes also showed significant differential expression. HvAGO18, orthologous to AGO18 which is a suggested 24-nt phasiRNAs binding partner (Zhai et al., 2015; Fei et al., 2016), is significantly enriched in four out of six contrast groups—most significantly (P = 0.0000108; log2FC + 4.17) in anthers at leptotene–zygotene compared to pre-meiosis. HvAGO18 transcript increases significantly in abundance in anthers again at pachytene–diplotene (P = 0.005; log2FC + 1.5) and declines slightly, although not significantly at metaphase–tetrad. HvAGO18—and several other AGO orthologs—are enriched in meiocytes compared to anthers in both staged prophase I sample groups. phasiRNAs have been shown to be enriched in anthers at different meiotic stages in maize (Zhai et al., 2015) and rice (Komiya et al., 2014; Fei et al., 2016) and shown to be present in most angiosperms (Xia et al., 2019). Deep sequencing technologies have enabled global analysis of meiotic sRNA in sunflower showing sequence similarity to 40% of lncRNAs (Flórez-Zapata et al., 2016). In Arabidopsis, Huang et al. (2019) revealed different meiotic and mitotic sRNA landscapes and found that meiocyte-specific sRNAs (ms-sRNAs) are significantly enriched in genic regions contrary to somatic small interfering RNAs that are enriched in intergenic regions. A high proportion of these ms-sRNAs (69%) were found to be DSB-dependent (Huang et al., 2019). A detailed analysis of genome-wide distribution of PHAS loci and their expression throughout anther development in wheat and barley has uncovered a new 24 nt pre-meiotic class (Bélanger et al., 2020). Further analyses will be needed to fully elucidate the relationship of lncRNAs described in our study and sRNAs.
Methylation, Demethylation, and AGO Genes
sRNAs biogenesis and function require several players including RNA polymerases (Pol), Dicer-like (DCL) proteins, double-stranded RNA-binding (DRB) proteins, RNA-directed RNA polymerases (RDRs) and AGO proteins (Borges and Martienssen, 2015; Yu et al., 2018). We checked the expression dynamics of sRNA biogenesis and AGO genes, for which protein products were detected in our recent anther proteomic study (Lewandowska et al., 2019), and found the expression of 10 out of 21 AGO genes changed significantly during barley anther and meiocyte development (Supplementary Figure 7). This gene expression profile highlights the importance of sRNA biogenesis and their guiding AGO proteins during meiosis. A new study characterized the components of phasiRNA biogenesis pathway in wheat and barley and found that most AGO genes were expressed in anthers of both species (Bélanger et al., 2020). Expression profiles of different sRNA biogenesis compounds, including AGOs we found in this study, were similar to those reported by Bélanger et al. (2020). These authors correlated the expression pattern of different gene family members with those of pre-meiotic 21 and 24 nt and meiotic 24 nt phasiRNA classes delineating their potential role in different anther developmental stages. However, further analysis of these pathways and their target loci is needed to fully describe their role in meiotic development. One such function has been demonstrated in a recent study by Walker et al. (2018) showing that de novo RNA-directed DNA methylation (RdDM) induces cell-lineage-specific epigenetic marks regulating meiotic gene expression in Arabidopsis.
mRNA methylation is emerging as another important level of meiotic gene regulation. The N6-methyladenosine (m6A) modification was shown to be associated with mRNA translatability in yeast (Bodi et al., 2015) and Xenopus (Qi et al., 2016). More recently, Bushkin et al. (2019) showed that m6A methylation at the 3′ UTR of the mRNA encoding Rme1p, a transcriptional repressor of meiosis, results in its degradation allowing meiosis initiation in yeast. Current epitranscriptome methods necessitate a large amount of mRNA making it difficult to implement in plant meiocytes. Instead, we determined transcript levels of known plant m6A pathway genes (Yue et al., 2019), for which protein products were detected in our anther proteomic study (Lewandowska et al., 2019). The expression of these genes was relatively low and stable throughout anther and meiocyte development. However, the translation status and potential function in these tissues remains to be elucidated. It will be interesting to compare the expression profiles of these enzymes in different plant species. Indeed, Liu J. et al. (2020) profiled m6A and N6,2′-O-dimethyladenosine (m6Am) across human and mouse tissues and suggested that the difference in these epitranscriptomic marks is greater between species than tissue types.
Significant Post-transcriptional and Post-translational Regulation During Prophase I
Comparison of meiocytes and anthers at leptotene–zygotene presented the greatest contrast in protein-coding gene expression. GO enrichment analysis of DEGs in this contrast group suggests heightened importance of ubiquitination, RNA modification, and histone demethylation in prophase I regulation.
Post-transcriptional RNA Modification Through PPR Proteins
Significant enrichment of the RNA modification GO term in meiocytes compared to anthers at leptotene–zygotene reflected several up-regulated DEGs annotated as PPR proteins. PPR proteins are found in all eukaryotes but are remarkably abundant in land plants, indicative of massive expansion during land plant evolution (Barkan and Small, 2014; Gutmann et al., 2020). Through sequence-specific binding, PPR proteins can mediate RNA folding, splicing, degradation, cleavage, and editing (Barkan and Small, 2014). PPR proteins are associated with organelles, deleterious mutations in which lead to defects in photosynthesis or oxidative phosphorylation (Barkan and Small, 2014). Dukowic-Schulze et al. (2014b) reported up-regulation of 24 mitochondrial genes in maize meiocytes compared to whole anther tissues and seedlings, arguing that this indicated high energy demand concomitant with chromosome movement in early prophase I. Further, Yuan et al. (2018) and Nelms and Walbot (2019) also reported significant GO enrichment of various mitochondrial processes in developing germinal cells. Up-regulation of PPR proteins in early barley prophase I may support this, indicating a general increase in mitochondrial activity in step with a detectable spike in PPR protein mediated RNA regulation.
Ubiquitin Ligases and Deubiquitinating Enzymes
E3 ubiquitin ligases are an important regulatory component during meiosis in many organisms (Okamato et al., 2012; Mohammad et al., 2018). E3 ligases interact with target proteins to facilitate, directly or indirectly, modification of the substrate with ubiquitin—conferring substrate specificity to the ubiquitination cascade (Iconomou and Saunders, 2016). Several E3 ligases, and multi-subunit E3 ligase components, have been described with significant roles in meiosis in Arabidopsis (Yang et al., 1999, 2006; Wang and Yang, 2006), rice (He et al., 2016; Zhang F. et al., 2017), and wheat (Li et al., 2006; Hong et al., 2013). Further, enrichment of E3 ligases in pre-meiotic pollen mother cells has been reported in both rice (Tang et al., 2010) and maize (Yuan et al., 2018). Our findings suggest that ubiquitination through E3 ligase activity is of continued importance throughout prophase I. E3 ligase genes up-regulated in meiocytes include an ASK1 ortholog, an Arabidopsis S-phase kinase-associated protein 1 (SKP1) which interacts with Cullin and F-box proteins to form the SKP-Cullin-F-box (SCF) E3 ligase complex (Yang et al., 1999). SKP1-like proteins are functionally conserved in plants, with the wheat ASK1 equivalent, TSK1, able to partially rescue fertility in Arabidopsis ask1 mutants (Li et al., 2006). ASK1 was identified as a negative regulator of male recombination, essential for the proper release of chromatin from the nuclear membrane (Wang and Yang, 2006; Yang et al., 2006). Several F-box proteins, putatively interacting in the SCF E3 ubiquitin ligase complex, also displayed differential expression throughout prophase I. In rice, mutations in two F-box proteins interacting with the ASK1 equivalent OSK1—MOF and ZYGO1—have also been shown to result in male sterility (He et al., 2016; Zhang F. et al., 2017). Several genes annotated as SINA E3 ligases also displayed differential expression during prophase I. A recent study in Drosophila melanogaster females identified a SINA E3 ligase which regulated both assembly and disassembly of the SC, preventing aberrant polymerization and polycomplex formation of SC components (Hughes et al., 2019). Although no such function of SINA E3 ligase is so far reported in plants, the E3 ligases highlighted in this work are attractive candidates to investigate a similar role for these proteins in regulating SC formation in barley. DUBs, particularly those targeting lys48 ubiquitin chains, are also differentially expressed throughout prophase I. DUBs facilitate removal of ubiquitin from proteins; in this way they reverse ubiquitination mediated by E3 ligases. Lys48 is one of seven lysine residues through which ubiquitin can form covalent C-terminally linked chains (Kulathu and Komander, 2012; López-Mosqueda and Dikic, 2014). Although targeting for proteasomal degradation is the canonical function of ubiquitination, first described by Ciehanover et al. (1978), many other effects of this modification have been observed such as recruitment of binding partners (Huang and D’Andrea, 2010), activation (Xu et al., 2009), or nuclear uptake (Plafker et al., 2004). The diversity of substrate fates is derived from diversity in chain topology (Komander and Rape, 2012). Simultaneous up-regulation of E3 ubiquitin ligase and DUB expression may represent both very tight control of substrate modification, where ligase and DUB build and remove the same chain topology, respectively, and specific reduction in lys48 ubiquitination while alternative chain linkage types are enriched, where they do not.
Expression of Meiotic Genes Does Not Always Reflect the Timing of Their Functional Roles in Meiocytes
Using phylogenomic analysis, we compiled a large inventory of barley orthologs of known meiotic genes. We determined their relative expression levels and found them to be distributed in 11 different expression modules. In total, 24% of selected meiotic genes were differentially expressed during anther and meiocyte development.
Like previous studies, our transcriptomic data provided evidence that many known meiotic genes are transcribed prior to prophase I. However, the timing of their translation is still not fully resolved due to limited availability of antibodies (this work), issues related to sample staging, and proteomic resolution (Zhang et al., 2014; Yuan et al., 2018). Genome-wide ribosome profiling could be achieved by using transgenic tools like RiboTag/RNA-seq (Lesiak et al., 2015). Such analysis was performed in mice revealing alteration of maternal mRNA translation in oocytes at meiotic re-entry (Luong et al., 2020).
Expression of Transcription Factors Changes Throughout Anther Stages
To decipher the mechanisms of meiotic transcriptome regulation, we analyzed our transcriptomic data for transcriptional regulation. The involvement of many transcription factors (TFs) in transcriptome reprogramming during pre-meiotic anther differentiation has been shown in maize (Zhang et al., 2014). This study was recently expanded by analyzing TFs expression at tissue level using microdissection (Yuan et al., 2018). However, both studies used microarrays, a closed technology, that may underestimate the number of TFs. We found a total of 1,353 annotated TFs (63 families) of which 28.2% were differentially expressed in at least one stage or tissue comparison. In agreement with previous maize studies, our data also reveal major TFs expression changes occurring at and beyond mitotic–meiotic transition in whole anthers. However, no such changes were observed in meiocytes. Beside these global analyses, functional studies of TFs are still lacking to fully understand meiotic gene regulation. More recently, TBP-ASSOCIATED FACTOR 4b (TAF4b), encoding a subunit of the RNA polymerase II general transcription factor TFIID, was shown to be enriched in meiocytes and to control transcription of genes involved in the meiotic cell cycle and recombination in Arabidopsis (Lawrence et al., 2019). In our study, a single barley TAF4b homolog gene (BAnTr.GP.6HG009004) was found in the green module showing increasing expression (P < 0.01) throughout anther and meiocyte stages, though below our threshold of log2FC > 1 in all comparisons.
Although alternative splicing (AS) is known to increase the genome coding capacity, its activity during meiosis has been studied as far as we are aware only in mouse testis (Schmid et al., 2013) and yeast (Kuang et al., 2017). Intron retention events were found to be enriched in both mouse (Naro et al., 2017) and yeast meiocytes (Kuang et al., 2017). They were also enriched during Arabidopsis flower development (Wang et al., 2014). Here, GO enrichment analysis showed enrichment of terms that are connected to alternative splicing. This contained genes involved in forming the spliceosome or spliceosome-associated non-snRNP proteins. Contrary to our expectations we only identified two meiotic genes as significantly differentially alternatively spliced. In future we would like to extend our potential meiotic gene list as the 121 genes are only a small subset and need to be extended to potential orthologous genes identified in studies from mice, yeast, and Arabidopsis. We also might be able to extend the list and add barley specific meiotic genes through this dataset.
Conclusion and Perspectives
Our analysis of barley anther and meiocyte transcriptomes provides evidence of a multi-faceted regulatory network orchestrating meiotic gene regulation. Sequencing of sRNA and DNA methylation, as done in maize by Dukowic-Schulze et al. (2016), in anther and meiocyte at pre-meiosis to leptotene–zygotene stages with the strongest transcriptomic changes will enhance the full picture of these interconnecting gene regulation mechanisms. Newly developed technologies, like single-cell RNA-seq (Denyer et al., 2019) and spatial transcriptomics (Giacomello et al., 2017), could increase the resolution of anther and meiocyte transcriptomes even at sub-cellular level. Nascent transcript sequencing by NET-seq (Churchman and Weissman, 2011), GRO-seq (Hetzel et al., 2016) or Neu-seq (Szabo et al., 2020) and ribosome profiling by RiboTag-seq (Lesiak et al., 2015), for example, could be used to directly evaluate transcriptional and translational activities throughout anther and meiocyte development.
The combination of immuno-cytological staging and RNA-seq data has allowed us to build a robust time-resolved barley anther and meiocyte transcriptomic dataset. We detected large scale down-regulation of lncRNAs at meiosis entry in anthers and enrichment of DAS at early prophase I. In addition to changes in expression, we also revealed the diversity of transcription factors accompanied by several other post-transcriptional and post-translational regulatory networks. Our data will be informative for research in anther and meiocyte development in other plant species. In addition, it has already contributed to building a first barley reference transcript dataset (BaRTv1) (Rapazote-Flores et al., 2019), to constructing a reference database for proteomic studies of staged anthers (D. Lewandowska et al., unpublished) and for annotation of barley reference genome and pangenome sequences (Monat et al., 2019; Jayakodi et al., unpublished).
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA558196; https://www.ncbi.nlm.nih.gov/, PRJNA593943; https://figshare.com/, https://doi.org/10.6084/m9.figshare.11974182, https://figshare.com/, http://doi.org/10.6084/m9.figshare.12136773, https://figshare.com/, http://doi.org/10.6084/m9.figshare.12136773, https://figshare.com/, http://doi.org/10.6084/m9.figshare.12136773.
AB, DL, IC, MS, LR, and RW: experimental design. NM: plant growth. AB, DL, IC, MA, MM, MS, and NM: collection of anthers. IC: collection of meiocytes and immuno-cytology. AB, PH, and JM: RNA extraction, library prep, and QC. MS: transcriptome assembly. AB: lead manuscript preparation. AB, DL, JO, IC, MS, and RW: data analysis and manuscript preparation. All authors contributed to the article and approved the submitted version.
This work was funded by the European Research Council (ERC Shuffle Project ID: 669182) and supported by the Scottish Government’s Rural and Environment Science and Analytical Services Division work program Theme 2 WP2.1 RD1 and RD2.
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.
We would like to thank Drs. Paulo Rapazote-Flores and Runxuan Zhang for their Bioinformatic help and suggestions and Mr. Philip Smith for proofreading. An earlier version of this manuscript was released as a pre-print at bioRxiv (Barakate et al., 2020).
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2020.619404/full#supplementary-material
Supplementary Figure 1 | 3D MDS plot of sample transcript expression.
Supplementary Figure 2 | 3D RNAseq overview and long non-coding RNA (lncRNA) categories.
Supplementary Figure 3 | WGNCA co-expression modules.
Supplementary Figure 4 | Eigengene network analysis.
Supplementary Figure 5 | WGCNA GO enrichment analysis.
Supplementary Figure 6 | Heat map of meiotic gene expression.
Supplementary Figure 7 | Heat map sRNA gene silencing protein expression.
Supplementary Table 1 | Raw read counts.
Supplementary Table 2 | Number of genes per WGCNA modules.
Supplementary Table 3 | GO enrichment analysis results.
Supplementary Table 4 | Differential expression analysis.
Supplementary Table 5 | lncRNA database BLAST hits by classification.
Supplementary Table 6 | Differentially Expressed E3 Ubiquitin ligases.
Supplementary Table 7 | Expression of sRNA gene silencing orthologs.
Supplementary Table 8 | Differential expression analysis results for Meiotic genes.
Supplementary Table 9 | Distribution of transcription factor families in different WGCNA modules.
Supplementary Table 10 | Transcription factor families in different comparisons.
Supplementary Method 1 | Immunocytology.
Supplementary Method 2 | Isolation of meiocytes.
Supplementary Method 3 | RNA extraction and sequencing.
Supplementary Method 4 | Exact parameters of Illumina read mapping.
Supplementary Video 1 | Isolated fresh meiocytes bag at Leptotene/zygotene.
Supplementary Video 2 | Isolated fresh meiocytes bag at pachytene/diplotene.
Supplementary File 1 | Assembled BAnTr Transcriptome fasta file. AntherTranscriptomeBAnTr.fasta at http://doi.org/10.6084/m9.figshare.12136773.
Supplementary File 2 | Padded BAnTr Transcriptome fasta file. AntherTranscriptomeBAnTrPadded.fasta at http://doi.org/10.6084/m9.figshare.12136773.
Supplementary File 3 | BAnTr Proteome fasta file. AntherProteomeBAnTr.fasta at http://doi.org/10.6084/m9.figshare.12136773.
Alamancos, G. P., Pagès, A., Trincado, J. L., Bellora, N., and Eyras, E. (2015). Leveraging transcript quantification for fast computation of alternative splicing profiles. RNA 21, 1521–1531. doi: 10.1261/rna.051557.115
Alexa, A., Rahnenführer, J., and Lengauer, T. (2006). Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics 22, 1600–1607. doi: 10.1093/bioinformatics/btl140
Ariel, F., Jegu, T., Latrasse, D., Romero-Barrios, N., Christ, A., Benhamed, M., et al. (2014). Noncoding transcription by alternative RNA polymerases dynamically regulates an auxin-driven chromatin loop. Mol. Cell. 55, 383–396. doi: 10.1016/j.molcel.2014.06.011
Arrieta, M., Colas, I., Macaulay, M., Waugh, R., and Ramsay, L. (2020). “A modular tray growth system for barley,” in Plant Meiosis, Methods in Molecular Biology, Vol. 2061, eds M. Pradillo and S. Heckmann (New York, NY: Humana Press), 367–379. doi: 10.1007/978-1-4939-9818-0_26
Barakate, A., Higgins, J. D., Vivera, S., Stephens, J., Perry, R. M., Ramsay, L., et al. (2014). The synaptonemal complex protein ZYP1 is required for imposition of meiotic crossovers in barley. Plant Cell 26, 729–740. doi: 10.1105/tpc.113.121269
Barakate, A., Orr, J., Schreiber, M., Colas, I., Lewandowska, D., McCallum, N., et al. (2020). Time-resolved transcriptome of barley anthers and meiocytes reveals robust and largely stable gene expression changes at meiosis entry. bioRxiv [Preprint], doi: 10.1101/2020.04.20.051425
Bélanger, S., Pokhrel, S., Czymmek, K., and Meyers, B. C. (2020). Pre-meiotic, 24-nt reproductive phasiRNAs are abundant in anthers of wheat and barley but not rice and maize. bioRxiv [Preprint], doi: 10.1101/2020.06.18.160440
Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. Roy. Stat. Soc. B Stat. Methodol. 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x
Bodi, Z., Bottley, A., Archer, N., May, S. T., and Fray, R. G. (2015). Yeast m6A methylated mRNAs are enriched on translating ribosomes during meiosis, and under rapamycin treatment. PLoS One 10:e0132090. doi: 10.1371/journal.pone.0132090
Braynen, J., Yang, Y., Wei, F., Cao, G., Shi, G., Tian, B., et al. (2017). Transcriptome analysis of floral buds deciphered an irregular course of meiosis in polyploid Brassica rapa. Front. Plant Sci. 8:768. doi: 10.3389/fpls.2017.00768
Bushkin, G. G., Pincus, D., Morgan, J. T., Richardson, K., Lewis, C., Chan, S. H., et al. (2019). m6A modification of a 3’ UTR site reduces RME1 mRNA levels to promote meiosis. Nat. Commun. 10:3414. doi: 10.1038/s41467-019-11232-7
Chen, C., Farmer, A. D., Langley, R. J., Mudge, J., Crow, J. A., May, G. D., et al. (2010). Meiosis-specific gene discovery in plants: RNA-Seq applied to isolated Arabidopsis male meiocytes. BMC Plant Biol. 10:280. doi: 10.1186/1471-2229-10-280
Chen, J., Su, P., Chen, P., Li, Q., Yuan, X., and Liu, Z. (2018). Insights into the cotton anther development through association analysis of transcriptomic and small RNA sequencing. BMC Plant Biol. 18:154. doi: 10.1186/s12870-018-1376-4
Ciehanover, A., Hod, Y., and Hershko, A. (1978). A heat-stable polypeptide component of an ATP-dependent proteolytic system from reticulocytes. Biochem. Biophys. Res. Commun. 81, 1100–1105. doi: 10.1016/0006-291X(78)91249-4
Colas, I., Darrier, B., Arrieta, M., Mittmann, S. U., Ramsay, L., Sourdille, P., et al. (2017). Observation of extensive chromosome axis remodeling during the “diffuse-phase” of meiosis in large genome cereals. Front. Plant Sci. 8:1235. doi: 10.3389/fpls.2017.01235
Denyer, T., Ma, X., Klesen, S., Scacchi, E., Nieselt, K., and Timmermans, M. C. P. (2019). Spatiotemporal developmental trajectories in the Arabidopsis root revealed using high-throughput single-cell RNA sequencing. Dev. Cell 48, 840–852. doi: 10.1016/j.devcel.2019.02.022
Ding, D. Q., Okamasa, K., Yamane, M., Tsutsumi, C., Haraguchi, T., Yamamoto, M., et al. (2012). Meiosis-specific noncoding RNA mediates robust pairing of homologous chromosomes in meiosis. Science 336, 732–736. doi: 10.1126/science.1219518
Dukowic-Schulze, S., Harris, A., Li, J., Sundararajan, A., Mudge, J., Retzel, E. F., et al. (2014a). Comparative transcriptomics of early meiosis in Arabidopsis and maize. J. Genet. Genom. 41, 139–152. doi: 10.1016/j.jgg.2013.11.007
Dukowic-Schulze, S., Sundararajan, A., Mudge, J., Ramaraj, T., Farmer, A. D., Wang, M., et al. (2014b). The transcriptome landscape of early maize meiosis. Plant Biol. 14:118. doi: 10.1186/1471-2229-14-118
Dukowic-Schulze, S., Sundararajan, A., Ramaraj, T., Kianian, S., Pawlowski, W. P., Mudge, J., et al. (2016). Novel Meiotic miRNAs and indications for a role of PhasiRNAs in meiosis. Front. Plant Sci. 7:762. doi: 10.3389/fpls.2016.00762
Edelmann, M. J., Iphöfer, A., Akutsu, M., Altun, M., di Gleria, K., Kramer, H. B., et al. (2009). Structural basis and specificity of human otubain 1-mediated deubiquitination. Biochem. J. 418, 379–390. doi: 10.1042/BJ20081318
Fei, Q., Yang, L., Liang, W., Zhang, D., and Meyers, B. C. (2016). Dynamic changes of small RNAs in rice spikelet development reveal specialized reproductive phasiRNA pathways. J. Exp. Bot. 67, 6037–6049. doi: 10.1093/jxb/erw361
Flórez-Zapata, N. M. V., Reyes-Valdés, M. H., Hernandez-Godínez, F., and Martínez, O. (2014). Transcriptomic landscape of prophase I sunflower male meiocytes. Front. Plant Sci. 5:277. doi: 10.3389/fpls.2014.00277
Flórez-Zapata, N. M. V., Reyes-Valdés, M. H., and Martínez, O. (2016). Long non-coding RNAs are major contributors to transcriptome changes in sunflower meiocytes with different recombination rates. BMC Genom. 17:490. doi: 10.1186/s12864-016-2776-1
Garcia, J. F., Dumesic, P. A., Hartley, P. D., El-Samad, H., and Madhani, H. D. (2010). Combinatorial, site-specific requirement for heterochromatic silencing factors in the elimination of nucleosome-free regions. Genes Dev. 24, 1758–1771. doi: 10.1101/gad.1946410
Giacomello, S., Salmén, F., Terebieniec, B. K., Vickovic, S., Navarro, J. F., Alexeyenko, A., et al. (2017). Spatially resolved transcriptome profiling in model plant species. Nat. Plants 3:17061. doi: 10.1038/nplants.2017.61
Gómez, J. F., and Wilson, Z. A. (2012). Non-destructive staging of barley reproductive development for molecular analysis based upon external morphology. J. Exp. Bot. 63, 4085–4094. doi: 10.1093/jxb/ers092
Guo, W., Tzioutziou, N., Stephen, G., Milne, I., Calixto, C., Waugh, R., et al. (2019). 3D RNA-seq - a powerful and flexible tool for rapid and accurate differential expression and alternative splicing analysis of RNA-seq data for biologists. bioRxiv [Preprint], doi: 10.1101/656686
Gutmann, B., Royan, S., Schallenberg-Rüdinger, M., Lenz, H., Castleden, I. R., McDowell, R., et al. (2020). The expansion and diversification of pentatricopeptide repeat RNA-editing factors in plants. Mol. Plant 13, 215–230. doi: 10.1016/j.molp.2019.11.002
Haas, B. J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P. D., Bowden, J., et al. (2013). De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8, 1494–1512. doi: 10.1038/nprot.2013.084
Hackl, T., Hedrich, R., Schultz, J., and Förster, F. (2014). proovread: large-scale high-accuracy PacBio correction through iterative short read consensus. Bioinformatics 30, 3004–3011. doi: 10.1093/bioinformatics/btu392
He, Y., Wang, C., Higgins, J. D., Yu, J., Zong, J., Lu, P., et al. (2016). MEIOTIC F-BOX is essential for male meiotic DNA double-strand break repair in rice. Plant Cell 28, 1879–1893. doi: 10.1105/tpc.16.00108
Hetzel, J., Duttke, S. H., Benner, C., and Chory, J. (2016). Nascent RNA sequencing reveals distinct features in plant transcription. Proc. Natl. Acad. Sci. U.S.A. 113, 12316–12321. doi: 10.1073/pnas.1603217113
Higgins, J. D., Perry, R. M., Barakate, A., Ramsay, L., Waugh, R., Halpin, C., et al. (2012). Spatiotemporal asymmetry of the meiotic program underlies the predominantly distal distribution of meiotic crossovers in barley. Plant Cell 24, 4096–4109. doi: 10.1105/tpc.112.102483
Hobo, T., Suwabe, K., Aya, K., Suzuki, S., Yano, K., Ishimizu, T., et al. (2008). Various spatiotemporal expression profiles of anther-expressed genes in rice. Plant Cell Physiol. 49, 1417–1428. doi: 10.1093/pcp/pcn128
Hong, M. J., Kim, D. Y., and Seo, Y. W. (2013). SKP1-like-related genes interact with various F-box proteins and may form SCF complexes with Cullin-F-box proteins in wheat. Mol. Biol. Rep. 40, 969–981. doi: 10.1007/s11033-012-2139-1
Hord, C. L., and Ma, H. (2007). “Genetic control of anther cell division and differentiation,” in Cell Division Control in Plants, Plant Cell Monographs 9, eds D. P. S. Verma and Z. Hong (Berlin: Springer-Verlag), 361–376. doi: 10.1007/7089_2007_136
Huang, J., Wang, C., Wang, H., Lu, P., Zheng, B., Ma, H., et al. (2019). Meiocyte-specific and AtSPO11-1-dependent small RNAs and their association with meiotic gene expression and recombination. Plant Cell 31, 444–464. doi: 10.1105/tpc.18.00511
Huerta-Cepas, J., Forslund, K., Coelho, L. P., Szklarczyk, D., Jensen, L. J., von Mering, C., et al. (2017). Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol. Biol. Evol. 34, 2115–2122. doi: 10.1093/molbev/msx148
Huerta-Cepas, J., Szklarczyk, D., Heller, D., Hernández-Plaza, A., Forslund, S. K., Cook, H., et al. (2019). eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 47, D309–D314. doi: 10.1093/nar/gky1085
Hughes, S. E., Hemenway, E., Guo, F., Yi, K., Yu, Z., and Hawley, R. S. (2019). The E3 ubiquitin ligase Sina regulates the assembly and disassembly of the synaptonemal complex in Drosophila females. PLoS Genet. 15:e1008161. doi: 10.1371/journal.pgen.1008161
Jalali, S., Jayaraj, G. G., and Scaria, V. (2012). Integrative transcriptome analysis suggest processing of a subset of long non-coding RNAs to small RNAs. Biol. Direct. 7:25. doi: 10.1186/1745-6150-7-25
Jan, S. Z., Vormer, T. L., Jongejan, A., Röling, M. D., Silber, S. J., de Rooij, D. G., et al. (2017). Unraveling transcriptome dynamics in human spermatogenesis. Development 144, 3659–3673. doi: 10.1242/dev.152413
Johnson, C., Kasprzewska, A., Tennessen, K., Fernandes, J., Nan, G. L., Walbot, V., et al. (2009). Clusters and superclusters of phased small RNAs in the developing inflorescence of rice. Genome Res. 19, 1429–1440. doi: 10.1101/gr.089854.108
Jones, P., Binns, D., Chang, H. Y., Fraser, M., Li, W., McAnulla, C., et al. (2014). InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240. doi: 10.1093/bioinformatics/btu031
Komiya, R., Ohyanagi, H., Niihama, M., Watanabe, T., Nakano, M., Kurata, N., et al. (2014). Rice germline-specific Argonaute MEL1 protein binds to phasiRNAs generated from more than 700 lincRNAs. Plant J. 78, 385–397. doi: 10.1111/tpj.12483
Kubo, T., Fujita, M., Takahashi, H., Nakazono, M., Tsutsumi, N., and Kurata, N. (2013). Transcriptome analysis of developing ovules in rice isolated by laser microdissection. Plant Cell Physiol. 54, 750–765. doi: 10.1093/pcp/pct029
Lawrence, E. J., Gao, H., Tock, A. J., Lambing, C., Blackwell, A. R., Feng, X., et al. (2019). Natural variation in TBP-ASSOCIATED FACTOR 4b controls meiotic crossover and germline transcription in Arabidopsis. Curr. Biol. 29, 2676–2686. doi: 10.1016/j.cub.2019.06.084
Lesiak, A. J., Brodsky, M., and Neumaier, J. F. (2015). RiboTag is a flexible tool for measuring the translational state of targeted cells in heterogeneous cell cultures. Biotechniques 58, 308–317. doi: 10.2144/000114299
Lewandowska, D., Zhang, R., Colas, I., Uzrek, N., and Waugh, R. (2019). Application of a sensitive and reproducible label-free proteomic approach to explore the proteome of individual meiotic phase barley anthers. Front. Plant Sci. 10:393. doi: 10.3389/fpls.2019.00393
Li, C., Liang, Y., Chen, C., Li, J., Xu, Y., Xu, Z., et al. (2006). Cloning and expression analysis of TSK1, a wheat SKP1 homologue, and functional comparison with Arabidopsis ASK1 in male meiosis and auxin signalling. Funct. Plant Biol. 33, 381–390. doi: 10.1071/FP06026
Li, L., Eichten, S. R., Shimizu, R., Petsch, K., Yeh, C.-T., Wu, W., et al. (2014). Genome-wide discovery and characterization of maize long non-coding RNAs. Genome Biol. 15:R40. doi: 10.1186/gb-2014-15-2-r40
Li, X., Shahid, M. Q., Wen, M., Chen, S., Yu, H., Jiao, Y., et al. (2020). Global identification and analysis revealed differentially expressed lncRNAs associated with meiosis and low fertility in autotetraploid rice. BMC Plant Biol. 20:82. doi: 10.1186/s12870-020-2290-0
Liu, C., Shen, Y., Qin, B., Wen, H., Cheng, J., Mao, F., et al. (2020). Oryza sativa RNA-dependent RNA polymerase 6 contributes to double-strand break formation in meiosis. Plant Cell 32, 3273–3289. doi: 10.1105/tpc.20.00213
Liu, J., Li, K., Cai, J., Zhang, M., Zhang, X., Xiong, X., et al. (2020). Landscape and regulation of m6A and m6Am methylome across human and mouse tissues. Mol. Cell 77, 426–440. doi: 10.1016/j.molcel.2019.09.032
Luong, X. G., Daldello, E. M., Rajkovic, G., Yang, C. R., and Conti, M. (2020). Genome-wide analysis reveals a switch in the translational program upon oocyte meiotic resumption. Nucleic Acids Res. 48, 3257–3276. doi: 10.1093/nar/gkaa010
Ma, X., Shao, C., Jin, Y., Wang, H., and Meng, Y. (2014). Long non-coding RNAs: a novel endogenous source for the generation of Dicer-like 1-dependent small RNAs in Arabidopsis thaliana. RNA Biol. 11, 373–390. doi: 10.4161/rna.28725
Mapleson, D., Venturini, L., Kaithakottil, G., and Swarbreck, D. (2018). Efficient and accurate detection of splice junctions from RNAseq with portcullis. Gigascience 7:giy131. doi: 10.1093/gigascience/giy131
McCarthy, D. J., Chen, Y., and Smyth, G. K. (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40, 4288–4297. doi: 10.1093/nar/gks042
Mohammad, A., Vanden Broek, K., Wang, C., Daryabeigi, A., Jantsch, V., Hansen, D., et al. (2018). Initiation of meiotic development is controlled by three post-transcriptional pathways in Caenorhabditis elegans. Genetics 209, 1197–1224. doi: 10.1534/genetics.118.300985
Naro, C., Jolly, A., Di Persio, S., Bielli, P., Setterblad, N., Alberdi, A. J., et al. (2017). An orchestrated intron retention program in meiosis controls timely usage of transcripts during germ cell differentiation. Dev. Cell 41, 82–93.e84. doi: 10.1016/j.devcel.2017.03.003
Okamato, S.-Y., Sato, M., Toda, T., and Yamamoto, M. (2012). SCF ensures meiotic chromosome segregation through a resolution of meiotic recombination intermediates. PLoS One 7:e30622. doi: 10.1371/journal.pone.0030622
Osman, K., Higgins, J. D., Sanchez-Moran, E., Armstrong, S. J., and Franklin, F. C. H. (2011). Pathways to meiotic recombination in Arabidopsis thaliana. New Phytol. 190, 523–544. doi: 10.1111/j.1469-8137.2011.03665.x
Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., and Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419. doi: 10.1038/nmeth.4197
Paytuví Gallart, A., Hermoso Pulido, A., Anzar Martínez de Lagrán, I., Sanseverino, W., and Aiese Cigliano, R. (2016). GREENC: a Wiki-based database of plant lncRNAs. Nucleic Acids Res. 44, D1161–D1166. doi: 10.1093/nar/gkv1215
Pérez-Rodríguez, P., Riaño-Pachón, D. M., Corrêa, L. G., Rensing, S. A., Kersten, B., and Mueller-Roeber, B. (2010). PlnTFDB: updated content and new features of the plant transcription factor database. Nucleic Acids Res. 38, D822–D827. doi: 10.1093/nar/gkp805
Plafker, S. M., Plafker, K. S., Weissman, A. M., and Macara, I. G. (2004). Ubiquitin charging of human class III ubiquitin-conjugating enzymes triggers their nuclear import. J. Cell Biol. 167, 649–659. doi: 10.1083/jcb.200406001
Qi, S. T., Ma, J. Y., Wang, Z. B., Guo, L., Hou, Y., and Sun, Q. Y. (2016). N6-Methyladenosine sequencing highlights the involvement of mRNA methylation in oocyte meiotic maturation and embryo development by regulating translation in Xenopus laevis. J. Biol. Chem. 291, 23020–23026. doi: 10.1074/jbc.M116.748889
Rapazote-Flores, P., Bayer, M., Milne, L., Mayer, C.-D., Fuller, J., Guo, W., et al. (2019). BaRTv1.0: an improved barley reference transcript dataset to determine accurate changes in the barley transcriptome using RNA-seq. BMC Genomics 20:968. doi: 10.1186/s12864-019-6243-7
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Schmid, R., Grellscheid, S. N., Ehrmann, I., Dalgliesh, C., Danilenko, M., Paronetto, M. P., et al. (2013). The splicing landscape is globally reprogrammed during male meiosis. Nucleic Acids Res. 41, 10170–10184. doi: 10.1093/nar/gkt811
Schreiber, M., Barakate, A., Uzrek, N., Macaulay, M., Sourdille, A., Morris, J., et al. (2019). A highly mutagenised barley (cv. Golden Promise) TILLING population coupled with strategies for screening-by-sequencing. Plant Methods 15:99. doi: 10.1186/s13007-019-0486-9
Schreiber, M., Mascher, M., Wright, J., Padmarasu, S., Himmelbach, A., Heavens, D., et al. (2020). A genome assembly of the barley ‘transformation reference’ cultivar golden promise. G3 10, 1823–1827. doi: 10.1534/g3.119.401010
Schwacke, R., Ponce-Soto, G. Y., Krause, K., Bolger, A. M., Arsova, B., Hallab, A., et al. (2019). MapMan4: a refined protein classification and annotation framework applicable to multi-omics data analysis. Mol. Plant 12, 879–892. doi: 10.1016/j.molp.2019.01.003
Soneson, C., Love, M., and Robinson, M. (2016). Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences [version 2; referees: 2 approved]. F1000Research 4:1521. doi: 10.12688/f1000research.7563.2
Song, X., Wang, D., Ma, L., Chen, Z., Li, P., Cui, X., et al. (2012b). Rice RNA-dependent RNA polymerase 6 acts in small RNA biogenesis and spikelet development. Plant J. 71, 378–389. doi: 10.1111/j.1365-313X.2012.05001.x
Szabo, E. X., Reichert, P., Lehniger, M. K., Ohmer, M., de Francisco Amorim, M., Gowik, U., et al. (2020). Metabolic labeling of RNAs uncovers hidden features and dynamics of the Arabidopsis transcriptome. Plant Cell 32, 871–887. doi: 10.1105/tpc.19.00214
Szcześniak, M. W., Bryzghalov, O., Ciomborowska-Basheer, J., and Makałowska, I. (2019). “CANTATAdb 2.0: expanding the collection of plant long Noncoding RNAs,” in Plant Long Non-Coding RNAs, eds J. A. Chekanova and H. L. V. Wang (New York, NY: Humana Press), 415–429. doi: 10.1007/978-1-4939-9045-0_26
Tang, X., Zhang, Z. Y., Zhang, W. J., Zhao, X. M., Li, X., Zhang, D., et al. (2010). Global gene profiling of laser-captured pollen mother cells indicates molecular pathways and gene subfamilies involved in rice meiosis. Plant Physiol. 154, 1855–1870. doi: 10.1104/pp.110.161661
Teng, C., Zhang, H., Hammond, R., Huang, K., Meyers, B. C., and Walbot, V. (2020). Dicer-like 5 deficiency confers temperature-sensitive male sterility in maize. Nat. Commun. 11:2912. doi: 10.1038/s41467-020-16634-6
Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578. doi: 10.1038/nprot.2012.016
Trincado, J. L., Entizne, J. C., Hysenaj, G., Singh, B., Skalic, M., Elliott, D. J., et al. (2018). SUPPA2: fast, accurate, and uncertainty-aware differential splicing analysis across multiple conditions. Genome Biol. 19:40. doi: 10.1186/s13059-018-1417-1
Veeneman, B. A., Shukla, S., Dhanasekaran, S. M., Chinnaiyan, A. M., and Nesvizhskii, A. I. (2016). Two-pass alignment improves novel splice junction quantification. Bioinformatics 32, 43–49. doi: 10.1093/bioinformatics/btv642
Venturini, L., Caim, S., Kaithakottil, G. G., Mapleson, D. L., and Swarbreck, D. (2018). Leveraging multiple transcriptome assembly methods for improved gene structure annotation. Gigascience 7:giy093. doi: 10.1093/gigascience/giy093
Walker, J., Gao, H., Zhang, J., Aldridge, B., Vickers, M., Higgins, J. D., et al. (2018). Sexual-lineage-specific DNA methylation regulates meiosis in Arabidopsis. Nat. Genet. 50, 130–137. doi: 10.1038/s41588-017-0008-5
Walter, W., Sánchez-Cabo, F., and Ricote, M. (2015). GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics 31, 2912–2914. doi: 10.1093/bioinformatics/btv300
Wang, H., You, C., Chang, F., Wang, Y., Wang, L., Qi, J., et al. (2014). Alternative splicing during Arabidopsis flower development results in constitutive and stage-regulated isoforms. Front. Genet. 5:25. doi: 10.3389/fgene.2014.00025
Wang, Y., and Yang, M. (2006). The Arabidopsis SKP1-LIKE1 (ASK1) protein acts predominately from leptotene to pachytene and represses homologous recombination in male meiosis. Planta 223, 613–617. doi: 10.1007/s00425-005-0154-3
Wierzbicki, A. T., Haag, J. R., and Pikaard, C. S. (2008). Noncoding transcription by RNA polymerase Pol IVb/Pol V mediates transcriptional silencing of overlapping and adjacent genes. Cell 135, 635–648. doi: 10.1016/j.cell.2008.09.035
Xu, M., Skaug, B., Zeng, W., and Chen, Z. J. (2009). A ubiquitin replacement strategy in human cells reveals distinct mechanisms of IKK activation by TNFalpha and IL-1beta. Mol. Cell 36, 302–314. doi: 10.1016/j.molcel.2009.10.002
Yamada, M., Yamagishi, Y., Akaoka, M., Ito, H., and Kato, A. (2014). Genomic localization of AtRE1 and AtRE2, copia-type retrotransposons, in natural variants of Arabidopsis thaliana. Mol. Genet. Genom. 289, 821–835. doi: 10.1007/s00438-014-0855-z
Yang, H., Lu, P., Wang, Y., and Ma, H. (2011). The transcriptome landscape of Arabidopsis male meiocytes from high-throughput sequencing: the complexity and evolution of the meiotic process. Plant J. 65, 503–516. doi: 10.1111/j.1365-313X.2010.04439.x
Yang, M., Hu, Y., Lodhi, M., McCombie, W. R., and Ma, H. (1999). The Arabidopsis SKP1-LIKE1 gene is essential for male meiosis and may control homologue separation. Proc. Natl. Acad. Sci. U.S.A. 96, 11416–11421. doi: 10.1073/pnas.96.20.11416
Yang, X., Timofejeva, L., Ma, H., and Makaroff, C. A. (2006). The Arabidopsis SKP1 homolog ASK1 controls meiotic chromosome remodeling and release of chromatin from the nuclear membrane and nucleolus. J. Cell Sci. 119, 3754–3763. doi: 10.1242/jcs.03155
Yuan, T. L., Huang, W. J., He, J., Zhang, D., and Tang, W. H. (2018). Stage-specific gene profiling of germinal cells helps delineate the mitosis/meiosis transition. Plant Physiol. 176, 1610–1626. doi: 10.1104/pp.17.01483
Zhai, J., Zhang, H., Arikit, S., Huang, K., Nan, G.-L., Walbot, V., et al. (2015). Spatiotemporally dynamic, cell-type-dependent premeiotic and meiotic phasiRNAs in maize anthers. Proc. Natl. Acad. Sci. U.S.A. 112, 3146–3151. doi: 10.1073/pnas.1418918112
Zhang, F., Tang, D., Shen, Y., Xue, Z., Shi, W., Ren, L., et al. (2017). The F-box protein zygo1 mediates bouquet formation to promote homologous pairing, synapsis, and recombination in rice meiosis. Plant Cell 29, 2597–2609. doi: 10.1105/tpc.17.00287
Zhang, R., Calixto, C. P. G., Marquez, Y., Venhuizen, P., Tzioutziou, N. A., Guo, W., et al. (2017). A high quality Arabidopsis transcriptome for accurate transcript-level analysis of alternative splicing. Nucleic Acids Res. 45, 5061–5073. doi: 10.1093/nar/gkx267
Zhang, H., Egger, R. L., Kelliher, T., Morrow, D., Fernandes, J., Nan, G. L., et al. (2014). Transcriptomes and proteomes define gene expression progression in pre-meiotic maize anthers. G3 4, 993–1010. doi: 10.1534/g3.113.009738
Zhang, M., Ma, X., Wang, C., Li, Q., Meyers, B. C., Springer, N. M., et al. (2020). CHH DNA methylation increases at 24-PHAS loci depend on 24-nt phasiRNAs in maize meiotic anthers. New Phytol. [Online ahead of print]. doi: 10.1111/nph.17060
Zhao, Y., Li, H., Fang, S., Kang, Y., Wu, W., Hao, Y., et al. (2016). NONCODE 2016: an informative and valuable data source of long non-coding RNAs. Nucleic Acids Res. 44, D203–D208. doi: 10.1093/nar/gkv1252
Ziolkowski, P. A., Underwood, C. J., Lambing, C., Martinez-Garcia, M., Lawrence, E. J., Ziolkowska, L., et al. (2017). Natural variation and dosage of the HEI10 meiotic E3 ligase control Arabidopsis crossover recombination. Genes Dev. 31, 306–317. doi: 10.1101/gad.295501.116
HSP90.7 HEAT SHOCK PROTEIN 90.7
PTB1a POLYPYRIMIDINE TRACT BINDING PROTEIN 1a
HEI10 HUMAN ENHANCER OF INVASION 10
NEDD8 NEURAL PRECURSOR CELL EXPRESSED DEVELOPMENTALLY DOWN-REGULATED PROTEIN 8
KDM3 LYSINE SPECIFIC DEMETHYLASE 3
PKDM9 PUTATIVE LYSINE SPECIFIC DEMETHYLASE 9
DCL5 DICER-LIKE 5
RDR6 RNA DEPENDENT RNA POLYMERASE 6
ASY1 ASYNAPTIC 1
ZYP1 MOLECULAR ZIPPER 1-LIKE PROTEIN
DMC1 DISRUPTED MEIOTIC CDNA 1
PCNA PROLIFERATING CELL NUCLEAR ANTIGEN
MET1B DNA (CYTOSINE-5)-METHYLTRANSFERASE 1B
MRE11B MEIOTIC RECOMBINATION 11 HOMOLOG B
XRCC3 X-RAY REPAIR CROSS COMPLEMENTING 3
DDM1C DECREASED DNA METHYLATION 1C
BUB3B BUDDING UNINHIBITED BY BENZIMIDAZOLES 3B
SDS SOLO DANCERS
SWI SWITCH 1/DYAD
RPA2c REPLICATION PROTEIN A2C
RPA1c.1 REPLICATION PROTEIN A1C.1
Keywords: barley, anther, meiocyte, transcriptome, meiosis, argonaute, lncRNAs, ubiquitin
Citation: Barakate A, Orr J, Schreiber M, Colas I, Lewandowska D, McCallum N, Macaulay M, Morris J, Arrieta M, Hedley PE, Ramsay L and Waugh R (2021) Barley Anther and Meiocyte Transcriptome Dynamics in Meiotic Prophase I. Front. Plant Sci. 11:619404. doi: 10.3389/fpls.2020.619404
Received: 20 October 2020; Accepted: 15 December 2020;
Published: 12 January 2021.
Edited by:Raju Datla, Global Institute for Food Security (GIFS), Canada
Reviewed by:Stefanie Dukowic-Schulze, University of Regensburg, Germany
Shun Sakuma, Tottori University, Japan
Copyright © 2021 Barakate, Orr, Schreiber, Colas, Lewandowska, McCallum, Macaulay, Morris, Arrieta, Hedley, Ramsay and Waugh. 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: Robbie Waugh, email@example.com
†These authors have contributed equally to this work