To study the proteomes and the phosphoproteomes of cell populations in the early embryo, we used mutants in which all cells in the embryo represent only one subset of the cell types present in the wild type embryo. Because we were interested in cell behaviors that affect the first step of gastrulation, which is driven by differences in cell behavior along the DV axis, we used mutants for genes of the DV patterning pathway. The embryos were derived from mothers mutant for the genes gastrulation defective (gd), Toll (Tl), or Serpin27A (spn27A). We chose those alleles that cause the strongest dorsalization, lateralization and ventralization of the dorso-ventral axis as judged by cuticle phenotypes and changes in gene expression patterns. To generate dorsalized embryos we used the gd⁹ allele, reported to generate the strongest dorsalization without affecting the length of the embryo (Konrad et al., 1988; Ponomareff et al., 2001). Mothers transheterozygous for the hypomorphic mutations in Tl^rm9 and Tl^rm10 were used to produce lateralized embryos, in which the entire dorso-ventral axis forms neuroectoderm (Stathopoulos et al., 2002; Anderson et al., 1985; Cowden and Levine, 2003). Ventralized embryos were generated in two different ways using mutations in Tl and spn27A: one was to make mothers transheterozygous for the dominant Tl gain-of-function allele Tl^10B (Anderson et al., 1985; Schneider et al., 1991) and deficiency Df(3R)ro80b, which uncovers the Tl locus; the other was to use mothers that were transheterozygous for spn27A^ex, an amorphic mutation (complete excision) of spn27A (Ligoxygakis et al., 2003), in combination with deficiency Df(2L)BSC7, which uncovers the spn27A locus. To confirm that the embryos produced by these mothers represented the dorsal, lateral and ventral cell populations, we analyzed the expression patterns of D-V fate determining genes (Figure 1A, Figure 1—figure supplement 1B, Supplementary file 1). ‘Lateralized’ and ‘dorsalized’ embryos from Tl^rm10/Tl^rm9 and gd⁹/gd⁹ mothers expressed neither twist nor snail, whereas ventralized embryos from Toll^10B/def and spn27A^ex/def mothers expressed twist and snail around their entire circumference in the trunk region (Figure 1B, Supplementary file 1). In embryos from Tl^rm10/Tl^rm9 mothers, sog expression expanded dorsally and ventrally, whereas dpp expression expanded ventrally (Figure 1B, Supplementary file 1). These expression patterns showed some variation and were not entirely homogeneous: ventralized embryos often had a gap in snail expression in a small dorsal-anterior domain around the procephalic furrow. In this region, we detected sog expression instead, suggesting ventralized embryos retain some cells with a neuroectodermal fate in a restricted area of the embryo (Figure 1B, sog probe).

Figure 1 with 1 supplement see all

Download asset Open asset

Because we wanted to use these mutants to identify proteins that reflect or control differential cell behavior, it was important to ascertain that the cells in these mutants recapitulate faithfully the biological qualities of the corresponding cell populations in the wild type embryo (Rauzi et al., 2015; Kölsch et al., 2007), specifically of the localisation of the adherens junctions and the cortical actomyosin meshwork. We find that, as in the mesoderm of wildtype embryos, the adherens junctions (as visualized by immunostaining for Armadillo/β-Catenin; Figure 1C) relocalize apically in the ventralized mutants, but remain apico-lateral in lateralized and shift slightly more basally in dorsalized mutants, again mirroring the morphology of lateral and dorsal regions of the wildtype embryo (Figure 1C, Figure 1—figure supplement 1A). Similarly, the apical actomyosin network, which we characterized in living embryos expressing a fluorescently tagged myosin light chain (sqh-mCherry, Figure 1E and F) forms a pulsatile apical network in ventralized embryos, whereas myosin accumulates at cell junctions in lateralized embryos, and dorsalized embryos dissolve the loose apical actomyosin of the early blastoderm (Figure 1E and F).

In summary, in terms of marker gene expression and cell behavior, the cells in these mutants resemble the corresponding embryonic cell populations of a wild type embryo, showing that these mutant cell populations are good sources of material to analyze the proteomic and phosphoproteomic composition of the natural cell populations at the onset of gastrulation.

To study the proteins and the phosphosites that might be relevant for cell behavior during gastrulation, we focused on a narrow developmental time window for sample collection. We synchronized egg collections and manually collected embryos from wild type and mutant mothers aged for 165–180 minutes after egg deposition at 25 °C (Stage 6, see Methods and Figure 1—figure supplement 1A and B). We analyzed their peptides and phospho-peptides with unbiased label-free quantification (LFQ) and SILAC (Stable Isotope Labeling with Amino acids in Culture Xu et al., 2012; Sury et al., 2010; Nolte et al., 2015, Figure 1—figure supplement 1C–E).

In the proteomic analyses, we identified 6111 protein groups (to which, for the sake of simplicity, we will refer simply as proteins; Supplementary file 3) across all genotypes. A total of 5883 of these were detected in wild type embryos (Figure 2A), exceeding previously reported number identified by proteomic approaches in early Drosophila embryogenesis (Casas-Vila et al., 2017). Most were detected in all genotypes (Figure 2B). The small number (519/6111) with restricted detection included the DV fate determinants Doc1, Snail, Twist, and dMyc (Figure 1A and B). The phosphoproteomic analysis identified 6259 phosphosites distributed over 1847 proteins (Figure 2C, Supplementary file 4). Only 73% of phosphosites were found across all genotypes (Figure 2D). Twenty-eight percent of the proteins (1699/6111) and 9% of the phosphosites (573/6259) differed significantly in an ANOVA test across all five populations (wild type and four mutants; permutation-based FDR <0.1, s0=0.1).

Figure 2 with 1 supplement see all

Download asset Open asset

We determined the degree of experimental variability by generating correlation matrices both for the proteome and the phosphoproteome. For the proteome, the replicates from the same genotypes clustered together (Figure 2E). For the phosphoproteome, the first replicate of each genotype was separated from the other two replicates (Figure 2F). We nevertheless kept all replicates for further analyses because it was impossible to determine the experimental source for this variation.

The enrichment for proteins or phosphosites in the mutant genotypes over the wild-type ranged from near-zero to 100 fold (Figure 2G and H) with about half changing by less than 1.5-fold. The fold-changes for the mesodermal fate determinants Snail (7.7 fold) and Twist (14.6 fold) measured in spn27A^ex/def were the largest positive fold-changes among the DV fate determinants (Figure 2G).

To test if the recovered protein populations represented the cell populations in the embryo, we analyzed whether they contained known marker proteins. We first looked for the protein product of the gene that was mutated in each group of embryos. We detected both Toll and Spn27A, and each of them was reduced in abundance in the respective mutant embryos. (Figure 2—figure supplement 1A).

Proteins that are known to be expressed differentially along the DV axis (Figure 1A) were more abundant in the appropriate genotypes: Snail, Twist, Mdr49, Traf4, and CG4500 in ventralized embryos; the pro-neuroectodermal (lateral) factor Sog in lateralized embryos; pro-ectodermal (dorsal) factors Zen, Doc1, Dtg, Net, and Egr in dorsalized embryos (Figure 1D, Figure 2—figure supplement 1C, D; p values for all comparisons in all figures are summarized in Supplementary file 2). Known ventral-specific proteins (Snail, Twist, Mdr49, and CG4500) were more strongly upregulated in serpin27A embryos than in Toll^10B, and most dorsal-specific proteins (e.g. Egr, Zen, Sdt, Net, and Ptr) were more strongly downregulated.

We also recovered known phosphosites in proteins acting in the early embryo. This included the serine 871 phosphosite in Toll (Zhai et al., 2008), and serines 463, 467, and 468 in Cactus that have been shown to be phosphorylated by CKII (Liu et al., 1997; Figure 2—figure supplement 1E, Supplementary file 2). Toll, a known target of the Ser/Thr kinase Pelle (Shen and Manley, 1998), was phosphorylated on serine-871, and this phosphosite was more abundant in ventralized embryos (Figure 2—figure supplement 1B, Supplementary file 2). Phosphosites in proteins associated with the Rho pathway will be discussed below. In summary, the proteomic and phosphoproteomic screen correctly identified known and differentially expressed proteins and phosphosites.

