
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)01884-4
10.1016/j.isci.2024.110659
110659
Article
Multi-omics analysis of aggregative multicellularity
Edelbroek Bart bart.edelbroek@gmail.com
1∗
Westholm Jakub Orzechowski 2
Bergquist Jonas 3
Söderbom Fredrik fredrik.soderbom@icm.uu.se
14∗∗
1 Department of Cell and Molecular Biology, BMC, Uppsala University, 751 24 Uppsala, Sweden
2 Department of Biochemistry and Biophysics, National Bioinformatics Infrastructure Sweden, Science for Life Laboratory, Stockholm University, Stockholm, Sweden
3 Department of Chemistry-BMC, Analytical Chemistry and Neurochemistry, Uppsala University, Uppsala, Sweden
∗ Corresponding author bart.edelbroek@gmail.com
∗∗ Corresponding author fredrik.soderbom@icm.uu.se
4 Lead contact

03 8 2024
20 9 2024
03 8 2024
27 9 11065914 3 2024
14 6 2024
31 7 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by/4.0/ This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Summary

All organisms have to carefully regulate their gene expression, not least during development. mRNA levels are often used as proxy for protein output; however, this approach ignores post-transcriptional effects. In particular, mRNA-protein correlation remains elusive for organisms that exhibit aggregative rather than clonal multicellularity. We addressed this issue by generating a paired transcriptomics and proteomics time series during the transition from uni-to multicellular stage in the social ameba Dictyostelium discoideum. Our data reveals that mRNA and protein levels correlate highly during unicellular growth, but decrease when multicellular development is initiated. This accentuates that transcripts alone cannot accurately predict protein levels. The dataset provides a useful resource to study gene expression during aggregative multicellular development. Additionally, our study provides an example of how to analyze and visualize mRNA and protein levels, which should be broadly applicable to other organisms and conditions.

Graphical abstract

Highlights

• Extensive capture of transcriptome and proteome through multicellular development

• High correlation of transcript and protein abundance is reduced during development

• Dynamically regulated mRNA results in linear regulation of cognate protein

• A time lag of 2–4 h is observed between mRNA and protein regulation

Organizational aspects of cell biology; Developmental biology; Expression study; Omics; Proteomics; Transcriptomics

Subject areas

Organizational aspects of cell biology
Developmental biology
Expression study
Omics
Proteomics
Transcriptomics
Published: August 3, 2024
==== Body
pmcIntroduction

One of the pillars of biology is the central dogma which states that, exceptions aside, the information stored in DNA is transferred to RNA and further passed on to produce proteins.1 Messenger RNA (mRNA) is transcribed from DNA and translated into protein in a quantitative manner. It follows that the levels of mRNA and protein are connected to each other. Both are dynamically regulated and change due to environmental and developmental cues.

Transcriptomic and proteomic analyses have experienced incredible improvement in throughput and accuracy, however, sensitivity of proteomics with LC-MS/MS is still lagging behind the sensitivity of transcriptomics by next-generation sequencing.2 Also, in proteomics the signal cannot be amplified by, e.g., PCR. When comparing these two omics approaches, studies based on transcriptomics have increased dramatically relative to studies based on proteomics, reflected by the number of datasets available in different repositories.3,4

It is not uncommon that mRNA levels are used as a proxy for the levels of the effector molecules – the proteins. There is, however, not always a linear relationship between mRNA and protein, which can be attributed to, e.g., differences in translation rates or protein stability.5 Studying the correlation between mRNA and protein levels is important in order to understand to what extent transcriptomics data can be used to predict gene expression.6

Organisms across the Tree of Life have evolved distinct strategies which control the balance between mRNA and protein.7 Due to these differences, the mRNA-protein correlation can differ between species, especially those that are phylogenetically distantly related. Additionally, the correlation between mRNA and protein can be different depending on the biological context of the cell. In steady-state cells, at the population level, the mRNA levels and protein levels are expected to be relatively stable, and their correlation high.5 On the other hand, when the cells are undergoing changes, many genes will be differentially regulated and the correlation might be lower. One example of this is development, where cells undergo major differential gene expression. Here, cells transition to a different behavior or identity through intrinsic or extrinsic signals.

Previously, the relationship between mRNA and protein has been studied by paired transcriptomics and proteomics at specific developmental stages in organisms dedicated to clonal multicellularity, such as animals and plants.8,9,10,11 However, less is known about the correlation between mRNA and protein in other organisms. Such studies would provide insight into evolutionary aspects of the control of genetic output. Importantly, the mRNA-protein correlation during the transition from uni-to multicellular life has not yet been extensively explored. This transition can be studied in organisms with facultative multicellularity, also referred to as aggregative multicellularity, i.e., they can alter between uni- and multicellular lifestyles.12,13 This enables studies of processes involved in transforming single cells into a multicellular organism. Facultative multicellularity evolved independently at least seven times during the course of evolution and has been found in the majority of all eukaryotic linages.13 The best studied example of organisms that go through aggregative multicellularity belong to Dictyostelia social amebae.14 One of these is the well-established model organism Dictyostelium discoideum.15 When the amebae run out of food, a developmental program is initiated where, D. discoideum transitions from free-living to multicellular. First, the cells form aggregates of up to 100,000 cells, which then continue to develop into a fruiting body, where dead stalk cells support a ball of spores.16 This aggregative multicellular development has been studied thoroughly at the RNA level, characterizing the main processes involved in the developmental program as well as differentiation into specialized cell types at the single-cell level.17,18,19,20,21 Thus far, the developmental proteome has not been extensively studied, and it remains unknown how well the observed transcriptional changes are reflected at the protein level.

In this study, we performed transcriptomics and proteomics analyses at several time points during growth and early development of D. discoideum to elucidate the mRNA and protein levels throughout multicellular aggregation. This setup allowed for a unique opportunity to study the regulation of genetic output during the transition from uni-to multicellular life. We confirmed previous findings, which identified differentially regulated genes involved in processes essential for early development. Additionally, we detected many genes that are dynamically regulated, where specific mRNAs can be, for example, upregulated early during development and downregulated at the later stages, and vice versa. However, at the protein level many of the dynamically regulated mRNAs result in linearly regulated proteins. Another observation was that, in general, protein expression is delayed several hours as compared to mRNA expression. Levels of mRNA and protein correlate to a high degree during growth (across genes Spearman correlation = 0.65). The correlation decreases after the onset of development, mainly due to the time lag between mRNA transcription and protein translation. Hence, the data presented here show that the correlation between transcriptomics and proteomics is dependent on the conditions being studied and it is important to proceed with caution when using the transcriptome as a proxy for protein expression. The data presented in this study will also be a valuable resource for investigating D. discoideum development and are available in an interactive web app for ease of use: https://edelbroek2024.serve.scilifelab.se.

Results

Experimental setup

The multicellular development of Dictyostelium discoideum starts when unicellular amebae starve and embark on a developmental program. During this process, cells go through distinct morphological stages including streaming where cells move together in response to the secreted chemoattractant cAMP. This is followed by the formation of a loose aggregate of about 100,000 cells, which subsequently produce an acellular matrix that covers the aggregate to make up the mound/tight aggregate stage. At this time cells start to differentiate, mainly into cells that will become spore- and stalk-cells, and form a motile slug that culminates into a fruiting body or “sorocarp”22 (Figure 1A). Here, we aimed to investigate how the transcriptome and proteome are regulated and correlated during early multicellular development of D. discoideum. Cells were starved on non-nutrient agar plates to induce multicellular aggregation after which they were harvested at time increments during 0–10h post starvation (Figures 1A and 1B). In order to minimize biological and technical variations, we collected cells from both halves of each plate and processed the cells for proteomics and transcriptomics, respectively (Figure 1A).Figure 1 Experimental setup

(A) Axenically grown cells were washed and plated on non-nutrient agar plates to induce multicellular development. Samples were taken from the same plate for transcriptomics and proteomics at 0h, 2h, 4h, 8h, and 10h post initiation of starvation, for a total of 15 paired transcriptomics and proteomics libraries.

(B) Multicellular development of D. discoideum on non-nutrient agar plates, at 4h, 6h, 8h, and 10h post initiation of starvation.

Major reorganization of the transcriptome during multicellular aggregation

