In this study, we aim to employ snRNA-seq to elucidate the cellular and molecular responses of rat kidney to CIH. By characterizing the gene expression patterns of individual kidney cell types under CIH conditions, we seek to identify key pathways and cell populations involved in OSA-associated renal damage. Our findings will provide critical insights into the mechanisms of CIH-induced renal damage and may reveal potential therapeutic targets for mitigating kidney complications in OSA patients.

12 male Sprague-Dawley rats, aged 8 weeks, were obtained from Fujian Medical University's Laboratory Animal Center. Following one week of acclimatization under standard housing conditions (12:12 h light-dark cycle) with ad libitum access to food and water, the animals were randomly assigned to two experimental groups using a computer-generated randomization scheme: normoxic control group (NC, n = 6) and CIH group (CIH, n = 6). This study was approved by the Animal Ethics Committee of Fujian Medical University (IACUC FJMU 2022-0031). All methods were carried out in accordance with relevant guidelines and regulations. All methods were also reported in accordance with ARRIVE guidelines.

Rats in the CIH group were exposed to IH for 12 weeks. The IH protocol was established based on previous study. Each IH cycle consisted of four phases: (1) reducing oxygen concentration from 21% to 6% over 120 s, (2) maintaining hypoxia (6% O₂) for 30 s, (3) rapidly reoxygenating to 21% over the next 10 s, and (4) sustaining normoxia (21% O₂) for an additional 20 s. The CIH group underwent daily 8-hour IH exposure (08:00-16:00) during the light phase, while NC rats were maintained in room air under identical temporal conditions.

All kidney specimens were fixed in 10% neutral buffered formalin, followed by dehydration and paraffin embedding. Sections were hematoxylin-stained (20 min) followed by eosin-phloxine counterstaining. Furthermore, sections were mounted in xylene-based medium. Histopathological examination was performed under an Olympus BX50 light microscope (Tokyo, Japan) and evaluated independently by a single board-certified pathologist blinded to the experimental groups to ensure unbiased assessment. Renal injury was assessed using the modified Jablonski scale (0-5): 0 indicates normal renal cells; 1 represents degeneration and necrosis of individual cells; 2 denotes degeneration and necrosis affecting a single kidney tubule; 3 corresponds to degeneration and necrosis adjacent to the proximal convoluted tubule, with surviving cells surrounding the affected tubule; 4 signifies necrosis confined to one-third of the distal convoluted tubule, with a necrotic band extending into the inner cortex; and 5 indicates necrosis involving all segments of the proximal convoluted tubule.

Three renal tissue samples were randomly selected from each group using a computer-generated randomization scheme and subsequently subjected to sequencing. Renal tissues were harvested and washed in pre-cooled PBSE (PBS buffer containing 2 mM EGTA). Nuclei isolation was carried out using GEXSCOPE Nucleus Separation Solution (Singleron Biotechnologies, Nanjing, China), according to the manufacturer's product manual. Isolated nuclei were resuspended in PBSE to 10 nuclei per 400 µl, filtered through a 40 μm cell strainer, and counted with Trypan blue. To assess nuclei quality, the enriched nuclei were stained with DAPI (1:1,000) (Thermo Fisher Scientific). We used DAPI fluorescence to assess nuclear integrity, ensuring it was greater than 80%. We also confirmed that non-nuclear debris was less than 20% and the nuclei count for each sample exceeded 50,000 before proceeding. Nuclei were considered as DAPI-positive singlets.

The concentration of single nucleus suspension was adjusted to 3-4 × 10 nuclei/mL in PBS. Single nucleus suspension was then loaded onto a microfluidic chip (GEXSCOPE snRNA-seq Kit, Singleron Biotechnologies) and snRNA-seq libraries were constructed according to the manufacturer's instructions (Singleron Biotechnologies). The resulting snRNA-seq libraries were sequenced on an Illumina novaseq 6000 instrument with 150 bp paired end reads.

Raw reads were processed to generate gene expression profiles using CeleScope v2.0.7 (Singleron Biotechnologies) with default parameters. Briefly, Barcodes and UMIs were extracted from R1 reads and corrected. Adapter sequences and poly A tails were trimmed from R2 reads and the trimmed R2 reads were aligned against the mRatBN7.2 (rattus_enemebl_111) transcriptome using STARSolo in STAR (v2.7.11a). Key STARsolo parameters included --outFilterMatchNmin 50, while --outFilterMismatchNmax was left at its default value of 10. Successfully assigned reads with the same cell barcode, UMI and gene were grouped together to generate the gene expression matrix for further analysis.

Scanpy v1.8.2 was used for quality control, dimensionality reduction and clustering under Python 3.7. For each sample dataset, we filtered expression matrix by the following criteria: (1) cells with gene count less than 200 or with top 2% gene count were excluded; (2) cells with top 2% UMI count were excluded; (3) cells with mitochondrial content > 10% were excluded; (4) genes expressed in less than 5 cells were excluded. After filtering, 42,581 cells were retained for the downstream analyses, with on average 766 genes and 1,157 UMIs per cell. The raw count matrix was normalized by total counts per cell and logarithmically transformed into normalized data matrix. Top 2000 variable genes were selected by setting flavor = 'seurat'. Principle component analysis (PCA) was performed on the scaled variable gene matrix, and top 20 principle components were used for clustering and dimensional reduction. We used the Harmony algorithm to integrate the data and effectively remove batch effects. Cells were separated into 24 clusters by using Louvain algorithm and setting resolution parameter at 1.2. Cell clusters were visualized by using Uniform Manifold Approximation and Projection (UMAP). The final UMAP visualization was performed on the batch-corrected dataset, ensuring that any clustering observed reflects genuine biological differences.

To identify DEGs, we used the Seurat FindMarkers function based on Wilcox likelihood-ratio test with default parameters, and selected the genes expressed in more than 10% of the cells in a cluster and with an average log (Fold Change) value greater than 0.25 as DEGs. Benjamini-Hochberg correction was applied to control the false discovery rate (FDR), and genes with adjusted p-values < 0.05 were considered differentially expressed. For the cell type annotation of each cluster, we combined the expression of canonical markers found in the DEGs with knowledge from literatures, and displayed the expression of markers of each cell type with heatmaps/dot plots/violin plots that were generated with Seurat DoHeatmap/DotPlot/Vlnplot function. Doublet cells were identified as expressing markers for different cell types, and removed manually.

The cell type identity of each cluster was determined with the expression of canonical markers found in the DEGs using SynEcoSys database. Heatmaps/dot plots/violin plots displaying the expression of markers used to identify each cell type were generated by Seurat v3.1.2 DoHeatmap/DotPlot/Vlnplot.

To investigate the potential functions of DEGs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used with the "clusterProfiler" R package 3.16.1. Pathways with p_adj value less than 0.05 were considered as significantly enriched. GO gene sets including molecular function (MF), biological process (BP), and cellular component (CC) categories were used as reference.

Data are represented as the mean ± standard deviation. The Student's t test was used to compare the two groups. The statistical significance was considered when p-values less than 0.05. GraphPad Prism version 8.0 was used for the statistical analysis (GraphPad Software Inc., San Diego, CA, United States).