Implications of stress-induced gene expression for hematopoietic stem cell aging studies – Nature Aging

Mice

Young (2–6 months) and aged (20 months) C57BL/6NTac (strain no. B6, Taconic Bioscience), C57BL/6.SJL (bred in house) and Nr4a1-GFP (strain no. 016617, The Jackson Laboratory) male and female mice were used. Mice were housed in a controlled environment with 12-h light–dark cycles with chow and water provided ad libitum. All animal procedures were approved by the Lund University Ethics Committee (no. 16468-20).

Reagents

All commercially available reagents used in this study are listed in Supplementary Table 5.

Flow cytometry analyses and cell sorting

BM cells were isolated from the tibia, femur and pelvis into ice-cold FACS buffer (2% FCS/PBS) with or without 5 µM TP (Tocris)²⁸ and filtered through 70-µm cell strainers. Cells were cKIT-enriched by anti-cKIT-APC staining, followed by incubation with anti-APC MicroBeads and magnetic separation on LS columns (Miltenyi Biotec). Aliquots of cells were resuspended in culture medium (DMEM with high glucose, 10 mM HEPES, 2% FCS) ± TP (0.1, 0.5, 1 and 5 µM) and incubated on ice or at 37 °C for 90 min. After incubation, cells were stained with biotinylated antibodies against B220, Gr-1, TER119, CD3, NK1.1 and Sca-1 Pacific Blue (Sony Biotechnology) for 30 min on ice in the dark. For scRNA-seq, cells were additionally stained with oligo-conjugated hashing (HTO) antibodies (TotalSeq-A0301, -A0302, -A0303, and -A0304, BioLegend). Secondary staining was performed with streptavidin-Brilliant Violet 605 (Sony Biotechnology). After staining, LSK HPSCs were FACS-purified on a BD FACS Aria III instrument (BD Biosciences) and subjected to single-cell profiling (10x Genomics) or analyzed on an LSRFortessa X20 analyzer with the FACS Diva software v.9.0 (BD Biosciences). For PB analysis using flow cytometry, blood samples were sedimented with 1% Dextran T500 (Sigma-Aldrich) for 30 min at 37 °C; the remaining erythrocytes were lysed using ammonium chloride solution (STEMCELL Technologies) for 3 min at room temperature. Cells were stained with CD19-PE-Cy7, TER119-PerCP-Cy5.5, CD11B-APC, NK1.1-Pacific Blue, CD3-Alexa Fluor 700, CD45.1-Brilliant Violet 650 and CD45.2-Brilliant Violet 785 (Sony Biotechnology). For the analysis of in vivo-fixed HSPCs, unfractionated BM cells were stained with PE-Cy5-conjugated B220, Gr-1, TER119, CD3, NK1.1, Sca-1 Pacific Blue, CD48-FITC, CD150-PE-Cy7, CD135-PE (Sony Biotechnology) and cKIT-APC-eFluor 780, CD201-APC (eBioscience) antibodies. Before analysis or sorting, cells were stained with propidium iodide (1:1,000, Invitrogen) to exclude dead cells. Data analysis was done using FlowJo v.10.5.3 (FlowJo LLC).

Transplantation experiment

cKIT-enriched BM cells derived from C57BL/6.SJL-CD45.1 donor mice were isolated into ice-cold FACS buffer and incubated on ice or at 37 °C for 90 min. Competitor cKIT-enriched BM cells derived from C57BL/6N-CD45.2 mice were kept on ice throughout the whole experiment. After incubation, donor and competitor cells were mixed at a 1:4 ratio and transplanted into lethally irradiated (9.5 Gy) C57BL/6N-CD45.2 recipient mice via intravenous injections. Recipient mice received antibiotic prophylaxis (ciprofloxacin-supplemented water, 125 mg l⁻¹, Krka) for 2 weeks beginning on the day of irradiation. Donor-derived reconstitution was monitored in PB using flow cytometry.

Analysis of in vivo-fixed HSPCs