In order to acquire a more detailed understanding of the dynamically regulated transcriptome during multicellular development, we decided to include additional RNA-seq libraries. Beside the 15 paired libraries, each time point was supplemented with a fourth biological replicate. In addition, we added four biological replicates of an extra time-point (6-h time point), for a total of 24 RNA-seq libraries. The broad transcriptional changes over time were investigated by principal component analysis (PCA) (Figure 2A). Developmental time correlates highly with the first principal component, whereas the second principal component oscillates from 0h to 6h back to 10h. The biological replicates showed minimal variation, both in the PCA plot and by calculating their correlation (Figure S1). From our data, we could identify 8310 protein coding transcripts differentially expressed at any of the time points (likelihood ratio test23 comparing a model where gene expression is determined by developmental time vs. a model of constant expression, FDR-adjusted p-value <0.01). This suggests that the great majority of the in total 11,866 proteins are regulated at the transcript level during the first 10h of multicellular development (Tables 1 and S1).Figure 2 Major reorganization of transcriptome during multicellular aggregation

(A) Principal component analysis (PCA) of the developmental transcriptome based on the 500 transcripts that show the most variation. The first two principal components (PC1, PC2) are shown, which together explain 92.6% of the variance. Four biological replicates were analyzed per time point. Each replicate is plotted as a number, representing the time point of the replicate.

(B) Regulation of the milestone gene transcripts over time. Milestone genes are included, which characterize the transition from “no aggregation” to “streaming” (streaming), from “streaming” to “loose aggregate” (loose aggregate), from “loose aggregate” to “mound” (mound) and from “mound” to “tipped aggregate”,17 from top to bottom. Illustrations of the morphological structures are included on the left. Whether the milestone genes are defined as downregulated (blue) or upregulated (red) in the transition, is annotated on the left.

(C) Hierarchical clustering of protein coding transcripts based on log fold change (logFC) versus the 0h time point. Transcripts were grouped into four main clusters, with the general regulation of each transcript shown in the heatmap to the left, and the general regulation of the cluster shown with boxplots for each time point to the right. The dashed line indicates a logFC of 0 versus the 0h time point. On the far right, the four most significant GO-terms for each cluster are shown, with Fisher’s exact test p-value for each GO-term (the dashed line indicates p-value 0.01). The size of the filled circles represents the fold enrichment of the GO-term in the cluster. For the full set of significant GO-terms, see Table S2.

Table 1 Numbers of quantified transcripts and proteins, in the separate and combined analysis

