Senescence atlas reveals an aged-like inflamed niche that blunts muscle regeneration – Nature

Animal models

C57Bl/6 (WT), p16-3MR (donated by J. Campisi)¹⁸, dystrophic mdx (DBA/2-background) and mdx/p16-3MR (dystrophic mdx mice crossed with p16-3MR mice) were bred and aged at the animal facility of the Barcelona Biomedical Research Park (PRBB), housed in standard cages under 12 h–12 h light–dark cycles and fed ad libitum with a standard chow diet. All of the experiments followed the principle of the ‘three Rs’—replacement, reduction and refinement according to Directive 63/2010 and its implementation in the Member States. All of the procedures had authorization from the PRBB Animal Research Ethics Committee (PRBB-CEEA) and the local government (Generalitat de Catalunya) and were conducted according to the European Directive 2010/63/EU and Spanish regulations RD 53/2013. Both male and female mice were used in each experiment unless stated otherwise. Live colonies were maintained and genotyped according to the Jackson Laboratories’ guidelines and protocols. The mice were housed together, their health was monitored daily for sickness symptoms (not age-related weight loss and so on) and they were euthanized immediately at the clinical end point when recommended by veterinary and biological service staff members. Mice were randomly allocated to experimental or treatment groups. No blinding was used. No statistical methods were used to predetermine the sample size. For PCR genotyping, the following primers were used: p16-3MR-1, 5′-AACGCAAACGCATGATCACTG-3′; and p16-3MR-2, 5′-TCAGGGATGATGCATCTAGC-3′. Positive mice show a band at 202 bp.

Human biopsies

Human muscle biopsy samples from the vastus lateralis muscle of patients undergoing surgery were obtained from the biobank of the EU/FP7 Myoage Consortium, as previously reported^5,6. Ethical approval was received from the local ethics committees at each of the five research centres of the Consortium. All of the participants provided written informed consent and were medically screened before participation. The biobanked muscle tissue had been directly frozen in melting isopentane and stored at −80 °C until analysed. Damaged areas were identified by morphological criteria on the basis of the presence of infiltrating mononuclear cells. Data are from female patients aged 69, 82, 80, 89 or 85 years old; the average age was 81 ± 7.5 years.

In vivo treatments

Quercetin (USP, 1592409; 50 mg kg⁻¹) and dasatinib (LC Laboratories, D-3307; 5 mg kg⁻¹) were administered orally (by gavage). Control mice were administered with an equal volume of vehicle (10% ethanol, 30% polyethilenglicol and 60% phosal). GCV (Sigma-Aldrich, G2536-100MG; 25 mg kg⁻¹) was injected intraperitoneally (i.p.). Anti-CD36 antibodies (Cayman Chemical, 10009893; 10 µg or 20 µg in young mice and 20 µg in old mice) diluted in PBS was administered via i.p., and control mice received an equal dose of IgA control antibodies (Southern Biotech/Bionova, 0106-14). Treatments with GCV, senolytics and CD36 were administered daily for 4–7 consecutive days as indicated in the figure legends. NAC (Sigma-Aldrich, A9165; 0.01 g ml⁻¹) was added into drinking water (exchanged every 3 days) 1 week before muscle injury and was prolonged until euthanasia. SIS3 (Sigma-Aldrich, S0447-5MG; 10 mg kg⁻¹) in PBS and bortezomib (Teva, TH000345/LF2.3; 0.5 mg kg⁻¹) in 10% DMSO diluted in PBS were administered i.p. (from 1 to 5 d.p.i.). For long-term treatments, 3-month-old mdx and mdx/p16-3MR mice were administered with D+Q or GCV, respectively, twice a week for 2 months.

Muscle regeneration

Mice were anaesthetized with ketamine–xylazine (80 and 10 mg kg⁻¹ respectively; i.p.) or isoflurane. Regeneration of skeletal muscle was induced by intramuscular injection of CTX (Latoxan, L8102; 10 µM) as described previously⁵¹. At the indicated times after injury, the mice were euthanized and the muscles were dissected, frozen in liquid-nitrogen-cooled isopentane, and stored at −80 °C until analysis.

Heterografting

Heterografting experiments were performed according to the protocol described previously⁵². In brief, the EDL muscle was removed from the anatomical bed of either p16-3MR or WT mice and was transplanted onto the surface of the TA muscle of the p16-3MR or WT recipient mouse or vice versa. Muscle grafts were collected on day 7 after transplantation.