Young (2 months old) C57BL/6 female mice were anesthetized with ketamine (MDS Animal Health)-xylazine (Bayer Animal Health) (10% KX, 10 µl g⁻¹); in vivo fixation was performed using transcardiac perfusion with PBS followed by 4% paraformaldehyde (PFA). BM cells were isolated from the tibia, femur and pelvis into ice-cold FACS buffer (2% FCS/PBS), filtered through 70-µm cell strainers and processed as described in the ‘Flow cytometry analyses and cell sorting’ section.

Quantitative PCR with reverse transcription

cKIT-enriched BM cells were isolated from young (2 months) and aging (20 months) C57BL/6 male mice and incubated at 37 °C for up to 16 h in DMEM High Glucose (Gibco) supplemented with 10% FCS, 0.1 mM 2-mercaptoethanol (Invitrogen), 1× penicillin-streptomycin-glutamine (Invitrogen), 50 ng ml⁻¹ stem cell factor, 10 ng ml⁻¹ TPO and 10 ng ml⁻¹ Flt3L (all from PeproTech). After incubation, RNA was extracted using the Direct-zol RNA MicroPrep Kit with on-column DNase I treatment (Zymo Research). Reverse transcription was carried out using the qScript cDNA Supermix according to the manufacturer’s instruction (Quantabio). The resulting complementary DNA (cDNA) was used for quantitative PCR with reverse transcription (RT–qPCR) with SsoAdvanced SYBR Green Supermix (BioRad Laboratories) and the following primers: Actb: 5′-CCACAGCTGAGAGGGAAATC-3′ (forward), 5′-CTTCTCCAGGGAGGAAGAGG-3′ (reverse); Dusp1: 5′-ACCATCTGCCTTGCTTACCTC-3′ (forward), 5′-CTCCGCCTCTGCTTCACAAA-3′ (reverse); Egr1: 5′-CCTATGAGCACCTGACCACA-3′ (forward), 5′-GAAGCGGCCAGTATAGGTGA-3′ (reverse); Fos: 5′-ATGGGCTCTCCTGTCAACAC-3′ (forward), 5′-GCTGTAACCGTGGGGATAA-3′ (reverse); Jun: 5′-GGAAACGACCTTCTACGACGAT-3′ (forward), 5′-GGGTTACTGTAGCCGTAGGC-3′ (reverse); Nr4a1: 5′-TTGAGTTCGGCAAGCCTACC-3′ (forward), 5′-GTGTACCCGTCCATGAAGGTG-3′ (reverse); Zfp36: 5′-CCCTCACCTACTTCGCCTAC-3′ (forward), 5′-ACTTGTGGCAGAGTTCCGTTT-3′ (reverse). Analyses were performed in two independent experiments with n = 2 biological replicates per group.

Analysis of RNA integrity from fixed HSPCs

cKIT-enriched BM cells were isolated from young (2 months old) C57BL/6.SJL female mice and split into four groups (approximately 900,000 cells per sample). Untreated cells were kept on ice throughout the whole procedure. Cells from remaining samples were spun down, resuspended in 2% PFA in PBS and incubated at room temperature for 15 min. Cells were spun down at 400g for 5 min after adding wash buffer (PBS supplemented with 1% BSA and 40 U ml⁻¹ SUPERase•In, Invitrogen) and were subsequently resuspended in the wash buffer. For reverse crosslinking, cells were incubated either at room temperature for 5 min in wash buffer supplemented with 125 mM glycine (Serva) or at 56 °C for 60 min in wash buffer supplemented with 40 U ml⁻¹ proteinase K (Roche), followed by incubation on ice for at least 5 min. RNA was isolated using the Direct-zol RNA MicroPrep Kit with on-column DNase I treatment (Zymo Research). RNA integrity was measured using the RNA 6000 Pico Kit and a 2100 Bioanalyzer instrument (Agilent Technologies).

scRNA-seq and bioinformatics analysis

