The transgenic line expressing ACT2:BirA, kindly given by Prof. R. Deal (Emory University, United States), and the reference ecotype Col-0 of A. thaliana (NASC stock N60000) were used in this work. The transgenic line was transformed by Agrobacterium tumefaciens (GV3101) according to floral dip method (Clough and Bent, 1998) to obtain plants carrying DMC1:NTF along with ACT2:BirA. About 1,500 seeds obtained after transformation were sterilized with 70% ethanol for 1 min., and bleach solution (10% commercial bleach and 0.05% Tween 20) for 10 min., finally washed with sterile water three times. Sterilized seeds were sown on a selective modified MS (Murashige and Skoog, 1962) medium containing hygromycin (50 mg/l) and 0.8% agar, kept for 2–3 days at 4°C, and germinated under long-day conditions (16 h light/8 h darkness) at 24°C. Seedlings transferred to pots were grown in a controlled growth chamber under long-day conditions. The ADMC1:NTF/ACT2:BirA double homozygous transgenic line (D.1) was selected from ten T₂ independent transformant lines. About 200 plants from the D.1 line were grown as above reported for subsequent experiments.

To confirm transgene insertion, PCR was conducted on the hygromycin resistance gene used as a selectable marker. Genomic DNA was extracted from leaf tissue of transgenic plants using the DNeasy Plant Mini Kit (Qiagen)¹ according to the manufacturer’s instructions. PCR reactions were performed on genomic DNA using primers listed in Supplementary Table 1. RT-qPCR was performed to verify expression of target BirA gene under the constitutive promoter ACT2 in the floral buds. Total RNA from transgenic ACT2:BirA plants was extracted using RNeasy plant minikit (QIAGEN) following the manufacturer’s instructions. mRNA was retrotranscribed to cDNA using Superscript III reverse transcriptase (Thermo Fisher Scientific, United States) and oligo dT(20) following the manufacturer’s conditions, and relative expression was verified by real-time qPCR. The experiment was performed on a QIAGEN Rotorgene 6,000 qRT-PCR machine, using Power Sybr Green real time mix (Thermo Fisher Scientific, United States). The reaction conditions were as follows: one cycle at 95°C for 10 min, and 40 cycles of 95°C × 10″ denaturation and 60°C × 45″ annealing and extension. The melting curve was run to verify the specificity of the primers. The ADENINE PHOSPHORIBOSYL TRANSFERASE 1 (APT1) gene was used as a housekeeping internal control. Target and housekeeping primers are listed in Supplementary Table 2.

The fusion gene NTF carried in the vector ADF8-NTF, kindly given by Prof. R. Deal (Emory University, United States), and the AtDMC1 promoter carried in the vector SLJ7753, kindly given by Prof. J.D.G. Jones (The Sainsbury Laboratory, Norwich, United Kingdom), were amplified with specific primers carrying GATEWAY adapters. Subsequently, they were cloned into Gateway vectors pDONR/Zeo and pDONR/P4P1R, respectively (Thermo Fisher Scientific). Generated entry clones were recombined with multisite GATEWAY reaction into destination vector pH7m24GW,3², and the final expression clone, named pEXPR-DMC1-NTF, was generated. All gateway reactions were performed following the manufacturer’s standard protocols. Primers used for vector construction are listed in Supplementary Table 3.

Flower buds of the size corresponding to male meiotic stage (Smyth et al., 1990) were collected, fixed with 1% formaldehyde solution, and stored at −80°C. Nuclei were purified from 1.5 g of frozen and homogenized tissue as described previously by Wang and Deal (2015) with some minor modifications. In particular, the tissue immersed in 40 ml of Nuclei Purification Buffer solution (NPBf) containing 20 mM MOPS (pH 7), 40 mM NaCl, 90 mM KCl, 2 mM EDTA, 0.5 mM EGTA, 0.5 mM spermidine, 0.2 mM spermine, 1% (v/v) formaldehyde was incubated in 50 ml tube under vacuum glass desiccator for 10 min., and vacuum release for 2 min. The procedure was repeated once. Then, glycine was added to a final concentration of 0.125 M under vacuum for 7 min. After NPBf solution elimination, the tissue was washed three times with water. Isolated nuclei were resuspended into RNAlater buffer, and RNA extraction was performed immediately.