Muscle force measurement

Ex vivo force measurements of EDL muscles were assessed as previously described⁵³ using and 300B apparatus (Aurora Scientific). Force was normalized per muscle area, determined by dividing the muscle mass by the product of length and muscle density of (1.06 mg mm⁻³), to calculate the specific force (mN mm⁻²).

p16-3MR Renilla luciferase reporter assay

In vivo, Renilla luciferase activity was measured in the TA, quadriceps and gastrocnemius muscles of p16-3MR mice. Anaesthetized mice were injected intramuscularly with coelenterazine H (PerkinElmer, 760506) and luciferase activity was immediately measured using the IVIS Lumina III (PerkinElmer) system. In vitro, Renilla luciferase activity was measured from the cryopreserved diaphragm and TA muscles using the Dual-Luciferase Reporter Assay Kit (Promega, E1910). The signal was measured using the luminometer Centro LB 960 (Berthold Technologies) and values were normalized to the total protein extracted measured using Bradford method (Protein Assay, Bio-Rad, 500–0006), and the damaged area was measured after haematoxylin and eosin (H&E) staining.

Cell isolation by flow cytometry

Muscles were mechanically disaggregated and incubated in Dulbecco’s modified Eagle’s medium (DMEM) containing liberase (Roche, 177246) and dispase (Gibco, 17105-041) at 37 °C with agitation for 1–2 h. When required, SPiDER-β-gal reagent (Dojindo, SG02; 1 µM) was added during the second hour. The supernatant was then filtered and cells were incubated in lysis buffer (BD Pharm Lyse, 555899) for 10 min on ice, resuspended in PBS with 2.5% fetal bovine serum (FBS) and counted. BV711-conjugated anti-CD45 (BD, 563709; 1:200), APC-Cy7-conjugated anti-F4/80 (BioLegend, 123118; 1:200), PE-conjugated anti-α7-integrin (Ablab, AB10STMW215; 1:200), APC-conjugated anti-CD31 (eBioscience, 17-0311-82; 1:200) and PE-Cy7-conjugated anti-SCA1 (BioLegend, 108114; 1:200) antibodies were used to isolate MCs (CD45⁺F4/80⁺), SCs (α7-integrin⁺CD45⁻F4/80⁻CD31⁻) and FAPs (SCA1⁺CD45⁻F4/80⁻α7-integrin⁻CD31⁻). PE-Cy7-conjugated anti-CD45 antibodies (BioLegend, 103114) were used to isolate CD45-positive and CD45-negative populations (Extended Data Fig. 1h). SPiDER-β-gal (SPiDER) was used to isolate senescent cells (SPiDER⁺) from non-senescent cells (SPiDER⁻) of each cell type (Extended Data Figs. 1h and 4a) (Supplementary Table 1). Cells were sorted using the FACS Aria II (BD) system. Cell lineage was confirmed by specific-cell marker expression (Extended Data Fig. 4b–d). Isolated cells were used either for RNA extraction, cell cultures, engraftments, proliferation assays or plated onto glass slides (Thermo Fisher Scientific, 177402) for immunostaining and SA-β-gal analysis.

To isolate ROS^high and ROS^low populations, the digested muscle was stained with CellRox Green reagent (Invitrogen, C10444; 5 µM) according to the manufacturer’s protocol and PE-Cy7-conjugated anti-CD45 (BioLegend, 103114; 1:200), PE-Cy7-conjugated anti-CD31 (BioLegend, 102418; 1:200), PE-conjugated anti-α7-integrin (Ablab, AB10STMW215; 1:200) and APC-conjugated anti-SCA1 (BioLegend, 108111; 1:200) antibodies to separate SCs (α7-integrin⁺CD45⁻CD31⁻) and FAPs (SCA1⁺CD45⁻α7-integrin⁻CD31⁻) (Supplementary Table 1). CellRox^high and CellRox^low cells were sorted using the FACS Aria II (BD) system. Isolated cells were used for cell cultures and proliferation assays. The acquisition was performed using the BD FACS Diva software.

Senescent cell transplantation