scRNA-seq was performed in a single batch. BM cells were isolated from a pool of young (2–3-month-old) C57BL.6/SJL female mice (n = 5), split into four groups and processed simultaneously. cDNA libraries were generated using the Chromium Next GEM Single Cell 3′ Reagent Kit v.3.1 (10x Genomics). Briefly, LSK cells (40,000 cells for each group) were sorted into ice-cold FACS buffer, resuspended at 1,250 cells per microliter in FACS buffer and 54,000 cells were loaded onto the Chromium Controller (10x Genomics). The cDNA sequencing library was generated according to the manufacturer’s instructions (10x Genomics). The HTO-derived sequencing library was prepared according to CITE-seq_and_Hashing_protocol_190213 (https://cite-seq.com/protocols). Final cDNA-derived and HTO-derived libraries were pooled and sequenced on an Illumina NovaSeq 6000 instrument using the S2 Reagent Kit v.1.5 (100 cycles). Preparation of libraries and sequencing were done at the Center for Translational Genomics at Lund University.

After sequencing libraries were processed using Cell Ranger v.7.0.0. The FASTQ files were aligned to the mouse reference genome (mm10) to create unique molecular identifier count tables of gene expression. HTO demultiplexing was performed with cellhashr (https://github.com/BimberLab/cellhashR). Quality control and downstream analyses were performed using Seurat²⁹ v.4.1.1. The filtering thresholds used for cell barcode exclusion were as follows: a high number of mitochondrial transcripts (>6%); low-quality or empty droplets (the number of detected genes below 2,000 or above 9,000 and the number of detected transcripts below 5,000 or above 60,000); and the number of HTO transcripts below 20 or above 400. Motivated by the simultaneous processing of all cells and the near complete overlap between two of the sample groups (the two groups where cells were kept on ice; Fig. 2e), no additional sample integration was performed. Identification of highly variable genes, data normalization and dimensionality reduction were performed using Seurat²⁹ with default parameters.

Cell cycle analysis was carried out using the CellCycleScoring function from Seurat. For the analysis of gene signatures, the aggregated expression of the respective genes was calculated using the AddModuleScore function from Seurat; cells were classified based on the aggregated expression using a modified version of the CellCycleScoring function that output only two classifications. The aggregated values for genes were calculated using a bin normalization approach for each individual gene and visualized on the UMAP embeddings through the FeaturePlot function.

DEGs were analyzed using Seurat’s FindMarkers function with default parameters and testing performed only on highly variable genes (two-sided Wilcoxon rank-sum test). DEG testing between cells incubated on ice with or without TP yielded three DEGs with nonsignificant P values; therefore, these samples were treated as one (‘ice’) for subsequent comparisons. As the dataset showed clear separation based on cell cycle phase classification, the samples were first compared cell cycle phase-wise. The DEGs between 37 °C with TP versus ice, 37 °C without TP versus ice across three cell cycle phases (G₁/S/G₂-M) were next used for MSigDB pathway analysis using Enrichr (https://maayanlab.cloud/Enrichr/).

All analysis steps except for the BCL to FASTQ conversion and the Cell Ranger run can be found in a Snakemake pipeline in the NCBI repository, with included conda (https://zenodo.org/record/4774217) environment specifications and version numbers of all packages used.

Analysis of published microarray and bulk RNA-seq data

The analysis was performed using R v.4.1.0, v.4.1.3 and v.4.3.2. For the microarray data, raw.CEL files were retrieved from the Gene Expression Omnibus (GEO) and RMA-normalized using the affy v.1.76.0 or oligo v.1.62.2 packages. For bulk RNA-seq data analysis, count tables were downloaded from the GEO and data were normalized in R using the DESeq2 (ref. ³⁰) v.1.38.3. When raw gene counts were not available, the normalized count tables provided by the authors were used. Genes with average expression counts less than the total number of samples were filtered out. GSEA³¹ was performed using normalized data of the input files and the following gene sets: MSigDB Hallmark provided by the software⁵ (MSigDB v.2022.1.Mm); custom-generated gene sets for IER and cell cycle; and aging and young signatures retrieved from Flohr Svendsen et al.². All analyses were performed using the GSEA software v.4.3.2 (https://www.gsea-msigdb.org/gsea/index.jsp). The false discovery rate (FDR) method was applied to correct multiple hypotheses testing. Mining of gene lists was performed using the online tool Enrichr³² (https://maayanlab.cloud/Enrichr/) with default parameters. Graphs were generated using Prism v.9.5.1 (GraphPad Software).

For the analysis of stress-associated genes from van Velthoven et al.¹⁹ and Wu et al.²⁰, data were downloaded from the GEO or using the supplementary tables provided by the authors. The DEGs with an FDR < 0.1 and fold change > 2 were used for comparison between IERsig genes identified in this study. The odds ratios (ORs) were computed in R by calculating the hypergeometric distribution, where the number of IERsig genes was 78 and the total number of genes was 19,570 (defined based on the total number of annotated genes on Affymetrix 430.2 microarrays from which the IERsig gene list was curated). The calculations of P values for the hypergeometric distributions were performed in R using a Fisher’s exact test.

For the analysis of the IERsig in the reference datasets on murine HSPCs, the R packages dplyr v.1.1.3 and ggplot2 v.3.4.4 were used. Briefly, the fraction of the IERsig in each sample was calculated by dividing the sum of gene expression levels for the IERsig genes by the total sum of gene expression levels in each sample. Data were visualized as dot plots with the fraction of IERsig across samples. To assess how the IERsig compared to randomly sampled gene sets, bootstrapping was used with a total of 10,000 iterations for each dataset. For reproducibility, a random seed (123) was set. In each iteration, a set of genes (RNA-seq data) or probe sets (microarrays), corresponding to the size of the IERsig, were randomly sampled from the gene expression data. The fraction of the bootstrapped genes in each sample was calculated in the same manner as for the IERsig. Descriptive statistics were computed for both the bootstrapped results and the IER signature fractions. To assess how the IER signature compared to the bootstrapped values, the centile rank of the coefficient of variation (CV) of the IER signature within the distribution of bootstrap CVs was computed and converted to a percentage. To visualize the results, histograms were generated to display the distribution of bootstrap CVs and the CV of the IER signature.

Statistics and reproducibility

Data were analyzed using R v.4.1.0, v.4.1.3 and v.4.3.2, Microsoft Excel v.16.54 and Prism v.9.5.1 (GraphPad Software). All experiments were repeated as indicated; n indicates the number of independent biological replicates. No statistical methods were used to predetermine sample sizes but our sample sizes are similar to those reported in previous publications^3,4. For the comparisons presented in Fig. 2d,h and Extended Data Fig. 5a, normality and equal variances were formally tested and data met the assumptions of the statistical tests used. Data collection and analysis were not performed blinded to the conditions of the experiments. For the transplantation and in vivo fixation experiments, no randomization method was applied to allocate animals to the experimental groups. Animals (n = 3–5) in the same cage received the same treatment. For the analyses of Nr4a1-GFP transgenic mice and the RT–qPCR experiments, no randomization method was applied and mice per sample were assigned to the experimental groups based on genotype or age. No animals or data were excluded from the analyses. For the flow cytometry analyses and RT–qPCR experiments, the analyses were performed in two or three independent experiments; no inconsistent results were observed. All results were successfully reproduced. The specific statistical test used for each experiment is indicated in the corresponding figure legend.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

• SEO Powered Content & PR Distribution. Get Amplified Today.
• PlatoData.Network Vertical Generative Ai. Empower Yourself. Access Here.
• PlatoAiStream. Web3 Intelligence. Knowledge Amplified. Access Here.
• PlatoESG. Carbon, CleanTech, Energy, Environment, Solar, Waste Management. Access Here.
• PlatoHealth. Biotech and Clinical Trials Intelligence. Access Here.
• Source: https://www.nature.com/articles/s43587-023-00558-z