Tissue-Specific Transcriptome Analysis Reveals Candidate Transcripts Associated with the Process of Programmed B Chromosome Elimination in Aegilops speltoides

Some eukaryotes exhibit dramatic genome size differences between cells of different organs, resulting from programmed elimination of chromosomes. Here, we present the first transcriptome analysis of programmed chromosome elimination using laser capture microdissection (LCM)-based isolation of the central meristematic region of Aegilops speltoides embryos where B chromosome (B) elimination occurs. The comparative RNA-seq analysis of meristematic cells of embryos with (Bplus) and without Bs (B0) allowed the identification of 14,578 transcript isoforms (35% out of 41,615 analyzed transcript isoforms) that are differentially expressed during the elimination of Bs. A total of 2908 annotated unigenes were found to be up-regulated in Bplus condition. These genes are either associated with the process of B chromosome elimination or with the presence of B chromosomes themselves. GO enrichment analysis categorized up-regulated transcript isoforms into 27 overrepresented terms related to the biological process, nine terms of the molecular function aspect and three terms of the cellular component category. A total of 2726 annotated unigenes were down-regulated in Bplus condition. Based on strict filtering criteria, 341 B-unique transcript isoforms could be identified in central meristematic cells, of which 70 were functionally annotated. Beside others, genes associated with chromosome segregation, kinetochore function and spindle checkpoint activity were retrieved as promising candidates involved in the process of B chromosome elimination.


Introduction
The genetic information is believed to be uniform in somatic and gamete precursor cells in most organisms. However, some species demonstrate regular elimination of specific DNA fragments as a part of the developmental program. This occurs either as a loss of entire chromosomes, chromosome fragments or DNA sequences, often during the differentiation of somatic and germline cells. The so-called programmed DNA elimination takes place in the ontogenesis of some unicellular ciliates and diverse metazoa species (reviewed by [1,2]).
In plants, only B chromosomes are known to be associated with programmed DNA elimination [3]. B chromosomes (Bs) can be observed in some plant, animal and fungi species as a dispensable addition to the basic karyotype. Their presence may be neutral in regard to the phenotypic characteristics of the host or it may exert detrimental effects on the fitness and fertility, often depending on the number of Bs To evaluate the impact of B chromosome elimination on the embryonic transcriptome and to identify candidate transcripts associated with the process of B elimination, we conducted a comparative analysis of central meristematic regions of B-containing embryos and embryos without Bs. Tissue-specific RNA-seq revealed that the elimination of Bs has a strong effect on the transcriptome of the meristematic central zone of young embryos. In particular, 35% from 41,615 analyzed transcript isoforms showed statistically significant altered expression during the elimination process. These transcripts are either associated with the process of B chromosome elimination, transcribed by Bs irrespectively of the elimination process or transcribed by the standard A chromosomes in response to the presence of Bs. Genes associated with chromosome segregation stand out as candidates involved in the process of B chromosome elimination.