Cells transplants were performed as described previously⁵, following an adapted protocol⁵⁴. FACS-isolated SPiDER⁺ and SPiDER⁻ cells were collected, resuspended in 20% FBS DMEM medium, labelled with Vybrant Dil Cell Labelling solution (Invitrogen, V22889) according to manufacturer instructions and injected into the TA muscles of recipient mice that were either uninjured or previously injured using the freeze crush method 2 days before⁵⁵. The cell-type proportions of MCs, SCs and FAPs were controlled in the transplanted SPiDER⁺ and SPiDER⁻ populations. Each TA muscle was engrafted with 10,000 cells, except when each senescent cell type was transplanted separately (Fig. 5k), where 5,000 cells were engrafted. Engrafted muscles were collected and processed for muscle histology 4 days after cell transplantation.

RNA interference

Freshly sorted cells or C2C12 cells (ATCC, CRL-1772) were transfected with siRNA targeting Cd36 (On-Target plus SmartPool, Dharmacon, L-062017-00-0005; 5 nM) or unrelated sequence as control (On-Target plus non-targeting siRNA Pool, Dharmacon, D-001810-10-05; 5 nM) using the DharmaFect protocol (Dharmacon, T-2003-02). Target sequences for Cd36 siRNA were as follows: 5′-CCACAUAUCUACCAAAAUU-3′, 5′-GAAAGGAUAACAUAAGCAA-3′, 5’-AUACAGAGUUCGUUAUCUA-3’, 5’-GGAUUGGAGUGGUGAUGUU-3’. Freshly sorted cells after incubation with siRNAs for 3 hours were washed and engrafted.

Cytokine array

Cytokine antibody arrays (R&D Systems, ARY028; Abcam, ab193659) were used according to the manufacturer’s protocol. For cells, freshly sorted cells were cultured for 24 h in serum-free DMEM. Cell culture supernatants were collected, centrifuged and incubated with the membranes precoated with captured antibodies. For tissue interstitial fluid, skeletal muscles of mice were dissected and slowly injected with a PBS solution with a Complete Mini EDTA-free protease inhibitor cocktail (Roche, 11836170001). The PBS exudate was then recovered centrifuged and incubated with the membranes precoated with captured antibodies. The membranes were then incubated with detection antibodies, streptavidin–HRP and Chemi Reagent Mix. The immunoblot images were captured and visualized using the ChemiDoc MP Imaging System (Bio-Rad) and the intensity of each spot in the captured images was analysed using the publicly available ImageJ software.

Proliferation assays

To assess proliferation in vivo, muscles were injured by local CTX injection, and the mice were administered with ethynyl-labelled deoxyuridine (EdU, Invitrogen, A10044; 25.5 mg kg⁻¹; i.p.) 2 h before euthanasia at 4 d.p.i. The muscles were collected and processed for immunofluorescence staining in tissue slides or cell isolation by FACS. EdU-labelled cells were detected using the Click-iT EdU Imaging Kit (Invitrogen, C10086). EdU-positive cells were quantified as the percentage of the total number of cells analysed. In vitro proliferation was quantified on freshly sorted SCs, seeded in 20% FBS Ham’s F10 medium supplemented with b-FGF (Peprotech, 100-18B-250UG; 2.5 ng ml⁻¹) in collagen-coated plates. After 3 days of culture, SCs were pulse-labelled with bromodeoxyuridine (BrdU, Sigma-Aldrich, B9285-1G; 1.5 μg ml⁻¹) for 1 h. BrdU-labelled cells were detected by immunostaining using rat anti-BrdU antibodies (Abcam, AB6326, 1:500) and a specific secondary biotinylated donkey anti-rat antibody (Jackson Immunoresearch, 712-066-150, 1:250). Antibody binding was visualized using Vectastain Elite ABC reagent (Vector Laboratories, PK-6100) and 3,3′-diaminobenzidine. BrdU-positive cells were quantified as the percentage of the total number of cells analysed.

Transwell assay

SCs were freshly isolated from regenerating muscle tissue at 3 d.p.i. and plated onto 24-well plates (Falcon, 353047) in 20% FBS DMEM supplemented with b-FGF. Subsequently, medium or freshly sorted SPiDER⁺ and SPiDER⁻ cell populations (Fig. 5i) or etoposide-induced senescent C2C12 cells (Fig. 6i) were seeded on a 0.4-µm-pore-size cell culture insert (Falcon, 353495) using the same medium. After 3 days of culture, a proliferation assay was performed on SCs with BrdU labelling as described above.

In vitro treatments

