Single-cell RNA sequencing reveals the role of immune-related autophagy in spinal cord injury in rats

Spinal cord injury refers to damage to the spinal cord due to trauma, disease, or degeneration; and the number of new cases is increasing yearly. Significant cellular changes are known to occur in the area of spinal cord injury. However, changes in cellular composition, trajectory of cell development, and intercellular communication in the injured area remain unclear. Here, we used single-cell RNA sequencing to evaluate almost all the cell types that constitute the site of spinal cord injury in rats. In addition to mapping the cells of the injured area, we screened the expression of immune autophagy-related factors in cells and identified signaling pathways by the measuring the expression of the receptor−ligand pairs to regulate specific cell interactions during autophagy after spinal cord injury. Our data set is a valuable resource that provides new insights into the pathobiology of spinal cord injury and other traumatic diseases of the central nervous system.

The proportion of mitochondrial genes to all genetic material may indicate whether the cell is in a steady state. We generally think that when the proportion of mitochondrial genes in a cell is higher than all genes, it may be in a state of stress. As a result, we filtered cells with mitochondrial gene content more than 10%.
We also filter the twins according to the doublet ratio of the 10x genomics platform according to the standard of 0.077 * sample cell number, and filter outliers by 2 * (mean + /-standard deviation), that is, low-quality cells that specify the number of genes and the number of UMI. After the above steps, we obtained 56287 cells.
Standardize the data using the LogNormalize method. After controlling the relationship between average expression and dispersion, highly variable genes were identified in a single cell. Then, we use principal component analysis (Principal Component Analysis, PCA), using variable genes as input, and identify significant principal components according to the ElbowPlot function. The first 15 principal components were selected as the statistically significant input of t-Distributed Stochastic Neighbor Embedding (t-SNE).

Cell clustering, cell marker, and cell type identification
The FindClusters function was used to cluster the cells, and the resolution parameter was set to 0.4, which is divided into 22 clusters. We further looked for the differences between cell clusters, marker. In order to verify our cell type annotation, we first replaced the homologous genes of rats and mice and annotated the standard cell types by using the MouseRNAseqData data set in SingleR (version 1.8.1) R package (Aran et al., 2019) and related literature (Milich et al., 2021).

Single cell subgroup analysis and quasi-sequential analysis
Subset function was used to extract immune cell groups (macrophages/microglia, monocytes, neutrophils, T_NK cells, B cells), macrophages/microglia, T_NK cells for subpopulation analysis. The differences of macrophage/microglia subsets between mild injury group and control group, moderate injury group and control group, severe injury group and control group were analyzed, and the significant difference gene screening threshold FoldChange was set as 1pvalue less than 0.05.
Monocle (version 2.22.0) R packet (Qiu et al., 2017) was used to infer the cell differentiation of the above three subsets. The integrated gene expression matrix from each cell subgroup is first derived from the Seurat object to Monocle to construct the cell dataset. Use the variable gene defined by the dispersionTable function, and then use the setOrderingFilter function to sort the cells.
The DDRTree method is used to reduce the dimension and the orderCells function is used to estimate the arrangement of cells along the trajectory. Based on the clustering characteristics and marker gene analysis, the trajectory map of the differentiation time of immune cells in single cell data set was obtained. The analysis of each track adopts a standard scheme with default parameters.

GO and KEGG enrichment analysis and Gene set variation analysis
Gene Ontology (GO)enrichment analysis is a common method for large-scale functional enrichment of genes in different dimensions and different levels. Fisher accurate test is used to calculate the enrichment significance of each Term in biological processes (BP), molecular function (MF) and cellular components (CC), and count the number of differential genes included in each GO entry. The hypergeometric distribution algorithm (Mi et al., 2017) was used to calculate the significance of differential gene enrichment in each GO entry. The calculation results will return a P value of rich significance, and the smaller this value, the more significant it is statistically.
Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis uses hypergeometric distribution to calculate the correlation between each pathway in KEGG pathway and this differential gene. The smaller the P-value, the correlation degree between the corresponding pathway and this differential gene.
Gene set variation analysis (GSVA) Analysis of immune Cell Group in single Cell data set based on R language GSVA package (Hänzelmann et al., 2013). GSVA analysis carried out rank statistics calculation similar to K-S test for the gene set corresponding to each pathway. Then the pathways with significant differences were obtained by limma package (Ritchie et al., 2015), so as to evaluate the enrichment degree of different pathways among different groups. Statistical significance was defined as P < 0.05.

