Single-cell profiling of the developing embryonic heart in Drosophila

Understanding the genetic control of heart development is essential to understanding the mechanisms underlying congenital heart diseases, the most common birth defect (Pierpont et al., 2018). To discover the evolutionary conserved genes and genetic pathways required for heart development, a variety of model systems have been employed over the past 30 years, including fly (Drosophila), zebrafish (Danio rerio), frog (Xenopus) and mouse (Sperling, 2011). However, until recently, our knowledge about the cell types that make up the heart had been largely dependent on observations of cell morphology. The emergence of single-cell RNA sequencing (scRNA-seq) technology has made it possible to discover new cardiac cell types and their markers, and understand cross-species correlation of cardiac cell types at the transcriptional level.

For these reasons, scRNA-seq technology has been used to study heart development in many of the model systems including mice (DeLaughter et al., 2016; Farbehi et al., 2019; Gladka et al., 2018; Goodyer et al., 2019; Hu et al., 2018; Jia et al., 2018; Lescroart et al., 2018; Li et al., 2016, 2019a,b; Skelly et al., 2018), zebrafish (Burkhard and Bakkers, 2018; Honkoop et al., 2019; Weinberger et al., 2020; Yuan et al., 2018) and sea squirt (Ciona robusta) (Wang et al., 2019), as well as humans (Asp et al., 2019; Cui et al., 2019; Litvinukova et al., 2020; Tucker et al., 2020). Previously unreported cell types, including specialized subtypes of cardiomyocytes, fibroblasts, endothelial cells and epicardial cells, have been identified through these studies. In addition, these studies have led to the identification of many new unique markers for each cardiac cell type, as well as the characterization of cell lineage trajectories. Remarkably, to date, scRNA-seq has not been performed in a developmental context for the Drosophila heart, a widely used and valuable model system for studying heart development.

The Drosophila heart shares many similarities with the human heart at the molecular, cellular and genetic levels, and has been used extensively to model many types of heart disease, including congenital heart disease (Basu et al., 2017; Bier and Bodmer, 2004; Klassen et al., 2017; Olson, 2006; Wolf et al., 2006; Zhu et al., 2017a,b,c). The Drosophila heart is not just a simple contractile tube. It has distinct morphological features with newly recognized functional equivalence to the human heart. For example, the Drosophila heart can be separated into the anterior aorta and the posterior heart chambers. The anterior aorta is labeled by mthl5 (also called Gia), which encodes a G protein-coupled receptor (GPCR) required for aorta differentiation (Patel et al., 2016). The aorta and heart chambers are separated by an aortic valve (Rotstein and Paululat, 2016), which ensures posterior to anterior blood flow. Besides the cardioblasts that form the heart tube (Haack et al., 2014) several types of pericardial cells, which are required for heart function (Fujioka et al., 2005), surround the Drosophila heart tube, including the even-skipped (eve)-expressing pericardial cells (Han et al., 2002). Adjacent pericardial nephrocytes play essential roles in filtration and protein reabsorption, similar to mammalian podocytes (Zhang et al., 2013a,b). However, it remains unclear how many pericardial cell types are present during embryonic heart development, and their origin; with inconsistent literature on cell types, numbers, marker expression and location, as well as their naming, urging clarity.

To investigate the transcriptome of the developing Drosophila embryonic heart at the single-cell level, we collected and analyzed 34,203 green fluorescent protein (GFP)⁺ cells obtained from over 3000 Hand-GFP transgenic fly embryos. These embryos were stage-synchronized and collected at five consecutive time points (developmental stages 13-16) before the formation of the fused, beating heart tube. This dataset provides a high-resolution atlas of cardiac cell types of the Drosophila embryonic heart. We identified six unique cardiac cell types with hundreds of unique markers. The data further revealed that fly cardioblasts combine the molecular features of the mammalian first and second heart fields. Using antibodies against the cell-type-unique markers we confirmed and localized the cardioblasts and different pericardial cell types within the embryonic fly heart.

To ensure that our scRNA-seq data contained all Drosophila heart cells during multiple stages of embryonic development, we used the Hand-GFP transgenic line (Table S8). Hand is an evolutionarily conserved basic helix-loop-helix (bHLH) transcription factor required for cardiogenesis in fly (Han and Olson, 2005; Han et al., 2006). Hand-GFP is the only known fluorescent marker expressed in all heart cells from embryonic stage 13 through to adult flies (Han and Olson, 2005; Han et al., 2006). Hand-GFP flies were used for embryo collection over 5 consecutive hours from stage 13 to stage 16, monitored by tracking the migration of Hand-GFP cells from bilateral rows to a fused dorsal heart tube. More than 3000 stage-synchronized embryos were collected, dissociated into single cells and passed through flow cytometry to sort and collect the GFP⁺ cells (Fig. 1A). Around 4% of total embryonic cells were GFP⁺ (Fig. 1B) and these cells were further processed through a 10x chromium controller for scRNA-seq. To ensure the enrichment of heart cells in the sorted GFP⁺ cells, we verified the expression of cardiac genes, including Hand and tinman (tin), in the GFP⁺ and GFP⁻ cell populations using quantitative real-time (qRT)-PCR. GFP⁺ cells showed much higher expression of Hand and tin (Fig. 1C), thus demonstrating that heart cells were greatly enriched among the sorted GFP⁺ cells.