ROS^high and ROS^low SCs and FAPs were freshly isolated from regenerating muscle at 24 h after injury, seeded and cultured in the presence of NAC (10 mM) or vehicle for 3 days. After the treatment, cells were fixed and further processed for staining. C2C12 cells maintained in 10% FBS DMEM were treated with etoposide (Sigma-Aldrich, E1383, 1 µM) for 5 days to induce senescence and were collected for RNA extraction and RT–qPCR. Cells were stained with a β-galactosidase staining kit (as described below) to confirm their senescent state.

Cell staining

SA-β-galactosidase (SA-β-gal) activity was detected in freshly sorted cells and cell cultures using the senescence β-galactosidase staining kit (Cell signalling, 9860) according to the manufacturer’s instructions. Lipid droplets were stained with Oil Red O (Sigma-Aldrich, O0625) according to manufacturer instructions. ROS levels were measured by immunofluorescence using CellRox Green reagent (Invitrogen, C10444; 5 µM) according to the instructions. TUNEL assays were performed using the In Situ Cell Death Detection Kit, Fluorescein (Roche, 11684795910), cells treated with DNase were used as a positive control of the staining according to the manufacturer’s description.

Muscle histology, immunofluorescence and immuno-FISH

Muscles were embedded in OCT solution (TissueTek, 4583), frozen in isopentane cooled with liquid nitrogen and stored at −80 °C until analysis. Muscle cryosections (thickness, 10 μm) were collected and stained for SA-β-gal (AppliChem, A1007,0001), H&E (Sigma-Aldrich, HHS80 and 45235), MYH3 (DSHB, F1.652), Sirius Red (Sigma-Aldrich, 365548) or used for immunofluorescence (Supplementary Table 1). The CSA of H&E- and MYH3-antibody-stained sections, the percentage of muscle area positive for Sirius Red staining and the number of SA-β-gal⁺ cells were quantified using Image J. Double immunofluorescence was performed by the sequential addition of each primary and secondary antibody using positive and negative controls. The sections were air-dried, fixed, washed on PBS and incubated with primary antibodies according to the standard protocol after blocking with a high-protein-containing solution in PBS for 1 h at room temperature. Subsequently, the slides were washed with PBS and incubated with the appropriate secondary antibodies and labelling dyes. Telomere immuno-FISH was performed after γH2AX immunofluorescence staining with telomeric PNA probe (Panagene, F1002-5) as described previously⁵⁶.

Digital image acquisition

Digital images were acquired using an upright DMR6000B microscope (Leica) with a DFC550 camera for immunohistochemical colour pictures; a Thunder imager 3D live-cell microscope (Leica Microsystems) with hardware autofocus control and a Leica DFC9000 GTC sCMOS camera, using HC PL FLUOTAR ×10/0.32 PH1 ∞/0.17/ON257C and HC PL FLUOTAR ×20/0.4 CORR PH1 ∞/0-2/ON25/C objectives; a Zeiss Cell Observer HS with a ×20 and x40 air objective and a Zeiss AxioCam MrX camera; and a Leica SP5 confocal laser-scanning microscope with HCX PL Fluotar ×40/0.75 and ×63/0.75 objectives. The different fluorophores (three or four) were excited using the 405, 488, 568 and 633 nm excitation lines. The acquisition was performed using the Leica Application (v.3.0) or LAS X (v.1.0) software (Leica) or Zeiss LSM software Zen 2 Blue.

RNA isolation and RT–qPCR

Total RNA was isolated from snap-frozen muscles using the miRNAeasy Mini Kit (Qiagen, 1038703). PicoPure (Thermo Fisher Scientific, KIT0204) was used for RNA isolation from sorted cells. For RT–qPCR experiments, DNase digestion of 10 mg of RNA was performed using 2 U DNase (Qiagen, 1010395). cDNA was synthesized from total RNA using SuperScript III Reverse Transcriptase (Invitrogen, 18080-044). For gene expression analysis in freshly sorted SCs, FAPs and MCs, cDNA was pre-amplified using the SsoAdvanced PreAmp Supermix (Bio-Rad, 172-5160) according to the manufacturer’s instructions. qPCR reactions were performed as described previously⁵⁷. Reactions were run in triplicate, and automatically detected threshold cycle values were compared between samples. Transcripts of the Rpl7 housekeeping gene were used as the endogenous control, with each unknown sample normalized to Rpl7 content (a list of the primers used in this study is provided in Supplementary Table 2).

RNA-seq sample and library preparation

