Recent advances in single-cell RNA sequencing (scRNA-seq) have dramatically expanded our understanding of the cellular complexity and heterogeneity of human tissues, including endocrine glands (Zhang et al., 2020; Kameneva et al., 2022; Dong et al., 2020). However, further progress in this field struggled with the incomplete molecular characterization of the particular cell types being responsible for the functional complexity of human endocrine tissues. One of the most problematic issues in the scRNA-seq profiling of human tissues that significantly impacts the biological relevance of the ultimate data is the sample preparation step. In that, the most commonly used immediate processing of freshly isolated tissues is extremely poorly integrated in the clinical logistics of tedious surgical procedures and subsequent surgical pathology assessments. This issue is particularly important in the context of the inherent problems in obtaining the high-quality single cell suspensions from solid tissues, requiring complex disaggregation/dissociation steps while preserving the high cellular yield and viability, unbiased cellular contents, transcriptional profiles, cellular states, etc. The suboptimal procedures used for tissue collection, storage, and especially disaggregation/cell dissociation, which involve excessive mechanical stress, suboptimal temperature conditions, prolonged enzymatic digestion, and the loss of the original tissue context, have been shown to significantly distort the resulting scRNA-seq data and may lead to cell type misclassification (Slyper et al., 2020; van den Brink et al., 2017).

Since the late 1970s, when the first methods for the disaggregation of solid tissues were described (Waymouth, 1974; Cerra et al., 1990), a variety of protocols utilizing mechanical, enzymatic, and chemical methods of dissociation (and combinations thereof) have emerged (Cunningham, 2010; Denisenko et al., 2020). Probably the most popular approach for obtaining single cell suspensions for scRNA-seq is enzymatic digestion implying the incubation of gross tissue samples with various proteases at 37°C as a key step. However, a number of studies have demonstrated that the employment of these techniques is associated with the profound activation of the stress signaling pathways and increased cell death, resulting in a significant bias in the scRNA-seq profiles (van den Brink et al., 2017; Lake et al., 2016). Such undesirable effects may be largely diminished via employment of the so-called “cold” dissociation techniques with a cold active protease (6°C) (O’Flanagan et al., 2019), however, at the cost of less efficient target cell dissociation (Deleersnijder et al., 2021). An alternative strategy to prevent shifts in transcriptional profiles is the use of general transcription inhibitors, such as Actinomycin D (ActD) (Lai et al., 2019). ActD has been shown to effectively mitigate artifact formation and block transcriptional changes during the dissociation of neuronal tissue, thus preserving the integrity of the transcriptional profile (Wu et al., 2017). However, the routine use of such inhibitors is limited by their toxicity, which depends on both concentration and cell type.

Furthermore, the enzymatic dissociation of normal and diseased endocrine tissues may be particularly challenging due to a number of confounding structural issues, such as a high lipid content in the normal adrenal cortex and adrenocortical neoplasia, or extensive stromal/capsular fibrosis and calcifications in the well-differentiated thyroid tumors. Single-nucleus RNA sequencing (snRNA-seq) protocols being compatible with the use of the fresh-frozen tissue samples may be successfully employed to overcome these limitations (Lake et al., 2016; Schmøkel et al., 2023); however, at the cost of the loss of a significant amount of the mature cytoplasmic mRNA, resulting in a lower coverage and poor representation of the rare cell types (Thrupp et al., 2020).

Each of the abovementioned approaches contributes to the emergence of a variety of artifacts related to the distortion of the transcriptional profiles of individual cells (Deleersnijder et al., 2021), an issue that must be explicitly addressed in the process of the ultimate data analysis (Massoni-Badosa et al., 2020). For example, the immediate-early response genes (e.g., the members of the FOS and JUN gene families) are primary candidates for changing their expression during single-cell dissociation at 37°C. Artifactual changes in gene expression patterns were investigated by comparing the transcriptional profiles of cryopreserved and living cells or methanol-fixed and living cells obtained from tissue dissociation using cold-active protease and enzymatic digestion at 37°C. This study showed that cold-active proteases dramatically reduce the number of scRNA-seq artifacts in the mouse kidney (Adam et al., 2017). However, addressing of these issues remains rather fragmentary and limited to certain tissue types and protocols, so numerous artifacts still need to be confidently addressed.

