Single-nucleus RNA-sequencing in pre-cellularization Drosophila melanogaster embryos.

Our current understanding of the regulation of gene expression in the early Drosophila melanogaster embryo comes from observations of a few genes at a time, as with in situ hybridizations, or observation of gene expression levels without regards to patterning, as with RNA-sequencing. Single-nucleus RNA-sequencing however, has the potential to provide new insights into the regulation of gene expression for many genes at once while simultaneously retaining information regarding the position of each nucleus prior to dissociation based on patterned gene expression. In order to establish the use of single-nucleus RNA sequencing in Drosophila embryos prior to cellularization, here we look at gene expression in control and insulator protein, dCTCF, maternal null embryos during zygotic genome activation at nuclear cycle 14. We find that early embryonic nuclei can be grouped into distinct clusters according to gene expression. From both virtual and published in situ hybridizations, we also find that these clusters correspond to spatial regions of the embryo. Lastly, we provide a resource of candidate differentially expressed genes that might show local changes in gene expression between control and maternal dCTCF null nuclei with no detectable differential expression in bulk. These results highlight the potential for single-nucleus RNA-sequencing to reveal new insights into the regulation of gene expression in the early Drosophila melanogaster embryo.

Introduction

Early animal development is largely driven by maternally-deposited RNAs and proteins. In Drosophila melanogaster, zygotic gene expression is detected as early as the 10th nuclear cycle; however, zygotic genome activation primarily occurs during the 14th nuclear cycle while maternal RNAs are degraded and cellularization begins [1, 2]. Much of the difficulty in understanding the regulation of early embryonic gene expression lies in the challenge to simultaneously capture expression level and patterning. Classic examples of patterned gene expression and regulation originate from in situ hybridizations [3-6], however the nature of in situ hybridizations does not allow for the study of many genes at once. In order to fully understand the regulation of gene expression across the genome, it is imperative that we establish new methods to examine changes in spatially-patterned genes.

Prior work from our lab demonstrated the use of RNA sequencing in patterning mutants following cryosectioning embryos across the anterior-posterior axis [7]. This work benefits from knowledge of the origin of each slice during analysis; however, with this method we cannot truly resolve from where the RNAs originated as many nuclei will contribute to expression within each slice. Recent work from Karaiskos, Whale, et al (2017) demonstrated the use of single-cell RNA-sequencing in the early Drosophila embryo and the ability to construct virtual in situ hybridizations from prior knowledge of patterned gene expression [8]. Others have used single-cell RNA-sequencing in dorsoventral mutant embryos and showed depletion of an entire subset of cells [9]. These studies demonstrate the potential for single-cell RNA-sequencing to answer questions relating to pattern and body axis formation in the early Drosophila embryo; however, whether single-cell RNA-sequencing is sensitive enough to detect subtle changes in gene expression in mutant embryos lacking major defects remains unclear.

To establish the use of single-nucleus RNA-sequencing in the early Drosophila embryo, we decided to examine gene expression in control, as well as maternal null dCTCF embryonic nuclei which are subsequently referred to as dCTCF. Insulator elements were first described for their enhancer-blocking activity [10-12], and have since been shown to affect genome and chromosome structure as well [13-18]. Interestingly, mammalian CTCF serves as the only insulator protein in mammals; however, Drosophila and other arthropods have evolved several insulator proteins [19, 20]. The redundancy of Drosophila insulator proteins allows us to understand the many functions of insulators without causing cell lethality. Intriguingly however, dCTCF is not actually required for embryonic viability [21]. Previous reports indicate that loss of individual Drosophila insulator proteins yields minimal changes in gene expression [19, 22–25], but others show that dCTCF is required for correct expression of certain genes observed by in situ hybridizations in embryos and larvae [26, 27]. The observed changes are slight however, which may explain why large-scale defects in transcription are not observed with RNA-sequencing in flies lacking dCTCF.

Using 10x Genomics, we assayed gene expression across over 8,000 nuclei from control and dCTCF embryos. Overall, the nuclei tend to cluster according to expression of spatially-patterned genes, indicating that the nuclei retain information regarding their position in the embryo prior to dissociation. This allows us to understand genome-wide expression in spatial regions of embryos prior to cellularization by sequencing, which was previously only possible by slicing embryos [7, 13]. As expected considering the viability of dCTCF embryos, we found fewer differentially expressed genes in bulk than in individual clusters. We also found several candidate patterned genes that may be differentially expressed in certain clusters but not in bulk. Our analyses are available in a reproducible and usable format (see Code Availability) allowing others to explore our data analysis as well as analyze other genes of interest not explored here. Altogether this work establishes the use of single-nucleus RNA-sequencing in the early Drosophila embryo to detect subtle changes in gene expression and encompasses a resource to explore candidate locally differentially expressed genes upon loss of maternal dCTCF.

Methods

CRISPR