Our knowledge of the genetics of the dorso-ventral patterning system gives us a biological criterion that we can use to analyze the data in a stringent manner. We know that region-specific protein sets should change in concert in a well-controlled manner in all of the mutants. Rather than simply looking for individual pair-wise changes, we can, and must, therefore impose this as an additional criterion in determining any potential proteins of interest: each protein must change in a manner that ‘makes sense’ genetically.

The assumption that each mutant represents a defined region of the embryo makes a simple prediction for the expected outcome of the measurements: if one adds up the quantities of protein found in the mutants representing the ventral, lateral and dorsal region (normalized to the fraction of the embryo the corresponding region occupies), the sum should equal the amount of protein in the wildtype. For example, the transcription factor Snail is expressed only in the prospective mesoderm (ventral domain) in the wildtype embryo, but practically in all cells of ventralized embryos, and nowhere in lateralized and dorsalized embryos (Figure 3—figure supplement 1A). This is also reflected correctly in the proteomes: Snail is absent in the dorsalized and lateralized proteomes, and its level is higher in the proteomes from the ventralized embryos (Figure 1D, Supplementary file 2). Thus, Snail shows an ideal behavior in each of the DV mutant genotypes because it recapitulates the expression of Snail in the corresponding domains of a wild type embryo.

We developed a ‘linear model’ that is based on this additional genetic criterion, which we then used to evaluate simultaneously all mutant proteomes. We calculated for each protein the sum of its normalized quantities in the mutants and compared that sum to its abundance in the wild-type embryo. In the absence of experimental measurements for the sizes of each of the areas in the embryo (except for the mesoderm), we determined in an analytical manner (see Methods: Development of a linear model) the optimal values for the proportions occupied by the dorsal and lateral populations in the wildtype embryo, and used these to calculate the ‘theoretical’ wildtype value ^twt_ProtX for each protein:

where D, L, and V are the measured abundance in the three mutant populations.

The deviation for each protein from the experimentally measured wiltype amount ^mwt_ProtX is the ratio ^twt_ProtX/ ^mwt_ProtX. When we apply this analysis to one of the marker proteins, Snail, we arrive at a deviation value of 1.07 in the case where Toll^10B is used to represent ‘ventral’. This shows that for this protein, the mutants represent the regional distribution in the wildtype very well. If we do the calculation with spn27A as the ventral population the deviation value for Snail is 2.91, which indicates that this mutant genotype may over-represent the ventral population.

For our further calculations, we use the log2 of this ratio, i.e. Deviation_ProtX = log2 (^twt_ProtX/^mwt_ProtX). We found that the majority of the proteins had a deviation around zero, that is the calculated value corresponds to the measured value in the wildtype (Figure 3A). This would in fact be expected for any protein that is expressed ubiquitously in the wild type (such as the non-regulated maternal proteins) and should therefore be present in equal amounts in all genotypes. But even the proteins that show significant differences between at least two mutant conditions, that is the ANOVA significant subset, also fall into the range between –0.5 and +0.5, that is less than 1.4-fold deviation (Figure 3A, Supplementary file 6). This shows that the majority of proteins fit the linear model, which in turn indicates that the mutant values are good representations of protein abundance in the corresponding domains of a wild type embryo. The proteins with the most extreme deviations (more than twofold) did not come from any well-defined class of proteins, but represented a wide range of ontologies (Figure 3—figure supplement 1D, Supplementary file 13).

Figure 3 with 1 supplement see all

Download asset Open asset

To find the proteins that function in a tissue-specific manner during gastrulation we sorted the proteins into sets that change in concert in all of the mutants in the predicted, 'correct' manner, again using the assumptions that underlie this study, i.e. that the changes in the different mutants would be expected to correlate with each other in logical ways, as described above.

Rather than focusing only on the proteins that the ANOVA had shown as significantly modulated, we included in this analysis all proteins that were detectable in the wildtype (5883/6111), even if they were undetectable in one or more mutant populations. This allows us to include the important group of proteins that show a ‘perfect’ behavior, like Twist, Snail, or Doc1, in that they are undetectable in the mutants that correspond to the regions in the normal embryo where these genes are not expressed.

We used hierarchical clustering to identify the sets of proteins that change in the mutants in the same manner. For this analysis, we ignored the quantitative extent of the changes in the mutants versus the wildtype, and only focused on the direction of change if a threshold of |0.5 log2 fold change| is exceeded (see Methods). We clustered the set of 3398/6111 proteins which excluded those proteins for which the changes between the mutants and the wildtype were either all in the same direction or below the threshold.

Based on known gene expression patterns along the DV axis in the wildtype one would expect six clusters (Figure 3C): expression restricted to ventral (snail), lateral (sog), or dorsal (dpp), or expression across two domains, that is dorsal and lateral (grh or std), dorsal and ventral (ama) or lateral and ventral (neur). However, in addition to these clusters (marked as 1/D, 2 /L, 5 /V, 6 /DL, 9/DV and 12/LV in Figure 3B and C, Supplementary file 7) the clustering yielded a further eight clusters (Figure 3B and C). This results from the surprising difference between the two ventralising genotypes (Figure 3E). A large number of proteins change in abundance in one but not the other mutant.

Most of the marker proteins were found in their proper predicted classes (Figure 3D). Among those allocated to clusters where the two ventral mutants differed in their behavior, there was no general rule as to which of the two ventral mutants represented the correct value. For example, both Heartless and Net are expressed in the mesoderm and also on the dorsal side of the embryo, but Heartless was seen with increased abundance only in serpin27A embryos, and Net only in Toll^10B embryos (Figure 3D). Similarly, for genes that are excluded from the mesoderm, that is expressed in dorsal and lateral regions, some scored as present in lower abundance in serpin27A (e.g. crb), whereas others were reduced in Toll^10B (eg. numb). We will return to the difference between Toll^10B and spn27A below.

Protein levels can be regulated post-translationally, and RNA and protein expression levels do not necessarily correlate strongly during development (Becker et al., 2018). However, the regional distribution of proteins in the early Drosophila embryo is thought to be achieved mainly through transcriptional regulation (Lasko, 2020; Ing-Simmons et al., 2021). We therefore investigated how well the proteomes reflected known dorso-ventral modulation of gene expression.