Many of the limitations of the currently employed techniques may be potentially overcome via simultaneous tissue dissociation and cell fixation, the procedure being capable of maintaining a high RNA Integrity Number (RIN) while minimizing the sample preparation-related distortions in the transcriptional profiles. Recently, García-Castro et al. (2021) introduced such a protocol whose prototype may be tracked back to the end of the 19th century (Canmillo Schneider, 1890; David, 1973), when Schneider reported the so-called “maceration” technique. In its contemporary variant, named ACME (ACetic (acid)-MEthanol), this technique reportedly produces quality suspensions of fixed single cells from planarians, D. melanogaster, D. rerio, and M. musculus tissues, where suspensions maintaining the high-integrity RNA may be further successfully cryopreserved using DMSO (García-Castro et al., 2021).

Here we extensively optimized and successfully implemented the unique ACME High Salt (ACME HS, see Methods) protocol for the single-cell transcriptomics analysis of human endocrine tissues represented by tumors arising from the adrenal medulla, adrenal cortex, pituitary gland, and thyroid follicular cells. Additionally, the incorporation of the high-salt washing buffer (3xSSC*) in ACME HS protocol instead of the standard PBS allowed us to stabilize RNA and prevent RNases reactivation during rehydration. The increased ionic strength of 3xSSC* buffer affects protein solubility and interactions, creating an environment that effectively impedes RNases function. Higher salt concentrations can alter protein structures, thereby reducing enzymatic activity and making it difficult for RNases to function effectively (Baba et al., 2017). We also compared our modified ACME HS and enzymatic dissociation methods for scRNA and nuclei isolation for snRNA profiling in terms of the number of the cells/nuclei recovered, RNA integrity, aligning of the resulting scRNA/snRNA data with the reference organ-specific profiles, and representation of specific cell types.

We clearly demonstrated that scRNA profiling of single cell suspensions obtained using ACME HS and enxymatic dissociation methods significantly outperformed snRNA profiling in terms of marker genes expression analysis while demonstrating in-between comparable performances in virtually all implemented analyses. Additionally, the modified ACME HS protocol allows successful cryopreservation of dissociated/fixed cells without sacrificing the mRNA yield and integrity. To our knowledge, this is the first report on successful implementation of the ACME HS technique in primary human tissues, and we believe that this protocol should significantly promote the scRNA studies in humans that are to be explicitly compliant with the real-life infrastructure and logistics of the surgical care centers.

Twenty-eight human adrenal gland neoplasms (including 12 adrenocortical tumors and 16 adrenal medullary tumors), nine thyroid carcinomas, and 12 pituitary neuroendocrine tumors (PitNETs) were acquired from the Endocrinology Research Centre, Moscow, Russia (Supplementary Table S1). In all patients, tumor specimens were definitively diagnosed by imaging, surgery, and histopathological examination. Each study participant gave written informed consent. In addition, single cells were isolated from fresh and fresh-frozen adrenal medullary tumor, adrenocortical tumor, thyroid carcinoma, and PitNET samples. After sampling, the tissues were placed in a cold Tissue Storage Solution (Miltenyi Biotec) pending dissociation. The fresh-frozen samples were stored at −80°C until processed.

