Sage Journals: Discover world-class research

Abstract

Ribosome biogenesis involves a large inventory of proteinaceous and RNA cofactors. More than 250 ribosome biogenesis factors (RBFs) have been described in yeast. These factors are involved in multiple aspects like rRNA processing, folding, and modification as well as in ribosomal protein (RP) assembly. Considering the importance of RBFs for particular developmental processes, we examined the complexity of RBF and RP (co-)orthologs by bioinformatic assignment in 14 different plant species and expression profiling in the model crop Solanum lycopersicum. Assigning (co-)orthologs to each RBF revealed that at least 25% of all predicted RBFs are encoded by more than one gene. At first we realized that the occurrence of multiple RBF co-orthologs is not globally correlated to the existence of multiple RP co-orthologs. The transcript abundance of genes coding for predicted RBFs and RPs in leaves and anthers of S. lycopersicum was determined by next generation sequencing (NGS). In combination with existing expression profiles, we can conclude that co-orthologs of RBFs by large account for a preferential function in different tissue or at distinct developmental stages. This notion is supported by the differential expression of selected RBFs during male gametophyte development. In addition, co-regulated clusters of RBF and RP coding genes have been observed. The relevance of these results is discussed.

Keywords

next generation sequencing orthologous prediction ribosome biogenesis MACE qRT-PCR tomato

Introduction

Ribosome biogenesis is a complex biological process during which pre-ribosome particles are processed to mature functional ribosome subunits. The process is best described for Saccharomyces cerevisiae.^1–3 Initially, the ribosomal RNA transcript (pre-rRNA) containing three of the four mature rRNAs is transcribed in the nucleolus priming the assembly of 90S preribosomes. Subsequent cleavage of the 35S pre-rRNA leads to the separation of the small and the large ribosomal subunit maturation. This maturation of the ribosomal subunits occurs largely in the nucleus, and only the final rRNA processing steps take place in the cytoplasm.^3,4 Here, both ribosomal subunits undergo further modifications and quality control events before final assembly to translational active ribosomes.⁵

The maturation of the ribosomal particles requires a large number of RNA and proteinaceous molecules. Small nucleolar RNAs (snoRNAs) base pairing with the pre-rRNA are required for ribosome biogenesis. Seventy-five snoRNAs have been identified in yeast as part of box C/D or box H/ACA snoRNPs,⁶ which catalyze the methylation and pseudouridylation of the rRNA, respectively.⁷ Besides snoRNAs, more than 250 proteins have been identified and assigned as ribosome biogenesis factors (RBFs) in yeast. These factors belong to protein families like RNA-helicases, GTPases, ATPases, RNA-binding proteins, endo- and exonucleases.^2,3