Description	transcripts	proteins	transcripts and proteins	
All annotated protein coding genes	11866	11866	11866	
Expressed transcripts/quantified proteins	10714	3663	3604	
% expressed transcripts/quantified proteins of all protein coding genes	90.3%	30.9%	30.4%	
Differentially expressed transcripts/proteins	8310	672	589	
% differentially expressed of total expressed transcripts/quantified proteins	77.6%	18.3%	16.3%	
All annotated protein coding genes represent all genes which feature a Uniprot protein identifier, accessed from dictyBase24 (http://dictybase.org/). The expressed transcripts are those which feature more than one read in any of the replicates, the quantified proteins are those that could be quantified in all replicates, after imputation. Differentially expressed transcripts are defined with DEseq2 likelihood ratio tests,23 and differentially expressed proteins with limma F-test,25 with adjusted p-values below 0.01. For differentially expressed transcripts/proteins, under heading “transcripts and proteins”, genes are included for which both the transcript and protein were differentially expressed in their respective analyses.

A previous study on multicellular development in D. discoideum identified transcripts which undergo sharp changes during developmental transitions – so called “milestone genes”.17 We investigated the regulation of these milestone genes in our transcriptomics dataset (and in our proteomic data, see below). Sharp down- and upregulation could be observed for the milestone genes which characterize the transition from no aggregation to streaming at the 8h time point (Figure 2B). Those that characterize the transition from loose aggregates to mounds are strongly regulated at the 10-h time point (Figure 2B). The temporal regulation of the milestone genes matches the morphological stages observed during development (Figure 1B). It should be noted that the milestone genes defined by Katoh-Kurasawa et al. were identified by developing cells on nitrocellulose filters17 whereas we developed our cells on non-nutrient agar plates. From our data it appears that the regulation of the milestone genes is robust and independent of the method used for development. We could further verify that changes due to the developmental method are minimal by comparing our dataset to the 0–10h data by Rosengarten et al., developed on nitrocellulose.19 Reanalysis of their data revealed broadly the same regulation for protein coding transcripts identified as differentially expressed in both datasets (Figure S2).

The 8310 protein coding transcripts that were identified as differentially expressed were clustered based on their fold change relative to 0h at the different time points. Subsequently, the genes were split into four groups based on hierarchical clustering, i.e., highly upregulated in cluster 1, genes moderately upregulated in cluster 2, genes moderately downregulated from 4 h of development in cluster 3, and genes strongly downregulated in cluster 4 (Figure 2C). In order to classify the differentially expressed genes in each cluster, we performed Gene Ontology-terms (GO-terms) enrichment analysis. This showed that the highly upregulated genes in general are associated with processes connected to development of multicellularity, such as cell-cell recognition, culmination involved in fruiting body development, and assembly of the spore wall (ultimately leading to formation of spores) (Figure 2C). Among the moderately upregulated genes, the formation of the multicellular aggregates is also represented, with terms corresponding to cAMP dependent chemotaxis and response to differentiation-inducting factor 1 (DIF-1) (Table S2). cAMP is secreted by starving cells and used as a chemoattractant to signal to cells to aggregate and DIF-1 is involved in induction of pre-stalk cells.26 Also, the cell cycle is represented, with an enrichment of genes involved in cell division and DNA replication (Figure 2C; Table S2). Even though genes for these processes are upregulated during development, the question whether this leads to cell division and chromosomal DNA replication is being debated (see Katoh-Kurasawa et al. and references within17). Terms in the clusters with downregulated genes include translation and metabolism. This is expected, since starvation induces growth arrest of the cells (Figure 2C). In conclusion, the transcriptomics dataset describes multicellular aggregation in detail, and the GO-terms associated with differently regulated clusters aptly represent biological processes regulated during early development.

The developmentally regulated proteome

In order to understand how well the transcriptomic data correlate with protein expression, we performed proteomic analysis using cells from the same plates from which RNA was isolated for RNA-seq. By performing mass spectrometry (LC-MS/MS; label free quantification) on cells collected from several time points during early development (Figure 1), 2478 proteins could be directly detected and quantified across all biological replicates and time points (Figure 3A). For 1185 proteins, quantification was possible in all biological replicates at one or more time points, but was missing in replicates at other time points. We decided to analyze why there is a discrepancy between the large number of transcripts identified and the lower number of proteins in our study. We hypothesized that the lack of quantification in these replicates was mostly due to a lack of- or low expression of the protein at these time points. This is supported by the fact that proteins which lack quantification in some replicates also have a lower maximum expression level in the replicates with quantification (Figure S3A). We therefore imputed missing values for these 1185 proteins, which were consistently expressed at a given time point, to avoid discarding them from analysis. The missing values were imputed using a probabilistic minimum, which accommodates for values missing due to low expression.27 Addition of the imputed proteins brought the total number of proteins that we could analyze to 3663, about a third of all protein coding genes. For 7061 proteins no quantified peptides could be detected in the dataset. Many of these proteins are likely of low abundance or not expressed in accordance with the mRNA levels of these genes (Figures 3A and S3B). Furthermore, the majority of unidentified proteins have a low annotation score, and low annotation quality may explain why some of the proteins were not identified (Figure S3C).Figure 3 The proteome and its regulation during multicellular development

(A) Donut plot representing the total number of proteins, i.e., protein coding genes (outer circle) and the full transcriptome (inner circle). The inner ring shows the expression level from the transcriptomics analysis from high to low, and is correlated to each group of proteins.

(B) Principal component analysis (PCA) of the proteomics dataset based on the 300 proteins that show most variation. The first two principal components (PC1, PC2) explain 60% of the variance in the dataset. The three replicates for each time point are plotted as separate numbers (time points in hours).

(C) Hierarchical clustering of proteins based on log fold change (logFC) versus the 0h time point. The proteins were grouped into three main clusters, with the general regulation of each protein shown in the heatmap on the left, and the general regulation of the cluster shown with boxplots for each time point next to the heatmap. The dashed line indicates a logFC of 0 versus the 0h time point. On the far right, the four most significant GO-terms for each cluster are shown, with the Fisher’s exact test p-value for each GO-term (the dashed line indicates p-value 0.01) and the size of the filled circle represents the enrichment of the GO-term in the cluster. For the full set of significant GO-terms, see Table S4.

Principal component analysis of the proteomics dataset resembles that of the transcriptomics and shows the same main trends, with the first principal component correlating with developmental time (Figure 3B). Here, too, variation between the biological replicates was low (Figure S1). By analyzing the regulation of the quantified 3663 proteins, 672 were identified as differentially expressed during development (FDR-adjusted p-value <0.01) (Tables 1 and S3). In a study by Kelly et al., the D. discoideum proteome was analyzed during early multicellular development at 0.5h and 8h after initiating development in tissue-culture-treated plates.28 In order to compare their findings with our proteomic analysis, we first reanalyzed their dataset. Subsequent comparison showed that the differentially expressed proteins from either dataset were regulated in similar manner over time (Figure S4A). When we examined the protein levels for the milestone genes (see above), we found that their corresponding proteins are regulated in the same manner as their transcripts (Figures 2B and S4B). In particular, we could detect sharp regulation at the 8h time point for genes which characterize the transition from no aggregation to streaming and this is even more pronounced at the 10h time point (Figure S4B).

In the same way as for the transcriptomic analysis, we grouped the differentially expressed proteins from our study based on their fold change at different time points relative to the 0h time point (Figure 3C). Interestingly, a distinct change can be observed at 8h of development where proteins are either upregulated (clusters 1 and 2) or down regulated (cluster 3). The first two clusters, which are made up of highly or moderately upregulated proteins, are associated with GO-terms linked to the development of multicellular aggregates and fruiting bodies, in line with what was observed in the transcriptomics analyses (Figures 2C and 3C), but additionally protein ubiquitination and proteolysis appear to be upregulated (Table S4). This is not unexpected since D. discoideum development is initiated by starvation and develops without support from external nutrients. Therefore, degradation and recycling of intracellular material, such as proteins, is essential to sustain the different processes required for development.29,30

The proteins that are downregulated from 8h are linked to growth arrest due to the lack of nutrients, where ribosomes and biosynthetic processes are broadly downregulated (Figure 3C). Taken together, the proteomic dataset presented in this study describes a significant fraction of the total D. discoideum proteome during development, and is the first study to follow the regulation of the proteome during aggregative multicellularity in detail.

High steady-state correlation of mRNA and protein levels

From the proteomics dataset it was possible to identify about one-third of the total D. discoideum proteome across replicates and time points. 589 protein-coding genes are differentially expressed in both the transcriptomics and proteomics datasets. These were clustered according to their fold change relative to the 0h time point, and they appear to be largely regulated in the same manner during development (Figures 4A and S5A; Table S5), illustrating that the mRNA and protein expression are correlated. Previously, a number of genes have been established as important regulators of multicellular development.31 We studied their mRNA and protein expression and found these to be generally upregulated over time in our data, in line with their function (Figure S5B).Figure 4 Correlation between mRNA and protein levels

(A) Regulation in log fold change (logFC) versus 0h time point for genes differentially expressed in both transcriptomics and proteomics datasets. The genes in the heatmap are hierarchically clustered based on their regulation. For an expanded plot with all genes differentially expressed in either dataset, see Figure S5A.

(B) Correlation of the mean mRNA and protein levels across the 0h time point. Each dot represents the mean protein and mRNA expression from a single gene. The dashed black line indicates the linear regression of the data, the red line is the y = x diagonal.

(C and D) Positive linear Pearson correlation of gp130, and negative correlation of DDB_G0281185, respectively. The expression in each sample is shown by fraction of total protein and fraction of total mRNA. Each biological replicate is plotted with a number, signifying the time point in hours, and letter, signifying the biological replicate. Only samples for which both transcriptomics and proteomics data was generated, are included. Linear regression is shown with a black dashed line, with the Pearson correlation above the plot.

(E) Distribution of per gene Pearson correlations for all genes with differentially expressed mRNA and protein.

In order to study the correlation of the transcriptomic and proteomic datasets in more detail, we normalized the expression of each gene by the quantification of all genes, such that all mRNA values or protein values at a given time point sum up to 1. Across all genes at the 0-h time point, the mRNA levels and protein levels correlated well (Spearman correlation 0.65, Figure 4B). At later time points however, the correlation dropped to 0.56 (Figure S6). It should also be noted that at all time points, the regression of the data had a slope greater than 1 (Figures 4B and S6). This is indicative of the fact that distribution of the data is different between the transcriptomics and proteomics datasets. We see, relative to the mRNA abundance, a wider dynamic range as well as a more uneven distribution of protein expression levels (Figure S7). For example, for genes quantified in both datasets, the top two proteins together encompass 10% of the total protein abundance, whereas the top two most abundant transcripts encompass a more modest 2.3%. A slope greater than 1 in the mRNA-protein correlation has previously been observed for human tissues as well as in yeast.32,33

Although the mRNA and protein levels correlated well for each time point, we wondered if this is also the case for the regulation of each of these molecules over all samples. For instance, if an mRNA is upregulated at a specific time point during development, does this also hold true for its cognate protein? As an example, the gp130 mRNA and its associated glycoprotein 130 are both downregulated during development, resulting in a high correlation (Pearson’s r = 0.93, Figure 4C). For other genes, a decrease in mRNA abundance between samples coincided with increased protein abundance resulting in a negative correlation (Figure 4D). When considering the Pearson correlation of all genes that were differentially expressed in both the transcriptomics and proteomics datasets, we observed a positive median correlation (median Pearson’s r = 0.54, Figure 4E; Table S6). Hence, in most cases an upregulation of mRNA corresponds to an increase in protein, and vice versa, for all differentially expressed genes (Table S5). Notably, the median correlation is much lower when considering genes which are only differentially expressed in one of the datasets, or not differentially expressed at all (median Pearson’s r = 0.16, Figure S8A). By calculating the correlation per gene, it is also possible to show the advantage of our sampling setup – to isolate mRNA and protein from the same plate (biological replicate) (Figure 1). When we compare to mismatched samples, i.e., when mRNA from one replicate (plate) is compared to protein from another replicate, for the same time point, the median correlation was significantly lower as compared to matching samples from the same plate (Figures S8B and S8C).

In sum, mRNA and protein expression are in general well correlated. Across genes, correlation is highest during steady state growth. The median per gene correlation is high for differentially expressed genes, but not for genes which lack regulation in the transcriptome or proteome during development.

Differences between mRNA and protein regulation

To better understand the details behind the expression patterns in the transcriptomics and proteomics data, we performed an integrated unsupervised analysis using MEFISTO.34 MEFISTO is a factor analysis method that reduces a multi-omics dataset into a few latent factors that explain most of the variance in the full dataset. Running MEFISTO on all genes for which we have both protein and mRNA data (Table 1), resulted in three factors that explain most of the variance in the transcriptomics data, and out of these three factors, Factor 1 also explained a large fraction of the variance in the proteomics data (Figure 5A). Factor 1, which makes up 25% of the variance in the transcriptomics data and 34% in the proteomics data, represents steadily increasing or decreasing expression over the developmental time course (Figure 5B). Factor 2, explaining 27% and 3% of the variance in the transcriptomics and proteomics data, respectively, represents a pattern where expression decreases between 0 and 4 h, after which it plateaus and then increases at 8 and 10 h, or vice versa. Factor 3, which explains 17% and 6% of the variance respectively, shows a dramatic increase or decrease in expression between 0 and 2 h, after which the expression gradually returns (Figure 5B). This shows that mRNA is more dynamically regulated, with more variable expression patterns, whereas proteins mostly show steadily increasing or decreasing levels over the developmental time course.Figure 5 Multi-omics factor analysis of mRNAs and proteins

(A) Percentage of variance explained by the first three factors of the multi-omics factor analysis for the transcriptomics and proteomics datasets.

(B) Trajectory of the first three factors over time. The different replicates are shown, and the dashed line is the loess (locally estimated scatterplot smoothing) regression through the replicates, with 95% confidence interval in grey.

(C) Analysis of genes with high Factor 1 values at both mRNA and protein modalities. The proportion of genes with different regulations shown in the pie-chart. Genes with opposing mRNA and protein regulation are represented in the heatmap. The factor 1 values for the full set of 280 genes are detailed in Table S7.

(D) Gene set enrichment analysis (GSEA) of Factor 2 mRNA values, with GO-terms classified under broad terms. The x axis shows the average Factor 2: mRNA values while the y axis shows the GSEA p-values, after negative log transformation.

(E) Regulation of members of the Arp2/3 complex. For each gene, z-scores were calculated from the mean expression per time point. The gray zone denotes the 95% confidence interval trajectory of all z-scores.

For 280 genes, both the mRNA and protein were highly associated with Factor 1, i.e., linearly up- or downregulated during development. The majority of these were regulated in the same direction in the mRNA and protein modalities (Figure 5C; Table S7). For 11 genes however, mRNA was linearly downregulated and protein upregulated, and vice-versa for 6 other genes (Figure 5C).

We investigated which GO-terms are associated with the three first MEFISTO factors using Gene Set Enrichment Analysis (GSEA).35 Full results are presented in Table S8. In accordance with Figure 2C, negative values for Factor 1, i.e., a steadily decreasing expression pattern, are associated with terms related to translation (e.g., translation, ribosome, and structural constituent of ribosome), for both mRNAs and proteins. Positive values for Factor 1 are associated with protein degradation (e.g., proteasome complex and endopeptidase activity) and also with actin binding/actin filament binding. Actin and actin regulating proteins play multiple roles in D. discoideum such as polarization of ameboid cells during movement but appears to be important also during aggregation and other stages of multicellular development.36 No GO-terms for mRNAs were associated with positive values for Factor 1. Since the majority of the mRNA regulation could be explained not by a steady increase or decrease but by more dynamic patterns, we focused on GO-terms associated with the dynamically regulated mRNA Factor 2 values (Figure 5D). Genes with high Factor 2 values are downregulated between 0h and 6h and then followed by upregulation at 6h–10h. Among the mRNAs matching this expression pattern, GO-terms related to translation and actin are enriched (Figure 5D). This is reflected in regulation of the Arp2/3 complex, a major regulator of the actin cytoskeleton. At the mRNA level, the regulation is similar to Factor 2 (Figure 5E). On the other hand, the complex appears to be linearly upregulated at the protein level, matching Factor 1 (Figure 5E). The proteasome complex on the other hand, is associated with negative Factor 2 values (mRNA), but is also linearly upregulated (protein) (Figures S9A and S9B). Finally, Factor 3 shows a positive association with similar terms (e.g., proteasome complex, endopeptidase activity, and proteolysis) for mRNAs, but no associations for proteins.

In conclusion, the factor analysis revealed several major trajectories of the dynamic mRNA regulation, which appear to be paired with linear up- or downregulation at the protein level.

Protein expression is generally delayed several hours compared to mRNA expression

From the previous analyses, we identified that for a number of genes, the mRNA and protein regulation opposed one another (Figures 5C and 5D). We hypothesized that for some of these genes, the difference may be explained by a time lag between the mRNA transcription and protein translation. To investigate this, we calculated the Spearman correlations of mRNA and protein levels, as before (Figures 4B and S6), but matched all of the transcriptomics time points with all of the proteomics time points. Matching mRNA and protein expression from the same time points shows Spearman correlations of 0.65 to 0.56, however, for all time points, we found higher correlations when matching the mRNA expression with protein expression 2–4 h later (Figure 6A). This observed time lag is in agreement with previous studies in yeast37 and Drosophila.38 When exclusively considering genes that are differentially expressed during development according to the protein and mRNA expression, the trend is largely the same (Figure 6B). Here, however, the maximum correlation is higher, similar to what we previously observed (Figure 4B). Since these are genes which are highly affected throughout development, there is a larger difference between highly correlated pairs of time points (e.g., mRNA at 0h vs. proteins at 4h, Spearman correlation = 0.71, Figure 6B) and those that show low correlation (e.g., mRNA at 10h vs. proteins at 0h protein, Spearman correlation = 0.02, Figure 6B).Figure 6 Time lag between mRNA and protein

(A and B) Spearman correlations of mean protein values and mean mRNA values across genes for time points from 0h to 10h. In A, all genes quantified in both datasets are included in the analysis. In B, only differentially expressed (DE) genes are included.

(C–E) Ratio of protein to mRNA for different time points. Boxplots are based on: (C) all genes identified in both datasets; (D) genes for which the mRNA was upregulated at 10h; (E) genes for which the mRNA was downregulated at 10h. Above the boxplots: Dunnett contrasts p-values relative to 0h time point reported for time points with p-values lower than 0.1. Dashed lines indicate the median protein to mRNA ratio for all genes.

Next, we investigated if the ratios of protein to mRNA are affected during multicellular development. By dividing the protein quantification by the mRNA quantification, the protein to mRNA ratio could be calculated for each gene, at each time point. It should be noted that these ratios are in no way indicative of the absolute numbers of protein or mRNA molecules, and are purely relative values. When considering all genes expressed in both omics datasets, there are no significant differences in protein to mRNA ratios at the different time points (Figure 6C). There are however significant differences between time points when considering genes for which the mRNA is either up- or downregulated during multicellular development (Figures 6D, 6E, and S10). 10h after onset of multicellular development, the ratio of protein to mRNA is significantly decreased for genes which are upregulated (mRNA) (Figure 6D). In contrast, downregulated genes show the opposite effect, with ratios of protein to mRNA significantly increasing over time (Figure 6E). We suspect that this is largely another effect of the time lag between mRNA transcription and protein translation. Those genes that are upregulated at the transcriptional level have a relatively lower level of protein until translation catches up or the mRNA is eventually downregulated again. The opposite is true for genes which are downregulated, here the protein needs to be turned over for the levels to agree with the mRNA.

Taken together, these results show a modest correlation between protein and mRNA levels analyzed at the same developmental time points, however, the correlation increases when considering protein samples taken 2–4 h after the mRNA samples.

Discussion

The evolution of multicellularity is thought to have occurred several times, through both clonal mechanisms, such as in animals and plants, and through aggregative mechanisms, where cells stream together to form multicellular structures.39,40 Aggregative multicellularity has been studied using the social amebae, where processes behind the transition from uni-to multicellular life has been investigated.41,42,43 Numerous studies that focused on individual genes, have identified some of the key players involved in this transition19,44 (Figure S5B), but in recent years next-generation sequencing methodology have paved the way to acquire a complete understanding of this process at the transcriptional level.17,18,19,20,21

Here, we report transcriptomic and proteomic analyses of D. discoideum during early multicellular development. Using LC/MS-MS, we were able to capture the regulation of roughly a third of the proteome during development. This was combined with transcriptomic analyses, covering the great majority of the genes. By analyzing the two datasets separately, but also combining them with factor analysis, we identified several biological processes, vital for D. discoideum multicellular development. Together, this provides us with a detailed picture of gene expression and regulation, from mRNA to protein, during early development of the social ameba.

Included in the upregulated processes, we find chemotaxis and development of the sorocarp, whereas ribosome biogenesis and general metabolism are among the most downregulated (Figures 2C and 3C), in line with what has been observed previously.19 Interestingly, some processes were dynamically up- and downregulated at the mRNA level, but linearly regulated at the protein level. This is similar to what has been observed in vertebrates, where a spike or dip in mRNA causes a switch in protein levels that are either down- or upregulated.10,45 Examples of this kind of regulation in D. discoideum are protein degradation processes and actin related genes, illustrating that relying on mRNA levels only for insights into the temporal impact of these processes may be misleading. Notably, for the majority of regulated mRNAs in our dataset, the protein response is delayed, and the proteins emerge first about 2h–4h after mRNA appearance (Figure 6). This result can likely explain the previous observation that early major morphological changes do not coincide with transcriptomic data, i.e., the phenotypes connected to the expressed mRNAs are delayed.17,19 This result demonstrates that the proteome more accurately describes the functional gene expression and resulting phenotype.6 Hence, it may be preferable to rely on proteomics, and not transcriptomics, for assigning genes to a specific temporal phenotype or morphological stage.

To what extent mRNA and protein expression correlate in different organisms remains largely unknown. For a solid comparison, the data should be generated from the same original sample, and should contain minimal technical variability.5 In our study, we could verify that technical variation was very small and observed a significant increase in correlation due to the sampling approach (Figures S1 and S8C). At the 0h time point, prior to development of the cells, we observed a Spearman correlation across genes of 0.648 (Figure 4B). This is somewhat lower than what has been recently reported for bacteria46 (Spearman = 0.80) and more in line with reports of mammalian cells and other eukaryotes.6 Maybe this is reflective of a more linear relationship between mRNA and protein in bacteria than in eukaryotes. Additionally, some of the dissimilarity might be due to different methods or genes selected for comparison. For example, we observed Spearman correlations as high as 0.70 when considering only differentially expressed genes (Figure 6B).

Besides the genome wide correlation between mRNA and protein levels at distinct time points, we also investigated the correlation between mRNA and protein for individual genes across all time points (Figure S8A). Here, however, we found a relatively low median Pearson correlation of 0.16. One reason for this low correlation is that while the majority of the transcriptome is regulated during development (Table S1), only a fraction of the proteome was clearly developmentally regulated. Thus, if we restrict our analysis to genes that were differentially expressed in both transcriptomic and proteomic datasets, we observe a drastically increased median correlation of 0.54 (Figure S8A). This is similar to what was observed in a xenograft model.47 Another factor contributing to the low correlation between individual mRNAs and their corresponding proteins over time, is the time-lag discussed above.

To conclude, the data presented here enable in-depth study of aggregative multicellularity at both transcript and protein levels, and can constitute a significant resource for comparative studies of other members of Amoebozoa. Notably, we show that the overall correlation between mRNA and protein in D. discoideum at steady-state is rather high, but correlations of individual genes vary, and care should be taken when inferring the presence of proteins from transcriptomic data. We are pleased to refer anyone interested to explore the RNA-protein expression during early development in D. discoideum to the easy-to-use interactive web application https://edelbroek2024.serve.scilifelab.se.

Limitations of the study

The transcriptomic and proteomic study described was performed on one ameba species, D. discoideum, going through aggregative multicellularity. Future work could expand this also to other species with different types of development. In this study, we were interested in early development, the transition from growth to multicellular life. However, to understand whether the same trends are present throughout the full developmental cycle, additional samples can be collected also at later developmental time points. The proteomic data were somewhat limited, which in large can be explained by, e.g., low gene expression (mRNA) for genes where proteins could not be detected. However, other more sensitive techniques for analyzing presence and levels of proteins may give a more complete picture.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Chemicals, peptides, and recombinant proteins	
	
HL5-C	Formedium	Cat#HLC0101	
Phenol stabilized: Chloroform : Isoamyl Alcohol 25:24:1	PanReacAppliChem	Cat#A2279	
TRIzol Reagent	Invitrogen	Cat#15596026	
TURBO DNase	Invitrogen	Cat#AM2238	
DC Protein Assay	BioRad	Cat#5000111	
3kDa centrifugal spin filter	Millipore	Cat#UFC5003	
	
Critical commercial assays	
	
TruSeq stranded mRNA library preparation kit	Illumina	Cat# 20020594/5	
	
Deposited data	
	
Proteomics data	This study	MassIVE: MSV000093620 https://doi.org/10.25345/C5H12VJ75, also available on ProteomeXchange: PXD047669	
Transcriptomics data	This study	GEO: GSE249880https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE249880	
Code for downstream analysis of transcriptomics and proteomics data	This study	Figshare: https://doi.org/10.6084/m9.figshare.25365283	
Proteome at 30min and 8h D. discoideum development	Kelly et al.28	ProteomeXchange: PXD023404	
0-10h D. discoideum developmental transcriptome	Rosengarten et al.19	GEO: GSE61914	
List of transcriptional milestones during D. discoideum development	Katoh-Kurasawa et al.17	https://genome.cshlp.org/content/suppl/2021/07/20/gr.275496.121.DC1/Supplemental_File_S7_JDedit.xlsx	
D. discoideum gene annotation	Singh et al.48	https://doi.org/10.6084/m9.figshare.3384364	
	
Experimental models: Organisms/strains	
	
D. discoideum AX4 wildtype cells	DictyStockCenter	DBS0237637	
	
Software and algorithms	
	
cutadapt v2.10	Martin49	https://github.com/marcelm/cutadapt	
STAR v2.7.5	Dobin et al.50	https://github.com/alexdobin/STAR	
samtools v1.10	Danacek et al.51	https://github.com/samtools/samtools	
featureCounts (subread v2.0.1)	Liao et al.52	https://github.com/ShiLab-Bioinformatics/subread	
FragPipe v20.0	Kong et al.53	https://github.com/Nesvilab/FragPipe	
DESeq2 v.1.41.12	Love et al.23	https://github.com/thelovelab/DESeq2	
topGO v2.54	Alexa, A., and Rahnenfuhrer, J. (2024). topGO: Enrichment Analysis for Gene Ontology. Version 2.56.0.	https://doi.org/10.18129/B9.bioc.topGO	
imputeLCMD v.2.1	Lazar et al.27	https://CRAN.R-project.org/package=imputeLCMD	
Limma v.3.57.11	Ritchie et al.25	https://doi.org/10.18129/B9.bioc.limma	
lmodel2 v.1.7.13	Legendre, P., and Oksanen, J. (2018). lmodel2: Model II Regression.	https://cran.r-project.org/web/packages/lmodel2/	
MOFA v.1.11.0	Velten et al.34	https://www.bioconductor.org/packages/release/bioc/html/MOFA2.html	
piano v.2.17.0	Väremo et al.54	https://doi.org/10.18129/B9.bioc.piano	

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Fredrik Söderbom (fredrik.soderbom@icm.uu.se).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• Complete proteomics data submitted to MassIVE, with accession number MSV000093620, and is linked to ProteomeXchange: PXD047669. The transcriptomics dataset of all 24 sequencing libraries has been submitted to GEO with accession number GSE249880. The links to the repositories are included in the key resources table.

• All code for downstream analysis of the transcriptomics and proteomics datasets can be accessed on GitHub (https://github.com/Bart-Edelbroek/multi-omics-dicty) and also on figshare55 (https://doi.org/10.6084/m9.figshare.25365283) together with the generated figures and tables.

• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Experimental model and study participant details

D. discoideum AX4 wildtype cells (DictyStockCenter ID: DBS0237637) were grown axenically in HL5-C (Formedium) to exponential phase. 3x108 cells were harvested at 400 x g, 5 min, and washed twice in 50 ml KK2 (2.2 g/l KH2PO4, 0.7 g/l K2HPO4). For the 0h time point (not developed), half the cells were harvested as before and stored at -80°C for subsequent processing for transcriptomics library preparation; the other half was harvested and stored at -80°C and later used for proteomics sample preparation. For the other time points, cells were plated on 92mm NN-Agar plates (1.2 g/l KH2PO4, 0.48 g/l Na2HPO4·2H2O, 15 g/l agar) and harvested at the defined time points using Nunc Cell Scrapers (Thermo Fisher), into KK2 buffer; half the plate for transcriptomics and half the plate for proteomics. The cells were treated as described above and the cell pellets were frozen at -80°C until further processed for transcriptomics or proteomics sample preparation.

Method details

Transcriptomics library preparation and sequencing

The frozen cell pellets were dissolved in 1 ml TRIzol Reagent (Invitrogen) and total RNA was prepared according to the user guide, except for an additional 75% EtOH wash of the RNA pellet. Following RNA extraction, 15ug total RNA samples were DNase treated using TURBO DNase (Invitrogen) according to manufacturer’s protocol and purified by phenol/chloroform extraction. 75 μl Phenol stabilized: Chloroform : Isoamyl Alcohol (25:24:1, PanReacAppliChem) was added to 75 μl DNase treated RNA, shaken for 20 s and centrifuged 5 min, 16 000 x g. The upper phase was transferred to new tubes with 187.5 μl EtOH (99%), 7.5 μl 3M Sodium Acetate, 5 μg glycogen, and the RNA was precipitated at -20°C overnight. The RNA was harvested (16 000 x g, 30 min, 4°C), washed with 150 μl 75% EtOH (16 000 x g, 10 min, 4°C), and resuspended in 50 μl RNase free H2O. Sequencing libraries were prepared from 700 ng total RNA using the TruSeq stranded mRNA library preparation kit (Cat# 20020594/5, Illumina Inc.) including polyA selection. The library preparation was performed according to the manufacturers’ protocol (#1000000040498). Libraries were sequenced on the NovaSeq 6000 System (Illumina) on two SP Flowcells, with single reads, 100bp read length (v1 chemistry).

Proteomics sample preparation and LC-MS/MS analysis

The cell pellets were lysed in 150 μL of 1% β-octyl glucopyranoside and 6M urea containing lysis buffer using a sonication probe for 60 seconds (3 mm probe, pulse 1 s, amplitude 30%) according to a standard operating procedure. After homogenization, the samples were incubated for 90 min at 4°C during mild agitation. The lysates were clarified by centrifugation for 10 min (16 000 × g at 4°C). The supernatant containing extracted proteins was collected and further processed. The total protein concentration in the samples was measured using the DC Protein Assay (BioRad) with bovine serum albumin as standard. Aliquots corresponding to 35 μg of proteins were withdrawn for digestion. The proteins were reduced, alkylated, and on-filter digested by trypsin using 3kDa centrifugal spin filter (Millipore, Ireland). The collected peptide filtrate was vacuum centrifuged to dryness using a Speedvac system. The samples were dissolved in 100 μL 0.1% formic acid and further diluted 4 times. For LC-MS/MS analysis, the peptides were separated in reversed-phase on a C18-column with 150 min gradient and electrosprayed on-line to a Q Exactive Plus Orbitrap LC-MS/MS system (Thermo Scientific). Tandem mass spectrometry was performed applying Higher-energy collisional dissociation.

Quantification and statistical analysis

Quantification of transcriptomics

To enable mapping of the sequencing reads, adapters were trimmed using cutadapt v2.10.49 Trimmed reads from different sequencing lanes were pooled and mapped using STAR v2.7.5, allowing a maximum intron size of 2000 bases.50 Mapped reads from both Flowcells were merged with samtools v1.10.51 Reads were assigned to genes with featureCounts, part of the subread v2.0.1 package.52 For both read mapping and counting, the improved D. discoideum gene annotation was used.48 mRNA library preparation and sequencing were performed at SciLifeLab Uppsala.

Quantification of proteomics

Label free quantification (LFQ) of the raw data was performed using FragPipe v20.0 (https://fragpipe.nesvilab.org/), which is powered by MSFragger.53 Analysis was performed with oxidation and lysine ubiquitination specified as variable modifications. Up to 3 missed cleavages were allowed. PSM validation performed with Percolator,56 and protein inference with ProteinProphet.57 Data is filtered at 1% FDR at the PSM, ion, peptide, and protein levels. Site localization with PTM-Prophet. For quantification, a minimum of 1 ion was required for MaxLFQ determination with IonQuant, using match between runs.58

Transcriptomics statistical analysis

Counts from transcripts encoding the same protein were summed, and transcripts not encoding proteins were discarded, to allow for analysis of protein-coding transcripts and enable downstream comparison to protein data. Differentially expressed genes over time from mRNA-seq were identified with DESeq2 v.1.41.12, using a likelihood ratio test to compare a model where gene expression is explained by developmental time to a null model of constant expression.23 Genes with an FDR-adjusted p-value below 0.01 were designated as differentially expressed. Normalized, transformed count data was extracted using variance stabilizing transformations, and shrunken log fold changes of differentially expressed genes were calculated with Approximate Posterior Estimation for generalized linear model.59 Counts from 0h to 10h time points from Rosengarten et al., were processed in the same manner.19,60 Genes plotted in heatmaps were hierarchically clustered based on their log fold changes or z-scores. Gene set enrichment for GO-terms was performed using topGO v2.54 with the weight01 algorithm and using Fisher's exact test to determine statistical significance (https://doi.org/10.18129/B9.bioc.topGO).

Proteomics statistical analysis

Protein quantification is based on MaxLFQ values from FragPipe. Values were imputed for proteins that were quantified in all biological replicates at a given time point, but where values were missing at other time points. Imputation was performed using a probabilistic minimum from the imputeLCMD v.2.1 R package.27 Differentially expressed proteins (FDR-adjusted p-value 0.01) were identified with Limma v.3.57.11 by fitting linear models, with empirical Bayes smoothing.25,61 The data by Kelly et al.,28,62 was imputed and processed identically for comparison of differentially expressed genes. Clustering and GO-term analysis was performed as for the mRNA data.

Integrative analysis

Genes which were quantified in both the transcriptomics and proteomics datasets, were utilized for integrative analysis. To enable comparison of the mRNA and protein levels, the values of each replicate, for each dataset, were divided by the total sum of values for that replicate such that the scaled values sum to 1. For across genes correlation at a single time point, the mean mRNA and protein levels were calculated from the biological replicates. Linear regression was calculated with ranged major axes using lmodel2 v.1.7.13 (https://cran.r-project.org/web/packages/lmodel2/). For per-gene correlations, Pearson correlations were calculated for each gene with all replicates available from both transcriptomics and proteomics.

Multi omics factor analysis was performed with MEFISTO, using MOFA v.1.11.0 on using data from all time points.34 For analysis based on Factor 1, genes were selected with a Factor 1 loading at both mRNA and protein modalities above 0.45 or below -0.45. Gene set enrichment analysis of Factor 2 mRNA was performed based on the Factor 2 gene weights using piano v.2.17.0.54

For time lag analysis, the Spearman across-genes correlation was calculated for each transcriptomics time point with each proteomics time point, either with all common quantified genes, or those that were differentially expressed in both datasets. To calculate the ratios of protein levels to mRNA levels, the normalized protein level was divided by the normalized mRNA level for each gene. Differentially expressed genes at the mRNA modality, which have a fold change above 2 at the 10h time point compared to the 0h time point, were identified as upregulated, and those with a fold change below 0.5 as downregulated.

Additional resources

The multi-omics data presented in this study are available in an interactive web app for ease of use: https://edelbroek2024.serve.scilifelab.se.

Supplemental information

Document S1. Figures S1–S10

Table S1. Transcriptomics: differential expression and quantification, related to Figure 2 and Table 1

Log2 fold changes (logFC) are relative to the 0h time point and shown with adjusted p-value. The logFC calculations are based on all replicates (a-d), with their individual quantification shown on the right in reads per million (rpm). All gene quantifications for each replicate (i.e., each column) therefore sum to 1e6. Synonyms are alternative gene names, if available.

Table S2. Significant GO-terms for differentially expressed transcript clusters, related to Figure 2C

Full set of significantly enriched GO-terms (Fisher’s p-value <0.01) for each of the clusters in Figure 2C. For each GO-term, the total number of genes in the set with the annotated term is shown, as well as the number of genes with that annotation in the cluster. The number of expected genes, given the size of the cluster, and the relative enrichment are shown right.

Table S3. Proteomics: differential expression and quantification, related to Figure 3 and Table 1

Log2 fold changes (logFC) are relative to the 0h time point and shown with adjusted p-value. The logFC calculations are based on all replicates (a-c), with their individual quantification shown on the right. For each replicate, the label free protein quantifications were normalized to a million (LFQpm). All gene quantifications for each replicate (i.e., each column) therefore sum to 1e6. Synonyms are alternative gene names, if available.

Table S4. Significant GO-terms for differentially expressed protein clusters, related to Figure 3C

Full set of significantly enriched GO-terms (Fisher’s p-value <0.01) for each of the clusters in Figure 3C. For each GO-term, the total number of genes in the set with the annotated term is shown, as well as the number of genes with that annotation in the cluster. The number of expected genes, given the size of the cluster, and the relative enrichment are shown right.

Table S5. Regulation over time of differentially expressed genes in both datasets, related to Figure 4A

Log2 fold changes (logFC) are relative to the 0h time point and shown with adjusted p-value for both the mRNA and protein data. Only genes are included which are differentially expressed over development in both datasets (adj. p-value <0.01). Synonyms are alternative gene names, if available.

Table S6. Per gene correlation of RNA and protein for all common genes, related to Figure 4E

The calculated Pearson correlation coefficient (r) is indicated for each gene, with p-value. The correlation coefficient is calculated on the log2 transformed, normalized RNA and protein quantification. Samples for which both transcriptomics and proteomics were performed are included, and only the 3604 genes that were quantified in both datasets (Table 1) For each included replicate, the mRNA and protein quantifications were normalized, so that all gene quantifications for each replicate (i.e., each column) sum to 1e6. Synonyms are alternative gene names, if available.

Table S7. Genes with factor 1 regulation (linear up/down) in both the transcriptomics and proteomics datasets, related to Figure 5C

All genes with an absolute factor 1 value above 0.45 in both datasets are included. The manner in which the genes are regulated in both datasets is classified in the regulation column.

Table S8. Gene Set Enrichment Analysis, associating MEFISTO latent factors 1–3 to GO-terms, related to Figure 5

Results are shown for each factor and separately for protein and RNA. For each GO-term the following numbers are given: a) the number of genes, b) the test statistic (enrichment score where 1 represents high accordance with the factor and −1 high accordance with the inverse pattern) and c) the p-value, adjusted for multiple hypothesis testing. Only GO-terms with at least 10 detected genes and with p-value <0.001 are included.

Acknowledgments

Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala. The facility is part of the National Genomics Infrastructure (NGI) Sweden and Science for Life Laboratory. The SNP&SEQ Platform is also supported by the 10.13039/501100004359 Swedish Research Council and the 10.13039/501100004063 Knut and Alice Wallenberg Foundation . Proteomics was performed at the MS-based proteomics facility platform at Uppsala University. Dr. Anna Widgren is acknowledged for her support with the analysis. We would like to thank Jonas Kjellin, Johan Reimegård, and Andrea Hinas for critical reading of the manuscript. The computations and data handling were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at Uppsala University, partially funded by the 10.13039/501100004359 Swedish Research Council through grant agreement no. 2022-06725 . J.O.W. is financially supported by the 10.13039/501100004063 Knut and Alice Wallenberg Foundation as part of the National Bioinformatics Infrastructure Sweden at SciLifeLab. This work was supported by grant no. 2021-05793 from 10.13039/501100004359 Swedish Research Council (Vetenskapsrådet) and grant no. CTS 18:381 from Carl Tryggers Stiftelse to F.S.

Author contributions

B.E., J.O.W., and F.S. designed the study. B.E. grew cell cultures and prepared samples for transcriptomics and proteomics. J.B. led the proteomics data generation. B.E. and J.O.W. analyzed the data and prepared figures. B.E., J.O.W., and F.S. drafted the manuscript. All authors read and approved the final manuscript.

Declaration of interests

The authors declare no competing interests.

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110659.
==== Refs
References

1 Crick F. Central Dogma of Molecular Biology Nature 227 1970 561 563 10.1038/227561a0 4913914
2 Timp W. Timp G. Beyond mass spectrometry, the next step in proteomics Sci. Adv. 6 2020 eaax8978 10.1126/sciadv.aax8978
3 Deutsch E.W. Bandeira N. Perez-Riverol Y. Sharma V. Carver J.J. Mendoza L. Kundu D.J. Wang S. Bandla C. Kamatchinathan S. The ProteomeXchange consortium at 10 years: 2023 update Nucleic Acids Res. 51 2023 D1539 D1548 10.1093/nar/gkac1040 36370099
4 Edgar R. Domrachev M. Lash A.E. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository Nucleic Acids Res. 30 2002 207 210 10.1093/nar/30.1.207 11752295
5 Liu Y. Beyer A. Aebersold R. On the Dependency of Cellular Protein Levels on mRNA Abundance Cell 165 2016 535 550 10.1016/j.cell.2016.03.014 27104977
6 Buccitelli C. Selbach M. mRNAs, proteins and the emerging principles of gene expression control Nat. Rev. Genet. 21 2020 630 644 10.1038/s41576-020-0258-4 32709985
7 Schaefke B. Sun W. Li Y.-S. Fang L. Chen W. The evolution of posttranscriptional regulation WIREs RNA 9 2018 e1485 10.1002/wrna.1485
8 Becker K. Bluhm A. Casas-Vila N. Dinges N. Dejung M. Sayols S. Kreutz C. Roignant J.-Y. Butter F. Legewie S. Quantifying post-transcriptional regulation in the development of Drosophila melanogaster Nat. Commun. 9 2018 4970 10.1038/s41467-018-07455-9 30478415
9 Grün D. Kirchner M. Thierfelder N. Stoeckius M. Selbach M. Rajewsky N. Conservation of mRNA and Protein Expression during Development of C. elegans Cell Rep. 6 2014 565 577 10.1016/j.celrep.2014.01.001 24462290
10 Peshkin L. Wühr M. Pearl E. Haas W. Freeman R.M. Gerhart J.C. Klein A.M. Horb M. Gygi S.P. Kirschner M.W. On the Relationship of Protein and mRNA Dynamics in Vertebrate Embryonic Development Dev. Cell 35 2015 383 394 10.1016/j.devcel.2015.10.010 26555057
11 Ponnala L. Wang Y. Sun Q. van Wijk K.J. Correlation of mRNA and protein abundance in the developing maize leaf Plant J. 78 2014 424 440 10.1111/tpj.12482 24547885
12 Kawabe Y. Du Q. Schilde C. Schaap P. Evolution of multicellularity in Dictyostelia Int. J. Dev. Biol. 63 2019 359 369 10.1387/ijdb.190108ps 31840775
13 Brown M.W. Silberman J.D. The Non-dictyostelid Sorocarpic Amoebae Romeralo M. Baldauf S. Escalante R. Dictyostelids: Evolution, Genomics and Cell Biology 2013 Springer 219 242 10.1007/978-3-642-38487-5_12
14 Glöckner G. Social Amoebae and Their Genomes: On the Brink to True Multicellularity Ruiz-Trillo I. Nedelcu A.M. Evolutionary Transitions to Multicellular Life: Principles and mechanisms 2015 Springer Netherlands 363 376 10.1007/978-94-017-9642-2_17
15 Bozzaro S. The past, present and future of Dictyostelium as a model system Int. J. Dev. Biol. 63 2019 321 331 10.1387/ijdb.190128sb 31840772
16 Devreotes P. Dictyostelium discoideum: a Model System for Cell-Cell Interactions in Development Science 245 1989 1054 1058 10.1126/science.2672337 2672337
17 Katoh-Kurasawa M. Hrovatin K. Hirose S. Webb A. Ho H.-I. Zupan B. Shaulsky G. Transcriptional milestones in Dictyostelium development Genome Res. 31 2021 1498 1511 10.1101/gr.275496.121 34183452
18 Nichols J.M. Antolović V. Reich J.D. Brameyer S. Paschke P. Chubb J.R. Cell and molecular transitions during efficient dedifferentiation Elife 9 2020 e55435 10.7554/eLife.55435
19 Rosengarten R.D. Santhanam B. Fuller D. Katoh-Kurasawa M. Loomis W.F. Zupan B. Shaulsky G. Leaps and lulls in the developmental transcriptome of Dictyostelium discoideum BMC Genom. 16 2015 294 10.1186/s12864-015-1491-7
20 Wang S.Y. Pollina E.A. Wang I.-H. Pino L.K. Bushnell H.L. Takashima K. Fritsche C. Sabin G. Garcia B.A. Greer P.L. Greer E.L. Role of epigenetics in unicellular to multicellular transition in Dictyostelium Genome Biol. 22 2021 134 10.1186/s13059-021-02360-9 33947439
21 Westbrook E.R. Lenn T. Chubb J.R. Antolović V. Collective signalling drives rapid jumping between cell states Development 150 2023 dev201946 10.1242/dev.201946
22 Kessin R.H. Dictyostelium: Evolution, Cell Biology, and the Development of Multicellularity 2001 Cambridge University Press
23 Love M.I. Huber W. Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 Genome Biol. 15 2014 550 10.1186/s13059-014-0550-8 25516281
24 Basu S. Fey P. Pandit Y. Dodson R. Kibbe W.A. Chisholm R.L. dictyBase 2013: integrating multiple Dictyostelid species Nucleic Acids Res. 41 2013 D676 D683 10.1093/nar/gks1064 23172289
25 Ritchie M.E. Phipson B. Wu D. Hu Y. Law C.W. Shi W. Smyth G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies Nucleic Acids Res. 43 2015 e47 10.1093/nar/gkv007 25605792
26 Söderbom F. Loomis W.F. Cell–cell signaling during Dictyostelium development Trends Microbiol. 6 1998 402 406 10.1016/S0966-842X(98)01348-1 9807784
27 Lazar C. Gatto L. Ferro M. Bruley C. Burger T. Accounting for the Multiple Natures of Missing Values in Label-Free Quantitative Proteomics Data Sets to Compare Imputation Strategies J. Proteome Res. 15 2016 1116 1125 10.1021/acs.jproteome.5b00981 26906401
28 Kelly B. Carrizo G.E. Edwards-Hicks J. Sanin D.E. Stanczak M.A. Priesnitz C. Flachsmann L.J. Curtis J.D. Mittler G. Musa Y. Sulfur sequestration promotes multicellularity during nutrient limitation Nature 591 2021 471 476 10.1038/s41586-021-03270-3 33627869
29 Mesquita A. Cardenal-Muñoz E. Dominguez E. Muñoz-Braceras S. Nuñez-Corcuera B. Phillips B.A. Tábara L.C. Xiong Q. Coria R. Eichinger L. Autophagy in Dictyostelium: Mechanisms, regulation and disease in a simple biomedical model Autophagy 13 2017 24 40 10.1080/15548627.2016.1226737 27715405
30 Fischer S. Eichinger L. Dictyostelium discoideum and autophagy – a perfect pair Int. J. Dev. Biol. 63 2019 485 495 10.1387/ijdb.190186LE 31840786
31 Loomis W.F. Genetic control of morphogenesis in Dictyostelium Dev. Biol. 402 2015 146 161 10.1016/j.ydbio.2015.03.016 25872182
32 Csárdi G. Franks A. Choi D.S. Airoldi E.M. Drummond D.A. Accounting for Experimental Noise Reveals That mRNA Levels, Amplified by Post-Transcriptional Processes, Largely Determine Steady-State Protein Levels in Yeast PLoS Genet. 11 2015 e1005206 10.1371/journal.pgen.1005206
33 Wang D. Eraslan B. Wieland T. Hallström B. Hopf T. Zolg D.P. Zecha J. Asplund A. Li L.-H. Meng C. A deep proteome and transcriptome abundance atlas of 29 healthy human tissues Mol. Syst. Biol. 15 2019 e8503 10.15252/msb.20188503
34 Velten B. Braunger J.M. Argelaguet R. Arnol D. Wirbel J. Bredikhin D. Zeller G. Stegle O. Identifying temporal and spatial patterns of variation from multimodal data using MEFISTO Nat. Methods 19 2022 179 186 10.1038/s41592-021-01343-9 35027765
35 Subramanian A. Tamayo P. Mootha V.K. Mukherjee S. Ebert B.L. Gillette M.A. Paulovich A. Pomeroy S.L. Golub T.R. Lander E.S. Mesirov J.P. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles Proc. Natl. Acad. Sci. USA 102 2005 15545 15550 10.1073/pnas.0506580102 16199517
36 Ishikawa-Ankerhold H.C. Müller-Taubenberger A. Actin assembly states in Dictyostelium discoideum at different stages of development and during cellular stress Int. J. Dev. Biol. 63 2019 417 427 10.1387/ijdb.190256am 31840780
37 Fournier M.L. Paulson A. Pavelka N. Mosley A.L. Gaudenz K. Bradford W.D. Glynn E. Li H. Sardiu M.E. Fleharty B. Delayed Correlation of mRNA and Protein Expression in Rapamycin-treated Cells and a Role for Ggc1 in Cellular Sensitivity to Rapamycin Mol. Cell. Proteomics 9 2010 271 284 10.1074/mcp.M900415-MCP200 19955083
38 Casas-Vila N. Bluhm A. Sayols S. Dinges N. Dejung M. Altenhein T. Kappei D. Altenhein B. Roignant J.-Y. Butter F. The developmental proteome of Drosophila melanogaster Genome Res. 27 2017 1273 1285 10.1101/gr.213694.116 28381612
39 Du Q. Kawabe Y. Schilde C. Chen Z.H. Schaap P. The Evolution of Aggregative Multicellularity and Cell–Cell Communication in the Dictyostelia J. Mol. Biol. 427 2015 3722 3733 10.1016/j.jmb.2015.08.008 26284972
40 Knoll A.H. The Multiple Origins of Complex Multicellularity Annu. Rev. Earth Planet Sci. 39 2011 217 239 10.1146/annurev.earth.031208.100209
41 Loomis W.F. Comparative Genomics of the Dictyostelids Eichinger L. Rivero F. Dictyostelium discoideum Protocols Methods in Molecular Biology 2013 Humana Press 39 58 10.1007/978-1-62703-302-2_3
42 Sucgang R. Kuo A. Tian X. Salerno W. Parikh A. Feasley C.L. Dalin E. Tu H. Huang E. Barry K. Comparative genomics of the social amoebae Dictyostelium discoideum and Dictyostelium purpureum Genome Biol. 12 2011 R20 10.1186/gb-2011-12-2-r20 21356102
43 Heidel A.J. Lawal H.M. Felder M. Schilde C. Helps N.R. Tunggal B. Rivero F. John U. Schleicher M. Eichinger L. Phylogeny-wide analysis of social amoeba genomes highlights ancient origins for complex intercellular communication Genome Res. 21 2011 1882 1891 10.1101/gr.121137.111 21757610
44 Gross J.D. Developmental decisions in Dictyostelium discoideum Microbiol. Rev. 58 1994 330 351 10.1128/mr.58.3.330-351.1994 7968918
45 Cheng Z. Teo G. Krueger S. Rock T.M. Koh H.W.L. Choi H. Vogel C. Differential dynamics of the mammalian mRNA and protein expression response to misfolding stress Mol. Syst. Biol. 12 2016 855 10.15252/msb.20156423 26792871
46 Balakrishnan R. Mori M. Segota I. Zhang Z. Aebersold R. Ludwig C. Hwa T. Principles of gene regulation quantitatively connect DNA to RNA and proteins in bacteria Science 378 2022 eabk2066 10.1126/science.abk2066
47 Koussounadis A. Langdon S.P. Um I.H. Harrison D.J. Smith V.A. Relationship between differentially expressed mRNA and mRNA-protein correlations in a xenograft model system Sci. Rep. 5 2015 10775 10.1038/srep10775
48 Singh R. Lawal H.M. Schilde C. Glöckner G. Barton G.J. Schaap P. Cole C. Improved annotation with de novo transcriptome assembly in four social amoeba species BMC Genom. 18 2017 120 10.1186/s12864-017-3505-0
49 Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads EMBnet. j. 17 2011 10 12 10.14806/ej.17.1.200
50 Dobin A. Davis C.A. Schlesinger F. Drenkow J. Zaleski C. Jha S. Batut P. Chaisson M. Gingeras T.R. STAR: ultrafast universal RNA-seq aligner Bioinformatics 29 2013 15 21 10.1093/bioinformatics/bts635 23104886
51 Danecek P. Bonfield J.K. Liddle J. Marshall J. Ohan V. Pollard M.O. Whitwham A. Keane T. McCarthy S.A. Davies R.M. Li H. Twelve years of SAMtools and BCFtools GigaScience 10 2021 giab008 10.1093/gigascience/giab008
52 Liao Y. Smyth G.K. Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features Bioinformatics 30 2014 923 930 10.1093/bioinformatics/btt656 24227677
53 Kong A.T. Leprevost F.V. Avtonomov D.M. Mellacheruvu D. Nesvizhskii A.I. MSFragger: ultrafast and comprehensive peptide identification in mass spectrometry–based proteomics Nat. Methods 14 2017 513 520 10.1038/nmeth.4256 28394336
54 Väremo L. Nielsen J. Nookaew I. Enriching the gene set analysis of genome-wide data by incorporating directionality of gene expression and combining statistical hypotheses and methods Nucleic Acids Res. 41 2013 4378 4391 10.1093/nar/gkt111 23444143
55 Edelbroek B. Westholm J.O. Source Code for: Multi-Omics Characterization of Aggregative Multicellularity 2024 figshare 10.6084/m9.figshare.25365283.v1
56 Käll L. Canterbury J.D. Weston J. Noble W.S. MacCoss M.J. Semi-supervised learning for peptide identification from shotgun proteomics datasets Nat. Methods 4 2007 923 925 10.1038/nmeth1113 17952086
57 da Veiga Leprevost F. Haynes S.E. Avtonomov D.M. Chang H.-Y. Shanmugam A.K. Mellacheruvu D. Kong A.T. Nesvizhskii A.I. Philosopher: a versatile toolkit for shotgun proteomics data analysis Nat. Methods 17 2020 869 870 10.1038/s41592-020-0912-y 32669682
58 Yu F. Haynes S.E. Nesvizhskii A.I. IonQuant Enables Accurate and Sensitive Label-Free Quantification With FDR-Controlled Match-Between-Runs Mol. Cell. Proteomics 20 2021 100077 10.1016/j.mcpro.2021.100077
59 Zhu A. Ibrahim J.G. Love M.I. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences Bioinformatics 35 2019 2084 2092 10.1093/bioinformatics/bty895 30395178
60 GEO Accession viewer https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE61914.
61 Phipson B. Lee S. Majewski I.J. Alexander W.S. Smyth G.K. Robust Hyperparameter Estimation Protects Against Hypervariable Genes And Improves Power To Detect Differential Expression Ann. Appl. Stat. 10 2016 946 963 10.1214/16-AOAS920 28367255
62 ProteomeXchange Dataset PXD023404 https://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD023404.