The raw sequenced reads were processed with 10X Cell Ranger (v6.1.1). Default Cell Ranger quality check measurements were used for further comparison methods through the Wilcoxon test. The expression matrixes for the filtered cells were submitted to Seurat (Hao et al., 2023) (v4.9.9 and v5.0.0) for basic analysis, including scaling and normalization. Cell filtering based on gene/molecule dependency was done by pagoda2 (pagoda2, 2023) (v1.0.11). For the pagoda2 package, we utilized a minimum cell size of 1,000 and a maximum cell size of 50,000. The p.level was set as the minimum between 0.001 and 1/(number of columns in countMatrix). Doublets and ambient RNA content were calculated with scrublet (Wolock et al., 2019) (v0.2.3), SoupX (Young and Behjati, 2020) (v1.6.2) and decontX (Yang et al., 2020) (v3.18), respectively, with default settings. Specifically, for the Scrublet package, the following hyperparameters were applied: the number of doublets simulated relative to the number of observed transcriptomes was set to 2.0. The estimated doublet rate for the experiment was 0.05, with a standard deviation of 0.02. The sampling rate for UMIs when creating synthetic doublets was set to 1.0. The number of principal components used to embed the transcriptomes before k-nearest-neighbor graph construction was fixed at 30. For the SoupX package, the term frequency–inverse document frequency (tf-idf) minimum was established at 1. We only included genes that were at or above the 0.9 expression quantile. The maximum number of markers was limited to 100, and the contamination range was set to c (0.01, 0.08). The maximal false discovery rate was fixed at 0.2, while the mode of the gamma distribution prior on the contamination fraction was set to 0.05, with a standard deviation of 0.1. For the decontX package, the maximum iterations of the expectation-maximization (EM) algorithm were set to 500. We defined concentration parameters for the Dirichlet prior for the contamination in each cell (native and contamination counts) as c (10, 10), updating these values with fit_dirichlet during each iteration. The convergence threshold for the EM algorithm was established at 0.001, and the log likelihood was set to 10. We selected 5,000 variable genes for dimensionality reduction prior to clustering, with the clustering resolution parameter used in ‘dbscan’ to estimate broad cell clusters set to 1. A seed value of 12,345 was provided to seed_with for reproducibility. The means for doublets and ambient RNA values per sample were compared between sample preparation methods (ACME HS, enzymatic, nuclei). Major cell types were identified by the label propagation function using Conos (conos, 2023) (v1.5.0) and reference datasets (Jansky et al., 2021; Kildisiute et al., 2021; Han et al., 2020; Torgersen et al., 1995; Zhang et al., 2022). Sample integration was conducted by applying RunHarmony on preprocessed Seurat objects (Korsunsky et al., 2019). Velocity analysis was performed with Velocyto (v0.17) on 1,000 cells subset per sample and visualized using Velocyto.R (La Manno et al., 2018) (v0.6) on integrated embeddings. In addition, we applied scVelo (v0.3.2) (Bergen et al., 2020) on whole samples in order to obtain confidence levels of velocity estimations and pseudotime. Differential expression analysis and cell cycle phase predictions were processed with Seurat (v5.0.0). Significance of phases enrichment over cell types calculated with chi-squared test. A functional enrichment test was performed for differentially expressed genes with clusterProfiler (Wu et al., 2021) and wiki pathways (Agrawal et al., 2023) as reference databases.

We used the PercentageFeatureSet function with default parameters from the Seurat package to evaluate the impact of various methods—ACME HS, enzymatic dissociation, and nuclei isolation—on cellular stress, focusing on their potential to induce stress or favor necrosis or apoptosis in cells. This function calculates the percentage of total counts assigned to a specified set of genes. For the apoptosis signature, we curated a gene signature encompassing CASP3, BAX, BAD, BID, APAF1, TP53, FAS, TNFRSF10B, CYCS, BCL2, and AIFM1(Kiraz et al., 2016). Meanwhile, for the necrosis signature, we selected HMGB1 (Nikoletopoulou et al., 2013), ATP5F1A, CALR, ARHGAP45, S100A8, S100A9, NAMPT, ANXA1, KRT18, TNF, and AGER genes (Kanehisa et al., 2021; Kanehisa, 2019). We assessed various modalities of cell stress, including oxidative stress, cellular senescence, DNA damage, heat shock, and the unfolded protein response. The oxidative stress signature was constructed using NFE2L2, KEAP1, SOD1, CAT, HMOX1, GCLC, GCLM, NQO1, and PRDX1 genes (Kanehisa et al., 2021; Kanehisa, 2019). The markers of cellular senescence included CDKN1A, CDKN2A, IGFBP3, GADD45A, CCND1, CDKN2B, IL1A, IL1B, IL6, IL10, HMGA1, HMGB2, and UBB(Kanehisa et al., 2021; Kanehisa, 2019). DNA damage signatures comprised TP53, BRCA1, CHEK2, ATM, RAD51, RPA1, MDM2, ATR, and XRCC5(Kanehisa et al., 2021; Kanehisa, 2019). For the heat shock signature, we considered the HSP family genes HSPB, HSPG2, HSPB11, HSPA6, HSPD1, HSPE1, HSPBAP1, HSPA4L, HSPB3, HSPA4, HSPA9, HSPA1L, HSPA1A, HSPA1B, HSP90AB1, HSPB1, HSPA5, HSPA14, HSPA14.1, HSPA12A, HSPB2, HSPA8, HSP90B1, HSPB8, HSPH1, HSPA2, HSP90AA1, HSPB9, HSPB6, HSPBP1, HSPA12B, and HSPA13(Yer et al., 2018; Kampinga et al., 2009). Lastly, the unfolded protein response signature included ATF4, ATF6, XBP1, HSPA5, DDIT3, HERPUD1, DNAJC3, ERN1, ERN2, and PDIA6 genes (Kanehisa et al., 2021; Kanehisa, 2019). Using a two-tailed Wilcoxon rank-sum test, we calculated statistically significant differences in the signature scores between the various dissociation methods.