Analysis of intercellular communication
CellChat (version 1.1.3) R packeage (Jin et al., 2021) is often used to analyze intercellular communication networks from single-cell rna sequencing data. CellChat was used to quantitatively infer and analyze the intercellular communication network from single-cell rna sequencing data, circle map was used to show the interaction of cell populations, and the bubble diagram was used to count all the important ligand pairs in intercellular signal transmission.

Analysis of cell score and Regulon regulator based on immune-related autophagy factors
We focus on autophagy genes related to immunity, so we download immune genes and autophagy genes from ImmPort database (Bhattacharya et al., 2014) (https://www.immport.org/) and HADb human autophagy database (http://www.autophagy.lu/index.html) for follow-up analysis. Based on our research species. Firstly, rat homologous genes were screened from the list of human immune genes and human autophagy genes by homologene (version 1.1.3) R package (https://oganm.github.io/homologene/). The differential immune genes were obtained in the difference analysis of immune cell subsets, and then through the intersection of these differential immune genes and autophagy-related homologous genes, the immunerelated autophagic factors (IRAFs) of rats were finally obtained.
Based on the above IRAFs dataset, we use the AddModuleScore function to calculate the score of a single cell in this gene set. Then all the cells were divided into high and low group cell groups according to the median score, and finally the diseased sample cell groups in the high score cell group were selected for follow-up SCENIC analysis.
Multicellular organisms contain a variety of cell types, each with its own morphology and function. Cell types are generally maintained by the coordination and interaction between transcription factors and their corresponding target genes. We can identify the co-expressed modules (regulon) between transcription factors (TFs) and potential target genes and the regulon activity scores of each cell (regulon activity score, RAS) through SCENIC (version 1.1.2) (Aibar et al., 2017) software.

Analysis of differences between groups based on RNA-seq data
We download the Counts matrix of the RNAseq dataset, convert Ensembl ID to GeneSymbol, and organize the matrix file. DESeq2 (version 1.34.0) R packet (Love et al., 2014) was used to analyze the difference between the uninjured group and the injured group. The significant difference gene screening threshold FoldChange was set as 1m p value less than 0.05p value. PCA principal component analysis is used to visualize the overall distribution of samples.

Transmission Electron Microscopy (TEM)
After the fresh tissue determines the material selection part, cut the required observation sample into small pieces of about 1mm³. Fix in 2.5% glutaraldehyde at 4℃ for 2-4h. The tissue block was rinsed three times with 0.1M phosphoric acid rinse solution for 15min each time, and then the sample was placed in a refrigerator with 1% osmic acid at 4℃ and fixed for 2h. The sample is added into 30%-50%-70%-80%-95%-100%-100% alcohol in turn for gradient dehydration, 15min each time, 100% propylene oxide for 3 times, 5 min each time. Propylene oxide + embedding solution (2: 1) at room temperature for 2h, Propylene oxide + embedding solution (1: 2) at room temperature for 3h, and pure embedding solution at room temperature overnight. Replace the pure embedding solution at room temperature for 3-4h. Then pick out the sample with toothpick and embed it on the embedding plate, Placing the embedding plate in an oven at 60℃ for 48h, and taking out the embedding block for later use after the resin is completely polymerized. After rough trimming, the resin block is semi-thin sliced in an ultra-thin slicer, and after trimming according to the positioning condition, the required position is ultra-thin sliced with a slice thickness of 70nm, and the copper net is used for fishing. Dyeing with 3% uranium acetate saturated alcohol solution for 8 min; 70% alcohol for 3 times and ultrapure water for 3 times; Staining with 2.7% lead citrate solution for 8 min; Clean with ultrapure water for 3 times, and slightly dry the filter paper.

Immunostaining
Immunofluorescence staining have been detailed in previous publications (Yan et al., 2022b). Briefly, The sections were incubated, rinsed, antigen-retrieval, permeabilized, blocked, and then incubated with primary antibody overnight at 4 ℃ in a wet box. The following primary antibodies were used: rabbit anti-CD68  table   Table S1, Summary of cell marker gene