We first looked for genes whose RNA expression patterns are reported in the BDGP in situ database (https://insitu.fruitfly.org/cgi-bin/ex/insitu.pl; Hammonds et al., 2013; Tomancak et al., 2002; Tomancak et al., 2007) to compare those with ventral expression in this data set and ours. We extracted all those genes that carry the labels ‘mesoderm’, ‘trunk mesoderm’, or ‘head mesoderm’ in BDGP (which are not mutually exclusive). A total of 107 of the resulting set of 109 genes had their proteins detected in our analyses, and 71 had been allocated to one of the DV clusters. Sixty were found in clusters that were fully or partially consistent with the reported RNA pattern (Supplementary file 8). Of the 11 proteins among these 71 that show consistent mesodermal upregulation in both ventralizing mutants (DV cluster 5), all are reported as ventrally expressed in BDGP.

There is also a database representing an atlas of differential gene expression at single cell resolution for precisely the time window of early gastrulation (Karaiskos et al., 2017) against which we compared the proteomes to regional RNA expression (Figure 4A). Filtering out ubiquitously expressed genes left 8924 differentially expressed genes of which 3086 coded for 3120 proteins in our clustered proteome dataset (Figure 4B).

Figure 4 with 1 supplement see all

Download asset Open asset

We first sorted these 3086 genes according to their expression patterns into the categories used above (D, V, L, DV, DL, LV) by virtue of similarity in their expression to six reference genes (Figure 4A, Figure 4—figure supplement 1A-C). In a second step, we excluded those that showed only spurious differences in expression along the DV axis, ending up with 155 genes with clear DV differences forming six DV RNA reference sets (Figure 4C, Figure 4—figure supplement 1B-D, Supplementary file 9).

We then compared the proteins in our 14 clusters against these six RNA reference sets. We asked for each protein which RNA reference set contained its corresponding gene. Theoretically, if both classifications, that is the RNA reference set and the proteomes, were perfectly correct, then genes from a protein cluster should be included only in the corresponding RNA reference set. We found that the majority of RNAs had proteins in partially or fully matching clusters of the proteomes (Figure 4D, Figure 4—figure supplement 1E). For example, nine of the thirteen proteins in cluster 5 (ventral-consistent) found their gene in the ‘twist’ similarity reference group (ventral; a perfect match: white in pie charts: Figure 4D, Figure 4—figure supplement 1E). The next best matches (e.g. ventral plus lateral, instead of only ventral; a partial match, gray) were often also highly represented: three of the four remaining cluster 5 proteins found their gene in the 'neur' similarity reference group (lateral +ventral).

Thus, the majority of proteins had perfect or partial matches with the RNA expression, showing that two independent measurements of regional expression patterns arrive at the same allocation. This confirms in an unbiased manner that the hierarchical clustering successfully sorted the proteomes in the correct manner, further supporting the initial assumption that the mutant populations were representative of specific regions in the embryo.

The difference between the results for the Toll^10B and spn27A embryos was an unexpected and potentially biologically interesting discovery. We investigated whether the matching of the protein distributions to their RNA expression patterns could give us further biological insights.

We find that for those clusters in which Toll^10B and spn27A agree, a larger proportion of proteins is allocated to the correct RNA reference set than in the clusters in which Toll^10B and spn27A differ (Figure 4D and E, Figure 4—figure supplement 1E). The ventral cluster 5, in which Toll^10B and spn27A agreed, included Snail,Twist and other genes expressed in the mesoderm (Figure 1D, Figure 2—figure supplement 1B), such as Mdr49 (Stathopoulos et al., 2002), CG4500 (Casal and Leptin, 1996), and Traf4 (Mathew et al., 2011).

For the proteins from the ‘ventral inconsistent’ clusters we found that the Toll^10B mutant differs from the spn27A mutant in a consistent manner. Proteins classified on the basis of being upregulated in Toll^10B (clusters 3, 7, 10, and 13) are often mismatched to genes with an ectodermal expression (dorsal and/or lateral RNA), whereas this does not occur for those classified based on their upregulation in spn27A (clusters 4, 8, 11, and 14, Figure 4—figure supplement 1E). This means that although Toll^10B mutants are strongly ventralized in terms of morphology and upregulation of ventral genes, the ectopic Toll signaling in the mutant fails to suppress all dorsal markers, which is consistent with our observation that spn27A mutants show a stronger reduction in dorsal-specific proteins. This confirms previous suggestions that spn27A mutants retain no or almost no DV polarity, whereas Toll^10B embryos retain residual polarity (Roth et al., 1989; Ligoxygakis et al., 2003). Determining the developmental source of these differences goes beyond the scope of this study, but will warrant further investigation.

We wondered whether there were consistent differences between those proteins that matched their RNA and those that did not. For example, a protein with large fold-changes may be more likely to match the correct RNA distribution. Because the clustering assigned proteins only on direction and not on extent of change, clusters also contain proteins with very small differences between the DV populations, even in cases where the RNA is known to show a clear difference (e.g. Traf4; Figure 2—figure supplement 1C).

To distinguish between strong and weak differential expression, we ranked proteins by comparing them to the most extreme protein in each cluster, that is the one that showed the greatest fold changes in the mutants over wildtype. We calculated the Euclidean distance (ED) between each protein and the most extreme (see Methods, Supplementary file 11). Thus, proteins with the lowest ED scores are those that are closest to the most extreme protein. We then analyzed if this score correlated with the degree to which a protein matched its RNA expression. We found that proteins from the ‘matching’ groups had ED-scores that were skewed towards lower values (Figure 4F) indicating that proteins with more extreme expression differences (low ED scores) are more likely to match the correct RNA expression pattern.

In summary, these approaches stratify our results in a useful manner: first, the DV clusters in which the two ventralized mutants behave consistently represent better the RNA expression patterns; second, proteins with strong fold-changes are more likely to represent the distribution of the corresponding RNA.

Changes in the abundance of phosphosites may occur for two reasons: either the protein itself varies in abundance, or the protein level is constant, but the protein is differentially phosphorylated. Combinations of these cases are possible, and protein abundance may be affected by phosphorylation itself. Since we know the changes in protein abundance, we can distinguish these cases by comparing the full proteome against the phospho-proteome (with the caveat that, for technical reasons, our measurements were done on parallel experiments rather than on the identical samples). 1765 of the phospho-proteins (96%) we identified were ones that we also found in the proteome, whereas 82 had not been detected in the proteome (Figure 5A). We found that most of the changes in phosphorylation were in proteins for which the level of the host protein was unchanged (black and white boxes in Figure 5B; 67 to 82% of the protein-phosphosite pairs). Among those for which the host protein showed differential abundance, 7–13% of their phosphosites changed in the same direction (both protein abundance and phosphorylation up, or both down), and 10–19% changed in the opposite direction.

Figure 5

Download asset Open asset

We tested if the phospho-proteomes fitted the 'linear model' (i.e. whether the sum of the weighted mutant values corresponded to the measured values in the wildtype) and found that the majority of the phosphosites did (Figure 5C, Supplementary file 6). Among the strongly deviating phosphoproteins, we find a number of kinases with known morphogenetic functions, such as Par-1, SRC42A and nucleoside-diphosphate kinase (awd) (Figure 3—figure supplement 1E, Supplementary file 13).

We clustered the phosphosites using the same procedure as for the proteome. After excluding sites that were unchanged or up- or down-regulated in the same direction in all mutants, clustering the remaining 3433 phosphosites again yielded 14 DV clusters (Figure 5D and H, Supplementary file 7). The two ventralising mutants now clustered together, and the dorsalized mutant showed the most distinct behavior (Figure 5G, compare to Figure 3E).

One aim of this study was to find cellular components that are differentially modified along the DV axis and that are candidates for regulating cell shape. Most likely, these cellular components are regulated by protein complexes or interacting protein networks, as already known for the regulation of actomyosin by the Rho pathway and some components of adherens junctions. Rho is activated and necessary for cell shape changes in the mesoderm, but we do not know the full set of the components of the pathway that are modulated in the mesoderm or elsewhere along the DV axis. We therefore looked at this pathway. Of 24 proteins associated with Rho signaling, we detected 21 in the wild type and at least one of the mutants (Figure 5E). Most, including the myosin light chain, occurred at similar levels in all genotypes, except Cofilin/Twinstar, Moesin, and Profilin/chickadee, which were more abundant in the ectodermal cell populations (clusters D and DL).

Fourteen of the 21 proteins were phosphorylated. These included the known phosphosites in myosin light chain (MLC) and Cofilin/Twinstar (Figure 5F, no statistical differences across genotypes for Sqh and Cofilin phosphosites, see Supplementary file 2), and the phosphorylation of the Cofilin/Twinstar kinase LIMK1 and phosphatase Slingshot (ssh), which were modulated in the D and DL clusters, as were RhoGEF2 and the MLC phosphatase Mbs (Figure 5E). In summary, we detected most of the elements of a well established pathway required for gastrulation and also identified new candidate regulation nodes within the Rho pathway.

To systematically find such networks, we used a diffusion-based algorithm (Giudice et al., 2024) on each of the DV clusters. The starting weight of each protein was based on either on its euclidean distance score ('ED', Supplementary file 11) or on the deviation from the linear model ('Dev', Supplementary file 6 and 10). Since these scores existed separately for the two ventralizing mutants, we also had to conduct the analyses twice in each case, that is once for each dataset. We focused our analyses only on the six DV clusters in which the ventralized mutants agree (D(1), L(2), V(5), DL(6), DV(9), and LV(12)). Overall, this resulted in 24 protein networks (ED score for each of the ventralizing mutants and deviation score for each mutant, each applied to the 6 clusters) for the proteome and 24 for the phospho-proteome. An ego network analysis (see Methods) yielded a set of 83 ontology terms in the proteome and 87 in the phospho-proteome that were significantly enriched in one or more networks. We concentrated our further analyses only on those ontology terms that were enriched in at least two of the four networks for each DV cluster and used a heatmap to represent them (Figure 6A and B, Figure 6—figure supplement 1A, B). The heatmaps illustrate that both experiments were highly enriched for cellular components associated with DNA and RNA metabolism or the regulation of gene expression. This is not unexpected for this developmental period of dynamic changes in gene expression. In agreement with this, the majority of the enriched proteins and phosphoproteins were characterized as nuclear ontology classes. Because of our interest in morphogenesis we focused on the cellular components that belong to cytoskeletal, cell adhesion, and vesicle trafficking categories. In the phosphoproteomes the networks enriched for cytoskeletal components were much more prevalent in the phosphoproteomes (14 of 62) than in the proteomes (3 of 63), with microtubules strongly represented (12 of 14 cellular components), in particular the alpha and beta tubulins and microtubule associated proteins (Supplementary file 12). Cytoskeletal proteins are often localized in the cell cortex, and we indeed find this association reflected in the results of the network analysis. The cell cortex is among the enriched components, and among the proteins in this category, we find cytoskeletal elements. For example, networks that include the actin-microtubule crosslinker Shot and the actin polymerase Profilin are enriched in the dorsal cluster; networks that include the apical polarity determinant Stardust or the Hippo pathway component Warts in the dorso-lateral cluster. A phosphoprotein network associated with adherens junctions and zonula adherens, one of which contains the junction-actin connectors Canoe and Girdin (Houssin et al., 2015; Sawyer et al., 2009) was enriched in the D cluster (Supplementary file 12).

Figure 6 with 1 supplement see all

Download asset Open asset

In summary, we can highlight two outcomes of the network propagation analysis. First, most networks, whether derived from the proteomes or the phosphoproteomes, are enriched for cellular components associated with regulation of gene expression (transcription, epigenetic regulation, translation, protein turnover). This is a useful validation of the approach, in that it reflects the main biological process that occurs at this stage of development: giving cells in the body different developmental fates, which is achieved through setting up different gene expression programmes. Secondly, the cytoskeleton emerges as a major target of regulation in the phosphoproteome, with the most prominent component being the microtubules. This is an interesting target for further exploration in the context of gastrulation and fits well with recent results that microtubules play a role in epithelial morphogenesis (Takeda et al., 2018; Ko et al., 2019; Booth et al., 2014; Pope and Harris, 2008).

We tested the biological relevance of the predicted phospho-regulation of microtubule networks. Before gastrulation, all cells have two subpopulations of MTs, which differ in their post-translational modifications: a disordered apical network of non-centrosomal MTs with short, non-aligned filaments, and an ‘inverted basket’ of basal-lateral MTs originating from the centrosomes and enclosing the nucleus (Takeda et al., 2018; Viswanathan et al., 2019; Figure 7A). The apical population contains only dynamic MTs, marked by tyrosinated α-tubulin, whereas the inverted basket also contains stable MTs, marked acetylated α-tubulin (Takeda et al., 2018; Wloga et al., 2017; Figure 7A, Figure 7—figure supplement 1A). During gastrulation MT acetylation patterns change. In the ectoderm, MTs become increasingly acetylated but retain their original organization whereas in central mesodermal cells, the basal-lateral MTs become less acetylated (Figure 7B, Figure 7—figure supplement 1A). Some MTs in non-constricting mesodermal cells align below the apical surfaces of these cells as they extend towards the ventral midline (Figure 7—figure supplement 1B, arrow). These MTs are non-acetylated, but partially tyrosinated (Figure 7—figure supplement 1A, blue arrowhead).

Figure 7 with 1 supplement see all

Download asset Open asset

We depolymerised microtubules with Colcemid and observed the ensuing cellular dynamics. Less than 1 min after the injection, most apical filamentous structures, astral MTs emanating from the centrosome, and the centrosomes themselves disappeared while the stable MTs associated with the nuclear envelope were partially retained (Figure 7—figure supplement 1C, D).

Colcemid treatment affected nuclear positioning and cell morphogenesis. Nuclei normally move basally for 1~2 μm in the last ~20 min of cellularization, and this failed in Colcemid-treated embryos, where the nuclei moved slightly further towards the apical cell surface (Figure 7C–K, Videos 1–3).

Video 1

Download asset

Video 2

Download asset

Video 3

Download asset

In normal embryos, nuclei in the constricting ventral furrow cells move a long way from the apical cell surface. In Colcemid-injected embryos, nuclear positioning was more random (Figure 7C and t0; 7L,7M). Ultimately, the ventral furrow failed to form (Figure 7C, t₀ +10’, Video 1).

Nuclei were also positioned apically in the neuroectoderm in Colcemid-injected embryos. The formation of the cephalic furrow was delayed by 5 min, but its progress was not affected by Colcemid-treatment (Figure 7F–H, Video 2).

Cells on the dorsal ectoderm form an apical dome with a characteristic, curved cell apex (Figure 7I, t₀-5’, insets) which is abolished in Colcemid-injected embryos, supporting the model that MT-dependent force is required for apical dome formation (Figure 7I, t₀-5’, insets). The dorsal epithelium forms folds which depend on the remodeling of apical MTs, but not on myosin contractility (Rauzi et al., 2015; Takeda et al., 2018; Wang et al., 2012) and involves the descent of the apical dome in initiating cells (Figure 7I, t₀ +5’). Dome descent does not occur in Colcemid-injected embryos, and dorsal fold formation eventually fails (Figure 7I, t₀ +5’, Video 3), supporting the current model that microtubule forces also engage on cell shortening during Dorsal Fold Formation (Takeda et al., 2018; Wang et al., 2012).

MT depolymerization also affected the apical plasma membrane dynamics. Blebs in the apical membrane of constricting mesodermal cells (Costa et al., 1994) were strikingly enlarged after Colcemid injection (Figure 7N). Lateral and dorsal cells lacked these constriction-dependent blebs. Nevertheless, after Colcemid injection, they accumulate excessive, tortuous subapical membrane (Figure 7O; Figard et al., 2013; Fabrowski et al., 2013). We also observed a distinct class of micron-scale membrane blebs in all dorsal cells, not limited to the dorsal fold initiating cells and unrelated to myosin-dependent apical constriction (Figure 7O). These blebs form during mid to late cellularization, exclude MTs and are stable for minutes.

In sum, and consistent with a role for microtubules predicted by diffused network analyses, MTs are required for correct nuclear positioning and cell shape homeostasis, and have distinct functional requirements in all three types of epithelial folds during Drosophila gastrulation. Distinct phenotypes of the apical membrane following Colcemid injection suggest differential functionality in the maintenance of membrane-cortex attachment or the dynamics of apical membrane retrieval for MT networks residing on different sides of the embryo.

Both Toll^10B and spn27A ventralising mutations produced embryos that recapitulated known biological qualities of the mesoderm along the entire DV axis, such as the expression of ventral fate determinants, or the apical localisation of the adherens junctions. However, these similarities were not fully mirrored in their proteomes. Curiously, the spn27A proteome seemed to be more similar to the dorsalized than to the Toll^10B proteome which would indicate that Toll^10B embryos are ‘more’ ventral than spn27A embryos. However, most of the mesodermal marker genes (snail, twist, mdr49, wntD, neur) make an exception are more abundant in spn27A embryos. Similarly ectodermal fate markers are more strongly downregulated in spn27A than in Toll^10B embryos. Specifically, Toll^10B mutants fail to repress the expression of ectodermal genes such as egr, zen and crb.

How can ventralizing mutations that act on the same pathway yield different proteomes? Spn27A is a serine protease inhibitor of the pathway that creates the active form of Spätzle (Spz), the ligand for Toll. Both mutations lead to constitutive activity of Toll, Toll^10B through a mutation in the receptor itself (Schneider et al., 1991; Erdélyi and Szabad, 1989), spn27A through enabling a homogeneously high level of Spz along the DV axis (rather than a peak on the ventral side; Ligoxygakis et al., 2003). Because Spz is highly abundant it should not be a limiting factor for the activation of Toll (Morisalo and Anderson, 1995; Schneider et al., 1994; Figure 2—figure supplement 1A) and loss of Spn27A should enable the full activation of Toll along the embryonic DV axis. Our results indicate that constitutively active Toll does not lead to the same level of signaling as the binding of the ligand to the receptor, and that these different levels lead to unexpected differences in the downstream targets of the signaling pathway.

The ‘linear model’ we formulated is based on the assumption that each mutant embryo faithfully represents one defined area of cells along the DV axis of the embryo, and that the full set of cell types in the embryo can therefore be reconstituted as the sum of the mutant cell types – weighted according to the area they occupy in the embryo. This should also be recapitulated for any individual protein expressed in the embryo. We found this to be true not only for the trivial cases of those proteins that occur at equal level in all genotypes, but also for most of the differentially modulated proteins. However, some proteins and phosphosites did not fit the model but showed strong deviations. One explanation could be that in the normal embryo the embryonic regions communicate with each other, and this communication is necessary for the expression or modification of certain proteins. These interactions cannot occur when the fates occur in isolation from each other in the mutants, and therefore some proteins would not be regulated properly and would not fit the model. Thus, wherever an interaction between the cell populations in the embryo is necessary for generating the correct expression or phosphorylation level, the linear model we proposed no longer applies; this means that strong deviations may indicate non-autonomous regulation. We do know some genes whose expression along the DV axis is determined by input from neighboring regions, such as Sog, Ind and single-minded. We indeed find that one of those proteins, Ind, is an outlier (deviation = 2.6) with higher than predicted expression in the dorsalized and ventralized mutants, consistent with repressive input from these regions in the wildtype.

Another case of proteins not following the linear model are those that are found in either decreased or increased abundance in all genotypes, a behavior we observed for a small percentage of proteins and phosphosites, perhaps as part of a general stress response related to the mutant situation. This is illustrated by the most extreme example, TM9SF4, which encodes an immune-related transporter that is present in all mutants at nearly 100-fold higher levels than in the wildtype. However, we did not find that this was a general rule either in the proteomes or in the phosphoproteomes: stress-related categories such as those from the 'chaperone' or 'immune response' ontology classes represented only a small percentage of proteins and phosphoproteins with the highest deviations (Figure 3—figure supplement 1D, E).

w¹¹¹⁸ (wildtype/control genotype in our studies, Bloomington stock 3605), gd⁹/FM6a (Konrad et al., 1988) (provided by S. Roth), Tl^rm9 and Tl^rm10 (Anderson et al., 1985) (provided by A. Stathopoulos), Toll^10B (Schneider et al., 1991; Erdélyi and Szabad, 1989), Df(3 R)ro80b/TM3 (Bloomington stock 2198), spn27A^ex/CyO (Ligoxygakis et al., 2003) (provided by S. Roth, Bloomington stock 6374), Df(2 L)BSC7/CyO (Ligoxygakis et al., 2003) (provided by S. Roth). To visualize non-muscle myosin in vivo, a sqh-sqh::mCherry transgene (Bloomington stock 59024) was used to construct the following stocks gd⁹;sqh-sqh::mCherry / CyO, sqh-sqh::mcherry /CyO;Tl^rm9 and sqh-sqh::mcherry / CyO;Df(3 R)ro80b/TM3.

Dorsalized embryos were derived from gd⁹ homozygous female mothers, lateralized embryos were derived from trans-heterozygous Tl^rm9/Tl^rm10 mothers, ventralized embryos we derived from Toll^10B/Df(3 R)ro80b and spn27A^ex/Df(2 L)BSC7 mothers (see Supplementary file 1). Female mutant mothers were crossed with w¹¹¹⁸ males, and the F1 from each of these crosses were collected and processed for mass-spectrometry analyses.

To visualize Myosin Light Chain, we generated the following mothers: Dorsalized: gd⁹;sqh-sqh::mCherry/+, Lateralized: sqh-sqh::mcherry/+;Tl^rm9/Tl^rm10, Ventralized: sqh-sqh::mcherry/+;Toll^10B/Df(3 R)ro80b. Female mutant mothers were then crossed with w¹¹¹⁸ males, and from F1 of these crosses, embryos in stage 5 a,b (Campos-Ortega and Hartenstein, 1997) were hand-selected under a dissecting microscope and mounted for live imaging (see below).

Embryos collected half an hour after egg-laying were allowed to develop for 2 hr 30' at 25 °C in a light and humidity-controlled incubator and then dechorionated in 50% bleach for 1' 30", washed with H₂0 and visually inspected under a dissecting microscope (Zeiss binocular) for 15'–20' at RT. To ensure younger embryos from each synchronized collection were in the target developmental stage (gastrulation stage, Stages 6 a,b Campos-Ortega and Hartenstein, 1997), we individually hand-selected the embryos on wet agar, which made the embryos semi-transparent, allowing the assessment of a range of morphological features, of which at least some are visible in each of the mutants:

• Yolk distance to embryonic surface: distinguishes between early (stage 5 a Campos-Ortega and Hartenstein, 1997) and late cellularization (stage 5b Campos-Ortega and Hartenstein, 1997).

• Yolk distribution within the embryo: identification of large embryonic movements of the germ band (eg.: Initiation of germ band extension, marking the initiation of stage 7 Campos-Ortega and Hartenstein, 1997). In DV patterning mutants, this is seen as twisting of the embryo.

• Change in the outline of the dorsal-posterior region: polar cell movement from the posterior most region of the embryo (stage 5 a/b Campos-Ortega and Hartenstein, 1997) to stage 6 a/b.

• Formation of the cephalic and dorsal folds: identification of stage 6 (Campos-Ortega and Hartenstein, 1997) (initiation of cephalic fold) and stage 7 (Campos-Ortega and Hartenstein, 1997) (dorsal folds).

The combined use of these morphological criteria, together with the synchronized egg collections allowed the accurate staging of wild type and mutant embryos. Any embryos that had developed beyond the initial stage of gastrulation (as judged by the abovementioned morphological criteria, Stage 7 Campos-Ortega and Hartenstein, 1997) were discarded and the remaining embryos were placed in 0.5 ml Eppendorf tubes and flash-frozen in liquid nitrogen. An 0.5 ml Eppendorf tube filled with embryos yields approximately 1 mg of protein.

Three transgenic lines were generated to visualize cell membranes and microtubules. For cell membranes, a single copy of EGFP or three copies of mScarlet interleaved with linkers (DELYKGGGGSGG) were trailed by a C-terminal CaaX sequence from human KRas4B (KKKKKKSKTKCVIM) for membrane targeting to yield EGFP-CaaX or 3xmScarlet-CaaX. For microtubules, EMTB-3xGFP from addgene #26741 was cloned into pBabr, a ΨC31 site-directed transformation vector, between the maternal tubulin promoter and the spaghetti-squash 3′ UTR (Takeda et al., 2018). These constructs were then integrated into the fly genome at attP2 or attP40 by Rainbow Transgenics Flies, USA, or WellGenetics, Taiwan.

SILAC metabolic labeling was performed using yeast transformed to produce heavy lysine (SILAC yeast with LysC13/6, Silantes: https://www.silantes.com/). To standardize each of the phosphoproteomic runs for each condition we labeled the proteome of w¹¹¹⁸ control flies (Figure 1—figure supplement 1D). Therefore, our analyses included a SILAC-labeled and an unlabelled w¹¹¹⁸ extract. To maximize the incorporation of LysC13/6 into the Drosophila proteome, unlabelled w¹¹¹⁸ adult flies were raised in bottles prepared with SILAC yeast as the only source for amino acids (2% SILAC yeast). The fly media were prepared in agreement with the recommendations of Silantes (https://www.silantes.com/). All emerging larvae were fed on SILAC fly medium from L1 until adult stage, in a temperature (25 °C), light and humidity-controlled incubator (Sanyo). The emerging labeled w¹¹¹⁸ adults were then transferred to cages for embryo collections, and fed with wet SILAC yeast until disposal of flies (after 2–3 weeks). SILAC w¹¹¹⁸ embryos were collected as described above.

Using this protocol, we labeled ~75% of the proteome of SILAC w¹¹¹⁸ embryos. SILAC labeling did not affect the phosphoproteome of wild type embryos, and had only a minor effect on phosphosite intensity distribution, indicating standardization with Lys 13/6 was a valid approach (Figure 1—figure supplement 1E).

For phosphopeptide enrichment, the High-Select TiO2 Phosphopeptide Enrichment Kit (#A32993) was utilized following the manufacturer’s instructions. In brief, desalted peptides were dried to complete dryness and resuspended in Binding Buffer (included in kit). The peptide solution was centrifuged (10 min, 12.500 x g, 22 °C) and the supernatant was transferred to TiO2 tips. Phosphopeptides were enriched and eluted using the provided elution buffer. The eluate was immediately dried in a SpeedVac concentrator (Eppendorf) and stored at –20 °C. Prior to LC-MS/MS measurement, peptides were solubilized in 10 µL of 2% formic acid and 2% acetonitrile. Three µL were injected per LC-MS/MS run. The phosphopeptide enrichment was performed in technical duplicates.

The LC-MS/MS instrumentation consisted of a nano-LC 1000 coupled to a QExactive Plus or of a nano-LC 1200 (Thermo Fisher) coupled to a QExactive HF-x instrument via electrospray ionization. The buffer system consisted of 0.1% formic acid in water (buffer A) and 0.1% formic acid in 80% acetonitrile. The column (75 µm inner diameter, 360 µm outer diameter) was packed with PoroShell C18 2.7 µm diameter beads. The column temperature was controlled to 50 °C using a custom-built oven. Throughout all measurements, MS1 spectra were acquired at a resolution of 60,000 at 200 m/z and a maximum injection time of 20ms was allowed. For whole proteome measurements, the mass spectrometer operated in a data-dependent acquisition mode using the Top10 (QExactive Plus) or Top22 (QExactive HF-x) most intense peaks. The MS/MS resolution was set to 17,500 (QE-Plus) or 15,000 (QE-HFx) and the maximum injection time was set to 60ms or 22ms, respectively. Samples of replicate one and three were measured on the QE-Plus system and replicate two was measured on the QE-HF-x system.

For phosphoproteome analysis, the MS2 resolutions were set to 30,000 (QEx-Plus) or 45,000 (QEx-HFx). Samples of all three replicates were measured on the QEx-HFx system. We added trial samples measured on the QEx-Plus system to increase the phosphosite coverage.

Raw files were processed using MaxQuant (v. 1.5.3.8; Cox and Mann, 2008) and the implemented Andromeda search engine (Cox et al., 2011). The Uniprot reference proteome for Drosophila melanogaster (downloaded: 07.2016, 44761 entries) was utilized. Phosphoproteome and proteome data were analyzed separately and the match-between-runs algorithm was enabled. For proteome analysis, the label-free quantification method (MaxLFQ) was enabled using the default settings. Default settings for the mass tolerances for FTMS analyser were used. The FDR was controlled using the implemented (Andromeda) reverse-decoy algorithm at the protein, peptide-spectrum-match and PTM site levels to 0.01.

For SILAC-based phosphopeptide quantification, a minimum ratio count of 2 was required; the minimum score for the modified peptide was set to 20. Technical duplicates were aggregated by using the log2 normalized SILAC ratio median.

The proteinGroups (proteome) and PhosphoSite(STY) tables were subjected to downstream analysis. Gene Ontology annotations were derived from the Uniprot database and annotated. The LFQ intensities were log2 transformed. Pairwise comparisons were performed using a two-sided unpaired t-test. One-way Analysis of Variance (ANOVA) was performed on genotypes and a FDR was calculated using a permutation based approach (s0=0.1, #permutations = 500) in the Perseus software (Tyanova et al., 2016).

Protein log2 fold change ratios were matched to the phosphosite table using the Uniprot identifiers. If the phosphorylation site was part of multiple protein groups, the average log2 fold change was utilized. The analysis of the correlation between the fold changes of phosphosites and their host proteins was performed as follows: for each proteome-matched phosphosite, a protein-phosphosite pair was assembled, yielding 6297 pairs of phosphosites and their host proteins taking into account all genotypes. We tested for each differentially phosphorylated site whether its respective protein was up- or down-regulated and made this comparison for each mutant genotype versus the wildtype. To consider a protein and a phosphosite regulated, we applied the same thresholds as used for clustering: ± 1.4 (0.5 in log2) fold change for proteins, and ± 1.3 (0.35 in log2) fold change for phosphosites (see below Hierarchical clustering analyses, Threshold determination). The protein-phosphosite pairs were placed in a scatter plot with 4 quadrants that were connected to 3 possible behaviors: correlation (fold change of host protein and phosphosite are both positive -green- or negative -blue-), anti-correlation (fold change of host protein and phosphosite have different signs -red/magenta-) or no-correlation, (fold change of host protein is within threshold range but phosphosite trespasses it -black/white-). Finally, we counted the number of protein-phosphosite pairs that displayed each of these described behaviors, and this quantification was used to assemble the bar plots.

Rehydrated embryos were blocked for 2 hr in 2% BSA (B9000, NEB) in PBS with 0.3% Triton X-100 (T9284 Sigma). Primary antibody incubations were done overnight at 4 °C. Primary antibodies used were: mouse anti α-tubulin 1:1000 (T6199, clone 6-11B-1, Sigma), mouse anti acetylated-α-tubulin FITC conjugated 1:250 (sc23950, Santa Cruz Biotechnology), rat anti tyrosinated-α-tubulin 1:250 (MAB1864-I, clone YL1/2, EMD Millipore/Merck), mouse anti-Armadillo/β-Catenin 1:50 (N27A1, Developmental Studies Hybridoma Bank), rabbit anti Snail (Tyanova et al., 2016) 1:500. Incubations with secondary antibodies were performed for 2 hr at RT. Alexa Fluor 488- and 594-coupled secondary antibodies were used at 1:600 (488 and 594 Abcam).

Immunostained embryos embedded in Fluoromount G (SouthernBiotech 0100–01) were visually inspected under a dissecting microscope (Zeiss binocular) to select the desired developmental stages. The embryos were sectioned manually with a 27 G injection needle at approximately 50% embryo length and slices were mounted for microscopy.

Images in Figure 1C were acquired with a Zeiss LSM880 Airyscan microscope, using a Plan-Apochromat 63 x oil (NA 1.4 DIC M27) objective at 22 °C, with a z-slice size = 0.3 μm. Acquired volumes were max-projected (along z axis) for a range of 1.5 μm (5 slices). Images in Figure 7A and Figure 7—figure supplement 1B were acquired using a super resolution Deltavision OMX 3D-SIM (3D-SIM) V3 BLAZE from Applied Precision (a GE Healthcare company). Deltavision OMX 3D-SIM System V3 BLAZE is equipped with 3 sCMOS cameras, 405, 488, and 592.5 nm diode laser illumination, an Olympus Plan Apo 60X1.42 numerical aperture (NA) oil objective, and standard excitation and emission filter sets. Imaging of each channel was done sequentially using three angles and five phase shifts of the illumination pattern. The refractive index of the immersion oil (Cargille) was 1.516. Acquired volumes were max-projected (along the z axis) for a range of 3 μm (10 slices). Images in Figure 7B and Figure 7—figure supplement 1A were acquired with a Leica SP8 microscope equipped with white laser. Gated detection on HyD detectors was used for each shown channel using a Plan-Apochromat 63 x oil (NA 1.4) objective at 22 °C, with a z-slice size = 0.3 μm. Acquired volumes were max-projected (along z axis) for a range of 1.5 μm (5 slices).

Dechorionated embryos expressing maternally provided sqh::mCherry were mounted in 35 mm glass-bottom petri dishes in two different ways with either the embryonic dorsal, lateral, or ventral surface facing the glass bottom or vertically glued (heptane glue) on their posterior end to the glass-bottom and embedded in 0.8% low melting point agarose in PBS that was previously cooled to 30 °C. For the superficial imaging of the sub-apical domain of embryos along their dorsal, lateral and ventral sides, we acquired volumes of 20 μm (z-slice size of 0.8 μm) using a PerkinElmer Ultraview ERS (microscope stand Zeiss Axiovert 200) with a Yokogawa CSU X1 spinning disk with a Plan-Apochromate 63 x (NA 1.4, oil) objective at 22 °C. The vertical mounting enabled the imaging of myosin light chain along the dorso-ventral cross-section of living embryos (Krueger et al., 2020) in the x-y plane. We acquired 1–3 slices (z-slice size = 1,1 μm) at 130–150 μm from the posterior end of the embryo using a Zeiss LSM780 NLO 2-photon microscope with a Plan-Apochromat 63 x objective (NA 1.4, oil, DIC M27) at 22 °C.

Correlation matrices were calculated using Matlab R2019b on the proteins and phosphosites that were detected in all genotypes and all the replicates. For this, we used the 'corr’ function to calculate the Pearson correlation between log2 intensities of proteins or phosphosites in the replicates of all genotypes. The resulting correlation matrix was plotted using the ‘clustergram’ function by applying the following settings: linkage: average; RowPDist & ColumnPDist: correlation.

We defined a ‘linear model’ that is based on the assumption that the three types of mutant embryos (dorsalized, lateralized, and ventralized) each represent one region along the DV axis of the embryo, and that protein abundance in the three regions should add up to the total protein abundance in the entire embryo, and therefore, the abundance in the three mutants should add up to the same value, if each is weighted by the region occupied in the wildtype embryo.

This linear model can be expressed for each protein ProtX as a sum, ^twt_ProtX, where D, L, and V are the abundance of a protein 'ProtX' in the proteomes of the three mutant genotypes (means of the log2 intensity values of the replicates, transformed to its linear value), and a, b, and c represent the proportion of each region along the dorsoventral axis:.

This model requires values a, b, and c for the three regions along the DV axis. For the mesoderm, this has been reported as 0.2 from measurements on cross-sections (Rahimi et al., 2016), but we wanted to determine the theoretical optimum for each of the values without any prior assumption about their real proportions in the embryo. The theoretical optimum would be one for which the proportions for the three regions when used in the sum yield a theoretical ('t') value ^twt_ProtX that is the closest to the experimentally measured ('m') value ^mwt_ProtX for that protein in the wildtype embryo.

To systematically explore the proportions for each region, we tested all possible combinations for a, b, and c at 0.05 steps in the range from 0 to 1 (i.e.: 0, 0.05, (...), 0.95,1). For each calculated ^twt_ProtX, we calculated the deviation from the measured abundance by subtracting the mean linear intensity measured for the same protein in a wild type embryo (^mwt_ProtX), and transforming this difference between the theoretical and measured wild type to log2 scale as follows:

A log2 value of 0 for the deviation therefore indicates a perfect match between the theoretical calculation and the actual measurement. For each possible combination of a, b, and c we obtained distributions of Deviation_ProtX values, for which we calculated the Interquartile Range (IQR) using the IQR function from Matlab 2019b. Based on the assumption that the best matching proportions should lead to the narrowest dispersion of the distribution of Deviation_ProtX, we sorted (smaller to largest) the combinations of proportion constants based on their calculated IQRs. This parameter screen yielded good fits for a range of combinations. Previous work indicated the mesoderm represents 20% of the circumference of the embryo Rahimi et al., 2016; however, for the two best, the area of the ventral region was slightly larger than the observed 20% of the circumference of the embryo in vivo (Figure 3—figure supplement 1B, C, Supplementary file 5). The third best was one for which the ventral domain corresponded to the experimentally measured value of 0.2 (20%, Figure 3—figure supplement 1B, C, Supplementary file 5), and for the lateral and dorsal domains the value was 0.4, which matched estimations based on the expression domains of lateral and dorsal genes (Stathopoulos et al., 2002; Jaźwińska et al., 1999). We therefore chose this set:

Because we used two different ventralising genotypes (Toll^10B and spn27A^ex), the ^twt_ProtX was calculated twice for each protein, one for for each mutant genotype combination (i.e.: D-L-V_Toll10B and D-L-V_spn27Aex).

We determined the threshold to be applied to the FCs between the DV mutant genotypes and the wild type of each proteomic experiment by analyzing the variability of the FCs within each biological replicate. We first took the mean log2 intensity per protein for the wildtype '^WTMean_ProtX', and next, calculated the log2 FC between the log2 intensity of each technical replicate with measurable intensity for each biological replicate and the ^WTMean_ProtX as follows:

where BP-G is the biological replicate of a particular genotype and TR-N a technical replicate (1<N < 3) for a biological replicate BP-G. Therefore, we obtained a set of ^BP-GFC(ProtX)_TR-N values per protein, and the number of ^BP-GFC(ProtX)_TR-N values per protein equals the number of technical replicates with measurable intensity for a particular biological replicate.

Next, we calculated ^BP-Gstdev(ProtX), which is the standard deviation of the set of ^BP-GFC(ProtX_TR-N values per protein and for each biological replicate. In this way, we obtained a distribution of the standard deviation of FCs per biological replicate. For each biological replicate distribution of FCs standard deviations we calculated the IQR (Prism Graphpad V8) and extracted the 3rd quartile value (^BR-GQ3) to capture up to 75% of the variability of each distribution:

Finally, we defined the FC threshold as the mean of the Q3 value across the biological replicates of each experiment (number of: biological replicates proteome = 6; biological replicates phosphoproteome = 4) as follows:

We excluded all genes that were scored as ubiquitously expressed in THE RNA ATLAS (Karaiskos et al., 2017). Of the remaining 8924 genes, we selected those that were also listed in our clustered proteomic dataset (3346 genes coding for 3383 proteins), yielding a list of 3086 non-ubiquitous genes that were present both in the RNA atlas and the clustered proteome. These are sorted into classes by comparing the expression pattern of each to that of six reference genes with restricted dorso-ventral expression that represent the six regulation categories D, L, V, DL, DV, and VL. We used as reference genes dpp for dorsal, the average between sog and soxN for lateral, twist for ventral, crb for dorsal +lateral, net for dorsal +ventral and neur for lateral +ventral (Figure 4A). To compare similarity for each of the 3086 genes we calculated their spatial Pearson correlation to each of the reference genes. Each gene was then classified as belonging to the category of the reference gene for which the Pearson correlation was the highest. We therefore obtained 6 clusters of genes, which we termed ‘DV RNA clusters’ each of them with their corresponding maximum Pearson spatial correlation value (Figure 4—figure supplement 1). To filter out false positives, we selected those genes with the largest similarity to the reference genes, for which we expected strong differential expression along the dorso-ventral axis. We did this by determining the Pearson spatial correlation value corresponding to the 95th Percentile of each dorso-ventral RNA cluster, using the 'prctile' function (Matlab R2019b). Finally, we used the 95th percentile value as a threshold to filter the 5th percentile of genes from each DV RNA cluster. We obtained a list of 155 genes that we used to compare against the proteome clusters (Supplementary file 9), which we termed ‘DV RNA Reference Sets’.

The filtered set of 155 genes codes for 157 proteins. We grouped the 157 proteins based on their DV cluster assignment and for each DV cluster, we classified its proteins based on the RNA reference set to which their genes had been allocated. Theoretically, if both classifications, i.e. the RNA reference set and the proteomes, were perfectly correct, then genes from a protein cluster should be included only in the corresponding RNA reference set. For DV clusters 1–12, we classified the results of this comparison as 'perfect match if the RNA expression pattern and the DV cluster belong to the same regulation category, as 'partial match if the RNA expression pattern and the DV cluster coincided only in one DV domain, or as mismatch if the RNA expression pattern and the DV cluster belong to mutually exclusive regulation categories.

We developed a score based on a calculated Euclidean distance to measure the proximity of each protein in a particular DV cluster to the most extreme fold changes measured in DV mutant vs. wild type comparisons. We used the same approach for the phosphosites.

We first calculated the mean log2 fold change (FC) between the DV mutants and the wild type (which meant we could not assess proteins nor phosphosites that were not detected in the wild type). Next, we rescaled each of the FC distributions [dorsalized (D) vs. WT, lateralized (L) vs. WT, ventralized Toll^10B/def (Vtl) vs. WT and ventralized spn27A^ex/def (Vsp) vs. WT] to transform the log2 FC = 0 in the original distributions as log2FC = 0.5 in the new, rescaled distribution. We first identified the upper and lower limits of each FC distribution, and next, transformed each pair (upper and lower) of limits to their absolute values. This enabled the identification of the largest absolute limit for each FC distribution, and depending whether the upper or the lower limit was the absolute largest, we used one of the following equations to rescale:

If the upper limit of a particular FC distribution is the absolute largest:

• i. log2FC_rescaled_i = (log2FC_original_i + |max_log2FC_original|) / (2* |max_log2FC_original|)

If the lower limit of a particular FC distribution is the absolute largest:

• ii. log2FC_rescaled_i = (log2FC_original_i + |min_log2FC_original|) / (2* |min_log2FC_original|)

where 'i' is a particular protein, and max_log2FC_original and min_log2FC_original are the upper and lower limit values for a particular FC distribution. Using this rescaling approach, we obtained a new set of rescaled FCs for each distribution (D vs. WT, L vs. WT, Vtl vs. WT and Vsp vs. WT). We considered those proteins that were undetected in a mutant genotype as being in decreased abundance in that genotype, and imputed the lower limit value of the rescaled FC distribution for that genotype.

For each protein, we assigned two vectors with their corresponding rescaled FCs (shown here for one ventralized genotype, but also calculate separately for the other):

Next, we assembled reference vectors representing the most extreme behaviors for each regulation category (Figure 3D, upper panel), using the rescaled FCs:

Where max^rFC and min^rFC are the maximum and minimum rescaled FCs for the corresponding distributions (D vs. WT, L vs. WT, Vtl/Vsp vs. WT). Finally, we calculated the Euclidean distance score (ED Score) between each protein in clusters 1–12, and the reference vectors that correspond to each DV class:

Where D_ref, L_ref and V_ref are the reference vector values* for the corresponding regulation category.

We applied the random-walk-with-restart-based algorithm (Giudice et al., 2024) for the following DV clusters: D (1), L (2), V (5), DL (6), DV (9) and LV (12), once each for the ED score and once for the deviation values, and each for the protein and phosphoproteomic datasets. In the case of the deviation values, we used only those proteins or phosphosites within the interquartile range. We assigned as seed value the reciprocal of the absolute value from both the ED Score or the log2 Fold Change Deviations. Note that if multiple phosphosites are assigned to the same protein we selected the median of the ED scores or log2 Deviations. We then partition the seed set in tyrosine kinases, remaining kinases and other proteins and perform the random walk with restart (RWR) from each of the three partitions separately. We also repeated the same procedure with the same set of initial nodes against 1000 random networks. We estimated the empirical p-value for each node of the network as the percent of its random scores that exceed the real score and selected only the nodes with a p-value <0.05 in at least one of the three partitions. The resulting subnetworks are further filtered using the ego decomposition (Giudice et al., 2024). Briefly, for each seed node we extracted ego networks with a maximum distance of 2 steps from the ego. We then filtered the ego networks by selecting only the most similar functional nodes to the ego using this formula:

where ${\displaystyle Resnik(ego,j)}$ is the semantic similarity between the ego and a node j in the ego network, ${\displaystyle \mathrm{m}\mathrm{e}\mathrm{a}{\mathrm{n}}_{\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{n}\mathrm{i}\mathrm{k}\mathrm{e}\mathrm{g}\mathrm{o}}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}\mathrm{t}{\mathrm{d}}_{\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{n}\mathrm{i}\mathrm{k}\mathrm{e}\mathrm{g}\mathrm{o}}}$ are the mean and the standard deviation of the ego against all the other nodes in the initial network. Nodes with a z-score>1.28 (equivalent to a p-value <0.1) are retained. After this filtering, the resulting ego networks with less than 5 nodes are discarded. Additionally, the weight of the ego networks are changed according to (3) to reflect the functional impact of the dysregulation of the ego on the neighboring nodes.

where Γ(ego) represents the nodes at distance 1 from the ego and Γ_Γ(ego) represents the nodes at distance 2 from the ego. The ego networks obtained are normalized again to correct for hubs using the Laplacian normalization using (1). For each ego network, we then calculate the topological distance vector and the functional distance vector as in Giudice et al. The topological distance vector is calculated using the following formula (4):

where $jsd$ refers to the Jensen-Shannon distance, representing the similarity between two probability distributions. The RWR_node refers to the RWR probability vector when one of the nodes of the ego network is selected as seed, and the RWR_ego refers to the RWR probability vector when the ego is the seed node. The functional vector is defined as the logarithm of the semantic similarity between the ego and any other nodes in the network (5).

where Resnik(ego, node) represents the semantic similarity measure between the ego and the node under consideration. To assess the most similar nodes to the ego, the Kernel Density Estimation (KDE) measure (with Gaussian kernel and bandwidth estimated using the Silvermann formula) to assess the most similar nodes to the ego, is employed. KDE estimates the probability density function (PDF) of the topological and semantic similarity vectors obtained at the previous step. For each ego network, we selected only those nodes within a 0.7≤PDF ≤ 1.0 of both topological and functional similarity. All the nodes overcoming this threshold are selected for the enrichment analysis against the cellular component domain of GO.

We filtered the proteins and phosphosites with an absolute log2 deviation value larger than the 95^th percentile (prctile Matlab function) of the distribution of absolute log2 deviation values. Because deviation values were calculated separately for each ventralized genotype (Toll^10B/df or spn27A^ex/df), we obtained two independent lists of proteins and phosphosites with deviation values exceeding the 95^th percentile. From these lists, we selected the proteins and phosphosites whose log2 deviations exceeded the 95th percentile threshold with both ventralized genotypes. For phosphosites, we used the host phosphoprotein for the ontology analyses. When two or more phosphosites with large deviations were hosted by the same phosphoprotein, we counted the phosphoprotein only once. We therefore obtained 206 proteins and 154 phosphoproteins (191 phosphosites) that were used in the ontology analyses.

The ontology analyses were performed using the PANTHER platform (http://www.pantherdb.org/, release PANTHER 17.0 dated February 23rd 2022). We queried the 'Functional classification gene list', based on 'Drosophila melanogaster' organism data. We used the FBgn (Flybase Gene Number) of the proteins and phosphoproteins to produce the query in Panther, and focused on the 'Protein Class' classification. The protein class query allocated 125/206 proteins and 110/154 phosphoproteins to protein class terms. Using gene ontologies from Flybase we manually allocated 48/206 proteins and 35/154 phosphoproteins to one or more of the 24 parental protein class categories. 33/206 proteins and 9/154 phosphoproteins could not be allocated to any parental class category, remained unassigned and were therefore excluded from the reported analyses. In summary, the reported protein class ontology analyses of extremely deviating proteins and phosphoproteins is based on 173/206 proteins and 145/154 phosphoproteins.

For quantitation of nuclear position, single z-slice H2Av-mRFP1 images were converted into tiff format, blurred with a Gaussian filter (σ=2), and segmented with CellPose (v2.2) in 2D using a custom-pretrained nuclear model tailored for each side of the embryo based on manual correction on segmentation generated by a default nuclei model with nuclear diameter set as 20 pixels. The stitch mode was used with a stitch threshold of 0.4 to generate 3D nuclear segments. Nuclear segments were filtered by size (1000–5000 pixels) and height (>10 μm), while those located at the edge of the imaging area were excluded for data processing. The regionprops function implemented in the Skimage Python library was used to define the bounding box of each nuclear segment, from which the middle Z coordinate of the bounding box was designated as the nuclear position. For time alignment, t₀ (the onset of gastrulation) was defined as the time point, at which apical constriction in ventral furrow produces a 2.5 μm gap between the cell apex and vitelline membrane for the datasets acquired on the ventral side. Using this t₀ designation (from water-injected embryos imaged on the ventral side), the pre-injection cellularization depths were fitted to a linear function based on the assumption that cellularization depth is linear with time during mid-cellularization. For datasets acquired on the lateral and dorsal sides, the pre-injection timing relative to the onset of gastrulation was derived by plugging in the pre-injection cellularization depth, from which the t₀ frame of the dataset was derived. Data processing and plotting were performed with custom-made Python codes using Numpy, Pandas, Matplotlib, and Seaborn libraries.

For en face membrane visualization, 3xmScarlet-CaaX images were deconvolved using the Huygens Software (Scientific Volume Imaging) with the deconvolution algorithm Classic MLE using custom parameter sets.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank NH Brown, N Bulgakova, S Huelsmann, N Perrimon, V Riechmann, A Stathopoulos, D Stein, S Roth and A Wodarz for reagents and fly stocks, N Lawrence and the Gurdon Institute Imaging Facility for help with 3D-SIM imaging, EMBL Advanced Light Microscopy Facility (ALMF) and CECAD Cologne Imaging Facility for continuous support. Flybase was used throughout the project and is gratefully acknowledged. We also thank Siegfried Roth for critical discussions, and all members of the Leptin lab for discussions throughout the work. This work was supported by funding from the European Molecular Biology Organisation (EMBO), the University of Cologne and the Deutsche Forschungsgemeinschaft (grant DFG LE 546/12).

© 2024, Gomez et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.