Single cells were isolated from tissue samples using both ACME HS and enzymatic dissociation methods, freshly and followed by cryopreservation in 3xSSC*10% DMSO or methanol cell fixation, respectively (Figure 1A; Supplementary Table S1). The key adaptations of the ACME HS method included supplementing the solution with 0.1M NAC and adding extra washing steps for isolated cells using cold high-salt 3xSSC* buffer. ACME dissociation was performed for approximately 1 hour on a rotator at room temperature, with periodic pipetting. Afterward, the ACME solution was removed, and the pellet was washed using a two-step procedure with cold 3xSSC* (see Methods for details).

Figure 1. Comparison of RNA integrity, morphology, and storage of ACME HS and enzyme-dissociated adrenocortical cells. (A) Schematic representation of a workflow for single-cell or single-nuclei processing and analysis of fresh and fresh-frozen tissues (created with BioRender.com). (B) Gel image of isolated total RNA from cryopreserved ACME HS-dissociated adrenocortical cells after 1,3,7,14, and 28°days of freezing at −80°C, and of freshly isolated adrenocortical cells kept at +4°C (0°days). (C) Gel image of isolated total RNA from adrenocortical cells obtained by enzymatic dissociation and fixed in 80% methanol after 1, 3, 7, 14, and 28°days of freezing at −80°C, and of freshly isolated cells kept at +4°C (0°days). (D) RNA integrity of fresh-frozen adrenocortical tumor, from which ACME HS-dissociated (B) and enzyme-dissociated (C) cells were obtained. (E). Bright field (BF) and confocal fluorescence microscopy images of freshly isolated ACME HS and enzyme-dissociated adrenocortical cells stained with Hoechst 33,342 (blue) and anti-CYP11B1 antibody (red), showing single cells, aggregates, and debris.

Total RNA extracted from freshly prepared cell suspensions (0 days) indicated that the RINs for cells obtained through the ACME HS dissociation method and enzymatic digestion were similar, at 7.6 and 8.3, respectively (Figures 1B, C). Similar RINs scores of RNA were found for cells obtained from adrenal medullary tumor, thyroid carcinoma and pituitary neuroendocrine tumor (PitNET) samples (Supplementary Figure S1). The obtained RIN scores were compared with the RIN score (9.8) of undissociated adrenocortical tumor (control) (Figure 1B).

Next, we visualized freshly prepared dissociated adrenocortical cells (0 days) using bright field and confocal microscopy. The cells maintained their morphology, displaying minimal aggregation and debris (Figure 1E). Additionally, microscopy was conducted for thyroid cells and one replicate of adrenocortical cells. (Supplementary Figure S2B). To identify adrenocortical cells, we stained the fixed cells with Hoechst 33342 and an anti-CYP11B1 (11β-hydroxylase) antibody conjugated with Alexa Fluor 594 (Figure 1E; Supplementary Figure S2B). CYP11B1 is localized in the inner mitochondrial membrane and is normally expressed in the zona fasciculata of the human adrenal cortex (Duparc et al., 2022). We observed intense immunofluorescence of CYP11B1 (red) in ACME HS-dissociated adrenocortical cells, in contrast to cells isolated using enzymatic digestion. This difference could be attributed to the extended permeabilization of cell membranes with methanol during the ACME HS protocol. The voids observed in the nuclei of adrenocortical stained cells are likely associated with their functional ability to efflux the DNA binding dye Hoechst 33342, resulting in the so-called side population (SP) (Crowley et al., 2016, p. 33). Thyroid follicular cells were also visualized by staining fixed cells with Hoechst 33342 and an anti-TSHR (TSH receptor) antibody conjugated with Alexa Fluor 594 (Supplementary Figure S2B).

A key disadvantage of commonly used enzymatic digestion methods is that freshly isolated tissues are extremely poorly integrated in the clinical logistics of surgical procedures. In contrast, the method we used enabled the freezing of fresh tissue samples before ACME HS dissociation. By combining tissue dissociation and cell fixation, the ACME HS method preserves cells within their context. Furthermore, the resulting suspensions can be cryopreserved for subsequent analysis (García-Castro et al., 2021). To confirm that cryopreservation of ACME HS-dissociated cells in 3xSSC* with 10% DMSO maintains RNA integrity, we evaluated the RNA quality in cells after different durations of cryopreservation.