The Transcriptome of Embryos Undergoing B Chromosome Elimination Is Strongly Affected
Elimination of A. speltodies B chromosomes occurs in the central meristematic zone between the developing apical meristem and embryonic root [3]. To evaluate whether the process of B chromosome elimination affects the transcriptome of developing embryos, the central meristematic regions of 17 to 20 days after anthesis (DAA) embryos of mother plants with (called Bplus) and without B chromosomes (called B0) were isolated by LCM (Figure 2), and RNA-seq was performed for transcriptome profiling. Laser capture microdissection of an Ae. speltoides embryo. (a) Schema is depicting an embryo with B chromosomes at the stage of 17 days after anthesis. Blue color indicates embryo parts containing Bs, white color shows the absence of Bs in the radicle. The red circle marks the region of ongoing Bs elimination between the apical meristem and developing radicle. SC-scutellum, AP-apical meristem, RA-radicle, CO-coleorhiza. Cryosection of a 17 DAA embryo (b) before and (c) after laser dissection of the region where Bs undergo elimination.
A de novo transcriptome assembly based on Trinity was applied since Ae. speltoides represents a non-model organism that contains in addition a high proportion of heterozygous loci. Trinity is highly effective for reconstructing transcripts and alternative spliced isoforms [12], and it allowed a better understanding of which transcript isoforms are truly expressed under a specific condition such as elimination of Bs. The de novo transcriptome assembly resulted in B0 and Bplus FASTA files containing 165,818 and 196,360 contigs, respectively. Subsequent filtering resulted in 117,441 contigs for downstream analysis (N50 contig size of 1026 bp). Detailed information about the assembly statistics is given in Table S2, Supplementary Materials. Transcriptome profiles of Ae. speltoides were evaluated at the transcript isoform expression level and quality was checked with the script PtR. A correlation matrix and PCA confirmed biological replications ( Figure 3, Figure S1, Supplementary Materials). Differential expression analysis identified 14,578 transcript isoforms that were differentially expressed during B elimination (FC of log 2 > 1, FDR-corrected p-value < 0.01, Table 1, Figure 4). Transcript isoforms significantly changing their expression during B chromosome elimination can either reflect the expression from B-located genic sequences or represent a response of A-located genes as a consequence of B elimination.  The high number of differentially regulated transcript isoforms (35% out of 41,615 analyzed transcript isoforms) indicates that the process of B chromosome elimination and the presence of B chromosomes have a strong effect on the transcriptome of the embryonic central meristematic region ( Figure 4). After stringent filtering, we could define 502 B0-specific (normalized mean read amount in Bplus is 0) and 341 Bplus-specific (normalized mean read amount in B0 is 0) transcript isoforms. In other cases, due to the high sequence similarity of A and B chromosome-located coding sequences [3], it was not possible to distinguish the origin of transcripts.
To verify RNA-seq data, reverse transcription quantitative PCR (RT-qPCR) was performed with a subset of transcripts known to be involved in cytokinesis and microtubule organization, similarly found to be B0-or Bplus-specific ( Figure 5). The kinetochore gene Nuf2 (nuclear filament-containing protein), Bplus-specific BREVIS RADIX (BRX), a regulator of cell proliferation and elongation in the root and shoot, SMC2, a component of the chromosome condensin complex, and TPX2, encoding a Xklp2-targeting protein, showed nearly no expression in the meristematic region of B0 embryos. Preferential expression in that cellular region of Bplus embryos was also shown for the chromosome segregation protein SPC25, whereas B0-specific expression was confirmed for a different isoform of BRX. Transcript levels and similar expression differences between candidate transcripts obtained by RT-qPCR validated the quality of the RNA-seq data analysis.