However, information on plant ribosome biogenesis in general as well as on plant snoRNAs is rather sparse. More than 100 binding sites of snoRNAs on plant pre-rRNA have been predicted (http://bioinf.scri.ac.uk/cgi-bin/plant_snorna/ home),^8,9 but their functional relevance has not been experimentally approached yet. Similarly, not much is known about plant RBFs. Recently, two bioinformatic studies have provided the first insights into the putative inventory of plant RBFs, one focusing on the family of RNA helicases in Zea mays, Glycine max, Arabidopsis thaliana, and Oryza sativa¹⁰ and the other describing the evolutionary development of the RBF set known in yeast.¹¹ The latter used 255 yeast RBFs to predict orthologs in plants including A. thaliana and Solanum lycopersicum.¹¹ This study realized a significant, but not complete overlap of about 75% between the RBF inventory in yeast and plants. However, the work focused on the discovery of the RBF family,¹¹ but did not evaluate the complexity caused by existence of co-orthologous genes in single species. The latter is important to explore whether plants have tissue- or developmental stage-specific RBFs or whether RBFs of a certain protein family might have taken on the function of factors not found in the plant genomes. For example, one of the two predicted A. thaliana orthologs of yeast nucleolin Nsr1 gene is abundantly expressed in all plant organs, while the second exhibits higher expression in flower buds.^12,13 Similarly, the two A. thaliana orthologs of yeast Lsg1 show a differential functional activity.¹⁴

We analyzed the complexity of RBF families and ribosomal proteins (RPs) in 14 plants. We further examined the expression of the predicted RBFs and RPs in tomato, using available RNA-seq data in combination with quantitative next generation sequencing (NGS) by Massive Analysis of 3'-cDNA Ends (MACE)¹⁵ on leaves and anthers, considering the importance of RBFs for sexual reproduction.^16–18 Additionally, the transcript profile of selected RBFs and RPs was investigated and compared by co-expression analysis. This allowed us to examine whether different (co-)orthologs exhibit a tissue- or developmental stage-specific expression, considering that several genetic studies have shown that mutations in specific genes involved in ribosome biogenesis are associated with developmental defects.¹⁹ We realized that most of the RBFs are expressed in different tissues, while for a subset of RBFs a tissue-specific expression can be concluded.

Materials and Methods

Ortholog Search and Clustering

The set of 255 yeast RBFs¹¹ and of 129 ribosomal proteins of the large subunit (RPL) and 89 ribosomal proteins of the small subunit (RPS) sequences from A. thaliana (http://ribosome.med.miyazakiu.ac.jp) ²⁰ were used for co-ortholog assignment followed by inspection of the domain architecture as described.²¹ For the ortholog search via OrthoMCL and InParanoid, we used 14 plant species including eudicots (A. thaliana, G. max, Lotus japonicus, Medicago truncatula, Populus trichocarpa, Solanum tuberosum, S. lycopersicum, Vitis vinifera), monocots (Brachy-podium distachyon, O. sativa, Sorghum bicolor, Z. mays), the Bryopsida Physcomitrella patens, the Chlorophyceae Chlamydomonas reinhardtii, and the fungi S. cerevisiae as outgroup and bait. We searched for RBFs in the orthologous groups for all plants (Supplementary Table 1) and additional information restricted to tomato and A. thaliana is provided (Supplementary Tables 2 and 3). The number of RPL and RPS orthologs is deposited in Supplementary Table 4.

Plant Material, RNA Extraction, cDNA Synthesis, and MACE

Total RNA was extracted from young leaves and anthers of eight-week-old tomato plants (cv Moneymaker) grown in a 16/8 hours day/night cycle (24-20 °C) in the greenhouse. Anthers were collected from flower buds of all developmental stages and pooled. Pollen cells from different developmental stages, namely tetrads, post-meiotic microspores, and mature pollen grains from open flowers, were isolated according to existing protocols.²² Total RNA was extracted using the NucleoSpin® RNA isolation kit (Macherey and Nagel) following manufacturer's instructions. One microgram of RNA was reverse-transcribed using Revert AidTM RvTranscriptase (Fermentas) and cDNA synthesis was performed using Oligo dT12-18 primer. The MACE libraries were prepared and sequenced by GenXPro GmbH.^15,23

NGS Analysis Pipeline

The analysis of the NGS data was performed with the flexible pipeline CRACPipe (Supplementary Method; Supplementary Figs. 1-3), which allows application to any species with available sequence information (genome or transcriptome) by screening NGS data with the BLAST-like Alignment Tool BLAT²⁴ for adapter sequence removal, quality score determination by Fast Quality Control FastQC,²⁵ and mapping with Sequence Search and Alignment by Hashing Algorithm 2 SSAHA2.²⁶ The pipeline parses genomic data provided in GenBank format to internally construct an annotated genome divided in Watson and Crick strand²⁷ containing information of the position, name and type of each feature. Supported feature types are mRNAs and non-protein coding RNAs, like tRNAs, rRNAs, snoR- NAs, snRNAs, miRNAs and other ncRNAs.²⁸ As output, mochi files are created for visualizing the results in the genome browser MochiView.²⁹ The pipeline can be used via web page front end (https://www.uni-frankfurt.de/fb/fb15/institute/inst-3-mol-biowiss/AK-Schleiff/cracpipe/). For details see Supplementary Method.

Figure 1

NGS and MACE comparison. Comparison of MACE and NGS data for leaf tissue. The line indicates the Spearman correlation for both leaf expression datasets using RNA-seq and MACE. The inset gives the number of genes assigned in the genome, detected by MACE and detected by NGS, and depicts the higher coverage achieved by MACE.

Figure 2

Prediction of RBF coding genes. (A) Phylogenetic relation of the 14 plant species used for the analysis. (B) Discovery rate of RBFs based on 255 yeast RBFs¹¹ in the genomes of 14 plant species in percentage of the number found in yeast given for the different families. The different colors indicate the phylogenetic relation (green: green algae; black: moss; light green: monocots; dark green: eudicots). (C) The percentage of plant RBFs encoded by the indicated number of co-orthologs. Atha: A. thaliana; Bdis: B. distachyon; crei: C. reinhardtii; Gmax: G. max; Ljap: L. japonicus; Mtru: M. truncatula; osat: O. sativa; Ppat: P. patens; Ptri: P. trichocarpa; Sbic: S. bicolor; Slyc: S. lycopersicum; Stub: S. tuberosum; Vvin: V. vinifera; Zmay: Z. mays.

Figure 3

Prediction of RP coding genes. The distribution of the number of co-orthologs identified for RPSs (A) and RPLs (B) in the 14 analyzed plants as shown in Figure 1

NGS Results

The MACE results are available at the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/; series: GSE56610). Twelve million reads (100 bp) with more than 14 nucleotides after linker removal have been obtained for both samples (Table 1). The cut-off resulted in the unique assignment of >97% of all reads. Assignment of reads to >22,000 genes and ~30,000 regions (non-feature) without current annotation has been achieved.

Table 1

Experimental overview for MACE analysis.

	LEAF	ANTHER
Total read number	11.762.204	12.687.296
Total mapped reads	11.534.512	12.411.606
Percentage of unique mapped reads	98%	97%
Hits on watson strand	6.235.294	6.386.089
Hits on crick strand	5.299.286	6.043.673
Genes^a	22.894	24.539
Non-feature	29.666	30.540

Analysis of properties of MACE libraries (column 1) in leaf (column 2) or anther (column 3) sample.

In total, 29.345 genes have been assigned.

MACE is a digital gene expression method in which each read represents a cDNA molecule.^15,30 By this, transcripts per million (TPM) is equivalent to reads per million, which is calculated by dividing the number of all reads mapped to an individual gene by all mapped reads normalized to a million reads. Furthermore, the cut-off for detection of a gene is the detection of a unique read for that gene.^15,30 In turn, a gene is annotated as not expressed if no unique read was mapped. The latter was justified by the estimation of the coverage obtained (Supplementary Fig. 4).^31,32 The analysis revealed that the call of “not expressed” has a reliability of 95% for the leaf dataset and of 99% for the dataset for anthers.

Figure 4

Prediction of co-orthologs to RBF coding genes. (A) Comparison of the total number of sequences in yeast (black bar), A. thaliana (red bar), and tomato (green bar). (B) The percentage of RBFs assigned to the specific complexes encoded by the number of co-orthologs shown on the right.

We used NoiSeq³¹ to simulate five replicates by a multinomial distribution and mapping probabilities were approximated from the available MACE library, where each feature had the chance to appear assigning a value bigger than 0. The five simulated replicates were used for differential expression analysis by the assumption that a high probability in NoiSeq correlates with a high probability that an observed difference in expression is significant.

Existing RNA-seq libraries containing expression values for the majority of tomato genes as well as Affymetrix microarray data were used to verify the quality of our MACE data. First, genes with low expression variance irrespective of developmental stages or experimental conditions were identified. We postulate that such genes are equally expressed under most experimental conditions and thus are best suited to judge the comparability of the two datasets generated. The variance of expression values for individual genes identified in Affymetrix studies of different tissues and developmental stages (91 samples; GEO, ArrayExpress) and after stress application (11 samples; GEO, TFGD; Supplementary Table 5) was determined. The 20 genes with the lowest variance were selected (Supplementary Table 6). The expression of these genes in our MACE datasets for leaves and anthers shows a linear relation of 0.97 (Supplementary Fig. 5). This confirms that the extracted TPM values for the two datasets are comparable. The R² value for linear regression shown in Supplementary Figure 5 is 0.91.

Figure 5

RBF co-orthologs with alterations in domain architecture. Alignments given in Supplemental Alignment 1 are shown as bar diagrams including the Pfam domains assigned. See Supplementary Table 3 for visualization of differences. Italics indicates the genes excluded from further analysis; italics and underlined indicates the gene found in unigene database, but not yet represented by a Solyc ID.

Secondly, the RNA-seq data for leaves and flowers were directly compared to our MACE results. While we detected expression of a higher number of genes in leaves compared to the published NGS dataset,³³ the comparison between our MACE data and NGS data (GEO: GSE33507) revealed a high correlation (r_s = 0.74, P = 2 X 10^-7) confirming the validity of the data used (Fig. 1).

cluster Analysis of Expression Profiles of Genes Encoding RPs and RBFs

The expression values of the RPs and RBFs from NGS experiments of tomato plants cv. Moneymaker (Supplementary Tables 7 and 8)³³ were used for cluster analysis. RPs and RBFs were independently clustered according to their expression profiles in different tomato tissues. For k-means clustering, Multi Experiment Viewer (MeV; http://www.tm4.org/mev.html) ³⁴ was used to determine six different clusters of each set of RPs and RBFs. The number of clusters required was determined as previously described (Supplementary Fig. 6).³⁵ For k-means clustering, 1000 iterations have been performed and the Pearson correlation was used as distance.³⁶ The accession numbers of the co-orthologs in the individual clusters are listed in Supplementary Table 9.

Figure 6

The overall expression profile of RBF and RP (co-)orthologs in tomato. (A) Tissues used for NGS RNA-seq analysis (left) and the number of RBFs for which expression is detected in a least one tissue (beside) or in the individual tissue (below). On the right, the tissues analyzed by MACE in this study are indicated. (B) Relation of the TPM expression value in leaves and anthers for all genes (black circles), for all RBF genes (red circles), and for all RP genes (yellow star). Gray line indicates identical expression in leaves and anthers, the long dashed gray line expression with two-fold change, and the short dashed gray line expression with four-fold change (45% genes with less than two-fold change; 64% with less than four-fold change of expression). Indicated are RBF genes not expressed (top), expressed only in leaves (left) or anther (right). Inset on the right shows the distribution of the expression difference between leaf and anther. The gray section indicates a pool of genes with significantly higher expression in anthers, and the red line shows the least square fit analysis to a Gaussian equation.

Quantitative Real-Time PCR

Gene expression of selected genes was determined using two biological replicas by real-time PCR on a Stratagene Mx3000P cycler (Agilent Technologies). The qPCR reaction (10 μL) consisted of gene-specific primers (Supplementary Table 10), PerfeCTa® SYBR® Green FastMix Low ROX™ (Quanta Biosciences), and the template. The cycling conditions were 95 °C/3 minutes followed by 95 °C/15 seconds, 60 °C/30 seconds, and 72 °C/ 30 seconds for 40 cycles. Primers were designed using PRIMER3 (www-genome.wi.mit.edu/cgi-bin/primer/primer3.cgi/). Data were analyzed using the 2^-ΔΔCt method³⁷ and presented as relative levels of gene expression, using ubiquitin (Solyc07g064130) and EF1α (Solyc06 g005060) genes as internal standards. Expression values are the mean of two biological replicates.

Results

The Global Assignment of Co-Orthologs to RBFs and RPs in Plant Genomes

We selected 14 plant genomes including the green algae C. reinhardtii, the moss P. patens, four monocots and eight dicots including the agriculturally relevant S. lycopersicum (tomato), and the model plant A. thaliana (Fig. 2A). Previously, orthologs to about 170 of255 RBFs in yeast have been assigned in these 14 plant species and the highest number was found in the genomes of S. lycopersicum and A. thaliana with 173 and 174 genes, respectively (Supplementary Tables 1-3).¹¹ Remarkably, only 116 RBFs have been found in all 14 plant genomes. The discovery rate of RBFs associated with different complexes is comparable. However, orthologs to RBFs of the exosome, the mitochondrial RNA processing (MRP) complex, and Trf4/Air2/Mtr4p polyadenylation (TRAMP) complex are under-represented in all plants compared to yeast (Fig. 2B). At least two co-orthologs are present in most of the species for approximately 20-40% of RBFs (Fig. 2C), while the same holds true for more than 50% of the RBFs in S. bicolor, G. max, and B. distachyon (Fig. 2C).

In parallel to the analysis of RBFs, we inspected the diversity of RP-coding genes based on the sequences assigned in A. thaliana.²⁰ In most of the 14 plant species, the majority of the RPs are encoded by two or three co-orthologs (Fig. 3; Supplementary Table 4). Only in the algae C. reinhardtii they are encoded by one gene, while in P. patens, P. trichocarpa, Z. mays, and G. max, the majority is encoded by three or more genes. Comparing the number of co-orthologs found for RBFs (Fig. 2) and RPs (Fig. 3), we did not observe a correlation between the occurrence of co-orthologs of RBFs and RPs.

Sequence Analysis of Identified RBF Co-Orthologs in Tomato

In S. lycopersicum and A. thaliana, more than one co-ortholog exists for about 24% and 27% of all RBFs, respectively. The total number of sequences assigned to the complexes is by large the same or even larger compared to yeast (Fig. 4A). In total, 249/260 sequences have been assigned as (co)-orthologs to RBFs in S. lycopersicum/A. thaliana, respectively. The extent of the occurrence of multiple orthologs varied among different complexes (Fig. 4B). With the exception of MRP and TRAMP in both species, and the exosome in S. lycopersicum, co-orthologs are found for RBFs of all complexes in all plants (Fig. 4B).

However, ortholog assignment is not necessarily proof for a function of the encoded proteins in the same process as the ancestral protein.³⁸ To provide additional evidence for comparable functions between yeast RBFs and the proteins encoded by the S. lycopersicum genes, we categorized the predicted tomato RBF sequences based on (i) euKaryotic Orthologous Groups (KOG), (ii) their functional description, as well as (iii) Pfam domain presence and order (Supplementary Table 2). For six RBFs (DBP2, DRS2, FAP7, MDN1, YAR1 and YVH1), the co-orthologous sequences differ in at least one of the analyzed categories. Manual inspection of the corresponding alignments (Fig. 5; Supplementary Alignment S1 Supplementary Table 3) showed that one co-ortholog of DRS2 (Solyc01g011100) has a shorter N-terminal domain; one FAP7 co-ortholog (Solyc06g075310) lacks the C-terminal region: and the co-orthologs of MDN1 (Solyc04g072270) and YAR1 (Solyc04g008580) represent only a short fragment. The latter holds true for two of the co-orthologs assigned to YVH1 as wel (Solyc12g021180, Solyc12g021190; Fig. 5). Solyc12g021180 encodes for a protein with a partial dual specificity phosphatase catalytic domain (DSPc) and a shorter C-terminal domain compared to the full-length protein, Solyc00g185750. However, Solyc12g021180 and Solyc12g021190 together match the unigene SGN-U584680 indicating that the two Solyc IDs correspond to a single gene (putative Solyc12g021185), which encodes for a protein with a full-length DSPc domain and C-terminus matching the protein model of Solyc00g185750 (Fig. 5). Collectively, this evaluation reveals that at least all sequences of DRS2 and FAP7 can be considered as co-orthologs most likely with similar function, while four sequences are probably non-functional RBFs. In addition, YVH1 might be represented by two co-orthologs considering Solyc12g021185 as functional gene. Summarizing, we assigned 245 putative functional RBF (co-)orthologs in the genome of S. lycopersicum to the 255 RBFs from yeast.

Table 2

Genes expressed highest in leaves and anthers.

	GENE	PROTEIN FUNCTION	RANK LEAVES	RANK ANTHER
Highest TPM values in leaves	Solyc02g085950	Ribulose bisphosphate carboxylase small chain 3B	1	539
	Solyc07g066310	psbR-homologue	2	69
	Solyc09g010800	metallothionein LEMT1	3	9
	Solyc07g047850	chlorophyll a-b binding protein 4	4	33
	Solyc10g006230	chlorophyll a-b binding protein 7	5	116
	Solyc03g034220	Ribulose bisphosphate carboxylase small chain 2A	6	21
	Solyc07g064160	chloroplastic thiamine thiazole synthase	7	514
	Solyc02g071030	chlorophyll a-b binding protein 1B	8	113
	Solyc06g063370	chlorophyll a-b binding protein 1A	9	96
	Solyc03g005760	chlorophyll a-b binding protein 3C	10	129
Highest TPM values in anther	Solyc01g009590	carotenoid cleavage dioxygenase 1	1113	1
	Solyc01g086830	seed storage/lipid transfer protein family protein	1454	2
	Solyc07g007250	metallocarboxypeptidase inhibitor	1511	3
	Solyc09g008670	threonine deaminase	491	4
	Solyc06g064480	protein TAP1	2196	5
	Solyc07g006380	flower-specific gamma-purothionin protein	2144	6
	Solyc01g099540	Anther specific NMT19-like	4224	7
	Solyc01g006400	Hop-interacting cysteine-rich extensin-like protein-4	117	8
	Solyc07g043420	2-oxoglutarate-dependent dioxygenase	16	10

Notes: The TPM expression values are listed in Supplementary Table 4. The upper part contains the 10 most abundant genes in leaves and the lower part the 10 most abundant genes in anthers. The rank of the expressed genes is assigned concerning the descending TPM expression values obtained by MACE analysis.

Expression of S. Lycopersicum RBF and RP (co-) Orthologs

The presence of multiple co-orthologs in plant genomes might indicate redundant or tissue- and/or developmental stage-specific functions for some co-orthologs. Thus, we examined the expression of all predicted RBFs and RPs in an existing expression dataset derived from RNA-seq analysis of S. lycopersicum cv. Moneymaker.³³ From the 245 predicted RBF and 193 RP (co-)orthologs in tomato, 242 and 188 were expressed in at least one of the examined tissues, respectively (Supplementary Tables 7 and 8). However, transcripts for only 129 RBF and 166 RP (co-)orthologs were detected in open flowers without sepals (Fig. 6A). The low number of RBFs and RPs expressed in flowers containing reproductive tissues of Moneymaker³³ was surprising, because rapidly developing tissues are expected to require massive ribosome biogenesis, particularly as mutants of RBFs have been found to be defective in pollen development.^16,17,39 Thus, we performed NGS of RNA in leaves and anthers to re-evaluate the expression of RBFs in reproductive tissue. We focused on anthers, which are part of the male reproductive tissue and contain sporophytic tissues and gametophytic cells, including pollen corresponding to different developmental stages, ranging from meiotic to mature stage. In addition, we analyzed the transcript abundance in young expanding leaves for comparison.

First, the NGS data for leaves and anther were evaluated (GEO: GSE56610). 21,509 genes have been identified by at least one read in both tissues. Here, we use the MACE approach, for which it was discussed that occurrence of one read is sufficient for the detection of a transcript (see Methods).^15,30 Thus, in the two samples the expression of 25,924 of the 29,345 assigned genes was detected, which accounts for 88.3% of the total annotated genes. This is comparable to the detection of 83.6% of all annotated genes in a previous NGS study including different tissues and developmental stages of S. lycopersicum cv. Moneymaker.³³ 45% of all genes have similar transcript abundance (<2-fold change); while 36% show a >4-fold difference of expression in the two samples (Fig. 6B). 1,385 and 3,030 genes are expressed exclusively in leaves and anthers, respectively.

The genes with the highest expression in leaves (normalized to TPM) predominantly encode for proteins involved in photosynthesis (Table 2; Supplementary Table 11). Many of these genes show a high expression in anthers as well, which documents that anthers are a mixture of reproductive and sporophytic tissues. The genes with the highest expression in anthers are rather moderately expressed in leaves and encode for metabolic or regulatory proteins (Table 2; Supplementary Table 11). Inspecting the profile of differences in expression, one realizes that besides the 3,030 genes exclusively expressed in anthers, 1,455 genes have a 10-fold higher abundance in anthers than in leaves indicating strong preferential expression (Fig. 6B, inset; Supplementary Information 1).

We examined the expression of all genes coding for the RP (co-)orthologs and found that with the exception of three genes, all RPs are expressed in at least one of the two tissues (Fig. 6B; Supplementary Tables 7 and 8). RP coding genes are in general highly expressed in both tissues. RBF genes which showed tissue-specific expression had a low number of reads (median of 0.09TPM; Supplementary Table 7), indicating a low transcript abundance.

The three RBF coding genes for which no transcript was detected both in leaves and anthers are Solyc00g185750 coding for a co-ortholog of YVH1, Solyc00g095450 coding for YAR1, and Solyc05g018780 coding for RRB1 (Fig. 6B). Transcripts for these genes were not detected in the samples analyzed before by RNAseq³³ or in a global analysis of the heat-induced transcriptome.³⁶ This raises the question whether these genes are only induced under very specific stress conditions or are particularly expressed in a not yet analyzed tissue or cell type. Worth mentioning, the expressed co-orthologs of YAR1 code only for fragments of the protein (Fig. 5), and thus, expression of this RBF could not be demonstrated. In contrast, the transcripts of the newly assigned co-ortholog of YVH1 and the two additional co-orthologs of RRB1 are detected in anther and leaf (Supplementary Table 7). In summary, we detected transcripts for 239 (98%) assigned RBFs and for all RPs in anthers. The transcript abundance of most RP genes is slightly higher in leaves, while the transcript abundance of RBFs is globally comparable in both tissues (Fig. 6B).

The Expression Profile of Predicted RPs and RBFs in Tomato Leaves and Anthers

Nine RBF coding genes and 120 RP coding genes are either 2-fold higher or exclusively expressed in leaves. The RBF coding genes encode for factors possessing multiple co-orthologs with the exception of AIR2 (Table 3). Forty RBF genes and only one RP gene are more abundant (anther to leaf ratio >2) or exclusively expressed in anthers (Table 3). Nineteen genes represent the full set of coorthologs coding for a certain RBF. Among the genes preferentially expressed in one of the tissues we found four factors represented by multiple co-orthologs, where one is highly expressed in leaves and another in anthers (NAP1, PTC3, TIF6, UBC9). Indeed, the NAP1 co-ortholog, which is highly expressed in anthers, was not detected in roots and leaves in a previous study (Supplementary Table 7).³³ For the UBC9 co-ortholog, a high expression in roots, stem, and ripe fruits was reported, while the expression in leaves was found to be the lowest,³³ which is consistent with our results (Supplementary Table 7).

Table 3

RBF genes differentially expressed in leaves and anthers.

	GENE	CO	ACC. NUMBER	Log₂(TPM_ANTHER/TPM_LEAF)
Expressed higher in leaves	AIR2	1	Solyc05g007090	-1,93
	HRR25	11	Solyc12g007280	-Inf
	HRR25	11	Solyc08g066870	-2,57
	NAP1	5	Solyc08g063000	-1,34
	NMD3	2	Solyc11g011780	-1,50
	PTC3	3	Solyc08g074230	-Inf
	REI1	2	Solyc08g006470	-3,19
	SMT3	4	Solyc09g059970	-1,84
	TIF6	3	Solyc12g010210	-Inf
	UBC9	5	Solyc02g093110	-1,18
Expressed higher in anthers	ASC1	3	Solyc12g040510	1,03
			Solyc06g069010	1,28
			Solyc03g119040	1,38
	BUD23	1	Solyc06g072250	1,06
	CBF5	1	Solyc02g081810	1,80
	CIA1	2	Solyc09g009410	1,31
	DIM1	1	Solyc01g100430	1,32
	EBP2	2	Solyc01g006090	2,08
	ERB1	1	Solyc01g103740	1,43
	ERP2	5	Solyc12g015740	1,17
			Solyc09g082770	1,26
			Solyc12g006090	+Inf
	FAL1	6	Solyc07g040750	3,65
	FAL1	6	Solyc06g062800	+Inf
	FPR3	1	Solyc01g091340	1,73
	FPR4	1	Solyc03g007170	1,59
	GAR1	1	Solyc10g078540	1,01
	GRC3	1	Solyc03g113610	1,76
	JJJ1	2	Solyc03g115140	1,13
	JJJ1	2	Solyc03g115120	4,01
	KRI1	1	Solyc07g064350	1,02
	LSG1	1	Solyc11g071910	5,91
	NAP1	5	Solyc06g062690	5,70
	N0B1	1	Solyc08g007990	1,18
	NOP4	1	Solyc01g109990	1,11
	NOP58	2	Solyc09g083080	1,12
	P0P1	1	Solyc11g007800	1,14
	PTC3	3	Solyc12g042570	1,41
	PWP2	1	Solyc03g120080	3,22
	RLI1	2	Solyc08g075360	+Inf
	RNA1	2	Solyc01g079680	1,10
	RNA1	2	Solyc09g065190	1,25
	RRB1	4	Solyc05g025630	1,82
	RRP45	1	Solyc05g047420	1,89
	SNU13	2	Solyc09g092080	1,03
	TIF6	3	Solyc12g015790	4,96
	UBC9	5	Solyc12g088680	1,82
	UTP22	1	Solyc02g085230	1,68
	YTM1	2	Solyc01g111030	+Inf
	YVH1	3	Solyc12g021180	1,48

To further support the observed differential expression of RBF coding genes based on MACE results, we randomly selected genes and analyzed their expression by qRT-PCR on mRNA isolated from anthers and leaves (Table 4). We included two factors represented by multiple co-orthologs into this analysis, namely RIL1 and ASC1. We obtained a qRT-PCR signal for each of the five genes coding for one of the two factors. This confirms the existence of a transcript for the multiple co-orthologs. For quantification, the expression values from qRT-PCR were normalized to either Solyc07g064130 coding for ubiquitin or Solyc06g005060 coding for EF1α. We observed by large the same ratios of expression levels in the two tissues as obtained by MACE analysis. Thus, we confirmed the leaf-specific expression of AIR2, which is only encoded by a single gene. Furthermore, we observed the anther-specific expression of one RIL1 coding gene. Remarkably, by qRT-PCR we also observed a higher expression for an additional gene coding for ASC1 in anthers (Table 4).

Table 4

RBF expression comparison of MACE and qRT-PCR analysis in anther and leaf.

NAME	GENE	qRT-PCR (NORMALIZED)				TPM-EXPRESSION VALUE
		NORM	LEAF	ANTHER	A/L	LEAF	ANTHER	A/L
RIL1	Solyc08g075360.1	EF1α	(6.4 ∓ 0.9)*10^-3	0.15 ∓ 0.05	32 ∓ 7	0.0	0.44	∞
	Solyc08g075360.1	Ubi	(3.06 ∓ 0.03)*10^-4	(1.2 ∓ 0.2)*10^-2	32 ∓ 7	0.0	0.44	∞
	Solyc07g008340.2	EF1α	(6.2 ∓ 0.9)*10^-2	(9 ∓ 2)*10^-2	2 ∓ 1	47.9	69.0	1.44
	Solyc07g008340.2	Ubi	(2.8 ∓ 0.1)*10^-3	(8.2 ∓ 0.6)*10^-3	2 ∓ 1	47.9	69.0	1.44
SSF1	Solyc01g096800.2	EF1α	0.30 ∓ 0.04	0.19 ∓ 0.08	0.8 ∓ 0.1	13.2	14.9	1.13
SSF1	Solyc01g096800.2	Ubi	(1.4 ∓ 0.4)*10^-2	(1.27 ∓ 0.01)*10^-2	0.8 ∓ 0.1	13.2	14.9	1.13
LSG1	Solyc11g071910.1	EF1α	0.2 ∓ 0.1	0.8 ∓ 0.2	8 ∓ 2	3.9	233.6	59.9
LSG1	Solyc11g071910.1	Ubi	(1.3 ∓ 0.2)*10^-2	0.14 ∓ 0.05	8 ∓ 2	3.9	233.6	59.9
PWP2	Solyc03g120080.2	EF1α	0.3 ∓ 0.1	1.0 ∓ 0.2	10 ∓ 4	0.09	0.8	9.0
PWP2	Solyc03g120080.2	Ubi	(1.1 ∓ 0.2)*10^-2	0.19 ∓ 0.06	10 ∓ 4	0.09	0.8	9.0
ASC1	Solyc12g040510.1	EF1α	0.3 ∓ 0.1	0.26 ∓ 0.06	0.9 ∓ 0.1	2.6	5.3	2.0
	Solyc12g040510.1	Ubi	(2 ∓ 1)*10^-2	(2 ∓ 1)*10^-2	0.9 ∓ 0.1	2.6	5.3	2.0
	Solyc03g119040.2	EF1α	28.8 ∓ 0.9	19 ∓ 4	0.8 ∓ 0.1	73.6	192.2	2.6
	Solyc03g119040.2	Ubi	1.02 ∓ 0.03	0.8 ∓ 0.1	0.8 ∓ 0.1	73.6	192.2	2.6
	Solyc06g069010.2	EF1α	(3.0 ∓ 0.6)*10^-2	0.14 ∓ 0.02	4.5 ∓ 0.5	140.5	342.2	2.4
	Solyc06g069010.2	Ubi	(5 ∓ 2)*10^-3	(2.1 ∓ 0.9)*10^-2	4.5 ∓ 0.5	140.5	342.2	2.4
RRP45	Solyc05g047420.2	EF1α	(4.5 ∓ 0.3)*10^-2	0.17 ∓ 0.01	3.3 ∓ 0.7	0.13	0.48	3.7
RRP45	Solyc05g047420.2	Ubi	(3.2 ∓ 0.2)*10^-3	(8.9 ∓ 0.1)*10^-3	3.3 ∓ 0.7	0.13	0.48	3.7
KRE33	Solyc04g051670.2	EF1α	0.84 ∓ 0.04	0.15 ∓ 0.03	0.4 ∓ 0.2	11.0	15.8	1.4
KRE33	Solyc04g051670.2	Ubi	0.05 ∓ 0.01	0.03 ∓ 0.01	0.4 ∓ 0.2	11.0	15.8	1.4
RRP5	Solyc03g051900.2	EF1α	0.2 ∓ 0.1	0.15 ∓ 0.03	0.9 ∓ 0.3	2.7	4.2	1.6
RRP5	Solyc03g051900.2	Ubi	0.02 ∓ 0.01	0.03 ∓ 0.02	0.9 ∓ 0.3	2.7	4.2	1.6
N0P12	Solyc08g076660.2	EF1α	0.8 ∓ 0.2	0.4 ∓ 0.2	0.7 ∓ 0.2	16.9	24.4	1.4
N0P12	Solyc08g076660.2	Ubi	(2.9 ∓ 0.7)*10^-2	(2.3 ∓ 0.6)*10^-2	0.7 ∓ 0.2	16.9	24.4	1.4
TSR1	Solyc01g096760.2	EF1α	0.46 ∓ 0.06	0.45 ∓ 0.06	1.4 ∓ 0.5	14.1	23.4	1.7
TSR1	Solyc01g096760.2	Ubi	0.11 ∓ 0.03	0.2 ∓ 0.1	1.4 ∓ 0.5	14.1	23.4	1.7
UBA2	Solyc01g109960.2	EF1α	0.2 ∓ 0.1	0.16 ∓0.01	1.2 ∓ 0.5	17.5	27.1	1.5
UBA2	Solyc01g109960.2	Ubi	(1.1 ∓ 0.2)*10^-2	(2 ∓ 1)*10^-2	1.2 ∓ 0.5	17.5	27.1	1.5
UTP20	Solyc09g092240.2	EF1α	0.2 ∓0.1	0.2 ∓0.1	1.1 ∓ 0.3	5.8	7.8	1.3
UTP20	Solyc09g092240.2	Ubi	(1.15 ∓ 0.06)*10^-2	(1.2 ∓ 0.1)*10^-2	1.1 ∓ 0.3	5.8	7.8	1.3
AIR2	Solyc05g007090.2	EF1α	0.28 ∓ 0.09	0.011 ∓ 0.009	0.05 ∓ 0.02	9.8	2.6	0.26
AIR2	Solyc05g007090.2	Ubi	(1.37 ∓ 0.09)*10^-2	(8 ∓ 2)*10^-4	0.05 ∓ 0.02	9.8	2.6	0.26
ENP1	Solyc07g056550.2	EF1α	0.3 ∓ 0.1	0.3 ∓ 0.2	1.6 ∓ 0.7	7.1	7.0	1.0
ENP1	Solyc07g056550.2	Ubi	(1.4 ∓ 0.7)*10^-2	(3 ∓ 1)*10^-2	1.6 ∓ 0.7	7.1	7.0	1.0
DED1	Solyc03g052980.2	EF1α	0.13 ∓ 0.02	0.10 ∓ 0.02	0.9 ∓ 0.2	76.8	66.8	0.9
DED1	Solyc03g052980.2	Ubi	(7 ∓ 3)*10^-3	(7 ∓ 2)*10^-2	0.9 ∓ 0.2	76.8	66.8	0.9
REI1	Solyc12g096450.1	EF1α	0.31 ∓ 0.01	(1.2 ∓0.7)*10^-2	0.08 ∓ 0.04	14.1	8.4	0.6
REI1	Solyc12g096450.1	Ubi	(3 ∓ 1)*10^-2	(4 ∓ 2)*10^-3	0.08 ∓ 0.04	14.1	8.4	0.6

We used the information on transcript abundance established previously³³ to analyze the expression profiles of RBF and RP coding genes. The expression profiles of the genes encoding for RPs and RBFs have been assigned to six clusters (see Methods). While each cluster contains RPs or RBFs with distinct profiles, one cluster for each set (21 rp and 68 rbf) represents a collection of genes that did not match (extremely high error bars) the profiles of any of the other clusters (Fig. 7A; 21 RP genes and 68 RBF genes). Subsequently, the correlations between the profiles of RP and RBF genes have been determined excluding the above-mentioned cluster (Fig. 7B). We realized that two profiles are highly correlative (Fig. 7B, black), with a correlation value of 0.97. The genes of these two clusters show a moderate expression in roots and stem tissues; a lower expression in leaves, anthers, and flowers; and a higher but by large similar expression in all analyzed fruit samples. These genes encode for 15 RPLs, 12 RPSs, 8 90S RBFs, 7 60S RBFs, 1 40S RBF, 1 Exonuclease RBF, 1 Exosome RBF, 1 TRAMP complex RBF, and one not assigned RBF. This might point toward a specific functional relation between these RPs and RFBs, which is discussed below.

Figure 7

Clustering of RP and RBF genes based on their expression. (A) Mean of the expression values of clustered genes normalized to the TPM determined for MACE leaf sample. Clusters were generated by k-means analysis. Error bars indicate the standard deviation. The clusters are ordered according to the number of factors present (first value) and the median of the cluster for the leaves (second value) is indicated. The two clusters on the top represent the set of genes that could not be related to a specific profile. (B) The Spearman correlation between the cluster profiles of RPs and RBFs was calculated, and the highest value for each column and each row is shown. Only values with strong correlation (above 0.5) are highlighted.

Expression of Predicted RBFs in Pollen

The role of ribosomes in development in eukaryotes has been shown in genetic studies with mutations in RPs, ribosome assembly, and biogenesis factors.¹⁹ Furthermore, the importance of ribosome biogenesis in physiological gametophyte development has been documented.^16,17,39 Our global analysis has focused on anthers of different developmental stages containing sporophytic tissues and gametophytic cells. Thus, the sample includes pollen ranging from meiotic to mature stage. The progression from proliferating microspores to terminally differentiated pollen grains is characterized by stage-specific gene transcriptional activation and repression as shown for A. thaliana pollen.⁴⁰ Therefore, we examined the levels of the genes selected for qRT-PCR analysis in developing tomato pollen (Table 5). Indeed, we identified differential expression for 12 of the analyzed genes, which accounts for 64% of the studied genes.

Table 5

Comparison of NGS and qRT-PCR analysis in pollen.

NAME	GENE	qRT-PCR (NORMALIZED)				NORM
		NORM	TETRADS	POST-MITOTIC	MATURE	PM/T	M/T
RIL1	Solyc08g075360.1	EF1α	(6.3 ∓ 0.8)*10^-2	(8 ∓ 1)*10^-2	(1.2 ∓ 0.3)*10^-2	1.26 ∓ 0.02	0.19 ∓ 0.01
	Solyc08g075360.1	Ubi	0.15 ∓ 0.05	0.19 ∓ 0.01	(2.7 ∓ 0.2)*10^-2	1.26 ∓ 0.02	0.19 ∓ 0.01
	Solyc07g008340.2	EF1α	(8.3 ∓ 0.4)*10^-3	(1.9 ∓ 0.3)*10^-2	(2.2 ∓ 0.2)*10^-2	2.4 ∓ 0.3	2.2 ∓ 0.5
	Solyc07g008340.2	Ubi	(1.9 ∓ 0.3)*10^-2	(4 ∓ 1)*10^-2	(3.4 ∓ 0.6)*10^-2	2.4 ∓ 0.3	2.2 ∓ 0.5
SSF1	Solyc01g096800.2	EF1α	(9 ∓ 1)*10^-2	0.14 ∓ 0.05	0.16 ∓ 0.04	1.6 ∓ 0.3	1.7 ∓ 0.1
SSF1	Solyc01g096800.2	Ubi	0.16 ∓ 0.03	0.28 ∓ 0.09	0.27 ∓ 0.03	1.6 ∓ 0.3	1.7 ∓ 0.1
LSG1	Solyc11g071910.1	EF1α	0.5 ∓ 0.2	0.4 ∓ 0.1	0.41 ∓ 0.06	0.8 ∓ 0.1	0.6 ∓ 0.1
LSG1	Solyc11g071910.1	Ubi	1.2 ∓ 0.4	0.9 ∓ 0.2	0.71 ∓ 0.03	0.8 ∓ 0.1	0.6 ∓ 0.1
PWP2	Solyc03g120080.2	EF1α	0.9 ∓0.3	1.3 ∓0.6	0.7 ∓ 0.3	1.3 ∓ 0.1	0.7 ∓ 0.1
PWP2	Solyc03g120080.2	Ubi	3 ∓ 1	4.0 ∓ 1	2.2 ∓0.5	1.3 ∓ 0.1	0.7 ∓ 0.1
ASC1	Solyc12g040510.1	EF1α	0.14 ∓0.03	0.57 ∓0.04	0.70 ∓ 0.01	0.33 ∓ 0.07	0.44 ∓0.05
	Solyc12g040510.1	Ubi	0.7 ∓0.2	0.2 ∓ 0.1	0.28 ∓ 0.05	0.33 ∓ 0.07	0.44 ∓0.05
	Solyc03g119040.2	EF1α	2.08 ∓0.07	3.9∓0.3	1.57 ∓ 0.04	2.0 ∓ 0.2	0.7 ∓ 0.1
	Solyc03g119040.2	Ubi	3.8 ∓ 0.5	8 ∓1	2.3 ∓0.2	2.0 ∓ 0.2	0.7 ∓ 0.1
	Solyc06g069010.2	EF1α	(2.8 ∓ 0.1)*10^-2	(7 ∓ 3)*10^-2	(6 ∓ 1)*10^-2	2.8 ∓ 0.2	2.1 ∓ 0.1
	Solyc06g069010.2	Ubi	0.117 ∓ 0.002	0.3 ∓ 0.1	0.24 ∓ 0.03	2.8 ∓ 0.2	2.1 ∓ 0.1
RRP45	Solyc05g047420.2	EF1α	(1.2 ∓ 0.4)*10^-2	(1.6 ∓ 0.5)*10^-2	(5.3 ∓ 0.7)*10^-2	1.3 ∓ 0.2	4.1 ∓ 0.7
RRP45	Solyc05g047420.2	Ubi	(5 ∓ 2)*10^-2	(6.2 ∓ 0.4)*10^-2	(19 ∓ 2)*10^-2	1.3 ∓ 0.2	4.1 ∓ 0.7
KRE33	Solyc04g051670.2	EF1α	(2 ∓ 1)*10^-2	(1.7 ∓ 0.7)*10^-2	(1.8 ∓ 0.3)*10^-2	0.7 ∓ 0.1	0.7 ∓ 0.1
KRE33	Solyc04g051670.2	Ubi	0.11 ∓ 0.05	(9 ∓2)*10^-2	(7.5 ∓ 0.2)*10^-2	0.7 ∓ 0.1	0.7 ∓ 0.1
RRP5	Solyc03g051900.2	EF1α	(3 ∓1)*10^-2	(6 ∓ 3)*10^-2	0.17 ∓ 0.03	1.6 ∓ 0.5	5.4 ∓ 0.9
RRP5	Solyc03g051900.2	Ubi	0.13 ∓ 0.06	0.18 ∓ 0.07	0.64 ∓ 0.02	1.6 ∓ 0.5	5.4 ∓ 0.9
NOP12	Solyc08g076660.2	EF1α	(9 ∓ 2)*10^-2	(5 ∓0.2)*10^-2	(1.2 ∓ 0.6)*10^-2	0.47 ∓ 0.08	0.15 ∓ 0.02
NOP12	Solyc08g076660.2	Ubi	0.4 ∓ 0.1	0.16 ∓ 0.08	0.06 ∓ 0.02	0.47 ∓ 0.08	0.15 ∓ 0.02
TSR1	Solyc01g096760.2	EF1α	(5.2 ∓ 0.7)*10^-2	0.12 ∓ 0.04	(6 ∓ 2)*10^-2	2.3 ∓ 0.1	1.0 ∓ 0.2
TSR1	Solyc01g096760.2	Ubi	0.11 ∓ 0.03	0.27 ∓ 0.09	0.10 ∓ 0.03	2.3 ∓ 0.1	1.0 ∓ 0.2
UBA2	Solyc01g109960.2	EF1α	(5 ∓ 2)*10^-2	(2.2 ∓0.6)*10^-2	(2.7 ∓ 0.2)*10^-2	0.39 ∓ 0.02	0.48 ∓0.01
UBA2	Solyc01g109960.2	Ubi	0.29 ∓ 0.09	0.109 ∓ 0.005	0.141 ∓ 0.005	0.39 ∓ 0.02	0.48 ∓0.01
UTP20	Solyc09g092240.2	EF1α	(2.7 ∓ 0.3)*10^-2	(1.8 ∓0.2)*10^-2	(2 ∓ 1)*10^-3	0.6 ∓ 0.1	0.10 ∓ 0.03
UTP20	Solyc09g092240.2	Ubi	0.16 ∓ 0.01	(8 ∓3)*10^-2	(2 ∓ 1)*10^-2	0.6 ∓ 0.1	0.10 ∓ 0.03
AIR2	Solyc05g007090.2	EF1α	(1.37 ∓ 0.07)*10^-2	(4 ∓ 1)*10^-3	(1.1 ∓ 0.2)*10^-2	0.37 ∓ 0.07	0.87 ∓ 0.09
AIR2	Solyc05g007090.2	Ubi	(2.1 ∓ 0.2)*10^-2	(9 ∓ 4)*10^-3	(2.0 ∓ 0.4)*10^-2	0.37 ∓ 0.07	0.87 ∓ 0.09
ENP1	Solyc07g056550.2	EF1α	(9 ∓3)*10^-2	0.12 ∓ 0.05	(6.8 ∓0.9)*10^-2	1.3 ∓0.1	0.5 ∓ 0.3
ENP1	Solyc07g056550.2	Ubi	0.5 ∓ 0.2	0.6 ∓ 0.2	0.12 ∓ 0.02	1.3 ∓0.1	0.5 ∓ 0.3
DED1	Solyc03g052980.2	EF1α	(8.6 ∓ 0.1)*10^-2	(5 ∓2)*10^-2	(2.8 ∓ 0.1)*10^-2	0.56 ∓ 0.02	0.52 ∓ 0.06
DED1	Solyc03g052980.2	Ubi	0.17 ∓ 0.01	0.10 ∓ 0.04	(8.4 ∓0.1)*10^-2	0.56 ∓ 0.02	0.52 ∓ 0.06
REI1	Solyc12g096450.1	EF1α	(2.2 ∓ 0.2)*10^-3	(1.0 ∓ 0.3)*10^-2	(1.2 ∓ 0.5)*10^-2	4.3 ∓ 0.5	4.9 ∓ 0.6
REI1	Solyc12g096450.1	Ubi	(1.3 ∓ 0.1)*10^-2	(5.2 ∓ 0.6)*10^-2	(6.0 ∓ 0.7)*10^-2	4.3 ∓ 0.5	4.9 ∓ 0.6

To evaluate the impact of expression in pollen on the transcript abundance observed for anthers, we calculated the ratio between the determined qRT-PCR values after normalization to the EF1 expression (Fig. 8). For two genes coding for ASC1 (Solyc12g0405l0) and RRP5 (Solyc03g05l900), we observed a significantly higher transcript level in one of the pollen stages than in anthers. This suggests that the level observed in anthers represents the expression in pollen. Remarkably, both of the genes are not significantly higher expressed in anthers compared to the transcript level in leaves (Table 4). For six of the analyzed genes, we observed a comparable expression in pollen tissues and anthers. However, the gene that shows the highest expression in anthers compared to leaves does not show an enhanced expression in pollen compared to anthers (RIL1; Solyc08g075360; Fig. 8, Tables 4 and 5). This documents that the enhanced expression in sporophytic tissues and gametophytic cells compared to leaves cannot directly be related to expression in pollen.

Figure 8

Comparison of expression in pollen and anthers. The ratio of the values listed in Tables 4 and 5 normalized to EF1α expression was calculated.

Discussion

RBF Co-Orthologs Show a Differential Expression in Tomato

We present the set of co-orthologs to RBFs and RPs in 14 plant genomes (Supplementary Tables 1-4). However, this analysis has two limits. On the one hand, additional information besides the assignment as co-ortholog is required to assign a gene to a specific functional pathway³⁸; on the other hand, RBFs without orthologs in yeast cannot be detected by this approach. Nevertheless, the number of co-orthologs to yeast RBFs assigned in plant genomes demonstrates the conservation of the process and the corresponding proteins can be used as markers to investigate plant-specific pre-ribosomal complex compositions eg by proteomic studies.

With focus on tomato we predicted 245 co-orthologs to yeast RBFs. They represent orthologs to 173 of 255 RBFs described in yeast (Figs. 2, 4, and 5). With the exception of three genes (Solyc00g095450, Solyc00g185750, Solyc05g018780) all are found to be expressed in leaves and/or anthers (Fig. 6). Remarkably, based on existing results (GEO: GSE33507) we could not provide evidence for expression of YAR1, while the proposed YVH1 (Solyc12g021185) is expressed. Thus, based on ortholog search and expression analysis we found 243 co-orthologs in tomato to 172 RBFs in yeast.

Previously, several factors have been experimentally confirmed to participate in ribosome biogenesis in plants. Most of them belong to the above-described orthologs to yeast factors. In many cases, the nomenclature used matches the nomenclature proposed for yeast (Pwp2, Rrp5, Noc4, Enp1, and Nob1¹⁷; Dim1A⁴¹; Rrp6L1 and Rrp6L2⁴²; Mtr4⁴³; Nuc1¹²; eIf6⁴⁴; RID2⁴⁵; OsNug2⁴⁶; AtLsg1¹⁴; Supplementary Table 12). In addition, AtFib1 and AtFib2⁴⁷ were assigned as co-orthologs of Nop1; AtSwa1¹⁸ as ortholog of Utp15; AtRtl2⁴⁸ as ortholog of Rnt1; AtSwa2³⁹ as ortholog of Mak21; AtRH36/AtSwa3¹⁶ as ortholog of Dbp8; AtXrn2⁴⁹ as one ortholog of Rat1/ Xrn1; AtGigantus1⁵⁰ as ortholog of YNL035C; AtRH57⁵¹ as ortholog of Rok1; and AtNuc2¹² as co-ortholog of Nuc1 (Supplementary Table 12).

To the best of our knowledge, only two proteins that might be involved in ribosome biogenesis in plants have been described, for which no ortholog was detected in yeast, namely AtNufIP⁹ and Domino1.⁵² In tomato, these two factors are encoded by Solyc01g096920 and Solyc06g069770, respectively (Supplementary Table 12). Co-orthologs for Domino1 are present in all analyzed plant species. In contrast, the C-terminus of AtNufIP shares a high sequence homology with the C-terminal portion of Rsa1, which can be identified in proteins found in the genome of monocots and dicots,⁹ while orthologs to AtNuflP can only be identified in dicots most likely because of the very variable N-terminus of the protein (Supplementary Table 12). Combining this information with our ortholog search and expression analysis, the 245 tomato proteins can be assigned to ribosome biogenesis with high likelihood. These proteins represent 172 RBFs described in yeast.

In general, about 25% of all RBFs predicted in A. thaliana and tomato are encoded by at least two different genes, which applies for all complexes involved in ribosome biogenesis (Figs. 2–4). The genes coding for RBFs are more strongly expressed in anther compared to leaf (Fig. 6; Table 3), which is not surprising, considering that anther is a male reproductive organ undergoing rapid metabolic changes. However, we also observed genes that show a stronger expression in leaves. These are mostly RBFs represented by multiple co-orthologs (Table 3), which is a hint for a tissue-specific set of assembly factors for ribosome biogenesis. In-depth analysis of the expression pattern of tomato RBFs with multiple co-orthologs in different tissues revealed differential expression for all (Fig. 9). Thus, we can conclude that the co-orthologs to individual RBFs found in plants account for a preferential expression in distinct tissues or at different developmental stages. This is consistent with the observed variation of the number of RBFs to which co-orthologs could be assigned (Fig. 2B). We conclude that plant ribosome biogenesis requires the action of tissue- or developmental stage-specific factors, which is in agreement with reports showing tissue-specific developmental defects in some RBF mutants.^16,17,19,39 On the other hand, considering the critical function of individual RBFs for a vital cellular process like ribosome biogenesis, we can assume that other proteins with similar functions might compensate for the factors not predicted in any plant, a hypothesis that needs to be experimentally challenged in the future. On the other hand, the absence of some factors in plants might be explained by differences in the processing pathway compared to yeast, since for example, not all cleavage sites on plant rRNA have been mapped so far.

Figure 9

The expression profile of RBF encoded by multiple genes. The expression values (Supplementary Table 4 for genes coding for RBFs with multiple co-orthologs have been normalized to the individual maximal value given next to the panel and the expression profile is shown (scale on the right). The order of samples is indicated on the right, legend in bold indicates the MACE results; the rest of the samples are derived from the NGS data.³³

Correlation of RP and RBF Gene Expression

In yeast, it has been shown that many RPS but only a few RPL proteins are present in the 90S pre-ribosomal particle.⁵³ In line, it is discussed that in eukaryotes the majority of the RPS of the body of the 40S subunit assemble co-transcriptionally.⁵⁴ Interestingly, the cluster with 30 RPs represents co-orthologs of 12 RPS genes, 9 of which are classified as very early assembling (RPS3A, RPS4, RPS7, RPS8, RPS11, RPS13, RPS23, RPS24, RPS27).⁵⁵ The other three RPS are considered as very late assembling (RPS17, RPS20, RPS26). The cluster with 20 RBFs contains genes encoding RBFs that are associated with the 90S particle (ENP2, FAL1, FYV7, HRR25, LCP5, NOC4, SAS10, TLR1/RRP36). Interestingly, this cluster of RBF genes has a similar expression profile to that of a cluster of 30 RPs (Fig. 7). Thus, these RBFs of the 90S particle and the RPSs might be functionally related. In line with such notion, yeast ENP2, FAL1, LCP5, and TLR1/RRP36 are required for processes leading to the 40S subunit.^56–59

In turn, RLI1, present in the cluster with 20 RBFs, associates with 40S subunits and is involved in the ribosome quality control and recycling.^60,61 This might link its function to the late assembly of RPS17, RPS20, and RPS26, representing the late assembling RPSs of the cluster with 30 RPs. Worth mentioning is that only one RIL co-ortholog shows the profile represented by the cluster with 20 RBFs (20 rbf).

Similarly, the RPLs of the cluster with 30 RPs suggest a regulation of the early processes of 60S maturation. In yeast, at least 10 of the 15 RPLs of this cluster have been identified in the pre-60S E₁, the first pre-60S complex formed during 60S biogenesis (RPL5, RPL7A, RPL13A, RPL17, RPL18A, RPL26, RPL27, RPL30, RPL35, RPL36A).⁶² However, only for three of the six assigned 60S associated RBFs of the cluster with 20 RBFs (NIP7, NOP16, SPB1), a direct association with the early pre-60S E₁ particle has been confirmed in yeast.^63,64 In contrast, ARB1 was described as an export factor and thus is a late factor.⁶⁵ Thus, the function of ARB1 might be related to the assembly of RPLs like RPL12, RPL23, and RPL37 present in the cluster with 30 RPs.

Summing up, it appears that co-orthologs of RBFs and RPs involved in early events of RP assembly are expressed in a similar manner. Furthermore, we observed that the highest populated clusters for each RP (Supplementary Table 13) and RBF (Supplementary Table 9) combine genes with maximum expression in young leaves (RP: 59 rp; RBF: 57 rbf). Interestingly, we also obtained one cluster of RBF factors, but not of RPs, which shows highest expression in anthers (19 rbf). This might suggest that a difference in tissue-specific ribosomes^66,67 is not always dependent on the composition of RPs.

Author Contributions

ES conceptualized the project. ES, SF, SS and KDS headed the project. SS, JE and MK were involved in co-ortholog prediction, pipeline CRACPipe design, and NGS data analysis. SF and PP performed the plant experiments. SS, SF and ES performed data analysis and manuscript drafting. All authors approved the manuscript.

Footnotes

Acknowledgments

We would like to thank Oliver Mirus for constant help and support and Markus T. Bohnsack for support during pipeline development. We acknowledge funding by the Deutsche Forschungsgemeinschaft and by the EU/Marie Curie.

Abbreviations

References

Dez

, Tollervey

Ribosome synthesis meets the cell cycle.

Curr Opin Microbiol. 2004; 7: 631–7.

Henras

A.K.

, Soudet

, Gerus

. The post-transcriptional steps of eukaryotic ribosome biogenesis. Cell Mol Life Sci. 2008; 65: 2334–59.

Woolford

Jr , Baserga

S.J.

Ribosome biogenesis in the yeast Saccharomyces cerevisiae.

Genetics. 2013; 195: 643–81.

Thomson

, Tollervey

The final step in 5.8S rRNA processing is cytoplasmic in Saccharomyces cerevisiae.

Mol Cell Biol. 2010; 30: 976–84.

Karbstein

Quality control mechanisms during ribosome maturation.

Trends Cell Biol. 2013; 23: 242–50.

Watkins

N.J.

, Bohnsack

M.T.

The box C/D and H/ACA snoRNPs: key players in the modification, processing and the dynamic folding of ribosomal RNA.

Wiley Interdiscip Rev RNA. 2012; 3: 397–414.

Bachellerie

J.P.

, Cavaille

, Huttenhofer

The expanding snoRNA world.

Biochimie. 2002; 84: 775–90.

Brown

J.W.

, Echeverria

, Qu

L.H.

. Plant snoRNA database. Nucleic Acids Res. 2003; 3: 432–5.

Rodor

, Letelier

, Holuigue

, Echeverria

Nucleolar RNPs: from genes to functional snoRNAs in plants.

Biochem Soc Trans. 2010; 38: 672–6.

10.

, Zhang

, Huang

, Zheng

Genome-wide comparative in silico analysis of the RNA helicase gene family in Zea mays and Glycine max: a comparison with Arabidopsis and Oryza sativa.

PLoS One. 2013; 8: e78982.

11.

Ebersberger

, Simm

, Leisegang

M.S.

. The evolution of the ribosome biogenesis pathway from a yeast perspective. Nucleic Acids Res. 2014; 42: 1509–23.

12.

Durut

, Abou-Ellail

, Pontvianne

. A duplicated NUCLEOLIN gene with antagonistic activity is required for chromatin organization of silent 45S rDNA in Arabidopsis. Plant Cell. 2014; 26: 1330–44.

13.

Pontvianne

, Matia

, Douet

. Characterization of AtNUC-L1 reveals a central role of nucleolin in nucleolus organization and silencing of AtNUC-L2 gene in Arabidopsis. Mol Biol Cell. 2007; 18: 369–79.

14.

Weis

B.L.

, Missbach

, Marzi

, Bohnsack

M.T.

, Schleiff

The 60S associated ribosome biogenesis factor LSG1-2 is required for 40S maturation in Arabidopsis thaliana.

Plant J. 2014; 80(6): 1043–56.

15.

Zawada

A.M.

, Rogacev

K.S.

, Müller

. Massive analysis of cDNA Ends (MACE) and miRNA expression profiling identifies proatherogenic pathways in chronic kidney disease MACE and miRNA profiling in CKD. Epigenetics. 2013; 9: 161–72.

16.

Huang

C.K.

, Huang

L.F.

, Huang

J.J.

, Wu

S.J.

, Yeh

C.H.

, Lu

C.A.

A DEAD-box protein, AtRH36, is essential for female gametophyte development and is involved in rRNA biogenesis in Arabidopsis.

Plant Cell Physiol. 2010; 51: 694–706.

17.

Missbach

, Weis

B.L.

, Martin

, Simm

, Bohnsack

M.T.

, Schleiff

40S ribosome biogenesis co-factors are essential for gametophyte and embryo development.

PLoS One. 2013; 8: e54084.

18.

Shi

D.Q.

, Liu

, Xiang

Y.H.

, Ye

, Sundaresan

, Yang

W.C.

SLOW WALKER1, essential for gametogenesis in Arabidopsis, encodes a WD40 protein involved in 18S ribosomal RNA biogenesis.

Plant Cell. 2005; 17: 2340–54.

19.

Byrne

M.E.

A role for the ribosome in development.

Trends Plant Sci. 2009; 14: 512–9.

20.

Nakao

, Yoshihama

, Kenmochi

RPG: the ribosomal protein gene database.

Nucleic Acids Res. 2004; 32: D168–70.

21.

Paul

, Simm

, Blaumeiser

. The protein translocation systems in plants - composition and variability on the example of Solanum lycopersicum. BMC Genomics. 2013; 14: 189.

22.

Chaturvedi

, Ischebeck

, Egelhofer

, Lichtscheidl

, Weckwerth

Cell-specific analysis of the tomato pollen proteome from pollen mother cell to mature pollen provides evidence for developmental priming.

J Proteome Res. 2013; 12: 4892–903.

23.

Yakovlev

I.A.

, Lee

, Rotter

. Temperature-dependent differential transcriptomes during formation of an epigenetic memory in Norway spruce embryogenesis. Tree Genet Genomes. 2014; 10: 355–66.

24.

Kent

W.J.

BLAT - the BLAST-like alignment tool.

Genome Res. 2002; 12: 656–64.

25.

James

L.R.

, Andrews

, Walker

. High-throughput analysis of calcium signalling kinetics in astrocytes stimulated with different neurotransmitters. PLoS One. 2011; 6: e26889.

26.

Ning

, Cox

A.J.

, Mullikin

J.C.

SSAHA: a fast search method for large DNA databases.

Genome Res. 2001; 11: 1725–9.

27.

Cartwright

R.A.

, Graur

The multiple personalities of Watson and Crick strands.

Biol Direct. 2011; 6: 7.

28.

Eddy

S.R.

Non-coding RNA genes and the modern RNA world.

Nat Rev Genet. 2001; 2: 919–29.

29.

Homann

O.R.

, Johnson

A.D.

MochiView: versatile software for genome browsing and DNA motif analysis.

BMC Biol. 2010; 8: 49.

30.

Asmann

Y.W.

, Klee

E.W.

, Thompson

E.A.

. 3’ tag digital gene expression profiling of human brain and universal reference RNA using Illumina Genome Analyzer. BMC Genomics. 2009; 10: 531.

31.

Tarazona

, Garcia-Alcalde

, Dopazo

, Ferrer

, Consea

Differential expression in RNA-seq: a matter of depth.

Genome Res. 2011; 21(12): 2213–23.

32.

Seesi

S.A.

, Tiagueu

Y.T.

, Zelikovsky

, Măndoiu

I.I.

Bootstrap-based differential gene expression analysis for RNA-Seq data with and without replicates.

BMC Genomics. 2014; 15(suppl 8): S2.

33.

Tomato Genome

The tomato genome sequence provides insights into fleshy fruit evolution.

Nature. 2012; 485: 635–41.

34.

Soukas

, Cohen

, Socci

N.D.

, Friedman

J.M.

Leptin-specific patterns of gene expression in white adipose tissue.

Genes Dev. 2000; 14: 963–80.

35.

Tibshirani

, Walther

, Hastie

Estimating the number of clusters in a data set via the gap statistic.

J R Stat Soc Series B Stat Methodol. 2001; 63: 411–23.

36.

Fragkostefanakis

, Simm

, Paul

, Bublak

, Scharf

K.D.

, Schleiff

Chaperone network composition in Solanum lycopersicum explored by transcriptome profiling and microarray meta-analysis.

Plant Cell Environ. 2014; doi: 10.1111/pce.12426.

37.

Livak

K.J.

, Schmittgen

T.D.

Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method.

Methods. 2001; 25: 402–8.

38.

Koonin

E.V.

Orthologs, paralogs, and evolutionary genomics.

Annu Rev Genet. 2005; 39: 309–38.

39.

, Yuan

, Liu

. SLOW WALKER2, a NOC1/MAK21 homologue, is essential for coordinated cell cycle progression during female gametophyte development in Arabidopsis. Plant Physiol. 2009; 151: 1486–97.

40.

Honys

, Twell

Transcriptome analysis of haploid male gametophyte development in Arabidopsis.

Genome Biol. 2004; 5: R85.

41.

Wieckowski

, Schiefelbein

Nuclear ribosome biogenesis mediated by the DIM1 A rRNA dimethylase is required for organized root growth and epidermal patterning in Arabidopsis.

Plant Cell. 2012; 24: 2839–56.

42.

Lange

, Holec

, Cognat

. Degradation of a polyadenylated rRNA maturation by-product involves one of the three RRP6-like proteins in Arabidopsis thaliana. Mol Cell Biol. 2008; 28: 3038–44.

43.

Lange

, Sement

F.M.

, Gagliardi

MTR4, a putative RNA helicase and exosome co-factor, is required for proper rRNA biogenesis and development in Arabidopsis thaliana.

Plant J. 2011; 68: 51–63.

44.

Kato

, Konishi

, Shigyo

, Yoneyama

, Yanagisawa

Characterization of plant eukaryotic translation initiation factor 6 (eIF6) genes: the essential role in embryogenesis and their differential expression in Arabidopsis and rice.

Biochem Biophys Res Comm. 2010; 397: 673–8.

45.

Ohbayashi

, Konishi

, Ebine

, Sugiyama

Genetic identification of Arabidopsis RID2 as an essential factor involved in pre-rRNA processing.

Plant J. 2011; 67: 49–60.

46.

C.H.

, Hwang

S.M.

, Son

Y.S.

. Nuclear/nucleolar GTPase 2 proteins as a subfamily of YlqF/YawG GTPases function in pre-60S ribosomal subunit maturation of mono- and dicotyledonous plants. J Biol Chem. 2011; 286: 8620–32.

47.

Saez-Vasquez

, Caparros-Ruiz

, Barneche

, Echeverria

Characterization of a crucifer plant pre-rRNA processing complex.

Biochem Soc Trans. 2004; 32: 578–80.

48.

Comella

, Pontvianne

, Lahmy

. Characterization of a ribonuclease III-like protein required for cleavage of the pre-rRNA in the 3'ETS in Arabidopsis. Nucleic Acids Res. 2008; 36: 1163–75.

49.

Zakrzewska-Placzek

, Souret

F.F.

, Sobczyk

G.J.

, Green

P.J.

, Kufel

Arabidopsis thaliana XRN2 is required for primary cleavage in the pre-ribosomal RNA.

Nucleic Acids Res. 2010; 38: 4487–502.

50.

Gachomo

E.W.

, Jimenez-Lopez

J.C.

, Baptiste

L.J.

, Kotchoni

S.O.

GIGANTUS1 (GTS1), a member of Transducin/WD40 protein superfamily, controls seed germination, growth and biomass accumulation through ribosome-biogenesis protein interactions in Arabidopsis thaliana.

BMC Plant Biol. 2014; 14: 37.

51.

Hsu

Y.F.

, Chen

Y.C.

, Hsiao

Y.C.

. AtRH57, a DEAD-box RNA helicase, is involved in feedback inhibition of glucose-mediated abscisic acid accumulation during seedling development and additively affects pre-ribosomal RNA processing with high glucose. Plant J. 2014; 77: 119–35.

52.

Lahmy

, Guilleminot

, Cheng

C.M.

. DOMINO1, a member of a small plant-specific gene family, encodes a protein essential for nuclear and nucleolar functions. Plant J. 2004; 39: 809–20.

53.

Grandi

, Rybin

, Bassler

. 90S pre-ribosomes include the 35S pre-rRNA, the U3 snoRNP, and 40S subunit processing factors but predominantly lack 60S synthesis factors. Mol Cell. 2002; 10: 105–15.

54.

Karbstein

Inside the 40S ribosome assembly machinery.

Curr Opin Chem Biol. 2011; 15: 657–63.

55.

O'Donohue

M.F.

, Choesmel

, Faubladier

, Fichant

, Gleizes

P.E.

Functional dichotomy of ribosomal proteins during the synthesis of mammalian 40S ribo-somal subunits.

J Cell Biol. 2010; 190: 853–66.

56.

Gerus

, Bonnart

, Caizergues-Ferrer

, Henry

, Henras

A.K.

Evolutionarily conserved function of RRP36 in early cleavages of the pre-rRNA and production of the 40S ribosomal subunit.

Mol Cell Biol. 2010; 30: 1130–44.

57.

Kressler

, de la Cruz

, Rojo

, Linder

Fal1p is an essential DEAD-box protein involved in 40S-ribosomal-subunit biogenesis in Saccharomyces cerevisiae.

Mol Cell Biol. 1997; 17: 7283–94.

58.

Soltanieh

, Lapensee

, Dragon

Nucleolar proteins Bfr2 and Enp2 interact with DEAD-box RNA helicase Dbp4 in two different complexes.

Nucleic Acids Res. 2014; 42: 3194–206.

59.

Wiederkehr

, Pretot

R.F.

, Minvielle-Sebastia

Synthetic lethal interactions with conditional poly(A) polymerase alleles identify LCP5, a gene involved in 18S rRNA maturation.

RNA. 1998; 4: 1357–72.

60.

van den Elzen

A.M.

, Schuller

, Green

, Seraphin

Dom34-Hbs1 mediated dissociation of inactive 80S ribosomes promotes restart of translation after stress.

EMBO J. 2014; 33: 265–76.

61.

Yarunin

, Panse

V.G.

, Petfalski

, Dez

, Tollervey

, Hurt

E.C.

Functional link between ribosome formation and biogenesis of iron-sulfur proteins.

EMBO J. 2005; 24: 580–8.

62.

Fatica

, Cronshaw

A.D.

, Dlakic

, Tollervey

Ssf1p prevents premature processing of an early pre-60S ribosomal particle.

Mol Cell. 2002; 9: 341–51.

63.

Lebreton

, Rousselle

J.C.

, Lenormand

. 60S ribosomal subunit assembly dynamics defined by semi-quantitative mass spectrometry of purified complexes. Nucleic Acids Res. 2008; 36: 4988–99.

64.

Talkish

, Zhang

, Jakovljevic

, Horsey

E.W.

, Woolford

Jr . Hierarchical recruitment into nascent ribosomes of assembly factors required for 27SB pre-rRNA processing in Saccharomyces cerevisiae. Nucleic Acids Res. 2012; 40: 8646–61.

65.

Dong

, Lai

, Jennings

J.L.

, Link

A.J.

, Hinnebusch

A.G.

The novel ATP-binding cassette protein ARB1 is a shuttling factor that stimulates 40S and 60S ribosome biogenesis.

Mol Cell Biol. 2005; 25: 9859–73.

66.

Wang

, Lan

, Gao

, Zheng

, Li

, Schmidt

Expression changes of ribo-somal proteins in phosphate- and iron-deficient Arabidopsis roots predict stressspecific alterations in ribosome composition.

BMC Genomics. 2013; 14: 783.

67.

Xue

, Barna

Specialized ribosomes: a new frontier in gene regulation and organismal biology.

Nat Rev Mol Cell Biol. 2012; 13: 355–69.