After removing doublets and debris (Fig. S1), a total of 34,203 Hand-GFP⁺ embryonic cells were obtained, sequenced and subjected to cluster analysis (Fig. 1D). First-level cluster analysis revealed nine different types of cells, including cardiogenic progenitors, muscle, nervous system, blood, fat body, trachea, gut, nephrocytes and salivary gland. The largest clusters of Hand-GFP⁺ cells from first-level clustering were the muscle cells (36.30%) and neuronal cells (27.21%). Around 12.6% of Hand-GFP⁺ cells (a total of 4315 cells) were found to be cardiogenic progenitor cells (Fig. 1D). This is consistent with previous observations that Hand is expressed at higher levels in heart cells, lymph gland, nephrocytes and visceral mesoderm muscle cells (Han and Olson, 2005), and shows lower expression in nervous system cells, blood cells, fat body, trachea, gut and salivary gland (Fig. 1E). Using the scRNA-seq data, we identified unique markers for each of the nine different Hand-GFP⁺ cell types (Fig. S2). Examples of genes used to identify the different cell clusters: cardiac progenitors, pericardin (prc) and lonely heart (loh); muscle, sallimus (sls) and bent (bt); nervous system, found in neurons (fne) and scratch (scrt); blood, Malic enzyme (Men) and GIP-like (Gip); fat body, apolipophorin (apolpp; also known as Rfabg) and viking (vkg); trachea, Osiris (Osi6 and Osi7) and Ecdysone-inducible gene L1 (ImpL1); gut, mesh and Tetraspanin 2A (Tsp2A); nephrocytes, Amnionless and Cubulin (Cubn); salivary gland tissue, CG8708 and CG14756 (Fig. S2).

Previous studies have shown that lymph gland progenitors, pericardial nephrocytes and cardioblasts arise from the cardiogenic mesoderm (Mandal et al., 2004); in Drosophila, bipotential progenitors have been identified that give rise to blood cells and nephrocytes (Morin-Poulard et al., 2022). Our scRNA-seq results support these findings as the cardiac and blood progenitor cells exhibited a close lineage relationship; they mapped to a single cell cluster on the uniform manifold approximation and projection (UMAP) feature plot (Fig. 1D). Next, we isolated these cells for secondary subcluster analysis. The subclusters revealed that of the cardiac cell types, the hematopoietic lymph gland cells were closest in lineage to the Odd skipped (Odd)+ pericardial cells (OPC) (Odd⁺ Tin⁻) (Fig. S3A). This is especially clear when visualizing the OPC and lymph gland cells by developmental stage (Fig. S3B). At the earlier stages (13-14) the lymph gland progenitor cells and OPC coalesce; over developmental time they gradually separate (stage 15), to eventually form two defined clusters (stage 16), indicative of two further differentiated cell types.

In addition to odd, both cell types share the expression of several other genes, including zfh1 and GILT1 (Fig. S3C). Zn finger homeodomain 1 (encoded by zfh1) is crucial for pericardial cell specification and lymph gland cell differentiation (Sellin et al., 2006); whereas, Gamma-interferon-inducible lysosomal thiol reductase 1 (encoded by GILT1) plays a role in the innate immune response to bacterial challenges in Drosophila hemocytes (fly equivalent of blood cells) (Ma et al., 2017). Indeed, in fly, lymph gland hematopoiesis and its associated cells play a vital role in the immune response (Vlisidou and Wood, 2015; Williams, 2007). The close lineage between lymph gland progenitors and OPC in our data suggests that the OPC as they mature – eventually into nephrocytes (fly equivalent of kidney cells) – might also play a role in the defensive response. This notion is supported by a recent study which showed that nephrocytes carry out immunoregulatory functions through the degradation of microbiota-derived peptidoglycan (Troha et al., 2019). These results together demonstrate developmental and functional similarities between hematopoietic cells (lymph gland) and nephrocytes (OPC).

We separated the heart cells from the cardiogenic progenitor cluster, then determined the heart cell-type heterogeneity during development. This revealed six distinct cardiac cell second-level clusters (Fig. 2A). All six expressed cardiac markers such as zfh1 consistent with their cardiac cell identity (Fig. 2B) (Su et al., 1999). Cluster 1, cardioblast, can be defined by the expression of neuromancer 1 (H15) (Fig. 2C) and many other known cardioblast markers (Table S1) (Miskolczi-McCallum et al., 2005; Qian et al., 2005; Reim et al., 2005). In total, we detected 247 cardioblasts in our assay with a transcriptional profile highly distinct from that in each of the five pericardial cell types (Fig. 2B,C).