We sequentially extracted total RNA from six aliquots of ACME HS (Figure 1B) and six aliquots of enzyme-dissociated (Figure 1C) cells obtained from the same adrenocortical tumor at various time intervals: immediately after preparing the single-cell suspension (0 days) and after 1, 3, 7, 14, and 28 days of cryopreservation or methanol fixation. The RIN scores obtained were compared with the control RIN scores (Figure 1B).

We found that the RNA integrity of ACME HS-dissociated cells was well-preserved during cryopreservation in 3xSSC* and 10% DMSO over the specified intervals (1, 3, 7, 14, and 28 days), with RIN scores averaging around 6.7 (Figure 1B). This preservation was consistent across adrenocortical tumors, adrenal medullary tumors, thyroid carcinoma, and PitNET, allowing for subsequent scRNA-seq analysis (Figure 1B; Supplementary Figure S1). Additionally, after 6 months of storage, the RNA integrity of cryopreserved ACME HS-dissociated adrenocortical cells was 5.9 (Supplementary Figure S2A). The RNA integrity of dissociated PitNET cells was assessed at only two time points (0 and 1 day) due to the small tissue sample size. The RIN of enzyme-dissociated adrenocortical cells fixed in 80% methanol decreased over time with significant degradation of ribosomal RNA evident in reduced or absent signals for the 18S and 28S peaks (Figure 1C; Supplementary Figure S1).

The RNA integrity of fresh-frozen adrenocortical tumor samples after 1, 3, 7, 14, and 28 days of freezing at −80°C ranged from 9.3 to 9.0, with no discernible patterns in RIN score changes over time (Figure 1D). These RIN scores indicate that the material is suitable for subsequent dissociation using the ACME HS method.

We used the flow cytometry method to assess the quality of ACME HS and enzyme-dissociated samples. We compared the amount of debris, aggregates, and single cells obtained by the two protocols of tissue dissociation–ACME HS and enzymatic across different tumor types. We evaluated the ACME-HS method by three independent repeats for each tumor type. It appeared more straightforward to calculate the amount of cellular debris and the number of singlets and aggregates at the FSC-area/FSC-height dot plot (Figure 2A1) than at the FSC/SSC dot plot (Figure 2A2), as the boundary between debris and singlets, and between singlets and their aggregates, is usually poor. We gated out the area of debris (green) as small events with the highest ratio of FSC-height to FSC-area signal. Singlets (red) were selected based on their well-correlated height versus area signal, while aggregates of cells (black) had an increased area signal compared to the height signal. (Figure 2A1; Supplementary Figure S3, A1–H1).

Figure 2. Representative flow cytometry data of the samples prepared by ACME HS and enzymatic dissociation methods. (A) Flow cytometry data of the ACME HS dissociated adrenocortical tumor sample (replicate 1). A1. FSC-height/FSC-area dot plot was used to calculate cellular debris, single cells, and cellular aggregates (green events - for debris, red events - for singlets, black events—for aggregates). A2. FSC/SSC dot plot demonstrating the distribution of cells, their aggregates, and cellular debris based on their light-scattering properties. A3. Pie diagram of the debris, singlets and aggregates distribution in the sample. A4. PI-height/PI-area dot plot of singlets used for additional gating of single cells (shown in blue) among nucleated cells and their aggregates (shown in red). A5. DNA histogram from single events showing the cell cycle distribution for all cells in the sample, with percentages of the cell cycle phases (G0/G1, S, G2/M) inserted. (B–E) Flow cytometry data for the different samples (replicate 1). B1–B4 – for the adrenocortical tumor; C1–C4 — for the adrenal medullary tumor; D1–D4 — for the thyroid carcinoma; E1–E4 — for PitNET (green events - for debris, red events - for singlets, black events — for aggregates). Index one stands for FSC/SSC dot plots for the samples obtained by the ACME HS protocol. Index two stands for CFSE/PI dot plots for the samples obtained by the ACME HS protocol. Index three stands for FSC/SSC dot plots for the samples obtained by enzymatic dissociation. Index four stands for CFSE/PI dot plots for the samples obtained by enzymatic dissociation. Created with BioRender.com.