Laser-scanning confocal microscopy imaging was performed using a confocal Zeiss LSM 510. Isolated nuclei were observed using a Florescence Microscope (Leitz Aristoplan).

Total RNA was extracted from three replicates of isolated nuclei using an RNeasy Extraction Mini Kit (Qiagen, see text footnote 1) according to the manufacturer’s instructions. Isolated nuclei derived from the same population of inflorescences where the buds different from stage 9 (Smyth et al., 1990) had been manually dissected. Quantification was performed on a QUBIT fluorometer.

RNA sequencing was performed by IGA Technology Services Srl³. The libraries were produced using retrotranscribed cDNA previously amplified by Ovation Ultralow Library System V2 (NuGEN Technologies, Inc.). Library size and integrity were assessed using the Agilent Bioanalyzer (Santa Clara, CA) or Caliper GX (PerkinElmer, MA). Sequencing was performed by Illumina HiSeq 2,500 (Illumina, San Diego, CA) and 30-M paired-end reads (2 × 125) per replicate were generated.

Raw sequencing data (FASTQ files) were quality checked using the software FASTQC v0.11.5⁴, then low-quality bases and adapter sequences were removed with the software BBDuk v35 (Bushnell et al., 2017) setting 35 bp as the minimum read length after trimming. The read length ranged from 35 to 125 bp after trimming. More than 98% of the final reads had a length of 125 bp since adapter contamination was very low. The trimmed reads were then mapped against the A. thaliana reference genome (Araport11) (Cheng et al., 2016) using STAR v2.7.3a with the option – alignEndsProtrude 100 DiscordantPair. The statistics from the mapping step were produced with Qualimap v2.2.1 (Okonechnikov et al., 2016). Kallisto v0.46.0 was used to obtain transcript expression quantification levels, as estimated counts and TPMs, which were then summarized as gene expression values using the R package “tximport” (Soneson et al., 2016).

The correlation analysis of the samples was performed using the R in-built function “cor,” and the results were plotted with the function “corrplot.” The classification of the gene classes and the saturation analysis were performed with the R package “NOISeq” (Tarazona et al., 2015). Genes were classified based on their average expression levels using the following criteria: “undetected” if the expression was 0 across all the replicates, “Low” if log2 average TPM ≤ 1.252, “Low_Mid” if log2 average TPM > 1.252 and ≤ 2.207, “Mid_High” if log2 average TPM > 2.207 and ≤ 4.169, and “High” if log2 average TPM > 4.169. The abovementioned values correspond to the 25th, 50th, and 75th percentiles of the log2 average TPM distributions, respectively.

In order to compare the expression profile of the meiocytes against other datasets, TPM values from the following datasets were downloaded from EBI: E-GEOD-38612, E-GEOD-55866, and E-MTAB-4202 (Supplementary Table 4). In addition, gene expression profiles from pollen and tapetum cells were obtained from Loraine et al. (2013) and Li et al. (2017). When multiple replicates of the same tissue were available, an average was calculated. Similarly, when multiple developmental stages of the same tissue were available, the average was calculated. The only exception consisted of the tapetum datasets kept separate because the two available stages were highly different. Finally, a complete expression matrix was obtained by merging all the datasets and the ARSyN function from the NOISeq package was used to remove the batch effect. A PCA analysis was then performed with the in-built “prcomp” function in R and the results plotted with “ggplot2.”

A gene expression network was generated using the expression matrix obtained after the ARSyN correction, following the Aracne algorithm implemented in the R package “parmigene” (Sales and Romualdi, 2011). Cytoscape v3.8.1 was then used to plot the network and select only the first neighbors of known meiotic genes in order to find new candidates. The gene ontology (GO) analysis of the obtained genes was then performed using the Cytoscape app named “BinGO” using the hypergeometric test as a statistical test and an FDR lower than 0.05 as threshold. All the heat maps were generated using Clustvis (Metsalu and Vilo, 2015).