Maternal dCTCF nulls were created by using CRISPR mutagenesis to insert a dsRed protein followed by two consecutive stop codons immediately upstream of the dCTCF open reading frame. The homologous replacement template plasmid was constructed using a pUC19 backbone and ~1 kb homology arms generated by PCR (5’ homology arm primers: CCACAAAGAAACGTTAGCTAGTTCC and TCCTATGGACAAATTGGATTTGTTTTGG, 3’ homology arm primers: CCAAGGAGGACAAAAAAGGACGAG and CGTGAGTGGCGCGTGATC). Repair template was coinjected into Cas9-expressing embryos (Rainbow Transgenic Flies, Camarillo, California), along with two guide RNAs (ATTTGTCCATAGGAATGCCA, TGTCCATAGGAATGCCAAGG) expressed from a U6:3 promoter on a modified version of the pCFD3 plasmid [28]). Resulting flies were crossed to flies containing chromosome 3 balancer chromosomes, and screened by genotyping PCR. Putative hits were further screened by PCR and sequencing of the entire locus using primers outside the homology arms (CATTAGAATTCAAGGGCCATCAG and CACTTGAAGGATGGCTCG). A successful insertion line was recombined with an FRT site on chromosome 3L at cytosite 80B1 (Bloomington stock # BL1997).

Fly husbandry

All stocks were fed standard Bloomington food from LabExpress and maintained at room temperature unless otherwise noted. We used the FLP-DFS (dominant female sterile) technique [29] to generate dCTCF embryos. First, we crossed virgin hsFLP, w*;; Gl*/TM3 females to w*;; ovo, FRT2A(mw)/TM3 males (Bloomington Drosophila Stock Center ID: 2139). From this cross, we selected hsFLP,w*;; ovo, FRT2A(mw)/TM3 males and crossed them to virgin CTCF*,FRT2A/TM3 females. Larvae from this cross were heat-shocked on days 4, 5, and 6 for at least two hours in a water bath at 37°C. Upon hatching, virgin hsFLP, w*/+; CTCF*, FRT2A(mw)/ovoD*, FRT2A(mw) females were placed into a small cage with their male siblings. Flies were fed every day with yeast paste (dry yeast pellets and water) spread onto apple juice agar plates. These crosses were conducted simultaneously with another insulator protein, and control embryos were collected from the ovo line used to generate those germline clones (Bloomington Drosophila stock Center ID: 2149).

Western blots

Flies laid on grape-agar plates for two hours and embryos were either aged two hours at room temperature or taken directly after collection. Embryos were dechorionated with bleach, rinsed, and frozen in aliquots of ~25 embryos at -80 C. Embryos were homogenized in 25 μl RIPA buffer (Sigma cat # R0278) supplemented with 1 mM DTT and protease inhibitors (Sigma cat # 4693116001) using a plastic pestle. After homogenization, samples were mixed with 25 μl 2x Laemmli buffer (Bio-Rad # 1610737EDU), boiled for 3 minutes, and spun at 21,000 x g for 1 minute. Samples were loaded onto Bio-Rad mini Protean TGX 4–20% gels (# 4561096) and run at 200V for 30 minutes. Protein was transferred at 350 mA for one hour to Immobilon PVDF membrane (Millipore-Sigma # IPVH00010). Blots were blocked for one hour in PBST (1x PBS with 0.1% Tween) with 5% nonfat milk, and stained with primary antibodies (courtesy of Maria Cristina Gambetta [27], 1:1000 in PBST with 3% BSA) for one hour. Blots were then washed 3 times for 3 minutes rotating in PBST and probed with an HRP-conjugated anti-Rabbit secondary antibody (Rockland Trueblot, # 18-8816-33, 1:1000 in PBST with 5% milk) for one hour. After extensive washing with PBST, blots were developed with Clarity ECL reagents (Bio-Rad # 1705060) and imaged. Validation of the loss of maternal dCTCF is shown in S1 Fig.

Nuclear isolation and sequencing

Nuclei were isolated from early to mid-nuclear cycle 14 embryos (stage 5) according to several previously published works [30-32]. First, the cages were cleared for 30 minutes to 1 hour to remove embryos retained by the mothers overnight, followed by a 2 hour collection and 2 hour aging. Then, the embryos were dechorionated in 100% bleach for 1 minute, or until most of the embryos were floating, with regular agitation by a paintbrush. The embryos were transferred to a collection basket made of a 50 mL conical and mesh. After the embryos were rinsed with water, the embryos were transferred into an eppendorf tube containing 0.5% PBS-Tween. From this point forward, samples were kept on ice to prevent further aging of embryos.

A minimum of 9 early to mid-nuclear cycle 14 embryos were sorted using an inverted compound light microscope and transferred to a 2 mL dounce containing 600 uL of lysis buffer (10 mM 10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 1% Bovine Serum Albumin, 1% RNase Inhibitor (Enzymatics, Part Num. Y9240L)) + 0.1% IGEPAL. The embryos were homogenized 20 times with a loose pestle and 10 times with a tight pestle. Pestles were rinsed with 100 uL lysis buffer + 0.1% IGEPAL after use. The resulting 800 uL of buffer and nuclei were transferred into an eppendorf tube, filtered with a 40 uM filter. Nuclei were pelleted by spinning for 5 minutes at 900 g and 4°C, washed in 500 uL lysis buffer (without 0.1% IGEPAL), and pelleted again before resuspending the nuclei in 20 uL lysis buffer (without 0.1% IGEPAL). Nuclei concentration was then adjusted to 1000 uL nuclei per uL, then nuclei were barcoded with the 10X Chromium Single Cell 3’ Gene Expression Kit (v3). Control and dCTCF nuclei were processed on separate days, then sequenced together with the Illumina NovaSeq (SP flow cell).

Data processing and analysis

We used kallisto-bustools [33] to generate a custom reference index and generate a nucleus x gene matrix. The data were analyzed in both Python and R, using primarily scVI via scvi-tools [34, 35], scanpy [36], and custom scripts for analysis.

Control and dCTCF nuclei were filtered separately as follows: (1) nuclei were ranked by the number of UMIs detected and nuclei ranked below the expected number of nuclei (10,000) were removed; (2) nuclei with fewer than 200 expressed genes were removed; (3) nuclei with greater than 5% mitochondrial expression were removed; (4) nuclei with greater than 50,000 UMI counts were removed; (5) genes expressed in fewer than 3 nuclei were removed.

Prior to batch correction, the data were subset to the 6000 most highly variable genes using scanpy’s dCTCF based on log1p normalized expression. We ran scVI with gene_likelihood = ‘nb’ to correct for batch effects.

The nuclei were clustered using the Leiden algorithm [37] within scanpy and visualized on a 2D UMAP [38]. Prior to batch correction, nuclei were clustered on log1p normalized gene expression. After batch correction, nuclei were clustered on the latent space derived from the scVI model. Marker genes representing each cluster were found using the sc.tl.rank_genes_groups function from scanpy with the Wilcoxon signed-rank test. In situ hybridizations of representative marker genes were obtained from the Berkeley Drosophila Genome Project[39-41]. Colors representing Leiden clusters were projected onto a virtual embryo using novoSpaRc [42, 43].

Log2 fold change and associated p-values were obtained for each gene using diffxpy (https://diffxpy.readthedocs.io/). Statistically significant differential expression was determined following Bonferroni correction of the p-values and filtered for adjusted p-value less than 0.05 and absolute value of log2 fold change greater than or equal to 1.5. Intersecting sets of differentially expressed genes were found and visualized with an UpSet plot [44, 45], following correction of adjusted p-values for the number of comparisons (multiplied by 11; 10 for the total number of clusters + 1 to include bulk differential expression).

Data and code availability

Raw sequencing data and.h5ad files are available on DataDryad: https://doi.org/10.6078/D13D9R. Much of our analysis originated from work by Booeshaghi and Pachter (2020) [46] and Chari et al (2021) [47], with the addition of custom scripts.

All of the code used in the analysis and in generating the figures is available here: https://github.com/aralbright/2021_AAMSME. Single-nucleus data pre-processing, batch correction and clustering, virtual in situ hybridization, and differential expression analyses are available in this GitHub repository as Google Colab notebooks. These notebooks are available for anyone to run from a web browser with the option to enter any genes of interest not discussed in this manuscript.

Results

To establish the use of single-nucleus RNA-sequencing for examining gene expression in the early Drosophila embryo prior to cellularization, we hand-sorted 10 to 20 early to mid-nuclear cycle 14 control and dCTCF embryos and isolated nuclei for single-nucleus RNA-sequencing using 10x Genomics 3’ Gene Expression. After filtering the data for high quality nuclei and correcting for non-biological variability (S2–S4 Figs), we used Leiden clustering [37] to detect distinct groups of nuclei from control and dCTCF embryos, altogether resulting in 8,400 nuclei across 10 clusters (Fig 1A) composed of both control and dCTCF nuclei (Fig 1B). We also removed yolk and pole cell nuclei as these nuclei are not informative for patterned gene expression; however, the fact that subsets of nuclei clustered on marker gene expression for yolk or pole cell nuclei provided us with confidence that our data accurately represent single-nucleus expression (S5 Fig). After removing groups of nuclei as indicated, the nuclei no longer cluster according to expression of yolk or pole cell markers, indicating that our data are of high quality (S6 Fig). Once we finalized the dataset, we then asked whether gene expression in the clusters determined by the Leiden algorithm are truly distinct.

Fig 1

Single-nucleus RNA-sequencing analysis of pre-cellularization Drosophila melanogaster embryos (A) Two-dimensional UMAP embedding of nuclei shows 10 transcriptionally distinct clusters, as determined using the Leiden algorithm. Associated colors are maintained throughout the manuscript. (B) Two-dimensional UMAP embedding of nuclei labeled by condition shows the overlap of control and dCTCF nuclei in a reduced dimensional space following correction for non-biological variability. Associated colors are maintained throughout the manuscript. (C) Heatmap of scaled gene expression for top four marker genes of each cluster, clusters are hierarchically ordered.

Given how well characterized patterned gene expression is in the early Drosophila embryo and that we found several distinct clusters of nuclei, we suspected that the clusters may represent different spatial regions within the embryo. Expression of the top marker genes representing each cluster is certainly distinct, and we noticed that many of these genes are expressed in patterns, namely fkh, tll, and htl (Fig 1C). To determine if single-nucleus RNA-sequencing in the early embryo can be spatially resolved, we examined in situ hybridizations of top marker genes for each cluster. We found that representative virtual and published in situ hybridizations of the top 20 marker genes (S1 Table) correspond to specific spatial regions within the embryo for clusters 0–7 (Fig 2A–2H and 2A’–2H’). We also found by projecting our nuclei onto a virtual embryo that the identities we assigned to each cluster correspond to the spatial identities of these clusters (Fig 2I).

Fig 2

Leiden clusters correspond to spatial regions within the embryo (A-H) Representative virtual in situ hybridizations (top, A-H) for top marker genes representing each cluster as labeled and the corresponding published in situ hybridizations (bottom, A’-H’) from the Berkeley Drosophila Genome Project [39–41]. (I) Projection of nuclei onto a virtual embryo labeled by the Leiden cluster as colored in Fig 1A. Virtual in situ hybridizations and projection of clusters onto a virtual embryo were generated using novoSpaRc [42, 43].

The anterior, posterior, and ventral clusters are the most defined based on a projection of our nuclei onto a virtual embryo, while clusters that represent the middle of the embryo in general had less well defined borders (Fig 2I). Even so, our data shows that single-nucleus RNA-sequencing in the early Drosophila embryo yields information related to the spatial position of nuclei prior to dissociation. We should note however, that the virtual in situ hybridizations shown above, as well as additional virtual in situ hybridizations that represent each cluster, contain genes present in the list of reference genes used to generate the virtual patterns (see Ilp4, htl, fkh in Fig 2, and Antp, NetA, disco in S7 Fig). As such, we considered that the virtual in situ hybridizations may be biased in those cases; however, we believe the presence of several other genes representing each cluster with similar patterning validated with both virtual and published in situ hybridizations indicates that reference bias is not an issue. In the end, we were unable to determine a spatial identity for clusters 8 and 9; however, we decided to include these clusters in subsequent analyses because the nuclei passed our quality control filters. Interestingly, cluster 9 appears to be absent in dCTCF embryos (Fig 1A and 1B). Without knowing the identity of cluster 9, we can only speculate why this may be the case; however, this raises the possibility that dCTCF may play a role in nuclear fate.

In an effort to establish the use of single-nucleus RNA-sequencing to detect local changes in gene expression in embryos prior to cellularization, we then asked whether we could detect potential differential expression of genes in individual clusters, but not in bulk. In most clusters and in bulk, gene expression appears to be up-regulated upon loss of maternal dCTCF (S8 Fig). We also found that differentially expressed genes shared between all clusters and in bulk represent one of the largest shared sets (Fig 3A, left most black bar). However, a substantial number of candidate differentially expressed genes appear differentially expressed in single clusters (Fig 3A, blue bars). Many other genes appear differentially expressed in groups of clusters, but not in bulk. Because we found many candidate differentially expressed genes, we considered that this may be due to low expression given the sparsity of single-nucleus RNA-sequencing; however, we found that the mean expression of candidate differentially expressed genes in single or multiple clusters overall does not have a substantially different pattern from that of non-differentially expressed genes (S9 Fig). Each of these curves are right-skewed, or most genes are expressed in low levels at less than 100 transcripts per million (TPM). Altogether, these results show that single-nucleus RNA-sequencing in the early embryo can be used to detect candidate differentially expressed genes that would not appear in bulk-sequencing data.

Fig 3

Differential expression of genes detected in one or more clusters, but not in bulk (A) UpSet plot for visualizing the top 40 shared sets of candidate differentially expressed genes between control and dCTCF nuclei within each cluster and in bulk. Horizontal bar plot (A, left) represents the total number of candidate differentially expressed genes within the cluster of the corresponding row. The vertical bar plot (A, top) represents the number of shared candidate differentially expressed genes for the conditions indicated below and is sorted from largest to smallest intersecting set, with each count representing a unique gene. Connected dots (black) represent the corresponding group of genes in the vertical bar plot above that might be differentially expressed in the clusters represented by rows with a filled in circle. Candidate genes differentially expressed in a single cluster are represented in blue. (B-D) log(scvi normalized expression) of (B) stumps, (C) bowl, (D) Esp in each cluster (left) and bulk (right) for control (teal) and dCTCF nuclei (pink). Asterisks indicate statistically significant differential expression (absolute value of expression > = 1.5 and Bonferroni corrected p-value < 0.05).

Upon the loss of an early developmental factor like dCTCF, we expect to observe potential differential expression of patterned genes in specific clusters with single-nucleus RNA-sequencing. Interestingly, stumps, a ventrally-expressed gene is differentially expressed in one ventral cluster, but not the other (Fig 3B). bowl, a gap gene primarily expressed in the anterior, appears to be up-regulated in the posterior-biased and one of the ventral clusters (Fig 3C). Finally, Esp, a posterior-striped gene is differentially expressed in several clusters. Intriguingly, we did not detect differential expression in bulk for stumps, bowl, or Esp. Because dCTCF is required for proper expression of Abd-B [ a Drosophila Hox gene, we also examined expression of several Drosophila Hox genes within each cluster and in bulk. Upon the loss of maternal dCTCF, Antp and abd-A are potentially differentially expressed in certain clusters, with potential differential expression of Antp in bulk data (S10 Fig). However despite the requirement of dCTCF for proper expression of Abd-B shown by in situ hybridization [26], we found no evidence of differential expression of Abd-B, in agreement with bulk RNA sequencing in larval CNS dCTCF mutants [27].

We cannot be certain whether or not specific genes are truly differentially expressed spatially without further investigation; however, our results demonstrate the use of single-nucleus RNA-sequencing to detect possible local changes in gene expression upon perturbation in the early embryo prior to cellularization.

Discussion

We conducted the above analyses in order to determine whether we could use single-nucleus RNA-sequencing as a means of understanding the regulation of gene expression in the early Drosophila embryo. First, we show that nuclei can be grouped into clusters represented by distinct gene expression. Then, we show that representative marker genes from the majority of the clusters recapitulate known patterns of expression. Importantly, we also present examples of potential differential expression of patterning genes in individual clusters upon loss of maternal dCTCF, but not in bulk.

Prior to this work, studies towards our understanding of the regulation of patterned gene expression in a spatial context included cytoplasmic RNAs in measures of expression. We must acknowledge the caveat that we do not know the extent to which maternal RNAs enter the nucleus and some of our results may reflect the presence of both maternal and zygotic RNAs. Nonetheless, we believe that single-nucleus RNA-sequencing is better suited as opposed to bulk RNA-sequencing to understand changes in gene expression in pre-cellularization embryos upon the loss of important developmental factors because of the ability to resolve local changes in expression. Supporting this notion, single-cell RNA-sequencing has already shown to resolve the loss of an entire cell fate in cellularized dorsoventral mutant embryos [9].

Whether or not the changes in gene expression that we observed have implications in embryonic development related to the loss of dCTCF is unclear without further investigation, such as single-molecule FISH to validate the observed changes in gene expression of particular RNAs. Ultimately, using single-nucleus RNA-sequencing to examine changes in gene expression upon the loss of important developmental factors has the potential to uncover perturbation responses previously undetected by bulk RNA-sequencing.

Top 20 marker genes representing each cluster as determined by sc.tl.rank_genes_groups.

(CSV)

Click here for additional data file.

(A) Western blotting of OreR, 0h and 2h dCTCF embryos using an antibody to Cp190, another insulator protein, as a control. (B) Western blotting of OreR, 0h and 2h dCTCF embryos using an antibody to dCTCF. The 2h embryos were aged for an additional 2 hours with the majority of the embryos representing nuclear cycle 14, the same time point at which we conducted single-nucleus RNA-sequencing. A cross-reactive band appears at approximately 75 kd, with the dCTCF band appearing at approximately 130 kd.

(TIF)

Click here for additional data file.

(A) Knee plot for barcodes ranked by the number of UMIs versus UMI counts for control (left) and dCTCF (right) experiments. Black line indicates the position of the 10,000th (expected number of cells) on each axis. (B) Percent mitochondrial expression per nucleus in control (left) and dCTCF(right) nuclei. Dashed line represents 5% mitochondrial expression, or the cutoff used for filtering the data. (C) Number of genes detected per nucleus by UMI counts in control (left) and dCTCF(right) nuclei.

(TIF)

Click here for additional data file.

Number of genes detected (left), UMI counts (middle), percent mitochondrial expression (right) per nucleus after filtering in (A) control and (B) dCTCF experiments.

(TIF)

Click here for additional data file.

Two-dimensional UMAP embedding of control (teal) and dCTCFmat-/- (pink) nuclei (A) before, (B) after batch correction using scVI, and (C) after removing low quality nuclei.

(TIF)

Click here for additional data file.

Two-dimensional UMAP embedding of nuclei before additional filtering colored by (A) number of genes detected, (B) UMI counts, (B) percent mitochondrial expression. (D-F) log(scvi normalized expression) of three genes with representative in situ hybridizations below for (D) cell cycle gene aurB, (E) yolk nucleus marker, sisA (F) and pole cell marker pgc.

(TIF)

Click here for additional data file.

Two-dimensional UMAP embedding of nuclei after removal of clusters with high percent mitochondrial expression, aurB expression, sisA expression, and pgc expression colored by (A) number of genes detected, (B) UMI counts, (C) percent mitochondrial expression. (D-F) log(scvi normalized expression) of three genes with representative in situ hybridizations below for (D) cell cycle gene aurB, (E) yolk nucleus marker, sisA (F) and pole cell marker pgc.

(TIF)

Click here for additional data file.

(A-P) Representative virtual (top, A-P) and Berkeley Drosophila Genome Project (bottom, A’-P’) in situ hybridizations for additional marker gene expression within each cluster as indicated. This supplemental figure accompanies Fig 2 in the main text.

(TIF)

Click here for additional data file.

Volcano plots of log2FC (log2(fold-change)) by the log of the adjusted p-value (p-adj) for differential expression calculated in bulk (top middle) and in individual clusters as indicated. Colored dots indicate genes with significant differential expression, an absolute value of log2FC > = 1.5 and p-adj < 0.05. Significantly down-regulated genes are indicated in green, significantly up-regulated genes in pink, and non-significantly differentially expressed genes in light gray.

(TIF)

Click here for additional data file.

Histogram of average gene expression in transcripts per million (TPM) of differentially expressed genes in one cluster (yellow), differentially expressed in multiple clusters and/or in bulk (blue), and non-differentially expressed genes (red). Each count on the y-axis represents a single gene.

(TIF)

Click here for additional data file.

Plots of log(scvi normalized expression) of select (A) Antennapedia complex genes: Dfd (top), Scr (middle), and Antp (bottom) and (B) Bithorax complex genes: abd-A (top) and Abd-B (bottom) in each cluster (left) and in bulk (right) for control (teal) and dCTCF nuclei (pink). Asterisks indicate statistically significant differential expression (absolute value of expression > = 1.5 and Bonferroni corrected p-value < 0.05).

(TIF)

Click here for additional data file.

16 Feb 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Anton Wutz

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for stating the following in the Competing Interests section:

"I have read the journal's policy and the authors of this manuscript have the following competing interests: MBE is a founder and former member of the board of directors of PLOS."

Please confirm that this does not alter your adherence to all PLOS ONE policies on sharing data and materials, by including the following statement: "This does not alter our adherence to  PLOS ONE policies on sharing data and materials.” (as detailed online in our guide for authors http://journals.plos.org/plosone/s/competing-interests).  If there are restrictions on sharing of data and/or materials, please state these. Please note that we cannot proceed with consideration of your article until this information has been declared.

Please include your updated Competing Interests statement in your cover letter; we will change the online submission form on your behalf.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: No

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: No

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This study presents results of single-nucleus RNA-sequencing (snRNA-seq) in early fly embryos undergoing zygotic genome activation. Wildtype (WT) and CTCF[mat-] mutants lacking maternal CTCF were both analyzed. Transcripts detected in single WT nuclei could be mapped onto a virtual reference embryo using known marker genes, and recapitulated known spatial expression patterns similarly to single-cell RNA-seq (Karaiskos et al. 2017). Differential gene expression between WT and CTCF[mat-] embryos was analyzed. Measuring transcript abundance differences between individual snRNA-seq clusters in WT and mutant embryos identified more differentially expressed genes than measuring transcript abundance differences between all cells of WT and mutant embryos in bulk.

These results are interesting because they show that snRNA-seq is sensitive enough to detect relatively subtle gene misexpression defects in mutant embryos lacking major defects in cell fate decisions. The data appears to be of high quality. But I have 2 major confusions about the study’s design that must be addressed prior to publication (see major comments).

Major comments:

1. Please clearly describe whether the CTCF[mat-] embryos generated in this study zygotically express CTCF, and discuss whether this confounds analyses of differential gene expression in CTCF[mat-] embryos if these embryos already initiated zygotic transcription.

2. I am confused by line 59. How can snRNA-seq help understand how gene expression is established prior to zygotic genome activation? If the zygotic genome is not transcribed, what would snRNA-seq detect?

3. Fig. S1 is missing.

Minor comments

4. Line 40 is not well written: “Much of the difficulty in understanding the regulation of early embryonic gene expression lies in our ability [in the challenge?] to simultaneously capture expression level and patterning”.

5. Line 74-75: Kaushal et al. 2021 Nat Comm report differentially expressed genes in CTCF mutants lacking maternal and zygotic CTCF relative to WT by RNA-sequencing. This seems contradictory with the statement that differential gene expression in CTCF mutants has “not been found via sequencing”.

6. Lines 80-81: “Differential expression in spatial regions by sequencing” was previously performed by the references cited in the introduction, and therefore stating that this “was previously only possible by mechanical manipulation” should be toned down.

7. Line 169 must be amended (“nuclei ranked below the to the expected number of nuclei”).

8. Line 206: Figure references have typos.

9. Lines 218 and 220: First words of sentences should be capitalized.

10. I did not understand the meaning of the sentence in lines 234-236.

11. Cluster 9 should be discussed further, even if its spatial identity could not be determined. Why is this cluster only detected in WT but not CTCF[mat-] embryos?

12. Fig. S4 panel C is not described in the figure legend.

Reviewer #2: The authors Albright et al classified embryonic nuclei by single-nucleus RNA-seq and examined CTCF-regulated gene expression in these nuclei by comparing wild type to CTCF maternal null mutants during zygotic genome activation. They identified more cluster-specific differential expression than in bulk RNA-seq and highlighted several examples of differential expression of spatial marker genes in specific clusters. The work should be of general interest. However, the data requires further analyses, and the conclusions were not clearly presented or fully supported. Major revisions are needed to both the analyses and writing.

1. Figure numbering is incorrect for all supplementary figures, and it is not possible to understand which figures the authors were calling for. The intended Figure S1 is missing. There is no Figure S9 in the submission.

2. The authors should indicate whether there is zygotic CTCF expression in this mutant. A diagram will be very helpful. The missing Figure S1 makes it more challenging to understand the properties of this KO.

3. Although the CTCF maternal knockout is known to be viable, dysregulation of Hox gene expression has been reported in embryos. The authors should characterize their new mutants and compare with previous data (Gambetta and Furlong, 2018) to report consistent or distinct phenotypes, e.g., viability and expression pattern of Hox genes.

4. Figure 1B indicates loss of cluster 9 in CTCF maternal KO. Is this because the cells are absent or because their gene expression changes and they are classified into other clusters? This can be determined by in situ of cluster 9-specific genes.

5. What is the accuracy of spatial prediction based on the RNA-seq? How many top marker genes were checked, and how many have consistent expression patterns with the in situ data?

6. The authors should verify their knockout by western blots, which are mentioned but not presented.

7. The authors stated that spatial identities cannot be assigned to clusters 8 and 9, but some quick searches with gene IDs in Figure 1C identified embryonic CNS for cluster 8 (maybe 9) and another clear distinct spatial pattern for cluster 9. The authors need to dig deeper into the spatial identity by searching more genes in those clusters using publicly available data.

http://www.flyexpress.net/search/genes/CG2233/images/BDGP/LDVO

http://www.flyexpress.net/search/genes/CG2225/images/BDGP/LDVO

http://www.flyexpress.net/search/genes/CG3408/images/BDGP/LDVO

8. The authors need to indicate how much overlap there is between differential expression of different clusters. For example, cluster 9 has ~170 DE, cluster 3 has ~80 DE, but the common DE between clusters 3 and 9 is ~20 genes. Does this mean that most DE genes are cluster-specific? The authors need to quantify and present the data clearly. It is extremely hard to decipher the information presented in current Figure 3A.

9. The authors concluded that they identified more cluster-specific DE than in bulk but did not present data beyond the information in Figure 3A, horizontal bars. Do these bars include primarily the same genes, or they are totally different genes? It is very unclear how many more genes are detected as cluster-specific than in bulk. The authors state on page 12, “Many other genes are also differentially expressed in groups of clusters, but not in bulk.” This is very important and quantifiable information, and the authors need to present the data clearly.

10. The authors need to verify their RNA-seq data with another type of assay, e.g., RNA FISH. The differential expression (e.g., Figure 3B and C) should be verified between wildtype vs. CTCF KO cells, and between wildtype clusters of cells. Without such verification, it is hard to conclude that the single nucleus RNA-seq provides useful information. The verification should be done at least for the top differences.

11. How do the authors explain that most DE is up-regulation? Does this agree with Kaushal et al., 2021? It seems likely that decreases in DE are still masked by the inability to further separate cell types using this approach.

12. Lines 269-272: “Because we found many differentially expressed genes, we considered that this may be due to low expression given the sparsity of single-nucleus RNA-sequencing; however, we found that the mean expression of differentially expressed genes in single or multiple clusters overall does not have a substantially different pattern from that of non-differentially expressed genes.” The authors should compare and present mean expression of differential genes versus non-differential genes. This analysis is essential to rule out the possibility that more differential genes in cluster data than bulk data result from inaccurate quantification due to insufficient sequencing and coverage. Figure S9 may contain such information but is currently missing. The authors should also clearly define what they mean by “a substantially different pattern”.

13. Figure 3B-D: the expression change of patterning markers may lead to morphogenesis defects – did the authors examine the tissue morphology and distribution of marker gene expression by in situ or RNA FISH in embryos/larvae? Is the difference caused by strong differences in a small group of cells or weak differences in all cells in one cluster?

14. What are the most affected factors and signaling pathways in CTCF KO? How many of these genes are CTCF binding targets based on published CTCF ChIP-seq?

15. Please finish the sentence in lines 141, 169 (“to” and “the expected number of *UMI*”?).

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

17 May 2022

In the spirit of PloS we wanted to provide an open access resource from a rigorous experiment for others to use in future work. Regarding the availability of our data and analysis, we realized after receiving this review that some of the links in our Google Colab notebooks were broken. We fixed this issue prior to the submission of revisions.

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This study presents results of single-nucleus RNA-sequencing (snRNA-seq) in early fly embryos undergoing zygotic genome activation. Wildtype (WT) and CTCF[mat-] mutants lacking maternal CTCF were both analyzed. Transcripts detected in single WT nuclei could be mapped onto a virtual reference embryo using known marker genes, and recapitulated known spatial expression patterns similarly to single-cell RNA-seq (Karaiskos et al. 2017). Differential gene expression between WT and CTCF[mat-] embryos was analyzed. Measuring transcript abundance differences between individual snRNA-seq clusters in WT and mutant embryos identified more differentially expressed genes than measuring transcript abundance differences between all cells of WT and mutant embryos in bulk.

These results are interesting because they show that snRNA-seq is sensitive enough to detect relatively subtle gene misexpression defects in mutant embryos lacking major defects in cell fate decisions. The data appears to be of high quality. But I have 2 major confusions about the study’s design that must be addressed prior to publication (see major comments).

Major comments:

1. Please clearly describe whether the CTCF[mat-] embryos generated in this study zygotically express CTCF, and discuss whether this confounds analyses of differential gene expression in CTCF[mat-] embryos if these embryos already initiated zygotic transcription.

We have addressed this comment by adding the following statement at the end of the first paragraph in the ‘Results’ section: “We confirmed that our dCTCFmat-/- embryos lack maternal dCTCF at 0h and 2h after laying via western blot (S1 Fig). These experiments were conducted agnostic of the zygotic genotype given that zygotic genome activation largely does not occur until nuclear cycle 14 [1,2].”

2. I am confused by line 59. How can snRNA-seq help understand how gene expression is established prior to zygotic genome activation? If the zygotic genome is not transcribed, what would snRNA-seq detect?

To address this comment, we removed the following text: ‘zygotic genome activation.’ While it’s true that large-scale zygotic genome activation does not occur until nuclear cycle 14 (with a subset of genes expressed before then)… the true need for snRNA-seq lies in the fact that the embryo is not composed of cells at this point.

3. Fig. S1 is missing.

We thank you for pointing out this mistake and have corrected the figure submissions.

Minor comments

4. Line 40 is not well written: “Much of the difficulty in understanding the regulation of early embryonic gene expression lies in our ability [in the challenge?] to simultaneously capture expression level and patterning”.

We replaced “in our ability” with “in the challenge” as suggested.

5. Line 74-75: Kaushal et al. 2021 Nat Comm report differentially expressed genes in CTCF mutants lacking maternal and zygotic CTCF relative to WT by RNA-sequencing. This seems contradictory with the statement that differential gene expression in CTCF mutants has “not been found via sequencing”.

Kaushal et al. 2021 found differential expression in 386 genes upon loss of CTCF, only 10% of which had a CTCF peak within 1 kb of the transcription start site. This enrichment is significant over a distribution of similar non-differentially expressed genes; however, we interpret this as CTCF binding having a small effect on gene expression directly. Kaushal et al 2021 also mention that they expected indirect transcriptional changes.

We changed the sentence to the following: “The observed changes are slight however, which may explain why large-scale defects in transcription are not observed with RNA-sequencing in flies lacking dCTCF.”

6. Lines 80-81: “Differential expression in spatial regions by sequencing” was previously performed by the references cited in the introduction, and therefore stating that this “was previously only possible by mechanical manipulation” should be toned down.

The references cited in the introduction conducted single-cell RNA-sequencing after cellularization, therefore we added the words, “of embryos prior to cellularization” as well as changed “mechanical manipulation” to “slicing embryos” to be more specific.

7. Line 169 must be amended (“nuclei ranked below the to the expected number of nuclei”).

8. Line 206: Figure references have typos.

9. Lines 218 and 220: First words of sentences should be capitalized.

We addressed points 7-9 by correcting the typos.

10. I did not understand the meaning of the sentence in lines 234-236.

The original intention of this sentence was to illustrate that we could not confidently assign specific tissue fates (i.e. neuroectoderm, mesoderm, etc) to the different clusters so we assigned the clusters more generic identities (i.e. anterior, posterior, etc) instead. In order to avoid confusion here, we ultimately decided to delete this sentence.

11. Cluster 9 should be discussed further, even if its spatial identity could not be determined. Why is this cluster only detected in WT but not CTCF[mat-] embryos?

We have speculated why cluster 9 is not present in CTCF[mat-] embryos; however we believe the true answer will require further investigation outside of the scope of this manuscript. We agree that it’s interesting to note at the very least, so we have added the following:

“Interestingly, cluster 9 appears to be absent in dCTCFmat-/- embryos (Fig 1A-B). Without knowing the identity of cluster 9, we can only speculate why this may be the case; however, this raises the possibility that dCTCF may play a role in nuclear fate.”

12. Fig. S4 panel C is not described in the figure legend.

We added “and (C) after removing low quality nuclei.”

Reviewer #2: The authors Albright et al classified embryonic nuclei by single-nucleus RNA-seq and examined CTCF-regulated gene expression in these nuclei by comparing wild type to CTCF maternal null mutants during zygotic genome activation. They identified more cluster-specific differential expression than in bulk RNA-seq and highlighted several examples of differential expression of spatial marker genes in specific clusters. The work should be of general interest. However, the data requires further analyses, and the conclusions were not clearly presented or fully supported. Major revisions are needed to both the analyses and writing.

1. Figure numbering is incorrect for all supplementary figures, and it is not possible to understand which figures the authors were calling for. The intended Figure S1 is missing. There is no Figure S9 in the submission.

We have corrected this in the figure submission.

2. The authors should indicate whether there is zygotic CTCF expression in this mutant. A diagram will be very helpful. The missing Figure S1 makes it more challenging to understand the properties of this KO.

We agree that the missing Figure S1 makes interpreting our mutant challenging. We have corrected that as well as further clarified the mutant in the text with the following:

”We confirmed that our dCTCFmat-/- embryos lack maternal dCTCF at 0h and 2h after laying (S1 Fig). These experiments were conducted agnostic of the zygotic genotype given that zygotic genome activation largely does not occur until nuclear cycle 14 [1,2].”

3. Although the CTCF maternal knockout is known to be viable, dysregulation of Hox gene expression has been reported in embryos. The authors should characterize their new mutants and compare with previous data (Gambetta and Furlong, 2018) to report consistent or distinct phenotypes, e.g., viability and expression pattern of Hox genes.

Gambetta and Furlong 2018 do show that dCTCF is required for proper Abd-B expression; however, Kaushal et al 2021 (of which Dr. Gambetta is the principal author) report that their results in bulk RNA sequencing embryonic CNS do not show that Abd-B is differentially expressed upon the loss of dCTCF. This comment raises an excellent point in that dCTCF is potentially implicated in the regulation of Hox gene expression or topologically associating domains surrounding the Hox genes, therefore we have added a new figure (S10Fig) depicting expression in a select few Hox genes.

We added the following text discussing this figure: “Because dCTCF is required for proper expression of Abd-B [26], a Drosophila Hox gene, we also examined expression of several Drosophila Hox genes within each cluster and in bulk. Upon the loss of maternal dCTCF, Antp and abd-A are differentially expressed in certain clusters, but not in bulk (S10 Fig). However, we found no evidence of differential expression of Abd-B, in agreement with bulk RNA sequencing in embryonic CNS dCTCF mutants [27].”

We limited the number of Hox genes in this figure out of space concerns; however, we found no evidence of differential expression in the other Hox genes (lab, pb, Ubx).

4. Figure 1B indicates loss of cluster 9 in CTCF maternal KO. Is this because the cells are absent or because their gene expression changes and they are classified into other clusters? This can be determined by in situ of cluster 9-specific genes.

Reviewer #1 had a similar comment above, and in response we added the following sentences:

“Interestingly, cluster 9 appears to be absent in dCTCFmat-/- embryos (Fig 1A-B). Without knowing the identity of cluster 9, we can only speculate why this may be the case; however, this raises the possibility that dCTCF may play a role in nuclear fate.”

We acknowledge that this is interesting and would like to follow up with in situ hybridization among other experiments relating to CTCF function in early development; however, we believe that the optimization and execution of many in situ hybridizations is outside of the scope of this manuscript. We believe the addition of the statements above are sufficient to point out this peculiarity without over-speculating on the cause.

5. What is the accuracy of spatial prediction based on the RNA-seq? How many top marker genes were checked, and how many have consistent expression patterns with the in situ data?

This is a very important question that we appreciate you raising. We show additional markers in S7 Fig. We did check the top 20 marker genes as shown in S1 Table, as well as a few genes in the top 50 that we recognized by name.

We do acknowledge that some of the marker genes in S1 Table do not exhibit patterned gene expression; however, we believe that this could be explained by a number of possibilities. For one example, small differences in expression of ubiquitously expressed genes can appear statistically significant enough for these genes to appear as markers. This particular example highlights the importance of double checking the expression patterns by in situ hybridization before calling spatial regions. As shown in the manuscript, we used existing in situ hybridizations from the Berkeley Drosophila Genome Project to confirm patterning.

As mentioned previously, some of the marker genes in S1 Table do not show patterned gene expression; however, clusters 0-7 contained multiple spatially patterned genes beyond what was shown in Fig 2 and S7 Fig. Without conducting additional experiments, such as single-molecule FISH in control and dCTCF[mat-] embryos, we cannot give an exact estimate of accuracy or confirm patterning defects. However, the fact that several genes show patterning according to the spatial regions we called each cluster provided us with enough confidence to determine where the nuclei originated from within the embryo.

6. The authors should verify their knockout by western blots, which are mentioned but not presented.

We corrected this with the addition of Figure S1.

7. The authors stated that spatial identities cannot be assigned to clusters 8 and 9, but some quick searches with gene IDs in Figure 1C identified embryonic CNS for cluster 8 (maybe 9) and another clear distinct spatial pattern for cluster 9. The authors need to dig deeper into the spatial identity by searching more genes in those clusters using publicly available data.

This reviewer raises an excellent point about potential ambiguity in the way that we refer to spatial patterning, particularly in clusters 8 and 9. When double checking the Berkeley Drosophila Genome Project in situ hybridization database however, we noticed that many genes in clusters 8 and 9 show no staining in Stages 4-6. Stage 5 corresponds to nuclear cycle 14, the time point at which we collected our embryos. Here are a few examples of cluster 8 and 9 marker genes:

CG2233 -

https://insitu.fruitfly.org/cgi-bin/ex/report.pl?ftype=3&ftext=GH20802-dg

CG2225 -

https://insitu.fruitfly.org/cgi-bin/ex/report.pl?ftype=1&ftext=FBgn0032957

CG3408 -

https://insitu.fruitfly.org/cgi-bin/ex/report.pl?ftype=1&ftext=FBgn0036008

The lack of staining in these images does not preclude the possibility that these genes are actually expressed in pre-cellularization nuclei at small levels or that these nuclei are fated to become the embryonic CNS; however, we cannot confidently identify these clusters without further investigation.

8. The authors need to indicate how much overlap there is between differential expression of different clusters. For example, cluster 9 has ~170 DE, cluster 3 has ~80 DE, but the common DE between clusters 3 and 9 is ~20 genes. Does this mean that most DE genes are cluster-specific? The authors need to quantify and present the data clearly. It is extremely hard to decipher the information presented in current Figure 3A.

9. The authors concluded that they identified more cluster-specific DE than in bulk but did not present data beyond the information in Figure 3A, horizontal bars. Do these bars include primarily the same genes, or they are totally different genes? It is very unclear how many more genes are detected as cluster-specific than in bulk. The authors state on page 12, “Many other genes are also differentially expressed in groups of clusters, but not in bulk.” This is very important and quantifiable information, and the authors need to present the data clearly.

We agree that we need to clarify the UpSet plot in Figure 3A, and hope that our edits (relevant changes italicized/underlined below for reference) to the figure legend satisfy the reviewer in resolving comments 8 and 9.

“(A) UpSet plot for visualizing the top 40 shared sets of candidate differentially expressed genes between control and dCTCFmat-/- nuclei within each cluster and in bulk. Horizontal bar plot (A, left) represents the total number of candidate differentially expressed genes within the cluster of the corresponding row. The vertical bar plot (A, top) represents the number of shared candidate differentially expressed genes for the conditions indicated below and is sorted from largest to smallest intersecting set, with each count representing a unique gene. Connected dots (black) represent the corresponding group of genes in the vertical bar plot above that might be differentially expressed in the clusters represented by rows with a filled in circle. Candidate genes differentially expressed in a single cluster are represented in blue.”

The horizontal bars do not count unique genes as they represent the number of differentially expressed genes in a single cluster, a single gene can be differentially expressed in any number of clusters. However, the vertical bar plot counts are of unique genes because each gene can only belong to one set (similar to how one would consider counts in a Venn diagram, UpSet plots being a more recent solution to data with many intersecting sets).

10. The authors need to verify their RNA-seq data with another type of assay, e.g., RNA FISH. The differential expression (e.g., Figure 3B and C) should be verified between wildtype vs. CTCF KO cells, and between wildtype clusters of cells. Without such verification, it is hard to conclude that the single nucleus RNA-seq provides useful information. The verification should be done at least for the top differences.

We completely agree that RNA FISH will be necessary to verify differential expression; however, we have edited our intention for our manuscript in the abstract, “In order to establish the use of single-nucleus RNA sequencing in Drosophila embryos prior to cellularization, here we look at gene expression in control and insulator protein, dCTCF, maternal null embryos during zygotic genome activation at nuclear cycle 14.”

While we acknowledge that CTCF is an interesting case of biology, from the beginning this manuscript was intended to establish the use of single-nucleus RNA sequencing in pre-cellularization embryos. Rather than address the many questions surrounding the role of CTCF in early development, we hope that our work serves as a resource for people to follow-up on.

At the end of our manuscript, we also state: “Whether or not the changes in gene expression that we observed have implications in embryonic development related to the loss of dCTCF is unclear without further investigation, such as single-molecule FISH to validate the observed changes in gene expression of particular RNAs. Ultimately, using single-nucleus RNA-sequencing to examine changes in gene expression upon the loss of important developmental factors has the potential to uncover perturbation responses previously undetected by bulk RNA-sequencing.”

We believe that our statements in the introduction and discussion, with additional changes to language throughout the manuscript deemphasizing true differential expression, are sufficient enough to address the reviewer’s concerns that we do not have FISH data in this manuscript as we highlight this work more for its potential to aid future studies.

11. How do the authors explain that most DE is up-regulation? Does this agree with Kaushal et al., 2021? It seems likely that decreases in DE are still masked by the inability to further separate cell types using this approach.

We believe that this is an interesting point to note; however, we can only speculate as to why most differentially expressed genes appear to be up-regulated and decided not to comment on this in the manuscript. Kaushal et al. 2021 did report almost twice as many up-regulated genes compared to down-regulated genes in larval CNS (Figure 3a in their manuscript).

We would argue that spatial differential expression is unmasked in our ability to cluster nuclei into distinct groups despite their belonging to a syncytium where no cells are present. Because we have prior knowledge of spatial expression in embryonic nuclei, we were able to examine the spatial identities of clusters over many iterations of our pipeline and are confident in our data. If we were to over-cluster, or further separate our nuclei into more clusters, we run the risk of highlighting non-biological variation as something biological.

Global effects on gene expression are apparent in our pseudo-bulking results from individual nuclei, as well as in Kaushal et al. 2021 in bulk extractions and sequencing of larval CNS RNA. We believe that our results in Figure 3 are sufficient to show that we have unmasked local changes in differential expression. We acknowledge that the examples we chose for Figure 3B-D are all up-regulated, but maintain that this is valid regardless of whether genes are up or down-regulated.

12. Lines 269-272: “Because we found many differentially expressed genes, we considered that this may be due to low expression given the sparsity of single-nucleus RNA-sequencing; however, we found that the mean expression of differentially expressed genes in single or multiple clusters overall does not have a substantially different pattern from that of non-differentially expressed genes.” The authors should compare and present mean expression of differential genes versus non-differential genes. This analysis is essential to rule out the possibility that more differential genes in cluster data than bulk data result from inaccurate quantification due to insufficient sequencing and coverage. Figure S9 may contain such information but is currently missing. The authors should also clearly define what they mean by “a substantially different pattern”.

The addition of missing S1 Fig corrected the misnumbering of our figure. S9 Fig shows a comparison of mean expression in non-differentially expressed genes (not DE), genes differentially expressed in multiple clusters (multi DE), and genes differentially expressed in a single cluster (single DE).

To clarify what we mean by ‘substantially different pattern,’ we added: “Each of these curves are right-skewed, or most genes are expressed in low levels at less than 100 transcripts per million (TPM).”

13. Figure 3B-D: the expression change of patterning markers may lead to morphogenesis defects – did the authors examine the tissue morphology and distribution of marker gene expression by in situ or RNA FISH in embryos/larvae? Is the difference caused by strong differences in a small group of cells or weak differences in all cells in one cluster?

14. What are the most affected factors and signaling pathways in CTCF KO? How many of these genes are CTCF binding targets based on published CTCF ChIP-seq?

We thank you for comments 13 and 14 which allude to future studies on CTCF function in early embryonic development. Kaushal et al. 2021 does contain examples of local changes in expression upon loss of CTCF, but later than our study during embryonic development. Gambetta and Furlong 2018 show a slight change in the domain of Abd-B expression in addition to morphological changes, also later in development.

When hand-sorting embryos prior to isolating nuclei for these experiments, we did not observe any morphological defects in the pre-cellularization embryos. While optimizing our experiments we did keep old plates containing CTCF[mat-] embryos, and concur that maternal CTCF is not required for embryogenesis. Beyond this observation, we did not conduct careful morphological analyses of our CTCF[mat-] mutant.

Our main objective for this study is to highlight the potential for single-nucleus RNA-sequencing to complement studies in patterning prior to cellularization. We agree that understanding the biology of CTCF, as well as other insulators, is important and inherently interesting; however, outside of the scope of this manuscript. We decided to use a dCTCF mutant in this manuscript due to general interest in the protein, but do not intend for this manuscript to answer all of our questions on its function in early Drosophila embryonic development. We hope that our data will be useful to those conducting more careful studies into CTCF function and effects on patterning in the early embryo.

15. Please finish the sentence in lines 141, 169 (“to” and “the expected number of *UMI*”?).

We have corrected these mistakes in the updated manuscript.

Submitted filename: 20220517_reviewer_response.docx

Click here for additional data file.

30 May 2022

7 Jun 2022

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: (No Response)

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

6. Review Comments to the Author

1. Lines 212-213: Related to my original points 1+2, this sentence could cause confusion if the reader interprets this as meaning that snRNA-seq was performed prior to zygotic genome activation. In this study, snRNA-seq was performed on embryos undergoing zygotic genome activation; therefore, dCTCF[mat-/-] embryos probably transcribe CTCF mRNA (this could be verified in the snRNA-seq data) but do not translate (at least high levels of) CTCF protein yet. This could be explained more clearly.

We agree that we could better explain the timing of our collection and the assumptions we were working under. We added language to the figure legend for S1 Fig, “The 2h embryos were aged for an additional 2 hours with the majority of the embryos representing nuclear cycle 14, the same time point at which we conducted single-nucleus RNA-sequencing,” in order to clarify that we did not detect CTCF protein at this time point.

As you have alluded to, we are working under the assumption that any zygotically expressed CTCF would not have time to be translated, folded, and function. For this reason, we never intended to genotype the embryos for zygotic dCTCF. In the figures pasted below (see attached document), we do see low CTCF expression (average 0.81 TPM); however, we only detected expression in 270 control and 467 nuclei (less than 10% of our total nuclei). Although expression appears higher in the dCTCF[mat-/-] embryos and in individual clusters by eye, these differences are not statistically significant. We decided not to include these panels in S1 Fig as the Western blot definitively shows the absence of dCTCF protein in the early embryo.

As far as what we originally had in 212-213, we decided to remove those sentences entirely in order to avoid the confusion. We believe that the Western blot is enough to show that dCTCF protein will have no effect on our results, because we did not detect it, without speculating on the effect of zygotic gene expression.

2. Lines 56-60: Related to my original point 6, I feel that the authors should tone down their argument that previous scRNA-seq studies may not have primarily measured zygotic expression due to the presence of maternal cytoplasmic RNAs in cells. The present study, Karaiskos et al. 2017, and Ing-Simmons et al. 2021 likely all handled nuclei instead of cells because embryos were dounce-homogenized in all three studies – resulting in nuclei, not intact cells. In my view, the information obtained by snRNA-seq and scRNA-seq in fly embryos is thus comparable (if the authors don’t agree, please explain why). The novelty of the current study is, in my view, rather the demonstration that snRNA-seq is sensitive enough to detect relatively subtle gene misexpression defects in mutant embryos lacking major defects in cell fate decisions.

This is a great point, we were relying on the assumption that the other works isolated cells rather than nuclei, but realize that this may not necessarily be the case. In order to tone down the language in our manuscript, we removed the statements about zygotic gene expression in lines 56-60 and replaced them with the underlined, “These studies demonstrate the potential for single-cell RNA-sequencing to answer questions relating to pattern and body axis formation in the early Drosophila embryo; however, whether single-cell RNA-sequencing is sensitive enough to detect subtle changes in gene expression in mutant embryos lacking major defects remains unclear.”

We believe this change still highlights that our work was conducted earlier in development, as well as detects subtle changes in expression as you mentioned, but removes the unclear language on zygotic gene expression.

3. Line 317: “Embryonic” should be changed to “larval”, as RNA-seq was performed on central nervous systems of third instar larvae (not embryos) in Kaushal et al. 2021.

Thank you for catching this, we have fixed this in the manuscript.

Submitted filename: 20220607_revision_response.docx

Click here for additional data file.

13 Jun 2022

Single-nucleus RNA-sequencing in pre-cellularization Drosophila melanogaster embryos

PONE-D-22-01112R2

Dear Dr. Albright,

thank you for sending your further revised manuscript. I have now read through your revision and answers to the remaining points raised by reviewer 1. Your revision has succeeded to address all remaining concerns in a satisfactory manner and your study is now scientifically suitable for publication. Your manuscript will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Anton Wutz

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

17 Jun 2022

PONE-D-22-01112R2

Single-nucleus RNA-sequencing in pre-cellularization Drosophila melanogaster embryos

Dear Dr. Albright:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Anton Wutz

Academic Editor

PLOS ONE