Biological Processes Potentially Related to the Elimination of B Chromosomes
The final data set containing non-redundant B0 and Bplus transcripts (117,441 contigs in total) was screened using Transdecoder v5.3.0 to predict open reading frames and resulted in 61,817 protein-coding transcripts. Such a high number of transcripts is in line with previous findings [13], where a global gene expression atlas of developing embryos of wheat and their putative diploid ancestors including Ae. speltoides was provided.
The completeness of the transcriptome assembly was assessed by BUSCO. Detailed results are displayed in Table S2, Supplementary Materials. The level of transcriptome completeness depends on the amount of tissues and developmental stages within the samples [14]. Therefore, we would expect a limited set of detected BUSCOs in our tissue-specific analysis of meristematic embryo cells. Nevertheless, we were able to recover a total fraction of~65% BUSCO-confirmed genes (including fragmented and complete). Comparing our findings to a recently published study in an oat embryonic transcriptome detecting 60% complete and fragmented orthologues [15], the BUSCO results for the tissue-specific Ae. speltoides transcriptome are highly similar.
The resulting protein-coding transcripts were annotated based on Interproscan Galaxy [16,17] and 26,722 transcript isoforms with known functions could be obtained, which are represented in 23,836 unigenes (threshold of 1.0 × 10 −3 ). The differentially expressed transcripts were aligned [18] against the hexaploid wheat genome assembly (IWGSC Ref.Seq.1.0, HighConf_Protein_2017), and transcripts were selected using the first best hit criterion with a minimum length of matching sequences of 100 amino acids. As a result, 8415 differentially expressed transcript sequences were annotated using proteins of the wheat genome.
To determine overrepresented biological processes associated with the elimination of the B chromosomes, the Gene Set Enrichment Analysis tool (GSEA, AgriGo v2. analytical toolkit) with wheat as a background reference was used. Separated analyses of up-and down-regulated transcript isoforms undergoing B chromosome elimination were performed. The results revealed different enriched GO terms associated with down-and up-regulated transcript isoforms.
The dataset of up-regulated transcripts in meristematic regions of B-containing embryos (Bplus) contains 2908 annotated wheat genes. Among them, 27 terms for biological processes, nine terms for molecular function and three terms for cellular component were overrepresented (FDR-corrected p-value < 0.05, Table S3, Supplementary Materials). In the category of biological processes, "cellular process" (989, 34.0% from dataset) and "nitrogen compound metabolic process" (467, 16.1%) were highly enriched. Other overrepresented terms were "gene expression", "microtubule-based process", "movement of cell or subcellular component", "microtubule-based movement" and "organelle fission". Terms in the molecular function category, particular "guanyl nucleotide binding" and "GTP binding", but also other enriched categories, such as "cytoskeletal protein binding", "microtubule binding", "microtubule motor activity" and "tubulin binding" coincide with the biological process category and additionally highlight that cytoskeleton/microtubule organization and vesicle trafficking are mostly affected in cells undergoing B elimination. In the cellular component category, only three terms, namely "cytoplasm", "chromosome, centromeric region" and "chromosomal region", were enriched. The results of GO enrichment are depicted in Figure 6, Figure S2 and Figure S3 (Supplementary Materials). Similar enriched GO terms ("cellular process", "microtubule-based movement", "cell division") were identified by [19] in the analysis of rye B transcripts of leaf and anther tissues. Some kinesin-like proteins are up-regulated in Bplus embryos and are deemed to be involved in microtubule binding, movement and motor activity (Table S4, Supplementary Materials). Furthermore, a transcript (Bplus_DN37268_c0_g1_i1) encoding Katanin p60 ATPase-containing subunit A-like 1 (p60 katanin-like 1) was highly up-regulated in the Bplus condition (log 2 FC of 6.16, FDR-corrected p-value of 0.0009). Katanin is the only defined microtubule severing protein affecting microtubule organization in higher plants [20].
Three transcript isoforms of the Nuf2 gene, which is crucial for chromosome segregation and spindle checkpoint activity, were significantly up-regulated (log 2 FC value > 1, FDR-corrected p-value < 0.05) in Bplus samples (Table S5, Supplementary Materials), of which isoform c_Bplus_DN38911_c1_g1_i9 shows the highest specificity for Bplus condition (see also Figure 5). Nuf2 is evolutionary conserved and known to play a role in spindle checkpoint control by regulating the bipolar attachment of microtubules of the sister chromatids before and after anaphase. Overexpression of Nuf2 in mammals also caused defects in chromosome segregation [21]. Hence, as deregulation of Nuf2 expression results in chromosomal segregation defects, the balanced activity of Nuf2 is essential for regular chromosome segregation. In the context of the observed B-specific overexpression of Nuf2 in embryos undergoing chromosome elimination, it is tempting to speculate that overexpression per se or expression of B-specific Nuf2 variants is part of the observed B chromosome segregation alterations. Further transcripts encoding the spindle checkpoint controlling component Shugoshin (Sgo1), a Mis12 protein that is involved in kinetochore formation (Table S5, Supplementary Materials), and the spindle and kinetochore-associated protein 2 (Ska2) are highly activated in cells undergoing B chromosome elimination.
The dataset of down-regulated transcripts in meristematic regions of B-containing embryos is represented with 2726 annotated genes. The category biological processes included 38 enriched terms, such as "DNA repair", "methylation" and "compound nitrogen metabolic process". In the category molecular function, only "oligosaccharyl transferase activity" was significantly enriched. The category cellular component includes 16 enriched terms with the most prominent terms "cell" and "cell part" (Table S6, Supplementary Materials). Gene set enrichment analysis demonstrated extensive modulation of genes involved in stress response (term "cellular response to stress", Table S7, Supplementary Materials) during the process of B chromosome elimination. For example, a gene encoding an important tumor suppressor, BRCA1, was highly down-regulated in the Bplus condition. BRCA1, known as a regulator of crucial cellular processes, protects cells from aneuploidy and genomic instability. Its down-regulation promotes aberrant mitoses and aneuploidy and thereby it is a hallmark of cancer [22]. Further, genes encoding proteins such as DNA repair protein RAD51 (RAD51D), DNA repair helicase, DNA mismatch repair protein and DNA damage checkpoint protein Rad9 were highly down-regulated in meristematic cells undergoing B chromosome elimination indicating that chromosomal repair mechanisms might be defected.