Gene ontology enrichment analysis (GOEA) were performed using in-house scripts based on the method described in Tian et al. (2017) and setting a minimum FDR threshold of 0.05.

In order to isolate meiocyte nuclei by the INTACT method for RNA-seq analysis in Arabidopsis, we generated transgenic material carrying a NTF protein under the meiosis-specific AtDMC1 promoter (Klimyuk and Jones, 1997) along with the biotinilase BirA under a constitutive ACTIN2 (ACT2) promoter. The promoter of ACT2 (At3g18780) was used by Deal and Henikoff (2011) to drive expression of BirA in root cell types (An et al., 1996). Following detection of ACT2:BirA expression in floral buds (not shown), we used the same transgenic line generated by Deal and Henikoff (2011) as starting material. In this line, we introduced the new construct carrying the AtDMC1:NTF expression cassette. Although AtDMC1 (At3g22880) was expected to drive NTF expression in both male and female meiocytes (Klimyuk and Jones, 1997), we restricted our confocal microscopy analysis to anthers and male meiocytes. Indeed, the flower buds were collected at a stage corresponding to male meiosis which occurs earlier than female meiosis (Schneitz et al., 1995). Confocal microscopic examination of the AtDMC1:NTF/ACT2:BirA line showed that NTF was expressed in the expected cell type. Indeed, we detected GFP-positive male meiocytes during prophase I stage, arguably before the nuclear envelope breakdown (NEBD; Figure 1). To purify labeled nuclei from meiocytes, we extracted total nuclei from whole inflorescences with buds at floral stage 9 (Smyth et al., 1990) in the AtDMC1:NTF/ACT2:BirA line. Then, the nuclei were incubated with streptavidin-coated magnetic beads according to INTACT protocols (Deal and Henikoff, 2011; Wang and Deal, 2015). After nucleus isolation, we used a fluorescence microscope to observe the complex DAPI-stained nuclei/beads (Figure 2).

To obtain the meiotic nuclear transcriptome, RNA sequencing (RNA-seq) was performed downstream of INTACT. RNA extracted from isolated nuclei of meiocytes in three replicates was sequenced by using the Illumina technology following an amplification step. Pearson’s correlation coefficient indicated the reliability of the experiment (Supplementary Figure 1). After trimming, an average of 85,927,325 reads per replicate were obtained (Supplementary Table 5). About 77% of these reads were mapped onto the Arabidopsis genome (Supplementary Table 6). The percentages (calculated as average of the replicates) of the uniquely mapped reads located in known exons, introns, and intergenic regions are 38, 26, and 36%, respectively (Supplementary Figure 2). Given the high percentage of reads that mapped to multiple positions (Supplementary Table 6), the gene expression, measured as TPM (Transcript Per Million), was calculated by a specific algorithm (Kallisto) which is able to process multiple mapping reads (Supplementary Table 7). A frequency ranging between 75.2 and 81.3% of the genes (total no. 32,833) was detected across the replicates (Supplementary Figure 3). Transcript types of each replicate are summarized in Figure 3 showing the distribution of expression profiles for the different gene classes. Ribosomal RNA appeared the most abundant class whereas lncRNA is the less expressed class. This result is consistent with a recent study in barley in which 65% of the downregulated DEGs were lncRNAs in leptotene/zygotene vs. pre-meiosis comparison (Barakate et al., 2020).

Transcripts were grouped into four classes according to their expression levels from low to high (Supplementary Figure 4). GOEA was performed to identify enriched (i.e., over-represented) GO terms associated with the different expression classes (Supplementary Figures 5A–D). The most enriched GO terms, such as “DNA demethylation” in the mid/high-expression class and “RNA stabilization” in the high expression class, suggest that these biological processes are particularly relevant to meiosis. To assess the reliability of our results, we surveyed the genes with documented functions in Arabidopsis meiosis. A list of 197 meiotic genes was implemented using GO terms GO:0051321 (meiotic cell cycle) and GO:0140013 (meiotic nuclear division) (Supplementary Table 8). The majority of meiotic genes were distributed in the classes of mid/high and high expression while a small number were in the low-expression class (Figure 4). By comparison, this class was the third most represented for the total of transcripts (Supplementary Figure 4).