Further, we used a backgating strategy to evaluate the data in other dot plots to correct our subsequent gating for the best separation and calculation of the amount of cellular debris and the number of singlets and aggregates. First, we applied color gating and visualized debris, singlets and aggregates in the FSC/SSC dot plots for all tissues (Figures 2A, 2B1–E1, B3–E3; Supplementary Figure S3, A2–H2). We made sure that debris was located in the lower left area of the dot plot, and the aggregates formed clusters in the upper right area of the FSC/SSC dot plot. Variations in the ratios of singlets and debris across the different tumors studied are expected due to the unique structures and components of these tumor tissues. The highest amounts of singlets were found in adrenal medullary tumor samples, ranging from 15.15% to 54.51%, while the lowest were in thyroid carcinoma samples, ranging from 5.14% to 36.12%. PitNET and adrenocortical tumor samples had average single cell ratios of 28.57% and 27.95%, respectively. Debris levels were consistent across all tissue types, averaging 53.41%. Adrenocortical tumor samples had the lowest aggregate levels, averaging 15.26% compared to other tissues.

Next, we stained the samples with a DNA-binding dye–propidium iodide (PI), to better discriminate the nucleus-contained cells from the nuclear-free debris. We studied the PI-height/PI-area dot plots and checked the positions of debris, singlets, and aggregates among all ungated events including debris and aggregates (Supplementary Figure S3, A3–H3). Then, we gated the single cells (blue) among the nucleated cells and their aggregates (red) (Figure 2A4; Supplementary Figure S3, A4–H4) and analysed DNA histograms from single cells (Figure 2A5; Supplementary Figure S3, A5–H5). In addition, we found that more than 60% of analyzed single cells are located in the G0/G1 cell cycle phase, in all tumor types. Much less cells are located in G2/M and S-phase in all tumor types.

Although our dissociation protocols differed from those that were specifically elaborated for cell cycle analysis and frequently used (Darzynkiewicz, 2011), in most cases, we could resolve various phases of the cell cycle (G0/G1, S, G2/M) by mean fluorescence intensity (MFI) per cell in DNA histograms: DNA content in G2/M phase was as expected, two times more than in G0/G1 as shown in Supplementary Figure S3, A5–H5.

To discriminate the nature of the debris, we stained the samples with 5,6-carboxyfluorescein diacetate succinimidyl ester (CFSE). This dye is frequently used in flow cytometric protocols for live cells labeling due to its ability to bind to intracellular molecules, primarily to amine groups. In addition to its role in viable cell staining, CFSE can trace dying cells in composite samples (Dumitriu et al., 2001). As shown in Figures 2B2–E2, B4–E4, most green events matching debris turned out to be CFSE-positive and PI-negative, which suggested that the debris was generally nuclear-free. By comparing the frequency of debris and aggregates (Supplementary Figure S3, A1–H1) and analyzing dot plots, we suggest that both ACME HS and enzymatic methods induced a relatively similar number of aggregates and debris. Despite the large amount of debris and aggregates, which was expected, we observed a sufficient number of single cells in our samples obtained by the ACME HS and enzymatic dissociation protocols (Figures 2A, B1–E1, B3–E3; Supplementary Figure S3, A2–H2).

For comparative analyses, we selected 41 human endocrine tumor samples (Supplementary Table S1). We obtained 107,875 cells and nuclei isolated from adrenocortical tumors (n = 12), 94,807 cells and nuclei from adrenal medullary tumors (n = 15), 41,418 cells and nuclei from PitNETs (n = 9), and 60.365 cells from thyroid carcinomas (n = 5).

First, we examined the summary statistics for the generated single cell gene libraries (Figure 3A; Supplementary Figure S4). We found that the quality of the ACME HS dataset was almost identical to the enzymatic and nuclei datasets, obtained from ACME HS-dissociated and enzyme-dissociated whole cells, and isolated nuclei, respectively. Some differences in ribosomal and mitochondrial gene expression levels between nuclei and whole cells obtained by ACME HS and enzymatic dissociation methods were confirmed directly, as well as the exon/intron alignment ratio (Figure 3A). Nuclei-based data showed clear advantages in terms of total genes detected, as expected. At the same time, there were no significant differences in quality parameters, namely, the total number of cells, reads in cells for all samples (Figure 3A). We observed the patterns mentioned earlier in the four tissues with different dissociation methods (Supplementary Figure S4).