Sequencing libraries were prepared directly from the lysed cells, without a previous RNA-extraction step. RNA reverse transcription and cDNA amplification were performed using the SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing from Clontech Takara. The Illumina Nextera XT kit was used for preparing the libraries from the amplified cDNA. Libraries were sequenced using the Illumina HiSeq 2500 sequencer (51 bp read length, single-end, around 20 million reads).

Bulk RNA-seq data preprocessing

Sequencing reads were preprocessed using the nf-core/rnaseq (v.1.2) pipeline⁵⁸. Read quality was assessed using FastQC (v.0.11.8)⁵⁹. Trim Galore (v.0.5.0)⁶⁰ was used to trim sequencing reads, eliminating the remains of Illumina adaptors and discarding reads that were shorter than 20 bp. The resulting reads were mapped onto the mouse genome (GRCm38, Ensembl⁶¹ release 81) using HiSAT2 (v.2.1.0)⁶² and quantified using featureCounts (v.1.6.2)⁶³. Reads per kilobase per million mapped reads (RPKM) and transcripts per million (TPM) gene expression values were calculated from the trimmed mean of M-values (TMM)-normalized counts per million (CPM) values using the Bioconductor package edgeR (v.3.30.0)⁶⁴ and R (v.4.0.0)⁶⁵. Differential gene expression analysis and PCA were performed using the Bioconductor package DESeq2 (v.1.28.1)⁶⁶. Variance-stabilizing transformation of count data was applied to visualize the sample-to-sample distances in PCA. Genes were considered to be differentially expressed if showed an adjusted P < 0.05.

Functional profiling of cell subpopulations

Functional enrichment analysis of the subsets of differentially expressed genes was performed using g:Profiler web server⁶⁷ with the g:SCS significance threshold, ‘Only annotated’ statistical domain scope, and canonical pathway KEGG⁶⁸, Reactome⁶⁹ and Wiki Pathways⁷⁰ sets. For each gene subset, the top five significant gene sets were selected for representation.

GSEA

The RPKM matrix after the removal of low-count genes (edgeR (v.3.30.0)⁶⁴) was used as an input for the GSEA (v.4.0.3) software⁷¹. We used the signal-to-noise metric to rank the genes, 1,000 permutations with the gene set permutation type and weighted enrichment statistics. Gene set sizes were chosen as 15–500 for MSigDB 7.0 GO:BP and 10–1,000 for MSigDB 7.0 canonical pathways (BioCarta, KEGG, PID, Reactome and WikiPathways)⁷². Gene sets passing the FDR < 0.25 threshold were processed for further analysis. Network representation and clustering of GSEA results were performed using EnrichmentMap (v.3.2.1)⁷³ and AutoAnnotate (v.1.3.2)⁷⁴ for Cytoscape (v.3.7.2)⁷⁵ with the Jaccard coefficient set to 0.25.

Functional profiling of SASP

We checked whether upregulated genes (DESeq2 adjusted P < 0.05 and log₂[fold change] > 0) from each Sen versus NSen comparison can be expressed in a form of secreted proteins by combining the evidence from multiple data sources: GO⁷⁶ cellular component (GO:CC), Uniprot⁷⁷, VerSeDa⁷⁸, Human Protein Atlas⁷⁹ and experimental data reporting SASP^10,80. The genes encoding extracellular (GO:CC) and/or secreted (other sources) products, with evidence from at least one source, were included in the final list of SASP genes (1,912 in total). Functional enrichment analysis was performed using the g:Profiler web server⁶⁷ with the g:SCS significance threshold, ‘Only annotated’ statistical domain scope, and canonical pathway sets from KEGG, Reactome and Wiki Pathways. Gene sets passing the FDR < 0.05 threshold were processed for further analysis. Network representation and clustering of the g:Profiler results were performed using EnrichmentMap (v.3.2.1) and AutoAnnotate (v.1.3.2) for Cytoscape (v.3.7.2) with the Jaccard coefficient set to 0.25.

Comparative enrichment analysis of senescent cells and previously published ageing datasets

We used the minimum hypergeometric test implemented in the R package mHG (v.1.1)⁸¹ for the comparative enrichment analysis of senescent cells and previously published ageing datasets: mouse³⁶, rat (Gene Expression Omnibus (GEO): GSE53960), African turquoise killifish (GEO: GSE69122), and human (GTEx⁸² v6p). Data processing and analysis were performed as described previously³⁶.

