Preparation of ARM animal models and embryo embedding

Wistar rats, aged 10–12 weeks and weighing 200–230 g (provided by Beijing HFK Bioscience Co., Ltd.), were raised in a specific pathogen-free (SPF) facility. Sexually mature and healthy female Wistar rats were selected, and they were co-housed with the male rats at a 3:1 ratio overnight. The presence of sperm in the vaginal smears of female rats on the following morning was recorded as GD0. Pregnant female rats were housed separately and randomly assigned to either the ARM group or the normal group. On GD10, pregnant rats in the ARM group received an oral gavage of 1% ETU (Sigma-Aldrich; Merck Millipore, Darmstadt, Germany) at a dose of 125 mg/kg, while pregnant rats in the normal group received an equivalent volume of normal saline. Anesthesia was induced using isoflurane inhalation followed by an intraperitoneal injection of lidocaine on GDs 14–16, and cesarean sections were performed to collect the embryos. A total of 10 pregnant rats were collected for the normal group, comprising 4 pregnant rats at N14 (45 embryos), 3 at N15 (39 embryos), and 3 at N16 (32 embryos). For the ARM group, a total of 13 pregnant rats were collected, with 6 pregnant rats at A14 (62 embryos), 4 at A15 (40 embryos), and 3 at A16 (29 embryos). The embryos were washed with pre-cooled sterile normal saline to eliminate surface blood. Sterile gauze was used to dry any remaining liquid. Subsequently, the embryos were submerged in an embedding box filled with pre-chilled optimal cutting temperature compound (OCT; SAKURA, Japan). They were oriented horizontally in a sagittal position and rapidly immersed in pre-chilled isopentane (Sigma-Aldrich) for 1 min. Once the OCT was fully solidified, the samples were stored in a -80 °C freezer.

Slide preparation