Figure 3. ACME HS demonstrates consistency with the enzymatic protocol and keeps advantages over the nuclei isolation. (A) Standard single-cell sample features (number of cells; read the cells%; exon/intron mapped; mitochondrial, ribosomal, nFeature_RNa; and total genes) of ACME HS, enzymatic and nuclei samples. Statistical differences estimated by the Wilcoxon rank-sum test: **** (0.0001 < p < 0.001), *** (p < 0.001), ns - not significant–p > 0.05. (B) Major cell type compositions among preparation methods, namely, adrenocortical, chromaffin, pituitary, and thyroid follicular cells. (C) Fractions of defined cells identified by different dissociation methods (ACME HS, enzymatic, and nuclei) for each tissue type–adrenocortical tumor, adrenal medullary tumor, thyroid carcinoma, and PitNET. The diagram does not indicate the number of cells representing less than 5% of the total number.

To compare the representation of distinct cell types and states between the enzyme, nuclei, and ACME HS datasets, we used Seurat (Hao et al., 2023) to generate integrated embeddings and annotations for each tissue type. ACME HS successfully integrated into enzymatic and nuclei datasets (Supplementary Figure S5). The major cell types were consistently defined and integrated across all three methods. In all four types of tissues, the ACME HS data retained tissue-specific cells, namely, adrenocortical, chromaffin, pituitary, thyroid follicular cells, and other nonspecific cells (Figures 3B,C). The heterogeneity of the major cell populations (adrenocortical, chromaffin, thyroid follicular, and pituitary cells) was assessed by further clustering of the integrated cells (Figure 4A; Supplementary Figure S6). Minor subclusters (<100 cells) were excluded from the analysis. We found that most subclusters (A-2, A-3, A-5, A-4, A-7, A-18, and A-19) were lost for the adrenocortical samples in the nuclei datasets, unlike in the ACME HS and enzyme datasets.

Figure 4. Heterogeneity and functional characterization of cell clusters specific for ACME HS, enzymatic and nuclei datasets (A) Visualization of the major cell subpopulations and states for adrenocortical and chromaffin cells. Adrenocortical and chromaffin cells were segregated by cell clustering applied on integrated data sets. Cells segregated into small clusters (<100 cells) were combined into separate minor groups and excluded from the analysis. (B) Differentially expressed genes (upregulated in selected clusters compared to all the rest cells, logFC >0.25 and FDR <0.05) were calculated for each cluster within major cell types and analyzed by using wiki pathways as a reference database for main clusters of adrenocortical subpopulations—A-3, A-4, A-5, and chromaffin cells—C-2, C-3, C-4.

Loss of numerous subclusters was also observed in the nuclei datasets of the chromaffin and pituitary samples (C-2, C-4, C-8, C-20, C-23, and C-25 and P-1, P-2, P-3, P-5, P-10, P-15, P-16, and P-17, respectively). However, some subclusters were more enriched in the nuclei datasets for chromaffin samples (С-0, С-1, С-5, С-6, С-14, and С-19) (Figure 4A; Supplementary Figure S6).

While some cell subpopulation variability is expected due to the individual tissue conditions, the most apparent difference is determined for nuclei-based samples. The main clusters of adrenocortical subpopulations–A-3, A-4, A-5, and chromaffin cells–C-2 C-3 C-4 enriched in oxidative phosphorylation, electron transport chain, ribosomal, mitochondrial, and mRNA processing genes are missing from nuclear datasets as opposed to ACME HS and enzymatic samples (Figure 4B). Because of the low cell numbers obtained from thyroid carcinoma samples with ACME HS-dissociation, comparisons of major thyroid follicular cell population clusters obtained with ACME HS and enzymatic dissociation methods were not possible for these samples.

We next examined whether ACME HS dissociation was able to preserve the expression of tissue-specific marker genes in all four tissue samples. In total, 107.875 cells were obtained from adrenocortical tumor samples (n = 12) and integrated. Adrenocortical cells accounted for 85.6% of the annotated cell types for (Figure 3C). These cells expressed key literature-derived marker genes that identify adrenal cortex zones: CYP11B2 (van de Wiel et al., 2022) for zona glomerulosa, CYP11B1 (Duparc et al., 2022) for zona fasciculata, CYP17A1, SULT2A1 (Janšáková et al., 2020), and CYB5A (Nakamura et al., 2011) for zona reticularis. These genes were detected in all sample groups regardless of the extraction method, except for CYP11B2, which was not detected in the enzymatic and nuclei samples (Figure 5A; Supplementary Figure S7A). Similar results were obtained for the adrenal medullary tumor (n = 15) and PitNET samples (n = 9) with 94.807 (89% chromaffin cells) and 41.418 (79.6% pituitary cells) integrated cells, respectively (Figure 3C). Correspondingly, these cells expressed key marker genes, such as CHGA, SYP (Mete et al., 2022), DBH, PNMT (Konosu-Fukaya et al., 2018) (Figure 5B; Supplementary Figure S7B) and POMC (only in the ACME HS datasets), with the exception of GH1 and POU1F1 in the ACME HS datasets (Figure 5C; Supplementary Figure S7C). Although thyroid follicular cells represented the majority (95.2%) of the 60.365 cells in thyroid gland samples (n = 5), ACME HS-dissociated cells exhibited almost no expression of key markers such as TG and TSHR (Lee et al., 2015; Liu et al., 2022) (Figure 5D; Supplementary Figure S7D).