To further characterize our meiotic transcriptome, we compared it with publicly available transcriptomic data from other specific tissues and cell types of Arabidopsis. Initially, we planned to compare our RNA-seq data with those from male meiocytes isolated by micromanipulation in Arabidopsis (Chen et al., 2010; Yang et al., 2011) but, unfortunately, these datasets are not available.

Based on the low expression of AtDMC1 (Supplementary Figure 6), we discarded replicate 1 in all the subsequent analyses. Principal component analysis (PCA) performed with meiocytes (current study), meristem, leaf, root, flower buds, tapetum (from two different developmental stages), pollen, and silique revealed that meiocytes clustered close to flower buds, as expected, but also to tapetum and pollen (Figure 5). We evaluated the expression profile of genes (list reported in Supplementary Table 9) considered to be specific for these two types of cells (Loraine et al., 2013; Li et al., 2017). The heat map revealed expression patterns specific for meiocytes, tapetum, and pollen (Figure 6). On the other hand, a hierarchical sample dendrogram showed that meiocyte sample groups together with tapetum 6–7 as well as pollen with tapetum 8–10 equivalent to a later stage.

Afterward, we assessed the overall gene expression of meiotic genes (list reported in Supplementary Table 8) and of the other transcripts in meiocytes and multiple tissues (Figure 7). The analysis revealed that the median expression value of meiotic genes was higher in isolated meiocytes when compared to the other tissues/cell types with the exception of floral bud and meristem samples. This result is not surprising since the meiocytes are included in floral buds and they share cell-cycle genes with the meristem. The median expression value of the other transcripts was basically steady across the different tissues (Figure 7). Considering the meiotic gene expression fold changes between meiocytes and the other tissues/cell types, the heat map revealed a pattern in agreement with the above reported analysis (Figure 8 and Supplementary Table 10).

To identify new candidate genes with a meiotic function, we generated a gene expression network based on the genes known to be involved in meiosis (Supplementary Table 8). By selecting the first neighbors of the known meiotic genes, we found 1,503 genes with a total of 7,607 connections (Figure 9 and Supplementary Table 11). The network is characterized by a very large cluster, thereby suggesting that the known meiotic genes and their neighbors are all highly connected as also indicated by an average degree of 4.95 connections. The most connected genes with 435 and 323 connections have a documented role in meiosis. In particular, FIDGETIN-LIKE-1 INTERACTING PROTEIN (FLIP, AT1G04650) forms a protein complex with FIDGETIN-LIKE-1 (FIGL1) that is conserved from Arabidopsis to human, and it regulates meiotic crossover formation via RAD51 and DMC1 (Fernandes et al., 2018). The other gene, AXR1 (AT1G05180), is involved in the neddylation/rubylation protein modification pathway and TE methylation in meiocytes (Jahns et al., 2014; Christophorou et al., 2020). AXR1 plays a significant role in DNA repair (Martinez-Garcia et al., 2020). Collectively, this finding reinforces the reliability of the analysis performed in this study. New candidates which could play a role in meiosis are two transcription factors (TFs), AT1G06070 and AT1G02220, with 64 and 30 connections, respectively. These two genes belong to the bZIP and NAC TF families, and their function in meiosis has not been documented. Although with less connections than the genes reported above, AT1G04200, described as Dymeclin (DYM, Dyggve–Melchior–Clausen syndrome protein), is an interesting candidate for a possible meiotic function. Indeed, the mouse DYM homolog has a high expression in testis (Yue et al., 2014) and has a proved interaction with FANCD2, a component of the Fanconi anemia DNA repair pathway (Zhang et al., 2017).

A GO analysis performed using the genes identified by the network analysis revealed GO targets associated with meiotic processes and, particularly, to Meiosis I such as recombination, synapsis, chiasma assembly, and meiotic chromosome segregation (Figure 10 and Supplementary Table 12).