Deconvolution of spatial sequencing provides accurate characterization of hESC-derived DA transplants in vivo

Cell therapy for Parkinson’s disease has experienced substantial growth in the past decades with several ongoing clinical trials. Despite increasing refinement of differentiation protocols and standardization of the transplanted neural precursors, the transcriptomic analysis of cells in the transplant after its full maturation in vivo has not been thoroughly investigated. Here, we present spatial transcriptomics analysis of fully differentiated grafts in their host tissue. Unlike earlier transcriptomics analyses using single-cell technologies, we observe that cells derived from human embryonic stem cells (hESCs) in the grafts adopt mature dopaminergic signatures. We show that the presence of phenotypic dopaminergic genes, which were found to be differentially expressed in the transplants, is concentrated toward the edges of the grafts, in agreement with the immunohistochemical analyses. Deconvolution shows dopamine neurons being the dominating cell type in many features beneath the graft area. These findings further support the preferred environmental niche of TH-positive cells and confirm their dopaminergic phenotype through the presence of multiple dopaminergic markers.


INTRODUCTION
The degeneration of dopaminergic neurons in Parkinson's disease (PD) is the leading cause of motor, and some non-motor, symptoms of the disease. 1 This has resulted in the development of therapies attempting to re-supply the lost dopamine (DA) signaling in the brain, ranging from pharmacological interventions 2,3 to cell-based therapies. [4][5][6] Cell-based therapies placing DA-producing grafts in the putamen 7,8 entered the scene several decades ago. [9][10][11][12] Initially, the dopaminergic neuroblasts were isolated from the early ventral mesencephalon (VM) derived from human fetal tissue. 7,13 However, because of the fundamental issues with this approach, 14 research has moved toward the differentiation of VM dopaminergic neurons from pluripotent stem cells. 15,16 Human embryonic stem cells (hESCs) are the most researched, replenishable source of dopaminergic transplants [17][18][19][20] although human induced pluripotent stem cells (hIPSCs) have also been used successfully. 21,22 Many studies have shown the efficiency of DA transplants in ameliorating motor function in animal models of PD, [23][24][25][26][27][28] leading to the approach entering clinical trials. 5,29,30 New clinical studies based on hESC-derived DA progenitors have recently been initiated in the United States and Sweden (clinicaltrials.gov: NCT04802733 and clinicaltrials.gov: NCT05635409). Long-term follow-up of fetal transplants has shown their potential for alleviating motor symptoms in PD patients for nearly two decades after transplantation 7,31 without additional treatment. 7 As mature dopaminergic neurons will not survive grafting, dopaminergic neural VM progenitors, which are not fully maturated, are used for transplantation. [32][33][34] This approach, however, leaves the transplant's final composition challenging to predict and control. 35 The final graft composition in vivo from hESC-derived dopaminergic transplants has been analyzed in bulk 36 and at single-cell resolution. 18 However, the quantitative results of the latter study were at odds with histological assessments since the sequencing study suggested that neurons make up the minority of the transplant. When assessed in situ with immunohistochemistry (IHC), cells positive for neuronal and dopaminergic markers represent over 50% of all cells in a transplant. 18,37 As the authors suggest, 36 this discrepancy could result from the cells' differential sensitivity to the processing as the cell proportions were skewed toward the more robust non-neuronal cells. Differences in cell composition across RNA sequencing technologies is a well-known phenomenon, 36,38 further supporting the need for minimal processing of the samples for unbiased results.
To overcome cell type specific biases and gain more insight into the transplants in the context of their environment, we used the spatial transcriptomics (ST) workflow. 39 We analyzed, in depth, three hESC-derived dopaminergic transplants maturated and integrated within the host tissue. Among the advantages of ST analysis are the preservation of the spatial relationships and minimal tissue preprocessing, which favors an accurate snapshot of the tissue state. Spatial methods are developing in several directions. Some sacrifice singlecell resolution in favor of unbiased transcriptome capture, [39][40][41] while others focus on high resolution and detection of predetermined genes. [42][43][44][45] Although spatial data can give crucial information about the analyzed tissue, most currently available methods are not standardized enough to be used in clinical development while retaining the ability to capture the entire transcriptome. With commercialization, ST has overcome much of this hurdle and it has also been shown to produce valuable data in clinical research. [46][47][48][49] As hESC-based transplantation is now entering clinical trials in PD, we have chosen to generate the data for those cells on the basis of ST as a foundation for future postmortem analysis.