B Chromosome-Specific Transcripts
Based on strict filtering criteria, 341 B-unique transcript isoforms were identified and 70 of them were annotated (Table S8, Supplementary Materials). Transcript isoforms detected during chromosome elimination represented several categories: transposons (such as transposase, retrotransposon); potential transposons (zinc finger); unknown; and genes with specified functions. Among genes with specified functions, there were two genes encoding kinesin-like proteins. One full-length protein showed 94% identity to the kinesin-like protein KIN-14C of Ae. tauschii subsp tauschii. In A. thaliana, KIN_14C is important for regular microtubule accumulation at the spindle poles during the prophase of mitosis. Interestingly, B chromosome-specific genes related to the microtubules and cell division were also identified in several species [5].
The reason why B chromosomes undergo elimination in roots is unknown. However, it was suggested that elimination of Bs counteracts negative effects potentially associated with the fitness and fertility of the plants [3]. Our tissue-specific study revealed that some of the B chromosome-specific transcript isoforms encode regulators of plant growth and development. One of the B-specific partial transcripts (c_Bplus_DN27718_c0_g1_i3) encodes a protein of the OVATE family. The plant-specific proteins of this family act as transcription repressors [23]. In A. thaliana, overexpression of AtOFP1 results in dwarf plants with reduced cell elongation in rapidly elongating aerial organs, including hypocotyl, leaf petiole and inflorescence stems [24]. Another example is the B-specific transcript (c_Bplus_DN30654_c0_g1_i2) encoding an IAA9/IAA20 orthologue of A. thaliana. The members of the Aux/IAA family are essential during plant development, such as root development, shoot growth, flower organ development and fruit ripening [25].
In summary, our tissue-specific RNA-seq analysis identified candidate transcripts which might be involved in B chromosome elimination in Ae. speltoides embryos. Tissue-specific sampling of the cellular region where B elimination occurs enabled the capturing of rare transcripts and cell-specific differences, which are otherwise hidden in whole-embryo samples. The study demonstrates that the process of B chromosome elimination and presence of B chromosomes has a surprisingly strong effect on the global transcriptome of meristematic cells in young embryos. These transcripts are either associated with the process of B chromosome elimination or activated in response to the presence of B chromosomes. The identification of B-specific candidate transcripts provides the basis for future work on the elucidation of molecular mechanisms underlying B chromosome elimination and studies validating the function of candidate genes in planta.

Plant Material
Aegilops speltoides Tausch (PI 487238; USDA-ARS, Aberdeen, ID, USA) plants with and without B chromosomes (further referred to as Bplus and B0) were grown under greenhouse conditions (16 h light, 25 • C/19 • C day/night) at IPK, Gatersleben (Germany). At the beginning of the tillering stage, plants were kept for one month at 12 • C to ensure synchronous flowering. Presence of Bs in plants was identified by PCR with primers for the B-specific repeat AesTR-183 [26] (Table S1, Supplementary Materials) and the exact number of Bs was determined by flow-cytometric analysis according to [27]. Only plants with 2, 3 or 4 Bs were used for the generation of embryos. To ensure the presence or absence of Bs in the embryos, plants with and without Bs were kept in separate greenhouse chambers to prevent cross-pollination. The time of anthesis was recorded for each flower in the spikes. Embryos were sampled at 17 to 20 DAA, immediately frozen in liquid nitrogen and stored at −80 • C.