scRNA-seq and analysis

scRNA-seq was performed using the Chromium Single Cell 3′ GEM, Library & Gel Bead Kit v3, 16 rxns (10x Genomics, PN-1000075) according to the manufacturer’s instructions and targeting a recovery of 5,000 cells per dataset. Each dataset was obtained with a sample size of two mouse biological replicates. The libraries were constructed as instructed in the manufacturer’s protocol and sequenced using the MGI DNBSEQ-T7 sequencer platform. The average read depth across the samples was 15,551 per cell. Sequencing reads were processed with STARsolo (v.2.7.3a)⁸³ using the mouse reference genome mm10 (GENCODE vM23 (ref. ⁸⁴)).

From the filtered barcode and count matricesm, downstream analysis was performed using R (v.4.0.3). Quality control, filtering, data clustering, visualization and differential expression analysis were performed using the Seurat (v.4.0.3) and DoubletFinder (v.2.0) R packages^85,86. Datasets were processed following Seurat standard integration protocol according to the tutorial instructions. Genes expressed in less than 3 cells and cells with fewer than 500 features, less than 2,000 transcripts and more than 20% reads mapping to mitochondrial genes as well as cells identified as doublets by DoubletFinder were removed. PCA was performed for dimensionality reduction and the first 30 components were used for UMAP embedding and clustering.

ATAC-seq sample and library preparation

Omni-ATAC-seq was performed in freshly sorted cells as described previously^87,88. After the transposition reaction and purification, the transposed fragments were amplified using 50 μl of PCR mix (20 µl of DNA, 2.5 µl of custom Nextera PCR primers 1 and 2, and 25 µl of KAPA HiFi HS Ready Mix for a total of 15 cycles). The PCR amplification conditions were as follows: 72 °C for 5 min; 95 °C for 30 s; 15 cycles of 95 °C for 10 s, 63 °C for 30 s and 72 °C for 60 s; and a final extension at 72 °C for 5 min. After PCR amplification, the libraries were purified, and the size was selected from 150 to 800 bp using AMPure XP beads. Paired-end sequencing was performed with 50 cycles on the Illumina NovaSeq 6000 platform.

Bulk ATAC-seq data preprocessing

Read quality was assessed using FastQC (v.0.11.8). All adaptors were removed using Fastp (v.0.21.0)⁸⁹. The clean reads were then aligned to mm10 mouse genome assembly using Bowtie2 (v.2.2.5)⁹⁰ with the settings ‘–very sensitive’. Low-mapping-quality reads were removed using samtools (v.1.3.1)⁹¹ with the settings ‘-q 30’. BigWig files were generated using deeptools (v.3.3.1)⁹² with the settings ‘-normalizeUsing CPM’. Peaks were called using Macs2 (v.2.1.0)⁹³ with the options ‘–nomodel –keep-dup -q 0.01’. For differential accessibility analysis, union peak sets were created using Bedtools (v.2.29.2)⁹⁴, reads corresponding to each region were assigned by FeatureCounts. Differentially accessible peaks were identified using DESeq2 (v.1.24.0) with the criteria of adjusted P < 0.1 and an absolute value of log₂[fold change] > 1. Differentially accessible peaks were further annotated by HOMER (v.4.10.4)⁹⁵, the associated motif enrichment analysis was performed by HOMER using the default settings.

Analysis of senescence-induced changes in promoter chromatin accessibility

An MA plot (log₂-transformed fold change versus mean average) was used to visualize changes in chromatin accessibility for all peaks. As a peak score, we used an average of TPM-normalized read counts: (1) reads per kilobase were calculated by division of the read counts by the length of each peak in kilobases; (2) the per million scaling factor was calculated as a sum of all reads per kilobase for each sample; (3) reads per kilobase were divided by the per million scaling factor; (4) peaks with the ‘promoter-TSS’ annotation TSS ± 1kb were selected and the average was calculated for each group. For MA plots, we included only those peaks with an average normalized signal > 5. The number of peaks with a log₂-transformed fold change of >1 or <−1 was calculated. Normalized ATAC-seq signal profiles of proximal promoters were visualized for key genes using the Integrative Genomic Viewer (v.2.8.13)⁹⁶.

Transcription factor analysis and activity prediction