Figure 5. Expression of tissue-specific genes in nuclei and whole cells obtained by ACME HS and enzymatic dissociation methods. (A) Key tissue-specific gene expression according to UMAP visualization, namely, CYP11B2, CYP11B1, and SULT2A1 for adrenocortical tumors; (B) CHGA for adrenal medullary tumors; (C) POMC for PitNETs; (D) TG for thyroid carcinoma samples. (E) The distribution of the enrichment scores of heat shock, apoptosis, necrosis, DNA damage, cell senescence, and unfolded protein response signatures across preparation methods of the datasets obtained from adrenocortical tumor, adrenal medullary tumor, thyroid carcinoma, and PitNET samples (combined dataset, n = 41).

Then, we analysed a panel of the top genes specific for adrenocortical, chromaffin, thyroid follicular, and pituitary cells (Supplementary Figure S7E). It turned out that each dissociation method allows to estimate the expression signatures of different genes that have a little overlap in ACME HS, enzymatic, and nuclear samples, with the exception of PitNETs.

Since dissociation and preservation techniques can induce cellular stress as evidenced by changes at the transcriptomic level, we examined the expression of key markers of stress and cell death. Specifically, markers associated with apoptosis, necrosis, cellular senescence, DNA damage, heat shock, and the unfolded protein response (UPR) were assessed in the ACME HS (n = 13), enzyme (n = 19), and nuclei (n = 9) datasets for all tissue samples (Figure 5E; Supplementary Figure S8; Supplementary Tables S1, 2).

We found no differences in the expression levels of genes related to apoptosis in any of the datasets. The necrosis gene signature showed no differences between the ACME HS and enzymatic samples but was lower in the nuclei samples. Nevertheless, the most significant expression of TNF (tumor necrosis factor) was observed in adrenocortical tumors and thyroid carcinomas in ACME HS and nuclei samples (Supplementary Table S2). The DNA damage gene signature exhibited minimal variations among all three methods. In addition, the heat shock protein signature showed minimal differences between the ACME HS and enzymatic samples but was significantly lower for nuclei samples. The UPR signature did not differ between ACME HS and enzymatic samples but was significantly lower for nuclei samples. In particular, ERN2, a UPR marker, was highly expressed in adrenal medullary tumor and PitNET datasets obtained by ACME HS and nuclei isolation methods, in adrenocortical tumors by enzymatic digestion, and in thyroid carcinomas by the ACME HS method. In addition, the cellular senescence signatures also showed minimal differences between the ACME HS and enzymatic samples. IL1B was highly expressed in adrenocortical and medullary tumor samples obtained from the ACME HS and nuclei datasets, while IL6 was highly expressed in all tissues compared with the PitNET datasets obtained by the enzymatic digestion method. Another important senescence marker, CDKN2B, was highly expressed in adrenal medullary tumor samples obtained by nuclei isolation as well as in thyroid carcinoma samples obtained by an enzymatic approach. Senescence markers HMGA1 and UBB were highly expressed in the nuclei datasets for adrenocortical and medullary tumor samples (Supplementary Table S2).

Next, we performed velocity analysis for individual samples to demonstrate the consistency between the ACME HS (n = 1) and enzymatic (n = 1) protocols together with their advantages over the nuclei isolation method (n = 1). For adrenocortical (Figure 6A) and chromaffin cells (Figure 6B), we identified similar velocity directionalities from ACME HS and enzymatic-specific clusters towards cell populations commonly shared between methods. This was confirmed with significance assessments (Supplementary Figure S9). The RNA velocity recapitulated the transcriptional dynamics within these datasets, including the general movement of the differentiating adrenocortical and chromaffin cells, as well as movement towards and away from the intermediate differentiation state. The velocity also captured the cell cycle dynamics involved in cell differentiation.