Laser Capture Microdissection (LCM) and RNA Extraction
Embryos were transferred to a cryostat (−20 • C) and then glued onto sample plates using the Tissue-Tek ® OCT™ compound (Sakura Finetek Europe BV, Alphen aan den Rijn, The Netherlands). Serial cross-sections of 20 µm thickness were prepared using a cryotome (Cryostar NX 70, Microm GmbH, Neuss, Germany), mounted on RNAse-free PEN membrane slides (MMI, Eching, Germany) and stored in the cryostat at −20 • C until complete dryness (5-7 days). Prior to microdissection, cryosections were adapted to room temperature for some minutes. Cells were isolated from the central meristematic root region of Bplus and B0 embryos (Figure 2), where the elimination of Bs was shown to occur [3]. LCM-based isolation using the MMI Cell Cut system, RNA isolation, mRNA amplification and RNA-seq library preparation has been performed as described in [28] with slight modifications.
Total RNA was extracted from 30-55 tissue sections for each of the three biological replicates from different genomic backgrounds using the Absolutely RNA Nanoprep Kit (Agilent Technologies, Santa Clara, CA, USA). From 3 to 4 embryos were pooled per one replicate. Number of tissue sections obtained from each embryos and total square area of the sections are given in Table S9 (Supplementary Materials). mRNA was amplified by one round of T7-based in vitro transcription using the MessageAmp TM II aRNA Kit (InvitrogenTM) to generate 1-2 µg antisense RNA.

RNA Sequencing and Data Preprocessing
After quality control of RNA samples, Illumina sequencing was performed by Novogene Co. Ltd. (Hong Kong, China) using the NEBnext Ultra RNA Kit for library preparation to produce 150 bp paired-end reads.

Sequence Assembly and Annotation
A de novo transcriptome assembly was performed with separated B0 and Bplus datasets using Trinity v2.6.5 [28] with default parameters which resulted in 214,060 and 254,246 contigs in B0 and in Bplus conditions, respectively. Trinity results provided a set of sequences (called "isoforms") grouped into clusters.
Quality and completeness of assemblies were determined using Transrate v1.0.3 [31], which resulted in optimized assemblies B0 and Bplus containing 165,818 and 196,360 contigs, respectively. In order to track the origin of each contig, the sequences were renamed with B0 and Bplus prefixes for each contig. In addition, B0 and Bplus files were merged into one dataset (B0&Bplus) and further processed by minimap2 v2.9 [32]. To find overlapping as well as unique contigs, a threshold of 80% identity within the aligned part of the shorter sequence was applied. The resulting FASTA files were processed by the CD-HIT-EST, v4.6.8 program (sequence identity threshold 0.90; [33,34]) to cluster highly homologous sequences and reduce redundancy. Finally, files were combined into a single file for usage in downstream applications. In particular, the final data set (117,441 contigs in total) was exploited for the detection of differentially expressed genes.

Quality Analysis and Differential Expression Analysis
For extraction of DEGs, the final file, including common and unique contigs in B0 and Bplus conditions (117,441 in total), was used for Kallisto indexing and read counting (v0.44.0; [35]). After that, log 2 -transformed counts-per-million (CPM) values were generated and "PtR" of the Trinity pipeline was employed to examine the data quality based on Pearson's correlation and principal component analysis (PCA). DEG analysis was conducted using DESeq2 [36] with an FDR of 0.05 and 0.01 [37] to identify significant DEGs.

GO Term Enrichment Analysis
Statistically enriched GO terms in DEGs were identified using the AgriGo v2.0 analysis toolkit [40]. Biological processes, molecular functions and cellular components were assessed using Triticum aestivum as a background reference and Fisher's exact test with FDR correction. Redundant GO terms were removed using REVIGO [41] with the following parameters: medium similarity (0.7); GO categories associated to: p-values; GO term sizes database: O. sativa; semantic similarity measure to use: SimRel.