For the analysis of transcription regulation, we combined the results of several methods: (1) motif enrichment analysis of differentially expressed genes with the TRANSFAC_and_JASPAR_PWMs and ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X libraries using the R package EnrichR (v.2.1)⁹⁷; (2) upstream regulator analysis of differentially expressed genes using the commercial Ingenuity Pathway Analysis (IPA, QIAGEN) software⁹⁸; (3) analysis of transcription factor differential expression using DESeq2 (v.1.28.1); (4) motif enrichment analysis of differentially accessible regions using HOMER (v.4.10.4).

Potential regulators from EnrichR and IPA results passing the threshold of P < 0.05 were used to build a union set of transcription factors, which was further filtered to retain only the molecules with DESeq2 baseMean value > 0. For further validation of the activity status, transcription factors were matched to the known HOMER motifs passing the Benjamini Q < 0.05 threshold.

A discrete scoring scale (inhibited, possibly inhibited, unknown/contradictory, possibly activated, activated) was used to evaluate transcription factor activity based on combined evidence from the EnrichR, IPA, DESeq2 and HOMER results. We used z-score statistics to define the activity status of transcription factors from the EnrichR analysis results by matching the differential expression of target genes with activatory and inhibitory interactions from the Bioconductor package DoRothEA (v.1.0.0)⁹⁹ and the web-based TRRUST v.2 database¹⁰⁰. To define the activity status of transcription factors from IPA upstream regulators analysis results, IPA-calculated z-score and analysis bias was taken into account. Activity predictions were further corrected by differential expression of transcription factors using DESeq2. The expression z-score statistical value was calculated to functionally classify transcription factors as activators or repressors on the basis of the proportion of upregulated and downregulated target genes. We further calculated the chromatin accessibility z-score to estimate the prevalence of HOMER motif enrichment in open versus closed regions that together with the predicted transcription factor function enabled us to validate the RNA-seq activity predictions using ATAC-seq data.

To estimate the level of confidence, for each enrichment result, we calculated a discrete ‘trust’ score, with each point assigned for: (1) EnrichR adjusted P < 0.05; (2) IPA P < 0.05; (3) activity status ‘activated’ or ‘inhibited’; (4) unidirectional absolute z-scores of >2 from both the EnrichR and IPA results; (5) concordance between transcription factor differential expression and the prediction of its activity score; (6) activity validated by the analysis of ATAC-seq data. Transcription factors with average trust > 1 were processed for further analysis.

Functional profiling of transcription factor target gene regulation

For each transcription factor, we merged the target genes from EnrichR and IPA results, split them into upregulated and downregulated and processed them to functional enrichment analysis of canonical pathways (KEGG, Reactome) and GO:BP using R package gprofiler2 (v.0.1.9)¹⁰¹ with the following parameters: correction method ‘FDR’, ‘custom_annotated’ domain score consisting of target genes for all studied transcription factors. Electronic GO annotations were excluded. Gene sets that passed the FDR < 0.05 threshold were processed for further analysis. For GO:BP, we selected the gene sets with a term size of >15 and <500 genes. Transcription factors were further mapped based on the matching terms from the gprofiler2 results to the main functional clusters of the gene sets created previously in GSEA/Cytoscape analysis.

Transcription factors mapped to the same functional cluster in ≥8 (out of 12) Sen versus NSen comparisons were processed for further filtering. We scored as 1 point in each case when any of the following attributes had a value above the upper quartile for a given cluster: number of comparisons, percentage of GSEA terms among all terms with enrichment, −log₁₀ of the average minimum FDR and average trust score. Moreover, we scored as 1 point if the transcription factor was associated with senescence in literature. For graphical representation, we selected examples of transcription factors and target genes based on literature research: 19 transcription factors (out of 29 with a score of ≥2) mapped to 9 clusters (matrix remodelling/fibrosis, interferon signalling, chemotaxis, lipid uptake, IGF regulation, detoxification, gene expression and protein translation, cell cycle, and DNA repair).

Functional profiling of transcriptional regulation of SASP

For each transcription factor upregulated, target genes from the EnrichR and IPA results were merged and intersected with the list of SASP genes. For SASP genes, we extracted GO:MF terms, clustered them into 12 categories (adhesion molecule, chemokine, complement component, cytokine, enzyme, enzyme regulator, extracellular matrix constituent, growth factor, hormone, ligand, proteinase and receptor) and estimated the enrichment of GO:MF clusters with a hypergeometric test using the R function phyper. Correction for multiple comparisons was performed using the Benjamini–Hochberg procedure.

