w¹¹¹⁸ , Beadex^MS1096-Gal4 (referred here as Bx: 8860/Bl), Beadex^MS1096-Gal4; EaterDsRed, w¹¹¹⁸ , w¹¹¹⁸;p35 (5072/Bl), w¹¹¹⁸;UAS-Ras^V12 (4847/Bl), w¹¹¹⁸;p35;UAS-Ras^V12,;Arc1^esm18; (37530/Bl), and Bx;Arc1^esm18;,;Arc1^esm18;Ras^V12 flies were used in the experiments. Flies were cultured in a 25°C, 12-h dark/light cycle room. Female virgins were collected for 5 days and crossed with respective males on day 7. Progeny larvae were kept as described in (20). Approximately 20–30 salivary gland pairs were fixed in 4% paraformaldehyde (PFA) for 20 min. Samples for extracellular staining were washed 3 × 10 min with 1× phosphate-buffered saline (PBS).

Antibodies against P1 (Plasmatocytes, 1:50), L1 (Lamellocytes, 1:50), and C4 (Crystal cells, 1:50) were incubated for 1 h at room temperature in PBS and subsequently washed 3 × 10 min in PBS. For intracellular proteins, the samples were incubated overnight at 4°C in anti-pJNK (1:250), anti-Idgf3 (1:50), or anti-CC3 (1:400) diluted in PBST (1% TritonX-100) and subsequently washed 3 × 10 min in PBS. The samples were incubated with secondary antibody anti-mouse-546 (1:500, Thermo Fisher #A11030), anti-mouse-488 (1:250, Thermo Fisher #A11001), anti-rabbit-568 (1:500, Thermo Fisher #A21063), and anti-rabbit-488 (1:500, Thermo Fisher #A11008), or 4’,6-Diamidino-2-phenylindole dihydrochloride (DAPI) (1:500, Sigma-Aldrich D9542) for 1 h at room temperature and washed 3 × 10 min with PBS before mounting in Fluoromount-G (Thermo Fisher, 00-4958-02).

Drosophila larvae at 120 h after egg deposition (AED) were bled into 10 µL of PBS, and the sample was incubated for 5 min at room temperature. The sample was fixed in 4% PFA for 20 min and washed 3 × 5 min with PBS. The sample was stained with antibody against pH3 (phosphorylated histone 3, 1:800, Sigma-Aldrich 06-570) and H2 (Hemese, 1:5 kindly provided by Istvan Ando, Szeged) overnight in 4°C. After incubation, the sample was washed 3 × 5 min with PBS and stained with secondary antibodies anti-rabbit-568 (1:500, Thermo Fisher A-21069), anti-mouse-488 (1:250, Thermo Fisher A-11001), and DAPI (1:500, Sigma-Aldrich D9542) for 2 h in room temperature. Subsequently, the sample was washed 3 × 5 min with PBS and mounted in Fluoromount-G (Thermo Fisher, 00-4958-02). A total of 15 samples per genotype were prepared and photographed using a tile scan at a confocal microscope Zeiss LSM780 (Zeiss, Germany), and the count of cells was evaluated in ImageJ (version 1.53t). Statistics and graph plotting were performed in Prism software (GraphPad Software, 9.3.0, USA).

Whole salivary glands were photographed using a Zeiss Axioscope II (Zeiss, Germany) microscope, and images were exported as.tiff files. The intensity was measured using ImageJ (version 1.53t), and representative pictures were obtained from a confocal Zeiss LSM780 (Zeiss, Germany) microscope. Statistics was performed in Prism software (GraphPad Software, 9.3.0, USA).

A total of 10 pairs of salivary glands per genotype were dissected into individual droplets of 200 µL PBS. A fluorescent microscope Leica DMi8 (Leica, Germany) was used to locate the attached hemocytes at the surface of the salivary gland. The hemocytes were collected using a micromanipulator TransferMan NK2 (Eppendorf, Germany) with ø 7 µm collection capillary and separated by a piezo-vibrator PiezoXpert (Eppendorf, Germany). Single-cell state hemocytes were individually transferred to 2.3 µL of lysis buffer (21), and each sample library was prepared as described below.

cDNA libraries of single attached hemocytes were generated using a modified version of the Smart-seq2 protocol (21). In short, cDNA synthesis was performed using universal primers, and PCR amplification was carried out over 24 cycles. cDNA products were subsequently purified using CA beads (Sigma; catalog no. 81260) for size selection using 8.8% polyethylene glycol 6000 (PEG 6000) to exclude primer dimers and nonspecific amplicons with sizes <150 bp. Combinatorial indexing via tagmentation was carried out in 96-well plates using 200 pg (measured in a Qubit fluorometer) of amplified cDNA, for a final volume of 10 mL/well. cDNA fragmentation using Tn5 transposase was carried out for 20 min on ice using the Illumina Nextra XT DNA sample preparation kit. Ligation and amplification of adaptors were carried out over 15 cycles in a final volume of 25 mL/well. Primer indices were used in the reaction from Illumina (Nextera index primers i7 and i5, catalog no. FC-131-1001). Tagmented and barcoded amplicons were then purified using CA beads for size selection. Quality control and fragment size distribution of the cDNA libraries were performed on a Bioanalyzer with the Agilent high-sensitivity DNA chip (catalog no. 5067-4626). Concentrations of each sample of cDNA libraries were measured on a PicoGreen 96-well plate NucleoScan fluorometer using a high-sensitivity double-stranded DNA (dsDNA) (HS assay kit; catalog no. Q32851). To perform library dilutions, the average fragment sizes of all cDNA libraries were measured for a final concentration of 2 nM in each sample. Finally, cDNA libraries were pooled and sequenced using Illumina NextSeq with 75-bp paired-end reads.

A total of 10 larvae per genotype were bled into 500 µL of PBS, and the cells were sorted. The sorting was performed with a MoFlo Astrios EQ (Beckman Coulter, USA) cell sorter using a 488- and 532-nm laser for excitation, 100 µm nozzle, sheath pressure of 25 psi, and 0.1 µm sterile-filtered 1× PBS as sheath fluid. Flow sorting data was interpreted and displayed using the associated software, Summit v 6.3.1.

To test the precision of the adjustments made to center the drop in each well, a colorimetric test mimicking the sort was done based on (22). A 1.5 mg/µL solution of horseradish peroxidase (HRP) (catalog no. 31490, Thermo Fisher Scientific) with one drop of flow check beads (Beckman Coulter, USA) was sorted into each well of an Eppendorf 384-well plate (Cat no 34028, Thermo Fisher Scientific). A color change after sorting indicated that the drop hit the sort buffer and that the precision was adequate.

Single hemocytes were sorted directly into a 384-well plate containing 2.3 µL of lysis buffer (Eppendorf twin.tecTM PCR plates) using a CyClone™ robotic arm and at highly stringent single-cell sort settings (single mode, 0.5 drop envelope), and cDNA libraries were generated by the Eukaryotic Single Cell Genomics Facility at SciLifeLab, Stockholm using a slightly modified version of Smart-seq2 as previously described (21), but where we used 20 cycles for cDNA amplification. The plate and sample holder were always kept at 4°C during the sort. After the sort, the plates were immediately spun down and put on dry ice.

Single-cell libraries were sequenced at the National Genomics Infrastructure, SciLifeLab Stockholm, using the HiSeq2500 platform (Illumina) for 56 bp single-end sequencing. We sequenced a total of 384 single cells from the circulating hemocytes (negative controls, n=2, per plate) and 79 from attached hemocytes.

A reference genome D. melanogaster (Dm6 v r6.37) was indexed, raw fastq files were used in mapping to the genome using STAR v 2.7.2 (23), and gene expression was measured using featureCounts v 2.0.0 (24), using default settings. The expression matrix was filtered to obtain high-quality cells using the following criteria: cells with >5% mitochondrial transcripts (stressed/dead/dying cells), <200 genes (low-quality cells), and those expressing more than 4,000 features (genes) (potential doublets or triplets) were removed in each replicate, and the remaining cells were subjected to subsequent computational analysis. A total of 380 and 61 single cells have passed the computational filtering from the circulating and the attached hemocyte, respectively, and are used for the data integration.

The main computational analysis of read-count matrices was performed using the Seurat package (v 4.0.3) (25) in R (v 4.1.0). The complete R workflow can be assessed and reproduced in the R markdown (see code availability section). We used the default processing pipeline, https://satijalab.org/seurat/v3.2/pbmc3k_tutorial.html. First, count matrices and metadata were loaded. A mitochondrial gene count above 10% was filtered. Quality filtering was performed, and cells with a minimum of 200 genes expressed were kept for further processing. Subsequently, reads were normalized for sequencing depth using the “NormalizeData” function in the Seurat toolkit, selecting the top 2,000 variable genes. Thereafter, dimensionality reduction was performed using Principal component analysis (PCA) computing the first 50 Principal components (PCs). The first 10 PCs from the analysis were then subjected to shared-nearest-neighbor (SNN) inspired graph-based clustering via the “FindNeighbors” and “FindClusters” functions. For modularity optimization, the Louvain algorithm was used, and clustering was performed at a resolution of 0.4 for clustering granularity, resulting in five clusters. After clustering, a UMAP dimensionality reduction was performed using the first 10 dimensions of the PCA.

Differential gene expression analysis (DGEA) of genes in identified clusters was performed using the function “FindAllMarkers” from the Seurat package (v. 4.0.3). Following the default option of the method, differentially expressed genes for each cluster were identified using a non-parametric Wilcoxon rank sum test. Differentially expressed genes in a cluster were defined by setting initial thresholds above a logarithmic fold-change of 0.5 and being present in at least 25% of the cells belonging to the same cluster. Representative marker genes with an adjusted p-value below 0.05 for each cluster were further selected. p-Values were adjusted using a Bonferroni correction including all genes in the dataset. To find representative marker genes with elevated expression in comparison to the remaining clusters, only positive log fold-changes were considered. For individual analyses such as gene enrichment analysis (see “Gene set enrichment analysis (gsea)”), threshold values for differential gene expressions were modified and will be described in detail in the respective sections of the Material and methods and Results. To identify DEGs between specific clusters of interest, the “FindMarkers” function in Seurat was used, and the identities were set to the respective clusters of interest. The same thresholds as stated above were used to define DEGs.

To track tests for top functional class enrichment among the global clusters, we used “ClusterProfiler” package v 3.16 (26) tool to conclude the enriched ontology terms as previously mentioned specifying the database “Org.Dm.db” to calculate the top 5 biological pathway enrichment. The gene set enrichment analysis (GSEA) was performed on top differentially expressed genes over the identified clusters in regard to gene expression profiles “Log2FC” as input in cluster profiler v 3.18.1 and ggupset package v 0.3.1 with a p-value cutoff of 0.05, minGSSsize of 3, maxGSSize of 800, and scoreType of “pos” to estimate for biological process ontology across clusters.

The read alignment and gene count matrix generation were performed as previously described in Material and methods (Section 2.9) (27). The single-cell gene count matrix cells with fewer than 250 UMIs, more than 10,000 UMIs, reads mapping to more than 7,000 genes, or more than 10% of read counts mapping to ribosomal genes were excluded. Each single-cell transcriptome was mapped to its original time window from which it was extracted by using the RT barcode. We performed standard processing of the data split by experiments as recommended by Seurat v4 documentation including NormalizeData, FindVariableFeatures (with the method set to “vst”), ScaleData, RunPCA, RunUMAP (with dims set to 1:50 and n.components to 2), FindNeighbors (with reduction set to UMAP and dims set to 1:2), and FindClusters. This version uses the first 10 principal components, along with pN of 0.2 and pK of 0.005. This version uses the first 10 principal components, along with pN of 0.2and pK of 0.005 for dimensionality reduction, clustering, and identification of cluster-specific marker genes. The standard Seurat processing pipeline, as described in the previous section, was performed on each non-overlapping inferred age window separately. We found that the default Seurat clustering resolution parameter did not capture the dynamics of the presence of different cell types. For each batch, we clustered the data with a variety of resolutions (from 0.1 to 1.5 in increments of 0.3) and then computed the Silhouette score for finding the best resolution for clustering, we found that 0.4 resolution provides the best Silhouette score for each single-cell fitting in cluster community, as increasing the resolution resulted in lower number of cells per clusters and insignificant results when performing differentially expressed genes text. Finally, we integrated and clustered the data, and from these clusters, we used Seurat’s FindMarkers to iteratively loop through all clusters and identify marker genes.

To interpret the global transcriptional progression of hemocytes and their cell fate decision, we established the cell continuum of cell differentiation and data layers of unspliced and spliced mRNA for the entire data generated with Velocyto CLI (v.0.17.17) according to the CLI usage guide (Velocyto run-smartseq2). The output loom files were combined using “loompy.” The merged loom file was imported into the scVelo package (v1.0.6) (28, 29). The unspliced and spliced mRNA counts of cells from clusters C0–C4 were extracted. We used the “merged.utilis” function in the scVelo pipeline, where cells with low pre-mRNA counts were removed as part of the filtering. In short, the gene-specific velocities are obtained by fitting a ratio between unspliced and spliced mRNA abundances and then computing how the observed abundances change from those observed in a steady state. The ratio of “spliced,” “unspliced,” and “ambiguous” transcripts were calculated, and data were pre-processed using functions for detection of minimum count number, filtering, and normalization using “scv.pp.filter_and_normalise” and followed by “scv.pp.moments” function. The gene-specific velocities were then calculated using “scv.tl.velocity” with mode set to “deterministic” and “scv.tl.velocity_graph” function to generate velocity graph and visualization using “scv.pl.velocity_graph” function. In addition, we used the “scv.tl.recover_latent_time” function to infer a shared latent time from splicing dynamics and plotted the genes along the time axis.

The raw processed data generated for this study have been deposited on Zenodo repository under accession number 7997643 and can be accessed via https://zenodo.org/deposit/7997643 and interactively on https://mubasher-mohammed.shinyapps.io/Sc-drosophila/. The custom code scripts used to analyze data for this study are available at https://github.com/ANKARKLEVLAB/Single-cell-P.f-gametocyte.

We initially characterized the population of TAHs that attach to Ras^V12 salivary glands (SGs) and found that they consist mainly of plasmatocytes and the occasional lamellocyte ( Figures 1A, B ). To identify genes that are differentially regulated in hemocytes (DEGs) from tumor larvae, we profiled the transcriptome of single circulating hemocytes from wild-type and tumor larvae and from TAHs, which were extracted manually from tumorous SGs. Additionally, to characterize the contribution of effector caspases to tumor progression, we analyzed larvae in which the effector caspase inhibitor p35 was co-expressed with Ras^V12: p35;Ras. We had previously shown that similar to Drosomycin (Drs), the expression of p35 in Ras-SGs restores the spherical nuclear shape that is disturbed upon sole Ras expression, although SGs from both combinations show the same size increase as Ras^V12 larvae (20) and ( Figures 2A–D , quantified in Figure 2E ). In contrast to Drs (20), co-expression of p35 still led to hemocyte recruitment ( Figures 2A’–D’ quantified in Figure 2F ). This allowed us to compare TAHs both between wild-type and tumor larvae and in settings with and without active caspases downstream of JNK signaling. TAHs and circulating hemocytes were selected for expression of the plasmatocyte-specific marker Eater ( Figure 3A , see Material and Methods).

Collectively, hemocyte transcriptomes displayed additional plasmatocyte markers indicating successful purification of these macrophage-related cells ( Figures 3B–H ). Plasmatocyte markers include Nimrod C1 (NimC1), serpent (srp), Hemolectin (Hml), and spätzle (spz), several of which have been used as pan-plasmatocyte markers and more recently as markers for specific plasmatocyte subpopulations (30). As expected, the crystal cells and lamellocytes markers lozenge (lz) and Atilla, respectively; were not detected ( Figures 3C, D ). Crystal cell prophenoloxidase 1 (PPO1) was detected in a few cells ( Figure 3H ) (see also below).

Clustering analysis revealed true variation between the single cell transcriptomes for the identified five clusters (at p<0.05 and a difference in expression >2 if not stated otherwise; Figures 4A–C , middle panel), three of which included tumor cells (clusters C0, C2, and C3 ( Figures 4C, D ). Cluster C2 overlaps largely with TAHs ( Figure 4A ). Gene Ontology analysis for biological processes of upregulated genes in the clusters highlighted several gene sets that overlapped between clusters but also some cluster-specific sets ( Figure 4E ). The latter category included genes involved in cell motility (cluster C1), translation (C0), protein folding (C4), and immune response genes for cluster 2 (attached cells) and several categories specific for cell division (C3, circulating cells from both tumor models). Taken together, despite the technically limited number of cells compared to other studies (30), we identified subpopulations of plasmatocytes whose signatures differed significantly depending on genotype (tumor versus wild type) and cell status (attached versus circulating).

To obtain more detailed insight into the gene categories that specified the clusters, we identified the genes whose expression significantly specified the clusters ( Supplementary Tables S1 - S5 ). In parallel, we compared expression strengths between genotypes focusing on the three populations that comprised cells from the tumor models (clusters 0, 2, and 3 in Figures 4D , 5A–C ). For attached hemocytes (C2), the list of the most strongly expressed genes varied between the two tumor models ( Figure 5A ). To follow this up, we scanned cluster 2-specific genes for gene set enrichment and identified two enriched categories (1): pupal adhesion and (2) immune response ( Supplementary Figures S1A, B ). Pupal adhesion was due to the presence of genes that code for salivary gland secretions ( Supplementary Figure S1B ). Since these genes are hardly expressed outside salivary glands (FlyAtlas (31)), we interpret their presence in hemocytes as being passenger transcripts (32), i.e., transcripts that have been taken up by TAHs by way of cellular fragments released from the glands. Alternatively, SG transcripts may have been co-purified together with TAHs during extraction. This latter explanation implies that expression of p35 partially prevented the release of SG fragments from Ras^V12 SGs in line with its anti-apoptotic function. The second category (immune genes) contains genes that are most strongly (although not mutually exclusively) expressed in either tumor model ( Supplementary Figure S1A ). Taken together, we show that hemocytes in both tumor models differ qualitatively from hemocytes in control larvae and that inhibition of effector caspases in tumors reveals additional quantitative differences.

Generally, TAHs in the Ras model show strong expression of AMPs including Dpt(Diptericin)A and B, Att(Attacin)C, and Cec(Cecropin)C. In contrast, p35;Ras^V12 TAHs show higher levels of all three phenoloxidases PPO1–3 ( Supplementary Figure S1A ). This is in line with the previously identified function of the TGF-like protein dawdle (daw), which promotes AMP activation (33) and, indeed, daw clusters with AMP expression in our hands ( Supplementary Figure S1A ). A second TGF member (dpp) with immune-regulatory function (33) is less expressed in tumors compared to wild-type hemocytes. Similarly, the M1 marker iNOS is more strongly expressed in p35;Ras^V12 TAHs, while the M2 marker Arginase (Arg) peaks in Ras TAHs.

For both clusters 0 and 3 ( Figures 5B, C ), the same set of genes was more strongly induced in circulating tumor compared with circulating wild-type hemocytes, although the majority showed stronger expression when effector caspases were repressed (p35;Ras^V12) compared to the tumor-only (Ras) model. Gene set enrichment analysis fully confirmed the enrichment for genes involved in cell division ( Figure 4E ) in cluster 3 (circulating cells from both tumor models) but failed to deliver significant returns for the category “Biological process” for cluster 0. In line, the quantification of hemocyte proliferation fully confirmed an increase in cell divisions in both tumor models including the additional increase in the p35;Ras^V12 model ( Figure 5D ). This may contribute more to the differentiation of TAHs rather than circulating cells, which show similar counts (18).

Finally, when searching for pathways compatible with expression in the five clusters using the Reactome database (https://reactome.org), an enrichment for genes involved in neutrophil degranulation was found for cluster 0 ( Supplementary Figure S1C ), which, when using a split analysis, also indicated that this pathway was more strongly expressed in p35;Ras^V12 hemocytes. Pseudotime analysis using RNA velocity ( Figures 6A–D ) indicated that TAHs and circulating hemocytes from both tumor models follow distinct pathways, although due to the restricted number of TAHs, a rarer intermediate population may have been missed.

Recently, we and others have characterized the activation of crystal cells and the subsequent release of leaderless PPO2 into the extracellular environment via cell rupture (34, 35). Functionally, crystal cell activation bears similarities to pyroptotic cell death in mammals, including its dependence on caspase activity, which can be inhibited by p35 (36). While our findings explained how PPO2 was released into the hemolymph, the mechanism for secretion of PPO1, which also lacks a signal peptide, remained obscure. Since we found that PPO2 was enriched in TAHs ( Figure 7A , Supplementary S1A ), we wondered which secretion mechanism led to their spread to other plasmatocytes. In line with crystal cell rupture and similar to passenger transcripts from SGs, the mRNA encoding PPO2 was detected in plasmatocytes both as a mature (spliced) form and as a non-spliced (nuclear) immature mRNA ( Figures 7C, D ). This indicates that phagocytic hemocytes had access to both the nuclear and the cytosolic fraction of ruptured crystal cells. In contrast, PPO1 transcripts were only detected in their mature form ( Figure 7B ). We conclude that PPO2 transcripts in TAHs originate either from pyroptotic crystal cells or TAHs, which have acquired crystal cell characteristics, while PPO1 mRNA originates from live crystal cells. Taken together, this explains why both PPO1 and PPO2 had been located on Ras SGs (19) despite the apparent absence of crystal cells ( Figure 1 ).

In addition to PPO1, two transcripts (CG 13962, and the transcription factor Relish) were detected only in their mature form, indicating a possible origin outside TAHs.