Spatial transcriptomics slides were designed with four identical capture areas, each measuring 6.5×6.5 mm, and housing 5,000 spots containing barcoded primers supplied by 10×Genomics. These primers were attached to the slide at their 5′ end and included a cleavage site, a T7 promoter region, a partial read1 Illumina handle, a unique spatial barcode specific to each spot, a unique molecular identifier (UMI), and a 20(T) VN sequence. The spots, each with a diameter of 55 μm, were organized in a regularly spaced hexagonal grid, ensuring that each spot was surrounded by six other spots, with a center-to-center distance of 110 μm. Frozen tissue sections, cut to a thickness of 10 μm using a pre-cooled cryostat, were precisely positioned onto pre-chilled Visium Tissue Optimization Slides (10×Genomics, Cat#3000394) and Visium Spatial Gene Expression Slides (10×Genomics, Cat#2000233). These prepared slides were subsequently stored at –80 °C until further use.

Tissue optimization

The spatial transcriptomics protocol was meticulously adapted to suit the tissue under investigation, following established guidelines (Berglund et al. 2018). In summary, several adjustments were made to the staining procedure. This involved the exclusion of isopropanol, reducing the incubation time for hematoxylin and bluing buffer, and augmenting the eosin concentration. Furthermore, we determined the ideal permeabilization incubation time through experimentation with the Visium Tissue Optimization Slides. To enhance tissue removal in the one-step protocol, a higher ratio of proteinase K to PKD buffer was implemented. Once the optimal conditions were ascertained, we proceeded by cutting three cryosections per sample at a thickness of 10 μm onto spatial slides and conducting immediate processing.

Fixation, staining, and imaging

The slides underwent a series of procedures: incubation at 37 °C for 1 min, fixation in paraformaldehyde for 30 min, and subsequent washing in 1×PBS. Staining involved incubation with Mayer’s hematoxylin (Dako, Agilent, Santa Clara, CA) for 4 min, bluing buffer for 30 s, and eosin (Sigma–Aldrich) diluted 1:5 in Tris-base for 30 s. RNase-free water washed the slides after each staining step. Following air-drying, mounting occurred with 85% glycerol (Merck Millipore, Burlington, MA) and coverslips (Menzel-Glaser). Bright-field images at 20× magnification were captured using the Metafer Slide Scanning platform. For pre-permeabilization, sections were incubated with collagenase and bovine serum protein at 37 ° C for 20 min and finally infiltrated with 0.1% pepsin for 7 min.

Reverse transcription (RT), spatial library preparation, and sequencing

RT was conducted as previously described (Berglund et al. 2018). Following RT, the wells were rinsed with 0.1× SSC and subjected to incubation at 56 °C with periodic shaking for 1.5 hours using a mixture of Proteinase K (QIAGEN) and PKD buffer (QIAGEN, pH 7.5) at a 1:1 ratio for tissue removal. The spatially barcoded cDNA was enzymatically released as outlined in previous descriptions. The resulting supernatants containing released cDNA were gathered and transferred to 96-well plates for the preparation of spatial transcriptomics libraries using an automated MBS 8000 system (Jemt et al. 2016). In short, the process involved the synthesis of the second-strand cDNA, followed by in vitro transcription, adaptor ligation, and a subsequent round of RT. Sequencing handles and indices were incorporated during an indexing PCR, and the finalized libraries were purified and quantified using the methods outlined in a prior study (Ståhl et al. 2016). Sequencing was performed on the Illumina NovaSeq 6000 with a sequencing depth of at least 50,000 reads per spot and 150 bp (PE150) paired-end reads (performed by Biomarker Technologies Corporation, Beijing, China).

Spot visualization and image alignment

The process of spot staining and imaging was previously detailed. In summary, fluorescently labeled probes were employed to stain primer spots through hybridization, and these spots were subsequently captured using the Metafer Slide Scanning platform. The resulting spot image, along with the earlier acquired bright-field tissue image of the corresponding region, was loaded into the web-based Spot Detector tool for spatial transcriptomics. (Wong et al. 2018). The alignment of the two images took place, and the integrated tissue recognition tool was utilized to isolate spots covered by tissue.

Weighted Gene Co-expression Network Analysis (WGCNA)

WGCNA analysis was conducted using the R packages, hdWGCNA (version 0.1.2.1; https://smorabit.github.io/hdWGCNA/index.html) and WGCNA (version 1.71; https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-9-559).

The SetupForWGCNA function was used to set the name of the hdWGCNA experiment, the MetacellsByGroups function was used to construct the metacell expression matrix of the single-cell dataset, and the NormalizeMetacells function was used to normalize the cell expression matrix. The SetDatExpr function was used to specify the expression matrix to be used for network analysis, with the default being the expression matrix of metacells. The TestSoftPowers function was used to guide the construction of the soft threshold for the co-expression network, simulating the similarity between the co-expression network and the scale-free graph under different soft power thresholds. The PlotSoftPowers function was used to visualize the result of the parameter sweep. A threshold of ≥ 0.8 was generally selected for fitting the scale-free topology model. The PlotDendrogram function was used for the visualization of the co-expression network, and the ModuleEigengenes function was used to calculate the module eigengenes in individual cells. The ME matrix had one row for each cell and one column for each module. The GetMEs function was used to extract this matrix from the Seurat object. The ModuleConnectivity function was used to calculate the kME value, and the PlotKMEs function was used for data visualization. The analysis of key module genes was based on the fact that genes in the same module usually have similar functions. The feature genes of the module were analyzed using dimensionality reduction, and the connectivity of each gene was calculated based on the feature genes to identify the highly connected hub genes in the module and to visualize them for display.

Immunofluorescence staining

Paraffin sections, providing a clear and comprehensive view of the urethra, hindgut, URS, and rectourethral fistulas, were chosen for staining. Subsequently, the sections underwent deparaffinization and antigen retrieval. Blocking was carried out with serum (ZSGB-BIO, Beijing, China) at 25 °C for 1 hour, followed by overnight incubation at 4 °C with primary antibodies. The primary antibodies used were as follows: Anti-Pcsk9 (1:100), anti-Hmgb2 (1:100), and anti-Sod1 (1:100), all sourced from Proteintech (Wuhan, China). After PBS washes, the sections were incubated at 25 °C in the dark for 1 hour with fluorescent secondary antibodies: goat anti-rabbit 488 (1:100; Proteintech). Following PBS washes, the sections were stained with DAPI solution (Servicebio Technology, Wuhan, China) for 10 minutes. Finally, the sections were mounted using an anti-fluorescence fading medium (Solarbio, Beijing, China), and images were captured using a laser-scanning confocal microscope (LSM 880, Zeiss, Oberkochen, Germany).

Pathway RespOnsive GENes (PROGENy) and non-negative matrix factorization (NNMF)

We conducted pathway activity analysis using the PROGENy package (version 1.12.0) in R software (version 4.0.3). To standardize spatial transcriptomic data, we set the PROGENy standardization parameter to FALSE and used the ScaleData function from the Seurat package for normalization. After enriching the spatial transcriptome data with PROGENy activity scores, we scaled the data to obtain spot-specific scores. These activity scores were then integrated with the clustering results, and we computed pathway-specific mean activity scores for each cluster. This dataset was used to create a heatmap with the R package heatmap. To reduce dimensionality, we employed the NNMF approach with assistance from the STutility R package. We assessed the correlation between PROGENy activity scores and dimensionality reduction outcomes using Pearson's method and visualized it with a correlation heatmap using the R package heatmap.

Whole transcriptome joint analysis

Using the R function, phyper, the P-value for ceRNA pairs was calculated using the hypergeometric distribution test. The p.adjust function (selecting the FDR method) was the applied to correct the P-values, with a filter set to remove some ceRNA pairs that did not meet the standard with hypergeometric test P-values < 0.01 and FDR values < 0.01. Co-expressed relationships were identified by calculating the correlation between two types of RNA, and positively correlated ceRNA pairs with a correlation coefficient greater than 0 were selected based on co-expression relationships. Seurat (version 4.1.0) was used to screen for differentially expressed genes in the spatial transcriptome, with the most relaxed parameter as P-value < 0.05 and FC > 1.5. The edgeR package (version 3.32.1) for the whole transcriptome was used for differential parameter modulation, with the most relaxed parameter as P-value < 0.05 and FC > 1.5. Finally, Cytoscape software was used to visualize and sort the connectivity of the ceRNA network. The data produced for the ceRNA network are available on NCBI GEO datasets (Track No. GSE159306). To access the data, please visit https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159306.