The pericardial cell types comprised the EPC (Eve⁺ Tin⁺), CTPC (Ct⁺ Tin⁺), wing heart pericardial cells (WHPC; Eve⁺ Tin⁻), end of the line pericardial cells (ELPC; Ct⁺ Tin⁻) and OPC (Odd⁺ Tin⁻) (Fig. 2A). The EPC and OPC have been most documented. Notably, although each displayed unique differentially expressed marker genes (Tables S2-S6), all the pericardial cells can be classified based on the expression of just four transcription factors: eve, tin, cut (ct) and odd. The expression patterns of ct and eve were strictly complementary in the pericardial cells (Fig. 2B). Eve is a homeobox-containing transcriptional factor required for specification in a subset of pericardial cells and for heart function (Fujioka et al., 2005; Zmojdzian et al., 2018). Ct has been implicated in developmental regulation of neurons and muscle cells (Blochlinger et al., 1990; Sudarsan et al., 2001), but its role in pericardial cell specification has not been investigated. Odd is a DNA-binding transcription factor that is expressed in a subpopulation of pericardial cells that will become nephrocytes (Ward and Coulter, 2000; Ward and Skeath, 2000). The different combinations of these transcription factors may play a key role in determining the identities of the different pericardial cell types.

To understand the dynamic transcriptome changes taking place in cardiac progenitor cells during cardiogenesis, embryos were stage-synchronized and collected at five consecutive time-points ranging from stages 13 to 16 (Fig. 2D). These stages reflect major events in cardiogenesis: (1) dorsal vessel cells derived from the dorsal-most mesoderm form two bilateral lines of cardiac progenitors at stage 13; (2) the progenitors migrate toward the dorsal midline at stages 14-early, 14-late and 15; and (3) the progenitors meet at the dorsal midline at stage 16, where they mature into cardioblast cells and pericardial cells, shortly after which the heart begins to beat (Fig. 2D).

The ratio across the different cardiac cell types was largely maintained from stage 13 to 16 (Fig. S4A,B). Stages 13-15 already showed clearly distinct cell clusters; yet, their gene expression-based cluster maps showed considerable spatial overlap, whereas the stage 16 clusters appeared to be shifted from the others (Fig. 2E; Fig. S4C). These data suggest that the fate of the heart cells is determined before stage 13, reflected in their largely unchanging transcriptional profiles; until the cells reach the dorsal midline and start fusion into a linear heart tube (stage 16), which evoked transcriptional changes.

Cardioblasts, the equivalent of human cardiomyocytes, make up the contractile heart tube in Drosophila. The cardioblast cluster showed specific and high levels of expression for H15, an established marker for cardiomyocytes, as well as for midline (mid) and Sulfonylurea receptor (Sur) (Fig. 3A). This is consistent with their role in cardioblast differentiation and heart development (Akasaka et al., 2006; Miskolczi-McCallum et al., 2005; Qian et al., 2005; Reim et al., 2005). We performed Gene Ontology (GO) analysis on the top 100 most differentially expressed genes in the cardioblasts. Muscle structure development, muscle cell development, cardioblast differentiation and cell fate commitment, cell junction assembly and organization were detected among the biological functions (Fig. S5A). Signaling pathways related to mammalian cardiac muscles including I band, sarcomere, Z disc, myofibril and cell junction were observed among the enriched molecular functions in the fly cardioblasts (Fig. S5B). In addition to these expected cardiac muscle markers, we also identified previously unreported marker genes for the cardioblasts, which could lead to the discovery of new pathways important for cardiac muscle formation and maturation. These genes include Anaplastic lymphoma kinase (Alk), Chondroitin polymerizing factor (Chpf), phurba tashi (phu), Syntrophin-like 2 (Syn2), Delta and CG3961 (Fig. 3B). All of them have highly conserved human homologs (Table 2), suggesting that they might have previously unknown and evolutionarily conserved roles in heart development.

Alk is a receptor tyrosine kinase (Table 2) and has been shown to regulate mesoderm and visceral muscle development, photoreceptor development and axon guidance (Bazigou et al., 2007; Englund et al., 2003; Lee et al., 2003). Immunostaining for Alk was consistent with the scRNA-seq findings and showed cardioblast-specific expression (Fig. 3B,C). Notably, the intensity of staining for Alk was generally higher in the Svp⁺ cardioblasts compared with the Tin⁺ cardioblasts (Fig. 3C). Future work to verify and study the function of the other cardioblast-specific marker genes will confirm and identify more genes with crucial roles during heart development.

In mammals, the first and second heart fields contribute to different regions of the heart along the anteroposterior axis, and are marked by distinct genetic markers (Buckingham, 2016). Most of the mammalian first and second heart field markers are evolutionarily conserved in flies (Table 1), but how these markers are distributed among the fly cardioblasts has not been examined systematically. We isolated the cardioblast cluster as an independent Seurat object for third-level cluster analysis (first-level: clustering of all the Hand-GFP⁺ cells, Fig. 1D; second-level: clustering of all Hand-GFP⁺ cardiac progenitors, Fig. 2A). The fly homologs, of all the mammalian first and second heart field markers that we examined, were expressed in cardioblasts distributed across the cluster (Fig. 4A,B). These markers include well-known first heart field markers such as heart and neural crest derivatives expressed 1 (HAND1), T-box transcription factor 5 (TBX5) (Paige et al., 2015), TBX20 and myosin light chain 7 (MYL7) (Singh et al., 2005), as well as recently discovered first heart field markers such as platelet-derived growth factor receptor A (Pdgfra) (Xiong et al., 2019) and keratin 8 (Krt8) (de Soysa et al., 2019) (Fig. 4A). Fly homologs of well-established second heart field markers, such as islet1 (ISL1) (Paige et al., 2015), fibroblast growth factor 10 (FGF10) (Paige et al., 2015) and GATA binding protein 4 (GATA4) (Liu et al., 2019; Zhou et al., 2017), were also expressed in cells across the cardioblast cluster (Fig. 4B). Together, our data demonstrate that the key markers differentiating the mammalian first and second heart fields are conserved in Drosophila but show even expression among the fly cardioblasts.