Generation of dopaminergic transplants matured in an animal model
In this study, dopaminergic progenitors, derived from hESCs and differentiated in vitro toward a VM phenotype, were grafted into the lesioned striatum of a toxicant-induced immunocompromised (nude) rat model of PD ( Figure 1A). 6-hydroxy-DA (6-OHDA) was injected into the right medial forebrain bundle resulting in a nearcomplete unilateral depletion of dopaminergic neurons in the substantia nigra pars compacta (SNpc).
The animals were sacrificed 6-15 months after grafting to allow functional maturation of the grafted cells. Brain tissue was flash-frozen and cryosectioned, coronal sections were placed on amine-binding slides (Surmodics) with ST oligo-dT array-printed spots (100mm in diameter). 39 After fixation, tissue was stained with H&E, revealing grafted areas in the injected striatal region ( Figures 1B-1B 00 ).
The mRNA from the striatal tissue and the SNpc control sections was then subjected to the ST protocol. 50 Briefly, the RNA is hybridized to the barcoded oligo-dT probes on the glass surface, the cDNA is generated in situ, and the tissue is digested. After that, the cDNA is cleaved away from the glass, linearly amplified, and prepared for Illumina sequencing. With the help of the coordinate-specific, barcoded probes printed in the array spots, the transcriptomic profile of each spot can then be projected back onto the tissue. The number of transcripts per feature for all samples is shown in Figure S1. As there are transcripts of two species being analyzed, the transcripts have been aligned to a chimeric rat and human genome. For gene expression analysis, the species identity was stripped away. However, transcripts aligning uniquely to a single species give a good glimpse into the species composition of features.

Batch effect detection and removal
Because of the complexity of the procedure, we first assessed the animal-to-animal and section-to-section variability to see the similarity between biological replicates ( Figures 1C and C 0 ). Without performing any batch correction, joint analysis of tissues from close anatomical sites showed considerable variation, both among sections and animals, with no meaningful overlap between sections. These differences among such biologically close tissues pointed to nested batch effects. To correct for these effects, we compared two emerging batch correction approaches, Scanorama 51 ( Figures 1D and D 0 ) and DESC 52 ( Figures 1E and E 0 ), to the commonly used batch correction algorithm within Seurat 53 ( Figures 1F and 1F 0 ). Although aimed primarily at single-cell RNA sequencing (scRNA-seq) datasets, most of these have previously been applied to ST data as well, with Seurat being specifically incorporated into the spatial data-specific module STUtility. 54 Although all three algorithms significantly diminished the distance among biological replicates, there were visible differences in how clusters based on the corrected gene expression were separated. All algorithms separated striatal, cortical, and septal regions, the white matter tracts, and the transplanted areas (Figures 1C 0 -1F 0 ) with varying degrees of modularity. Whereas DESC and Seurat managed to separate tissues on the basis of morphological areas, Scanorama showed a level of separation based on the host animal ( Figures 1D and 1D 0 ). Notably, unlike most anatomical structures, grafts were united throughout all animals, indicating a higher degree of similarity among them.
DESC, on the other hand, fully integrated all the sections regardless of the host. However, besides a handful of major clusters connected to morphological areas, about twice as many small communities appeared, not clearly related to anatomical regions (Figures 1E and 1E 0 ).
Seurat was best able to correct for batch effects both on the level of animals and sections while not sacrificing division within the tissues themselves or creating outlier clusters (Figures 1F, 1F 0 , and S2), as seen with DESC. The uniform manifold approximation and projection (UMAP) representation of the corrected dataset does not point to a bias at either slide, animal, or tissue section level. There is also an apparent separation of layers and functional regions within the cortex, in line with what is expected from the known anatomy of the rat brain. 55 Taken together, these findings support further continuation with the Seurat algorithm clustering.

Identification of human tissue on the transcriptome level
To properly delineate the grafted and host areas (Figures 2A-2A 00 ), we classified features on the basis of the ratio of transcripts aligning solely to a single location on either the human or the rat genome. In line with the expectations, areas with a high percentage of human transcripts were observed in regions with a histologically visible transplant (Figures 2B-2B 00 and S3A). Outside of these regions, the fraction of human reads tends to be negligible and can serve as a proxy for the quality of the well processing. High signal from human reads uniformly distributed throughout tissues can point to an ST sample with template switching caused by over-amplification and can be used to exclude sections with low signal-to-noise ratios.