Transcription factors that had target enrichment in the same GO:MF cluster with P < 0.05 in ≥8 (out of 12) Sen versus NSen comparisons were processed for further filtering. We scored 1 point in each case in which any of the following attributes had a value above the upper quartile for a given cluster: the percentage of secreted proteins among targets, the number of comparisons with P < 0.05, the number of comparisons with adjusted P < 0.05, −log₁₀ of the average P value, the average trust score. Moreover, we used ATAC-seq data analysis to score 1 point in cases in which the transcription factor motif was present in the promoter region of at least one SASP gene within the cluster. For graphical representation, we selected 17 transcription factors with a score of ≥2 and with ≥3 comparisons with adjusted P < 0.05. They were associated with five GO:MF categories (extracellular matrix constituent, cytokine, chemokine, complement component and growth factor), for which we selected the most common target genes.

Analysis of lipid metabolism gene set

For the analysis of lipid metabolism, we constructed a gene set using data from multiple sources: KEGG pathway maps (fatty acid degradation, cholesterol metabolism, regulation of lipolysis in adipocytes), WikiPathways (fatty acid oxidation, fatty acid beta oxidation, mitochondrial LC-fatty acid beta-oxidation, fatty acid omega oxidation, fatty acid biosynthesis, triacylglyceride synthesis, sphingolipid metabolism (general overview), sphingolipid metabolism (integrated pathway), cholesterol metabolism (includes both Bloch and Kandutsch–Russell pathways) and cholesterol biosynthesis) and literature research^102,103,104. We further estimated the expression of these genes by filtering DESeq2 results (adjusted P < 0.05 in at least 3 out of 12 comparisons) and extracted log₂-transformed fold change values to plot the difference in expression between senescent and non-senescent cells.

Reconstruction of ligand–receptor mediated cell–cell communication networks

For reconstructing cell–cell communication networks, we modified the single-cell-based method, FunRes, to account for bulk gene expression profiles³⁹. In brief, transcription factors with an expression value of more than 1 TPM were considered to be expressed. Receptors regulating these transcription factors were detected using a Markov chain model of signal transduction to detect high-probability intermediate signalling molecules¹⁰⁵. Ligand–receptor interactions between two cell populations were reconstructed if (1) the receptor is expressed and regulates any transcription factor, (2) the ligand is expressed and (3) the receptor–ligand interaction is contained in the cell–cell interaction scaffold. Finally, a score is assigned to every interaction by multiplying the average receptor and ligand expression in their respective cell populations. Significance was assessed by permuting cell population labels 100 times and recomputing the interaction scores in the permuted datasets. Interactions were considered to be significant if they were at least 2 s.d. greater than the mean of the permuted interaction scores. Only significant interactions were retained in the final network.

Downstream analysis of senescence-induced ligand–receptor interactions

For the functional profiling, we selected ligand–receptor interactions between three senescent cell populations (SCs, FAPs and MCs) and a non-senescent SC population in old mice at 3 d.p.i. We used the Bioconductor package SPIA (v.2.40.0)⁴⁰ with a reduced set of non-disease KEGG pathway maps to evaluate the activity of a pathway’s downstream ligand–receptor interactions. For each interaction, differentially expressed target transcription factors in non-senescent SCs were split into upregulated and downregulated in comparison to senescent SCs. As a reference set of genes, we took a list of target transcription factors from all of the interactions studied. SPIA analysis was performed with 2,000 permutations, and pPERT and pNDE were combined using the Fisher’s product method. Pathways passing the pGFdr < 0.05 threshold were considered to be significantly enriched. For each pathway, we calculated the ratio of ligand–receptor interactions that activate or inhibit the pathway to the total number of interactions analysed. For results representation, we selected eight activated and eight inhibited pathways with the highest ratio of interactions.

Statistical analysis

The sample size of each experimental group is described in the corresponding figure caption, and all of the experiments were conducted with at least three biological replicates unless otherwise indicated. GraphPad Prism was used for all statistical analyses except for sequencing-data analysis. Quantitative data displayed as histograms are expressed as mean ± s.e.m. (represented as error bars). Results from each group were averaged and used to calculate descriptive statistics. Mann–Whitney U-tests (independent samples, two-tailed) were used for comparisons between groups unless otherwise indicated. P < 0.05 was considered to be statistically significant. Experiments were not randomized.

Reporting summary

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