Even though the Drosophila cardioblast combines the molecular features of both mammalian heart fields, the presence of two different types of cardioblasts in fly is well-established. These are the Svp⁺ cardioblasts that function like the mammalian cardiac inflow tract and the Tin⁺ cardioblasts that are equivalent to mammalian cardiomyocytes (Molina and Cripps, 2001). To gain insight into the different types of cardioblasts, we isolated them for subcluster analysis. We could not obtain two well-separated clusters but could distinguish the Svp⁺ and Tin⁺ cell groups within the cardioblast cluster (Fig. 5A-C). The findings demonstrate the close lineage relationship between Svp⁺ and Tin⁺ cardioblasts, while recognizing their individual characteristics.

The Svp⁺ cardioblasts, exclusively, expressed known ostia markers, including T-box genes Dorsocross1 (Doc1) and Doc3 (Reim and Frasch, 2010), Wnt oncogene analog 4 (Wnt4) (Chen et al., 2016) and wingless (wg) (Reim and Frasch, 2010) (Fig. 5D; Figs S6,S7). This high-dimensional analysis revealed new markers for the Svp⁺ cardioblasts, such as Transglutaminase (Tg), pyramus (pyr), Glucose transporter 4 enhancer factor (Glut4EF), headcase (hdc) and Heat shock protein 23 (Hsp23) (Fig. 5F). Tg and pyr were exclusively expressed in the Svp⁺ cardioblasts (Fig. 5D,F; Figs S8,S9). Tg protein possesses peptide cross-linking capabilities and may be involved in hemolymph coagulation (Lindgren et al., 2008). Immunostaining showed that Tg protein was highly expressed in the Svp⁺, but not the Tin⁺, cardioblasts (Fig. 5G). These findings are in line with an in situ hybridization study in Drosophila that reported strong Tg expression in Svp⁺ cells by stage 16, but low expression in Tin⁺ cardiac cells (Ikle et al., 2008). The gene pyr encodes a fibroblast growth factor (FGF) that is a ligand for the FGF receptor; it is expressed in ectodermal stripes adjacent to pericardial cells in the dorsal mesoderm (Stathopoulos et al., 2004). FGF signaling has been implicated in migration and differentiation, suggesting that pyr could play a role in the migration of cardioblasts and the differentiation of Svp⁺ cardioblasts.

The high-dimensional analysis similarly revealed markers for the Tin⁺ cardioblasts, including faint sausage (fas; CG17716; tei), NHP2, Myc (dm), tenectin (tnc), Netrin-B (NetB), tincar (tinc) and hoi-polloi (hoip) (Fig. 5E,F; Fig. S7). The gene fas encodes an immunoglobulin domain-like cell adhesion molecule that is required for the specification and alignment of cardioblasts (Haag et al., 1999). NHP2 encodes a component of the nucleolar ribonucleoprotein complex that catalyzes pseudo-uridylation of rRNA to regulate the translation of specific mRNAs in Drosophila ovary (Morita et al., 2018), but its role in heart development remains unclear.

Notably, all the pericardial cells could be classified based on just four transcription factors: tin, eve, odd and ct (Figs 2B,C, 6A). To examine the spatial localization of each, we performed immunofluorescence staining in the Hand-GFP embryos. This revealed that the three major pericardial cell types (EPC, Eve⁺ Tin⁺; OPC, Odd⁺ Tin⁻; CTPC, Ct⁺ Tin⁺) were distributed along the entire heart tube, i.e. both the thoracic and the abdominal segments (Fig. 6B-D), whereas the WHPC (Eve⁺ Tin⁻) and ELPC (Ct⁺ Tin⁻) showed a polarized distribution. Eight WHPC were localized only in the thoracic segments (Fig. 6B; Fig. S10A). The ELPC form a pair of Ct⁺ Tin⁻ PC located at the rear end of the A7 segment, hence their tentative name as end of the line PC (Fig. 6D); an unexpected find. Typically, each hemisegment has two EPC, three CTPC and four OPC (Fig. S10B-D). Aside from differences in localization and cell number, the different PC displayed differences in cell shape: EPC (Eve⁺ Tin⁺) and WHPC (Eve⁺ Tin⁻) were typically round; CTPC (Ct⁺ Tin⁺) and ELPC (Ct⁺ Tin⁻) were typically spindle shaped.