Identification of anatomical areas
The Seurat corrected data were used as a basis for differential gene expression analysis aimed at distinguishing transcripts specific to the transplants (Figures 2A-2A 00 ). The transplant area was defined as a Seurat cluster with the highest human content. Because of the placements of the grafts, we have selected the surrounding striatum as the contrast region. The host striatum has been selected through feature clusters rich in medium spiny neuron (MSN) marker genes, 56 specifically PENK, ADORA2A, and PPP1R1B. This cluster outlined the striatal area as expected on the basis of the morphology of the tissue. In all cases, the expression of MSN marker genes strongly diminished in the transplant areas, confirming the differing composition of the transplant to its surrounding tissue (Figures 2C-2C 00 and S3B). Furthermore, SNpc, white matter tracts, and cortical areas were To overcome variability on the level of unique features for differential gene expression, all regions of interest within a tissue were treated as a pseudobulk to support the robustness of the analysis.
Differential gene expression was assessed by modeling the data on the basis of negative binomial distribution as implemented in the DESeq2 statistical module. We assessed transcriptomic differences in three scenarios: Between the transplant and the striatum, the transplant, and all the remaining host tissue in the grafted sections, and the transplant and the SNpc area ( Figures 2D-2F, S5A, and S5B 0 ). When looking at the genes most differentially expressed in the transplant compared with the striatum, we can see many well-established dopaminergic markers ( Figure 2D). Although these markers are very sparsely represented in the remainder of the grafted host tissues, they can be seen highly expressed also in the SNpc, pointing to a shared robust dopaminergic population. Many of these genes are also highly differentially expressed in the SNpc compared with the rest of the midbrain, cortical, and white matter regions, supporting the uniqueness of the dopaminergic population in both regions ( Figure S5C). These differences are also recapitulated in the analysis between the transplant areas Interestingly, there are very few or no significantly downregulated genes in the analysis of differentially expressed genes (DEGs) (Figures 2F, S5B, and S5C). This effect is most likely caused by a partial volume effect when features at the edges of unique areas are populated by cells from neighboring tissues. This is further supported by the results of a differential gene expression analysis between SNpc and the transplants, where classical MSN marker genes PENK and PPP1R1B show up as being expressed in the transplant. When assessing the expression of these genes on the ST assays, we can see a gradient of MSN markers from the striatum past the edges of the transplanted regions ( Figure S6B).

ST versus scRNA-seq composition analysis
Dopaminergic transplants, generated using this hESC cell line (RC17) and differentiation protocol 37 have previously been deeply characterized using scRNA-seq. 20 In the reanalysis of this previous study, three major clusters were observed based on the transcriptomic profiles generated with the scRNA-seq data ( Figures 3A and S13A). The smallest cluster carries neuronal characteristics, and the two larger clusters express glial and vascular leptomeningeal cell (VLMC) markers. The single-cell extraction protocol resulted in a very low recovery of neuronal cells, and only a small subset of these cells ex-pressed canonical markers of dopaminergic neurons. 36 The neuronal population was, however, very well differentiated from the larger population of captured glial cells.
The presence of marker genes of dopaminergic neurons, astrocytes, and VLMC has been shown across the scRNA-seq data as well as in the ST analysis of both the grafted tissue sections and the nonlesioned substantia nigra regions ( Figures 3B and 3C). Among dopaminergic neuron markers (TH and DAT [SLC6A3] shown), their expression concentrates toward the edges of the large transplant in the ST analysis while still covering a significant proportion of the transplants' features. Furthermore, species-specific analysis of the expressed DAergic marker genes in the grafted sections showed them to be of nearly universal human origin. In contrast, the opposite is true for the same marker gene transcripts in sections containing rat substantia nigra ( Figure S7). In scRNA-seq, only a minority of the cells in the neuronal cluster show an expression of TH or DAT. Glial cell markers (AQP4 and SOX9 for astrocytes and COL1A1 for VLMCs) are expressed mostly at the center of the graft in situ ( Figure 3B).
Graft-induced dyskinesia observed from fetal grafts has been proposed to be due to the presence of serotonergic neurons in the grafts 14,57 or through an extrinsic serotonergic stimulation of the grafted DA neurons. 58 In the grafts differentiated with the current protocol, no serotonergic neuron markers were found in the ST analysis compared with markers found in the midbrain regions, pointing to an absence of a serotonergic population in the transplants at this stage of the maturation ( Figure S6C).

Analysis of transplants' composition through scRNA-seq
To estimate cellular proportions in the dopaminergic grafts, we analyzed the scRNA-seq data using a reference dataset. 59 This dataset is composed of cell types from striatal, cortical, and ventral midbrain regions. Because of the imbalances in the dataset composition and the underrepresentation of cells positive for dopaminergic markers in the neuronal cluster, we chose SingleR, which assigns the probability of belonging to a cell type to each cell in the dataset separately. This approach avoids assigning cell identity on the cluster level and can resolve related populations such as the DA and other neurons. The SingleR assignment resulted in the cluster carrying neuronal characteristics being divided into multiple neuronal types, with dopaminergic neurons being the most prominent. Astrocytes, radial glial cells, VLMCs, and neuronal progenitors were distributed in a gradient throughout the two larger clusters ( Figures 4A and S8).