RT-qPCR
Residual aRNA from the RNA-seq approach was used for cDNA synthesis. First-strand cDNA was synthesized using SuperScript III (Invitrogen) with random priming according to the manufacturer's instructions. The Power SYBR Green PCR mastermix was used to perform reactions in an ABI 7900 HT Real-Time PCR system (Applied Biosystems, Foster City, CA, USA). Data were analyzed using SDS 2.2.1 software (Applied Biosystems). Four replications were conducted for each transcript. Expression values were normalized with primers for the internal reference gene TaGAPDH (Table S1, Supplementary Materials) and calculated as an arithmetic mean of the replicates. Dissociation curves confirmed the presence of a single amplicon in each PCR reaction. Efficiencies of PCR reactions were determined using LinRegPCR software (http://www.gene-quantification.de/download.html). Values were calculated according [42] and given as relative expression (1+E)−∆Ct. All primers are listed in Table S1, Supplementary Materials.

Data Availability
Raw sequence reads can be obtained from the European Nucleotide Archive (ENA) under study accession number PRJEB39517.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1422-0067/21/20/7596/s1, Figure S1: Assessing variability of biological replicates and relationship among samples without B chromosomes (B0) and with eliminating B chromosomes (Bplus). Principal component analysis (PCA) plots displaying variability within the 6 samples along PC1 and PC2 as well as PC2 and PC3; Figure S2: GO term enrichment analysis for genes with statistically significant up-regulated changes during B chromosome elimination. An SEA tool of the AgriGo v2.0 displays hierarchical tree graphs of statistically significant GO terms (adjusted p-value < 0.05). in molecular function category. The significant terms are colored, while non-significant terms are shown as white boxes. The type of lines (solid, dashed, and dotted) reflects two, one and zero enriched terms at both ends connected by the line, respectively. The information within each box of the GO term includes GO term, adjusted p-value, description of the term, item number corresponding to the GO term in the query list and background and the total number of query list and background; Figure S3: GO term enrichment analysis for genes with statistically significant up-regulated changes during B chromosome elimination. An SEA tool of the AgriGo v2.0 displays hierarchical tree graphs of statistically significant GO terms (adjusted p-value < 0.05) in cellular compartment category. The significant terms are colored, while non-significant terms are shown as white boxes. The type of lines (solid, dashed, and dotted) reflects two, one and zero enriched terms at both ends connected by the line, respectively. Arrow in green describes negative regulation. The information within each box of the GO term includes GO term, adjusted p-value, description of the term, item number corresponding to the GO term in the query list and background and total number of query list and background; Table S1: Primers used for PCR and RT-qPCR; Table S2: Summary statistics and quality assessment for the de novo assembly based on Trinity and Transrate; Table S3: Results of Gene Set Enrichment Analysis (GSEA, AgriGo v2. analytical toolkit) using transcript isoforms up-regulated during the B chromosome elimination in Ae. speltoides. Background reference for the enrichment analysis is the wheat genome; Table S4: Transcript isoforms up-regulated during the B chromosome elimination in Ae. speltoides and representing the enriched term "microtubule-based process" based on AgriGo v2. Background reference for the enrichment analysis is wheat genome; Table S5: Transcript isoforms up-regulated during the B chromosome elimination in Ae. speltoides and representing the enriched term "chromosomal region" based on AgriGo v2. Background reference for the enrichment analysis is wheat genome; Table S6: Results of Gene Set Enrichment Analysis (GSEA, AgriGo v2. analytical toolkit) using transcript isoforms down-regulated during the B chromosome elimination in Ae. speltoides. Background reference for the enrichment analysis is wheat genome; Table S7: Transcript isoforms down-regulated during the B chromosome elimination and representing the enriched term "cellular response to stress" based on AgriGo v2. Background reference for the enrichment analysis is wheat genome; Table S8: Results of differentially expression analysis in embryos of Ae. speltoides based on DESeq2. Annotated transcript isoforms uniquely expressed during the B chromosome elimination are shown; Table S9: Tissue sections isolated by laser capture microdissection.

Conflicts 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. The authors declare no conflict of interest.