The PC also showed distinct distribution patterns along the dorsal-ventral axis. Considering the migration of cardiac cells during development, we investigated spatial localization at both stage 15 and 16. At stage 15, the heart tube has not yet fused, evident in the distance between the cardioblasts and PCs on the dorsal and ventral sides; with cardioblasts closer to the dorsal midline forming the front line of the migrating cardiac progenitors (Fig. 7A,B). Most EPC (Eve⁺ Tin⁺) located to the lateral side of the cardioblasts and above the CTPC (Ct⁺ Tin⁺) and OPC (Odd⁺ Tin⁻). Both CTPC (Ct⁺ Tin⁺) and OPC (Odd⁺ Tin⁻) typically located under the EPC (Eve⁺ Tin⁺), with OPC (Odd⁺ Tin⁻) located farthest from the dorsal midline.

At stage 16, the heart tube closes as the cardioblasts connect (Fig. 7C,D). At this stage, the WHPC (Eve⁺ Tin⁻) located to the dorsal side of the embryo in the thoracic segments. EPC (Eve⁺ Tin⁺) located on the top (dorsal) or the lateral side of the cardioblasts, whereas CTPC (Ct⁺ Tin⁺) had migrated underneath (ventral to) the cardioblasts (lower side of the heart tube). At the rear end of segment A7, the two ELPC (Ct⁺ Tin⁻) located underneath the cardioblasts. The OPC (Odd⁺ Tin⁻) located alongside the CTPC (Ct⁺ Tin⁺), underneath the cardioblasts and EPC (Eve⁺ Tin⁺).