Integration of the scRNA-seq and the ST datasets
To explore the cell type composition of the tissues analyzed with ST, we leveraged the development of transcriptome deconvolution algorithms. Using a reference dataset, 59 deconvolution tools estimate the proportion of each cell type that contributes to the small bulk cDNA indexed on each ST feature. After comparing multiple software packages (Figure S9), we settled on cell2location, on the basis of a Bayesian model, and explicitly aimed at ST. 60 We supplied cell2location with the reference dataset used for the SingleR analysis described above, allowing us to estimate and compare the proportions of relevant cell types. Both sections originating from the grafted striatum and the host substantia nigra regions were assessed.
After deconvolution of the grafted striatal sections, known cell populations occupy their expected regions in the host tissues ( Figure 4F), and the transplant areas clearly stand out from their surroundings. When comparing the distribution of cell types in the transplant area and the surrounding striatal tissues, we find a clear enrichment of dopaminergic neurons, astrocytes, and neuronal progenitor cells ( Figures 4C, 4D, and S10A). Part of the transplant population has also been assigned to an inhibitory neuronal population closely related to both septal and ventral midbrain neurons ( Figure 4F).
Even though transplants stand out with their human transcript content, not all cells in the transplants are of human origin ( Figures 2B  and S3). To see how coupled various cell types are to the human origin, we correlated the cell type content with that of human transcripts in features. Dopaminergic neurons showed the highest correlation with human content of the feature, distantly followed by excitatory neurons and neuronal progenitor cells ( Figure S10B).
Compared with the scRNA-seq analysis of the transplants, we find large differences in the proportion of cell types (Figures 4E and  S10A). The ST deconvolution estimates the number of cells of a specific type in features on the basis of the whole-transcriptome profile matched to a scRNA-seq reference. This deconvolution resulted in the largest fraction of transplanted cells being categorized as dopaminergic neurons with a mean value of 11.3% for all three grafts (95% . Strikingly, the same cells only accounted for fewer than 1% of cells in the scRNA-seq dataset. VLMCs were enriched in the transplant area in both datasets, but their numbers were much smaller in the ST analysis than in the scRNA-seq. This discrepancy was even more pronounced for the radial glia-like cell population, limited mostly to the subventricular zone in the ST analysis. Neuronal progenitors was the only population that was present in similar proportion in the scRNA-seq dataset and in the ST. However, in the ST dataset, unlike the remaining populations, neuronal progenitor cells are also detectable in low amounts throughout the surrounding tissue. Comparison of the enriched cell types' proportions within the transplanted striatal to the VM sections also pointed out significant differences ( Figures 4B, S11, and S12). Whereas both harbor robust and specific dopaminergic populations, confined to SNpc and ventral tegmental area (VTA) in the VM sections and to the transplant in the striatal sections, VLMCs are enriched to a detectable level only in the choroid plexus epithelium in the midbrain sections and nearly completely missing in the SN and VTA. Thus, presenting a population unique to the transplants. Neuronal progenitor cells were observed in small numbers throughout both section types; however, they were enriched in the grafts, and radial glia-like cells were only sparsely represented in both ventral midbrain and grafted striatal sections.
The largest transplant also allowed us to look at the spatial distribution of different cell types. The dopaminergic neurons tended to be detected closer to the transplant's edges, and VLMCs and astrocytes toward the center of the transplant. This is in line with IHC analyses (Figure 5), which show the same trend in the spatial distribution of the graft. Here, the spatial division of cell types becomes even more evident, with TH-positive neurons being positioned at the edges of a graft staining positively for human nuclear antigen and VLMCs identified with COL1A1 concentrating toward the center.

DISCUSSION
Here, we present a spatial analysis of the composition of hESCderived dopaminergic transplants grafted in the striatum of a rat model of PD. We also show that despite the number of novel batch correction tools for scRNA-seq, 51-53,61 not all perform equally well on ST. In our dataset, only Seurat has demonstrated reliable performance in most metrics, supporting it as a suitable batch correction algorithm for the exponentially increasing number of ST datasets.
We show that the genes differentially expressed between the transplant and the host striatum or the host tissue in general, featured many dopaminergic markers. These were also observed in abundance in SNpc and did not appear in other morphology-associated clusters, fitting with the expectation. The transplanted DA neurons appear to robustly express all phenotypic markers we would consider to be compatible with a mature dopaminergic signature from a transcriptomic angle. This does not necessarily mean that the neuron has reached full maturity. This requires axonal outgrowth, synaptic www.moleculartherapy.org interaction, and the adoption of mature membrane potential and firing patterns, all of which may or may not be reflected in the transcriptomic profile. For human DA neurons in vivo, the graft can take up to four years to reach full maturity and function. 7 This is impossible to achieve in the rodent brain, but 15 month survival is one of the longest observed to date.
Marker genes for astrocytes and other glial cells were not found to be differentially expressed in the transplant relative to the surrounding rat tissue. However, as these cell types were also represented in host forebrain tissue, their absence is not surprising. This finding supports the uniqueness of the dopaminergic population in the grafted tissue.
Further investigation of the spatial gene distribution reveals the presence of the discovered dopaminergic markers, mainly at the edges of the largest graft. This is in accordance with the IHC analysis, which shows the enrichment of TH-positive cells at the edges of the transplants either due to selective differentiation, migration, or survival. 32,33,62 Relative to dopaminergic markers, VLMC, and astrocyte-associated genes tended to be expressed higher toward the center of the transplants.
Since the advent of deconvolution algorithms, many new software packages have been released. Here we show the difference between  stereoscope and cell2location. The results are similar to those previously reported, 61,63 where in most scenarios, cell2location results are more accurate on datasets with known compositions.
We show that ST introduces less bias in cell type composition compared to scRNA-seq. 35,36 This is most likely because there is no tissue dissociation before cDNA synthesis, leading to an equal likelihood of capture for most cell types. This is especially obvious when looking at the vulnerable cell types, such as neurons, which become dwarfed by the glial cell numbers in scRNA-seq. The presented method has a 100 mm resolution, which reports on a larger cell community. Through deconvolution with non-spatial single-cell datasets, the cellular composition of the features and, therefore, the tissue at large can be estimated. The remaining limitations are mainly the uncertainty regarding the cell's position within the feature and the dependency on an accurate and comprehensive single-cell reference dataset. The validity of cell type quantification, generated by deconvolution of the ST data, was confirmed using immunohistochemical analysis, which remains the gold standard for cell quantification but is limited in its multiplexing capacity. With the emergence of higher resolution implementations of this sequencing technique (e.g., with the announced Visium HD or BGI Stereo-seq single-cell resolution platform), 64 many remaining limitations will be removed. However, the analysis pipeline presented here should still be valid and directly applicable. For subcellular resolution, microscopy-based targeted RNA hybridization techniques remain the gold standard. 65 However, they are limited in the number of genes that can be simultaneously detected per cell because of RNA crowding and thus depend on a deep a priori knowledge of the assessed tissue.
After ST deconvolution, a large fraction of the transplant is categorized as dopaminergic. The DA cell type also correlated the most with the human transcript content of features. Transplants also harbored astrocyte, VLMC, and neuronal progenitor populations. To facilitate the ST deconvolution, we performed classification at the single-cell level of the scRNA-seq dataset. 36,66 This revealed an additional cell population, the peptidergic inhibitory neurons. These were also observed in low numbers in the transplant areas after ST deconvolution ( Figure S10A).
One of the larger observed differences between the scRNA-seq Sin-gleR assignment and the ST analysis was the relative abundance of glial cells. Although they represent the vast majority of the cells recovered with scRNA-seq, they are not as enriched in the transplant area of the ST analysis. The identity of the glial cells also differs between the two analyses. SingleR analysis shows that most of the cells resemble radial glia. However, the gene expression profiles are similar between astrocytes and radial glia in the reference dataset, pointing to the difficulties in dissecting these two populations ( Figure S13). The SingleR classification also reveals a rather fluid continuum with multiple cell types overlapping and co-existing in the same UMAP cluster. We believe that this is reflective of the shared origin of the cells and that they retain many transcriptomic similarities. When assessing the AQP4 gene, specific for astrocytes in the mature but not the developing state, 67,68 we find that even though this is absent in the radial glia of the reference dataset, 59 SingleR still classifies the majority of AQP4+ cells as radial glia. This points to a limitation in the SingleR algorithm. Incompatibilities between the differentiated hESCs and the mouse reference could exacerbate the limitations of SingleR. This limitation does not appear to be shared by cell2location, 60 and thus a higher fraction of the features gets classified as astrocytes which also fits with the IHC analysis.
Assessing DA transplants in situ also allows the mapping of post-synaptic changes in gene expression. DA depletion in the parkinsonian brain as well as DA replacement therapies, give rise to numerous aberrant changes in gene expression in the GABAergic MSNs. 69 Although in situ hybridization has yielded deep insights into individual transcripts and the effect of ectopic DA transplants, 70 ST can provide a holistic assessment of post-synaptic normalization. Interestingly, we here observe a near complete normalization of PENK in the dorsolateral striatum (the principal target of A9 neurons) while still being upregulated in the medial striatum (A10 DA target area) ( Figure S6). This does not directly follow the proximity to the graft core, which is here placed in the central and medial striatum. Such gradient has not previously been observed in studies using rat fetal DA transplants. 69,71,72 In conclusion, this study shows the benefit ST can bring to clinical samples in the future as it deals with several important biases, which can skew our interpretation. It is much easier to dissect pathophysiological changes and brain repair when quantitative transcriptomic data are complemented by spatial information. Furthermore, as the throughput, resolution, and accuracy of these methods have increased exponentially from year to year, there is much that ST methods can elucidate.

Research animals
All experimental procedures performed in this study were approved by the Ethical Committee to use Laboratory Animals in the Lund-Malmö region. Adult athymic nude female rats were purchased from Envigo (Hsd:RH-Foxn1 rnu ). The rats were housed in groups of 3 or 4 in individually ventilated cages under a 12 h light/dark cycle with ad libitum access to sterile food and water. Seven animals were used in this study. Three animals with histologically verified grafts and three control animals (intact SN) were used for ST. One transplanted animal was retained for IHC. Ten ST arrays were confirmed to contain tissue containing the transplants, five from the largest TX [L], four from the intermediate TX [M], and one from the smallest TX [S], all shown in Figure S2. Seven sections from the intact SN from three animals were used as an internal control for host DA neurons, five of which are shown in Figure S4. For animal details, see Table S1.

6-OHDA lesions
Adult athymic nude female rats were lesioned using a unilateral infusion of 6-OHDA, as previously described. 73 In brief, 3 mL of 6-OHDA (Sigma-Aldrich) at a concentration of 3.5 mg/mL (calculated from freebase weight, HBr-salt in 0.2 mg/ml ascorbic acid in 0.9% sterile saline) was injected using a 30G stainless steel cannula. The cannula was connected via polyethylene tubing to a 10 mL Hamilton syringe at a flow rate of 1 mL/min at the following coordinates (from bregma): anteroposterior (AP), À4.0 mm; mediolateral (ML), À1.3 mm; and dorsoventral (DV) (from dura), À7.0 mm, with the incisor bar set at À4.5 mm. The injection needle was left in place for an additional 3 min to allow diffusion of the toxin.

Drug-induced rotational locomotion
The 6-OHDA lesion was assessed using the drug-induced rotation tests. 74 Rats were placed in automated rotometer bowls modeled after the design of Ungerstedt et al. 75 and recorded using rotation counting software (AccuScan Instruments Inc.). The animals were recorded for 90 min following amphetamine injection (2.5 mg/kg intraperitoneal [i.p.]), and only animals rotating more than 5.5 rpm net clockwise rotations were included in the study.
Cell differentiation and grafting RC17 human pluripotent stem cells (hPSCs) (Roslin Cells, hPSCreg #RCe021-A) at passage 28 were differentiated and prepared for transplantation, according to. 37,76,77 Cultures of hPSCs were first differentiated into accurate VM progenitor patterning in N2 media with the GSK3 inhibitor CHIR99021 (0.9 mM) and Shh-C24II (300 ng/mL) to obtain the correct VM regionalization in the rostral-caudal axis and floor plate ventralization. At nine days after the start of the differentiation protocol, the cell VM culture was fine-tuned to the caudal early VM progenitors by the addition of FGF8b. At 11 days after the start of differentiation, the cells were replated and expanded in B27 media with FGF8b, AA, and BDNF and cultured for an additional 5 days in vitro (DIV) to give rise to late caudal VM progenitors with a high proportion positive for FOXA2/LMX1/OTX2. At 14 DIV, a subset of cells was prepared for mRNA extraction or ICC (validation shown in Aldrin-Kirk et al. 76 ). The cultured cells were then prepared for transplantation by dissociation into single-cell suspension in Hank's balanced salt solution (HBSS) plus DNase. The rats were engrafted with 150,000-225,000 differentiated hESCs into the striatum. 65 The cells were engrafted at a concentration of 75,000 cells/mL at the following coordinates (relative to bregma): AP, +0.5 mm; ML, À2.6 mm; and DV (from dura), À4.5 mm, with the incisor bar set at À2.4 mm. www.moleculartherapy.org

Perfusion and PFA fixation
Animals for IHC analysis were deeply anesthetized by sodium pentobarbital overdose (Apoteksbolaget) and transcardially perfused with a 50 mL physiological saline solution followed by 250 mL of freshly prepared, ice-cold, 4% paraformaldehyde (PFA) in 0.1 M phosphate buffer adjusted to pH 7.4. The brains were then removed and postfixed for 2 h in cold PFA before storing in 25% buffered sucrose for cryoprotection over at least 24 h until further processing. The remaining PFA fixed brains were cut into 35mm thick coronal sections using a freezing microtome (Leica SM2000R), collected into eight series, and stored in an anti-freeze solution (0.5 M sodium phosphate buffer, 30% glycerol, and 30% ethylene glycol) at À20 C. Furthermore, samples were stained with a solution of 1 mg/mL DAPI in PBS to stain for nuclei and washed twice with PBS. After this final step, sections were mounted on coated glass slides, covered with PVA-DABCO, and imaged on a confocal microscope as described below.

Laser scanning confocal microscopy
Confocal imaging was performed using a Leica SP8 microscope. Images were captured using a HyD detector and always with the lasers set to be activated in sequential mode to avoid serial excitation. Solid-state lasers at 405, 448, 552, and 650 nm wavelengths were used to excite their respective fluorophores. A pinhole of 1 AU was always retained during image acquisition. Leica objectives 5Â/0.15 and 20Â/0.75 were used during imaging acquisition.

Spatial Transcriptomics
All samples were processed according to the published ST protocol, 50 with all the reagents as they are listed in the protocol. In summary, sections were cut with 10 mm thickness and mounted on the active areas of the ST slides. All sections have been fixed with 4% formaldehyde solution and stained with hematoxylin and eosin. Afterward, slides were covered with 85% glycerol, and a cover glass was mounted prior to imaging. The slides were scanned on a slide scanner. Following imaging, sections were permeabilized with pepsin-HCl mixture for 8 min, as determined by the quality control slides ran prior to the main experiment, quickly followed by cDNA synthesis, which took place as described in Salmén et al. 50 overnight at 42 C. This step was followed by second-strand synthesis, after which cDNA was detached from the slide and collected as prescribed.
After the dissociation of cDNA from the slide, probes remaining on the slide were labeled by fluorescent probes and imaged to be able to correct the theoretical feature positions to the actual positions.
The further processing of the detached cDNA was done accordingly with the protocol mentioned above without changes. The indexed library was finally sequenced with NextSeq 500/550 High-Output v2 kit.
Image alignment, detection of spatial spots and tissue ST Spot Detector was deployed as previously described. 78 H&E tissue image, as well as Cy3-labeled spatial spots image, were supplied, aligned and the precise positions of spatial spots under the tissue were exported.
Sequencing, genome alignment, and data filtering All samples were sequenced with Illumina NextSeq 500 system using 75 cycle reagents from the NextSeq 500/550 High-Output v2 kit running paired-end protocol: R1 31 bp, i7 6 bp, R2 46 bp. The sequencing data has been processed with STARsolo aligner (version 2.7.9a_2021-06-25). 79 Sequences were aligned to a chimeric (combined) genome where the human (GRCh38 release 93) and rat (Rnor6 release 93) were fused into a double-sized genome. STARsolo was executed with the -solo-MultiMappers EM parameter, which uses the maximum likelihood estimation (MLE)-expectation maximization (EM) algorithm to accurately quantify reads with multi-gene alignment. Feature IDs were provided as a whitelist for cell barcode detection, and deduplication was achieved using the feature unique molecular identifiers (UMIs). The output is a matrix with per feature gene expression read counts. This results in the high fraction of reads being mapped to both human and rat gene homologs as the length of the sequence is frequently insufficient to distinguish between the species. All reads aligning to a single location were allocated to the matching species, remaining reads were proportionally distributed among the genes spanning the matching sites as previously described. 79,80 For species assessment in Figures 2, S3, and S7, genes aligned to more than one location were discarded and only uniquely mapping transcripts were kept. Full details of the STARsolo parameters can be found on the GitHub repository. 79,80 Ensembl IDs were translated to general symbols to remove the species identity from the transcripts, and duplicate entries and gene isoforms were merged. The resulting gene-feature matrices have been filtered for genes exceeding 2 UMIs in 4 features and features with more than 300 UMIs. The remaining genes and features were discarded.
Feature positions have been detected as described in Salmén et al., 50 imaged, and processed through ST Spot Detector software 81 without modifications. The resulting mask of corrected feature positions underneath the tissue was used for the final assessment.

Computational data analysis Visualization of gene expression
To visualize the presence of genes in the tissue, the tissue images were converted into masks, maximizing the contrast between the tissue and its surroundings. These tissue outlines were then transformed into a set number of features, onto which the gene expression levels were linearly interpolated.

Normalization and batch correction
To assess our baseline, sections from the same morphological areas (hereafter "related sections") were processed without any batch correction. Reads in all features were normalized and log transformed. Three thousand variable genes from all sections were selected for further processing via the Scanpy module. Nearest neighbors were computed (10 neighbors considered) for each feature, and the data were reduced into 2 UMAP dimensions. The division into clusters was computed through Leiden algorithm (resolution = 1).
Seurat was implemented as in Hao et al. 53 with SCT transform. 82 Seurat has integrated related sections on the basis of 3,000 integration features. PCs were kept to the point where further PCs did not explain more than 5% over the previous component, and the dataset was further dimensionally reduced with UMAP. Clusters were assigned shared nearest neighbors (SNNs) through a Leiden algorithm with resolution 1.
Scanorama was implemented in integration with Scanpy. 83 The data were normalized through a Scanpy preprocessing module by scaling the number of reads from each feature to a median number of transcript/features in the dataset and log transformed. Three thousand highly variable genes were selected from (scanpy.pp.highly_varia-ble_genes: flavor = "seurat"), and datasets in this form were corrected with Scanorama (correct_scanpy, default parameters). Nearest neighbors were found for each feature, the output dataset was dimensionally reduced with UMAP, and clusters were assigned by Leiden algorithm with a resolution 1. DESC was implemented as described in Li et al., 52 the reads per feature were normalized to 10,000 reads/feature as by the developers' tutorial, log-transformed, and 3,000 highly variable genes were selected. The batch effect was corrected at the level of unique tissues. DESC was trained with default parameters (dimension = [len(genes), 32,16], tolerance = 0.005, n_neighbours = 3, batch_size = 256, learning_rate = 300, tSNE_core_number = 4) and Louvain resolution of 1.

Transplant selection
Sequencing data were mapped onto human (GRCh38) and rat (Rnor6), with only a single mapping onto the two genomes kept. Counts of all genes belonging to one species within a feature were added up, and the ratio of human to rat genes was carried on into the analysis. Clusters were assigned to features by Seurat, as described previously. The trans-plant cluster was selected as the cluster with the highest average human transcript proportion in features. Wells with more than 20% of human reads average in clusters made up of host tissues were excluded from joint analyses. The transplant size annotation TX[S], TX[M], and TX [L] was based on the sum of all features belonging to the Seurat clusters with an average of above 15% human reads.

Region annotation
Seurat-assigned clusters were used for region annotation. Clusters were assigned into regions on the basis of the expression of MSN genes for striatum (PENK, ADORA2A, and PPP1R1B), dopaminergic markers for substantia nigra (TH, PBX1, SLC6A3, ALDH1A1, DDC, and RET), MBP for white matter tracts and CCK and KCNC2 for cortex. The cluster was categorized into the corresponding region if the average difference of the mean expressions for the genes was larger than 1 standard deviation.

Differential gene expression
For differential gene expression, features from related tissues were assorted into their respective clusters. Mitochondrial and ribosomal genes were removed from the analysis. All features with the same tissue of origin and cluster were merged into a pseudobulk because of reported better performance and lower false discovery rate. 84 All pseudobulk reads were then analyzed with DESeq2.
Deconvolution with a single-cell reference dataset A reference single-cell dataset for stereoscope and SingleR assessment was selected from a scRNA-seq dataset of the mouse nervous system. 59 The whole dataset was subset for corresponding areas/regions and cell types on the basis of previous knowledge. 77 To prevent low-quality cells from the reference dataset, we removed all cells with fewer than 1,000 transcripts per cell to filter out low-quality cells while preserving less transcriptionally active cell types. Cell labels were assigned based on the level of clusters and description provided by the authors.
The cells were gated at 50-500 cells per cell label for stereoscope deconvolution. In cell categories with more than 500 members, cells to be included in the reference dataset were selected by random choice.
Features were deconvolved with cell2location and stereoscope. 85 With cell2location, substantia nigra and grafted tissues were deconvolved separately, with an estimation of the reference dataset being run for 650 epochs and deconvolution of ST data for 30,000 epochs.
For stereoscope, tissues from each animal were analyzed separately. However, the same reference dataset was used for all. The algorithm was run for 75,000 epochs, as suggested by the authors.
Pearson correlation coefficient was used to correlate cell type content to the human content of a feature.
Previous scRNA-seq dataset of the hESC transplant Dataset was adapted from Tiklová et al. 36 and analyzed per the authors' instructions. For correlation-based analysis, only clusters www.moleculartherapy.org used for cell2location deconvolution were used from the reference dataset. 59 Cell type assignment of the scRNA-seq reference dataset To assess the expected prevalence of cell types in a single-cell dataset of hESC-derived dopaminergic transplants from Tiklová et al. 36 , the processed scRNA-seq data has been normalized with SCTransform, and the effect of samples has been regressed out. The data and the output have been analyzed with SingleR, 66 using the same reference dataset 59 as cell2location. The reads in this reference dataset have also been normalized with SCTransform with the ChipID variable regressed out.

DATA AVAILABILITY STATEMENT
Sequencing files are available at Sequence Read Archive under the submission number SRA: SUB12409880 with a BioProject identifier SRA: PRJNA914308. The code is available at https://github.com/ jana-rajova/ST_hESC_DA_analysis.