Here, we have presented a comprehensive analysis of the transcriptional dynamics during embryonic development of the fly heart. Owing to its small size, the fly embryo cannot be dissected to remove its heart. Therefore, we used fluorescent-labeled cells in a flow cytometry strategy. This method does have technical constraints, like sample handling, cell death and gating stringency, that reduce the final cell number obtained per embryo. To drive GFP expression, we used Hand, which is a key transcription factor for cardiogenesis and hematopoiesis in Drosophila (Han and Olson, 2005; Han et al., 2006). In the Drosophila heart, Hand is expressed from embryonic stage, at which cardiac cell fate has already been determined. To study previous stages, earlier markers can be used, such as tin, which is expressed in the early mesoderm that develops into the heart and other organs (Bodmer et al., 1990). The Hand protein sequence and expression pattern are highly conserved across all animals with a heart, including Drosophila, sea urchin, zebrafish, frog, chicken and all mammals (Han et al., 2006; McFadden et al., 2000; Srivastava et al., 1995, 1997; Yelon et al., 2000). Even though Hand is mostly known for its expression in the heart, its expression in several other tissues has been identified (Han and Olson, 2005; Lo et al., 2007; Popichenko et al., 2007). In line with these earlier findings, our scRNA-seq data demonstrate that Hand-GFP is expressed highest in cardiogenic progenitors and nephrocytes, followed by muscle cells, and the weakest expression is in neurons, blood, fat body, trachea, gut and salivary glands (Fig. 1D), with heart cells making up ∼8% of the total Hand-GFP⁺ cells. To facilitate further investigation into the role of Hand and the dynamic developmental transcription profiles in these tissues, these data have been made available (https://singlecell.broadinstitute.org/; all Hand-GFP⁺, SCP1469; GEO accession number GSE168774).

The heart cells labeled by Hand-GFP were defined by high expression levels of known cardiac markers such as tin, loh, zfh1 and tup. Within this heart cell cluster, we identified five types of pericardial cells and the cardioblasts (Fig. 2). The five types of pericardial cell clusters could be distinguished based on the expression of just four markers: Tin, Eve, Odd and Ct. We used these markers to locate each of the distinct PC cell types in the developing fly heart, thus providing spatial information and confirming a pair of unique cells located at the end of the linear heart tubes we tentatively named ELPC (Fig. 8). Reports to date have been inconsistent, and incomplete in understanding, on the cell types, their naming, cell numbers, marker expression and location during fly heart development. The EPC and OPC are well known; however, information on the CTPC and WHPC is much more limited, with no data on the newly identified ELPC. We are currently using the unique PC type markers identified in our scRNA-seq data to target each specific PC type using the Gal4-UAS and/or split-Gal4 system (Luan et al., 2020) to express fluorescent markers for visualization and lineage tracing, and to induce cell death to investigate their function. The data have uncovered many new marker genes expressed by specific cardiac cell types; however, the role of many of these markers in heart development remains unknown. For example, we detected a single cardioblast cluster marked by the expression of known markers H15, mid and Sur, as well as the previously unreported markers Alk, Chpf, CG3961 and Delta. Of these only DLL1, the human homolog of fly Delta, has been implicated in heart disease including dilated cardiomyopathy and heart failure (Norum et al., 2016, 2017). Therefore, we are currently using fly to investigate their functions in the context of heart development and, ultimately, disease. Furthermore, it would be interesting to see whether Alk and other heart-expressed receptor tyrosine kinases, like Htl, act in conjunction or carry out distinct roles. Notably, unlike Alk, htl expression is not cardioblast-specific in favor of specialized functionality.

Altogether, these fly single cell transcriptional profiles combined with temporal (five stages) and spatial (immunochemistry) information, provide a comprehensive resource for tracking gene expression in different heart cell types during embryonic heart development.

ScRNA-seq data of the (developing) heart has been available for vertebrate species such as mice (DeLaughter et al., 2016; Farbehi et al., 2019; Gladka et al., 2018; Goodyer et al., 2019; Hu et al., 2018; Jia et al., 2018; Lescroart et al., 2018; Li et al., 2016, 2019a,b; Skelly et al., 2018), zebrafish (Burkhard and Bakkers, 2018; Honkoop et al., 2019; Weinberger et al., 2020; Yuan et al., 2018) and human (Asp et al., 2019; Cui et al., 2019; Litvinukova et al., 2020; Tucker et al., 2020). However, heart size in these species limits the scRNA-seq to tissue sampling, which inherently introduces bias. Owing to the small size of Drosophila, the entire heart, including every single cell, can be sequenced. Our analysis confirmed the conservation of the genes and pathways involved in the earliest stages of heart development from fly to human. However, it was not clear whether the Drosophila heart shares the mammalian division into first and second heart fields. Cells that originate from the first heart field develop into cardiomyocytes of the left ventricle, whereas those of the second heart field turn into cardiomyocytes located in the right ventricle, atria, and the outflow and inflow tracts (Buckingham, 2016). Even sea urchin, which sit at the phylogenetic base of vertebrate origin, show evidence of these two heart lineages (Wang et al., 2019). Although the fly heart is made up of distinct structures (anterior aorta, posterior heart chambers, aortic valve and ostia), we did not find any indication of separate heart fields in fly. Fly homologs of key markers of the mammalian first and second heart fields are conserved but show expression in cells distributed across the cardioblast cluster (Fig. 4). This suggests that separation of heart fields might be distinct to the phylum Chordata and evolutionarily arose after the split from flies, thereby implicating progressive lineage restriction. However, we cannot rule out that the expression of homologs of the first and second heart field markers will show preferential expression at later stages of development (i.e. after stage 16).

One previous scRNA-seq study of the fly heart, the Fly Cell Atlas, has been published (Li et al., 2022). However, as this massive project covered all major organs in the fly in both sexes, this study was necessarily limited to the adult stage. Here, we specifically looked at the fly heart across multiple embryonic stages. This has revealed new insights into fly heart development, for example the combined first and second heart fields in fly cardioblasts, that Svp⁺ and Tin⁺ cardioblasts share very close transcription profiles but with unique markers, and the early specification of the heart cells. Notably, whereas some Svp⁺ cardioblasts expressed Tin⁺ cardioblast markers, such as tin and fas, cells of the Tin⁺ cardioblast subcluster never expressed Svp⁺ cardioblast markers, including svp, Doc2, Doc3 and Tg (Fig. 5F; Fig. S7). This is consistent with the transient co-expression of svp and tin in some cardioblasts during stages 11 to 13 (Lo and Frasch, 2001), in which tin acts as a direct activator of svp via a conserved cardiac enhancer (Ryan et al., 2005). We did not detect the Tin⁺ cardioblast subtypes, like Tin⁺ Lbe⁺ cardioblasts and Tin⁺ Lbe⁻ cardioblasts (Ahmad, 2017; Han and Bodmer, 2003), possibly due to the lower-than-expected total number of cardioblasts collected in the stage 16 embryos. At this stage the heart tube starts to fuse, resulting in increased cell-cell junctions and connective tissue among the cardioblasts, making it more difficult to collect single cells at later stages. The scRNA-seq findings of early heart cell specification appear to coincide with the migration patterns, which were confirmed by our spatial immunohistochemistry data. Cardiac progenitors originate from the cardiac mesoderm, then undergo a series of cell divisions and differentiate into two major classes of heart cells, the contractile cardiac cells and the non-contractile PC. Our data show that these major classes specify before stage 13, followed by extensive migration – segmental arrangement of cell types is evident at stage 15 (Fig. 8) – before formation of the functional heart tube at stage 16.

Through the combined scRNA-seq, immunostaining and confocal imaging, we comprehensively resolved the cell types of the developing embryonic fly heart at the transcriptome level, and then demonstrated the existence and exact spatial localizations for each of these cell types. These findings have increased our understanding of the processes driving fly heart development and have revealed differences and similarities with heart development in vertebrate models. The data can also inform future explorations into defining developmental processes, for example, to study the collective migration of cardiac progenitor cells. This process has shown to be oscillatory and dynamic (Balaghi et al., 2023; Zhang et al., 2018, 2020) and appears to be well-suited to single-cell transcriptomics. Recent studies into the mechanics of cardiac progenitor migration in Drosophila have revealed the roles of several cell-cell adhesion molecules and supracellular actin structures (Balaghi et al., 2023; Zhang et al., 2018, 2020). Our data indicated the presence of fas, which is required to obtain the correct cardioblast number and alignment (Haag et al., 1999); it would be interesting to combine these datasets.

Drosophila has a strong track record as a model system for heart disease research. Along with the Fly Cell Atlas scRNA-seq data for the adult fly heart (Li et al., 2022) and other embryonic fly heart data (Vogler et al., 2021 preprint), our findings strengthen the fly as a model for studying human congenital heart disease and other genetic cardiac diseases – through the identification of new heart cell types and hundreds of cardiac markers for each. The next key step will be to elucidate the developmental lineages for each cell type from embryo to adult, and the functions of these cell types during development and cardiac disease conditions. As the key genes regulating fly heart development identified with our data are likely candidate genes for disease causation, we anticipate this resource will inspire new avenues of research. These mechanistic studies are needed to fully understand the cellular and molecular processes underlying cardiac development and associated diseases. Furthermore, we are using the data from fly and datasets from other model systems to develop tools that facilitate cross-species interpretation of findings to benefit research into cardiac development and diseases.

Flies were maintained at 25°C and standard food (Meidi Laboratories). The Hand-GFP and 4xHand-Gal4 fly lines have been generated and described previously (Han and Olson, 2005; Table S8).

Drosophila melanogaster expressing Hand-GFP were used to obtain embryonic heart cells, with w¹¹¹⁸ flies as control for setup of flow cytometry. One-week-old flies were used for egg collection. To collect synchronized embryos, grape juice agar plates were changed three times, with 1 h intervals, before egg collection – at which point the eggs were collected every hour and kept at 25°C before transfer to 18°C for further development until desired stages.

Embryos (female and male, indistinguishable at this stage) were collected and dissociated into single-cell suspension according to previously published methods (Salmand et al., 2011) with some modifications. Briefly, collected embryos were washed thoroughly with ice-cold 1× phosphate buffered saline (PBS) and then dechorionated in fresh 4-5% bleach (sodium hypochlorite solution; Sigma-Aldrich) for 3 min at room temperature. Dechorionated embryos were extensively rinsed with water. Then embryos (∼150 mg) were drained on absorbent paper and, using a paint brush, transferred to a Dounce homogenizer (15 ml Tissue Grinder; Dounce). A 5 ml suspension of embryos in Schneider medium was subjected to mechanical disruption using the Dounce homogenizer; achieved by five strokes using a pestle. Homogenized samples were transferred to 15 ml centrifuge tubes. Collagenase (Type I) (20 U/ml; Worthington), 100 U total, was added to the Schneider medium with the homogenized embryos and the samples were incubated at room temperature for 10 min. Next, 10 μl of ethylenediaminetetraacetic acid (EDTA) (0.5 M; final concentration of 1 mM) and 0.5 ml 2.5% trypsin were added to 5 ml cell solution to digest the samples for another 20 min at room temperature. The trypsin reaction was stopped by adding 20% fetal bovine serum (FBS) and cells were collected using 70 μm mesh for filtration. Next, the cells were pelleted by centrifugation (1000 g, 5 min, 4°C) and resuspended in 2 ml pre-chilled artificial hemolymph. Cells were further filtered using a 40 μm mesh, resuspended in 0.5-1.0 ml, of which 10 μl was used for visual examination under a fluorescence microscope (Apotome; ZEISS).

Isolated embryonic cells were sorted by flow cytometry (BD Influx cell sorter). Collected cells were pelleted by centrifugation (1000 g, 5 min, 4°C) and resuspended in 50 μl pre-chilled artificial hemolymph. Cell suspensions were examined under a fluorescence microscope (Apotome; ZEISS) and cell numbers were counted using a Neubauer Hemocytometer; at which point the samples were ready for downstream experiments.

Single-cell libraries were generated according to the manual of Chromium Single Cell 3′ Reagent Kits v3 (10x Genomics). Single cells were partitioned into Gel Bead-In-Emulsions (GEMs) in a Chromium Single Cell Controller (10x Genomics) with a target cell number range of 5000-8000. Sequencing libraries were constructed according to the methods described previously (Fu et al., 2020). Briefly, GEM samples were immediately reverse transcribed into cDNA and proceeded to downstream cleanup. The cDNAs were enzymatically fragmented and went through end repair and 3′ A-tailing after further amplification. Adaptor ligation was performed after double-sided size selection by SPRIselect beads and followed by sample index PCR. Unique molecule identifier (UMI) sequences, 10x barcode sequences, sequencing primer P5 and P7 on both ends and sample index sequences were all added to the cDNA samples. Finally, the library was purified and 1 μl of the amplified cDNA libraries was quality assessed and quantified by BioAnalyzer High Sensitivity Chip (Agilent) and Qubit dsDNA HS (High Sensitivity) Assay (Invitrogen). The libraries were sequenced by Novagene using an Illumina HiSeq 2500.

Raw sequencing reads of Chromium single-cell 3′ RNA-seq output were processed and a data matrix was generated using the Cell Ranger pipeline from 10x Genomics (http://10xgenomics.com). Cell cluster analyses were performed using the Seurat package (version 3.1.2) as described in the tutorials (http://satijalab.org/seurat/) (Butler et al., 2018). Previous results have reported that mitochondrial transcript abundance is much higher in tissues with high-energy demands, such as heart and muscle. Especially in the heart, mitochondrial transcripts comprise almost 30% of total mRNA, in contrast to tissues with lower energy demands (such as ovary, prostate, testes, lung and white blood cells) in which ∼5% mitochondrial transcripts contribute to their total mRNA content (Mercer et al., 2011). Therefore, we filtered out cells with more than 30% expression of mitochondrial genes. For the detection of potential heterotypic doublet artifacts from our scRNA-seq result, we used DoubletFinder 2.0.3 (McGinnis et al., 2019) without specifying any ground-truth information. The analysis requires an estimated homotypic doublet proportion, for this we used 7.5% (Fig. S1). We performed a global-scaling normalization method (LogNormalize), which normalizes the gene expression of each cell by the total expression, multiplies this by a scale factor (10,000 by default) and log-transforms the result in Seurat. Unsupervised UMAP for dimension reduction (McInnes et al., 2018) was used for non-linear dimensional reduction and cluster results were visualized in two-dimensional plots. We defined the cell types by cross-referencing their marker genes with known markers, as well as checking the function and tissue distribution of those marker genes in FlyBase (https://flybase.org) (Larkin et al., 2021). The FindAllMarkers function from Seurat was used to identify differentially expressed marker genes across different cell types, and the top ten genes were plotted in heatmaps.

To gain obtained detailed information on the embryonic heart cells, the data for the ‘cardiogenic progenitors’ cluster (cluster 1, Fig. 1D) was further defined by secondary clustering using the subset function from Seurat. Cardioblast cells were identified by their expression of H15 and mid, whereas pericardial and lymph gland cells were identified by expression of tin, eve and odd. Unique marker genes for each type of heart cell were visualized by violin plot and/or feature plot. The top ten differentially expressed genes in each heart cell cluster were used for heatmap graphics.

RNA from sorted cells (obtained from female and male embryonic flies) was isolated using the TRIzol method (Invitrogen) as previously described (Huang et al., 2021; Zhao et al., 2012) and samples were transcribed into cDNA using All-In-One RT MasterMix (Applied Biological Materials). qRT-PCR was performed on three replicates by real-time analysis using Power SYBR Green PCR Master Mix (Applied Biosystems). We used the gene ribosomal protein L32 (rp49; also known as RpL32) as an internal control. Primer sequences have been provided in Table S7. The qRT-PCR products were quality checked by 1% agarose gel with dye SYBR Safe (Invitrogen).

At the desired developmental stage, embryos were collected and washed thoroughly with ice-cold 1× PBS and subsequently dechorionated in fresh 4-5% bleach (sodium hypochlorite solution; Sigma-Aldrich) for 3 min at room temperature. Dechorionated embryos were extensively rinsed with water, then mounted with VECTASHIELD antifade mounting medium (Vector Laboratories, H-1000) for visualization and image capture using an Apotome microscope with 10× objective (ZEISS). GFP was excited at 488 nm and images were collected at 500-550 nm. To visualize the embryo morphology, differential interference contrast (DIC) images were simultaneously collected. For each stage, 10-20 embryos (female and male) were imaged; representative images are shown in the figures.

For embryonic staining, 20- to 24-h-old eggs were collected on juice agar plates. The embryos were then processed for immunofluorescence staining as previously described (Han et al., 2002; Zhu et al., 2017b). Briefly, the embryos were chlorox dechorionated, fixed in heptane and 4% paraformaldehyde (Polysciences) in 0.1 M phosphate buffer (pH 7.4). The fixed specimens were permeabilized and blocked in 0.3% Triton X-100 in 1× PBS (PBST) with 2% bovine serum albumin for 2 h and incubated overnight with primary antibodies at 4°C. Anti-Tin (1:500) and anti-Eve (1:500) (both gifts from Dr Olson, University of Texas Southwestern Medical Center, USA), anti-Alk (1:1000) (gift from Dr Ruth H. Palmer, University of Gothenburg, Sweden) (Lorén et al., 2003), anti-Tg (1:500) (gift from Dr Shun Kawabata, Kyushu University, Japan) (Shibata et al., 2010), anti-Ct and anti-Svp (1:50) (Developmental Studies Hybridoma Bank) were used to detect cardioblasts and pericardial cells. Cy5- and Cy3-conjugated secondary antibodies (1:1000; Thermo Fisher Scientific) were used to recognize the primary antibodies. All embryonic fly hearts were imaged by confocal microscopy using a 20× Plan-Apochromat 0.8 N.A. air objective or 63× Plan-Apochromat 1.4 numerical aperture oil objective under Airyscan SR mode. (LSM900; ZEISS). For each stage, 10-20 embryos (female and male) were imaged; representative images are shown in the figures.

Digital image processing and analysis were performed with ImageJ software (Fiji, 1.49; NIH) (Schneider et al., 2012), and statistical analyses were performed by Prism 7 (GraphPad Software). Statistical data are expressed as mean±s.d., and an unpaired two-tailed Student's t-test was applied to determine the significance of the difference between two groups. A value of P<0.05 was considered statistically significant.

We thank Dr Christopher Lazarski for helping with cell sorting (Children's National Hospital, Washington, DC, USA), Dr Ruth H. Palmer (University of Gothenburg, Sweden) for providing the Alk antibody and Dr Shun Kawabata for providing the Tg antibody (Kyushu University, Japan). We thank the Developmental Studies Hybridoma Bank at the University of Iowa for the antibodies. We thank Dr Peng Zhang (University of Maryland School of Medicine) for initial data analysis for the single-cell RNA sequencing.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201936.reviewer-comments.pdf