Stem cells tightly regulate dead cell clearance to maintain tissue fitness

The following previously generated mouse lines were used in this study: Rxra^fl (ref. ⁴¹; Jax stock 013086), Sox9-creER (ref. ⁴²), Krt14-rtTA (ref. ⁴³; Jax stock 008099), Rosa26^{lox-STOP-lox-YFP} (ref. ⁴⁴; Jax stock 006148; referred to as R26^YFP), Rosa26^mTmG (ref. ⁴⁵; Jax stock 007576; referred to as R26^mTmG), Rosa26^Brainbow2.1 (ref. ⁴⁶; Jax stock 013731, referred to as R26^Brainbow2.1), Rosa26^{lox-STOP-lox-Cas9-EGFP} (ref. ⁴⁷; Jax stock 026175, referred to as R26^Cas9-EGFP), Rosa26^{lox-STOP-lox-DTA} (ref. ⁴⁸; Jax stock 010527, referred to as R26^DTA) and Mertk^{−/−} (full knockout; ref. ⁴⁹). The Mertk-knockout mice used in this study are referred to as Mertk^{−/−V2} in the originating paper. Wild-type CD1 or C57BL/6 mice were originally purchased from Charles River and The Jackson Laboratories, respectively, and maintained as in house colonies.

Mice were maintained and bred under specific-pathogen-free conditions at the Comparative Bioscience Center (CBC) at The Rockefeller University, an Association for Assessment and Accreditation of Laboratory Animal Care (AALAC)-accredited facility. Mertk-knockout mice and C57BL/6J wild-type controls (maintained as separate colonies) were bred and maintained in a specific-pathogen-free facility at Yale University. All mice were bred and maintained under a strict 12-h light cycle and fed with standard chow. The temperature of the animal rooms was 20–26 °C, and the humidity was 30–70%. Adult mice were housed in cage with a maximum of five mice. All mouse protocols were approved by the Institutional Animal Care and Use Committee (IACUC) at The Rockefeller University, or by the IACUC at Yale University.

For comparative assessments of phenotype between control and mutant mice, age and sex matched mice were used, with preference given to littermate controls wherever possible, and sample size greater than three mice per genotype or condition across multiple litters whenever possible. For our inducible overexpression studies, a Krt14-rtTA^+/− (heterozygous) male was mated with CD1 females and all offspring were transduced with lentivirus at E9.5 (see following sections). Offspring of both genotypes received doxycycline by intraperitoneal injection (0.5 mg per mouse) at P14 to activate Krt14-rtTA within 12 h, and expression was maintained by feeding the mother and pups doxycycline (2 mg kg^−1) chow (Bioserv). Krt14-rtTA^− mice were used as control, with Krt14-rtTA⁺ mice as the experimental group. To generate Rxra control and cKO mice for experiments, the Rxra^fl line was crossed with Sox9-creER⁺; R26^YFP mice. Sox9-creER^− mice with any Rxra^fl;R26^YFPgenotype, and Rxra^fl/+;Sox9-creER⁺;R26-YFP^fl/+mice were used as controls, while experimental mice were Rxra^fl/fl;Sox9-creER⁺;R26-YFP^fl/+. All mice received tamoxifen (2% in corn oil) (Sigma-Aldrich) to activate Sox9-creER, administered by intraperitoneal injection once a day for 3 days, as indicated. Sox9-creER was similarly activated when crossed to R26^mTmG (to label HFSCs prior to FACS-isolation and culture), R26^Brainbow2.1 (to stochastically label HFSCs and identify functional phagocytes), R26^Cas9-EGFP (to mosaically knockout Rxra), and upon transduction with the inducible Cyp26b1 expression construct (to degrade RA in HFSCs). To activate Sox9-creER sparsely when crossed to R26^DTA, 2% tamoxifen was intraperitoneally injected once early in second telogen.

Hair cycle staging

Male and female mice have different hair cycle lengths due to a longer telogen quiescence phase in females, but otherwise progress through the hair cycle similarly. In addition to sex, strain and individuals also affect hair cycle stages. Therefore, we always determine hair cycle stage by visual inspection, and morphological staging on sectioned tissue. Specifically, for C57BL/6 pure and mixed backgrounds, visual inspection was performed by trimming full-length telogen hairs with electric clippers to reveal dorsal skin. Hair follicle entry into anagen was determined by darkening of skin and reappearance of hair. Catagen progression was determined by lightening of the skin, which appears black at the end of anagen, to a near complete loss of pigmentation (greyish-pink skin) by late catagen. Entry into telogen was marked by the appearance of completely unpigmented (pink) skin. In unpigmented mice (CD1 strains), histological analysis of hair follicle morphology was used to confirm hair cycle staging based on relative age.

For all experiments, a small piece of midline dorsal skin was taken in parallel from each mouse, fixed and processed for sectioning and immunofluorescence to precisely stage the hair cycle. Hair cycle was staged based hair follicle morphology^50,51,52,53, as well as by immunofluorescence for markers of anagen (EdU-incorporation following a 2 h pulse chase or Ki67 staining), catagen (cleaved caspase-3 and/or TUNEL positivity) or telogen (pSMAD1/5/9 and/or LEF1). For samples obtained from anagen or catagen stage mice the dorsal back skin was subdivided in two along the anterior-posterior axis prior to the experiment, as precise hair cycle stage differs anterior to posterior (with anterior generally one substage ahead). For anagen samples, 5′-ethynyl-2′-deoxyuridine (EdU) was injected intraperitoneally (50 μg g^−1) (Sigma-Aldrich) and chased for 2 h prior to collection.

To assess hair cycling defects following ablation of HFSC-mediated apoptotic corpse clearance in catagen, cohorts of catagen-specific Rxra control and cKO mice were generated as described above. Similar cohorts of HX531 or annexin V intradermally injected wild-type CD1 mice were generated as described below. Both sets of mice were shaved and examined weekly over the course of second telogen for skin darkening (Rxra line) and hair regrowth (all experimental paradigms). Mice were collected upon initial signs of anagen re-entry (skin darkening and/or small hairs breaking skin surface). For the Rxra line, co-housed littermates were collected once one animal showed signs of anagen entry. Hair cycle was staged based on hair follicle morphology and immunofluorescence as described and compared across co-housed littermates (Rxra line) or between contralateral vehicle- and inhibitor-injected back skin within each mouse (intradermal manipulations).

Induction and knockout constructs

RXRα induction

To make TRE-RXRa-Myc; pGK-H2B-RFP, human RXRa cDNA was PCR-amplified from pSV-Sport-RXRα (a gift from B. Spiegelman; Addgene #8882)⁵⁴, and a NheI site was introduced at the 5′ end. This was then inserted into the NheI and EcoRI restriction sites of a pLKO vector modified to contain the inducible tetracycline response element at the 5′ end, as well as a 3′ MYC epitope tag. Prior to packaging this construct as a lentivirus, induction of RXRα was tested in culture using FACS-isolated Krt14-rtTA⁺ keratinocytes grown in E300 medium. In brief, cells were transiently transfected with the TRE-RXRα construct using Effectene into keratinocytes in a 6-well plate format, following the manufacturer’s protocol (Invitrogen). Forty-eight hours later, doxycycline (100 ng ml^−1) was added to induce RXRα expression for a further 24 h. Cells were fixed and stained for Myc-tag and RXRa, as described for cell culture immunofluorescence.

RXRα-mosaic knockout

To identify efficient CRISPR single guide RNAs (sgRNA) against mouse Rxra, we synthesized oligonucleotides targeting exon 4 with BsmBI restriction sites at 5′ and 3′ respectively (IDT). Sequences used are available in Supplementary Table 6. Oligonucleotides were subcloned into pLentiGuide-Puro (a gift from F. Zhang; Addgene #52963), following the Zhang laboratory protocol⁵⁵. To select a guide for in vivo use, we first tested the cutting efficiency in culture using K14Cre⁺; R26^Cas9-EGFP expressing keratinocytes. pLentiGuide-Puro constructs were transiently transfected using Effectene into keratinocytes in a 6-well plate format, following the manufacturer’s protocol (Invitrogen). After 72 h, genomic DNA was collected using QuickExtract DNA Extraction Solution (Lucigen), and guide DNA was prepared by heating to 65 °C for 10 min followed by heat inactivation at 95 °C for 2 min. Following PCR amplification of each guide target region, a T7 endonuclease I cutting assay (NEB) was used to identify the extent of insertions and/or deletions for each guide. The most efficient guide showed ~70% genome editing events in vitro, and was cloned with its U6 promoter into a modified pLKO vector containing a constitutive pGK-driven mScarlet fluorophore (3′ to the sgRNA) for lentiviral preparation.

RARE reporter

To identify cells responding to RA, we obtained pGL3-RARE-luciferase as a gift from T. M. Underhill (Addgene plasmid #13458; http://n2t.net/addgene:13458) and subcloned the RA response element (RARE) into a modified pLKO backbone behind a minimal SV40 promoter to drive RFP expression with pGK-driven H2B-GFP⁵⁶.

Cyp26b1 induction

To deplete active RA metabolites from cells inducibly, we over-expressed mouse Cyp26b1 cDNA (Origene, MC205286) from a CAG promoter, interrupted by a Lox-Stop-Lox (LSL) cassette in a modified pLKO backbone. As a viral control, we used pGK-driven RFP in the opposite orientation⁵⁶.

Lentiviral preparation and Injection

High-titre lentivirus was prepared and E9.5 embryos of indicated genotypes were infected with lentivirus delivered by ultrasound-guidance microinjection into the amniotic sac as previously described^57,58. At E9.5 the surface ectoderm exists as a single layer of unspecified skin progenitors, which can be efficiently, selectively and stably transduced by the viral DNA, without transduction of dermal cell types⁵⁸.

HFSC culture

All primary HFSC lines were grown on a layer of mitomycin C-inactivated 3T3/J2 feeder fibroblast cells, and maintained in E intermediate (300 μM) calcium media⁵⁹ supplemented with 10 μM Y-27632 (Selleckchem) (E300-Y medium)⁶⁰. The 3T3/J2 fibroblast cell line⁵⁹ was expanded in DMEM/F12 medium (Thermo Fisher Scientific) with 10% CFS (Gibco), 100 U ml^−1 streptomycin and 100 mg ml^−1 penicillin. Cells were grown at 37 °C, with 7.5% CO₂, and medium was routinely changed every 2–3 days. Cell lines were grown to confluency, then propagated by digesting with 0.25% Trypsin EDTA (Gibco) for 5–10 min at 37 °C and resuspended with culture medium for passaging. Experiments were conducted with cells at passages 8–10. For experiments, cells were switched to E intermediate calcium medium without Y-27632 (E300 medium) and cultured for 24–48 h prior to the experiment. All cell lines were maintained in a culture facility routinely testing negative for mycoplasma contamination.

Primary HFSCs were derived from the following mouse crosses at second telogen: R26-mTmG^fl/+; Sox9-creER^− (mTomato⁺ HFSCs to make apoptotic corpses and necrotic debris), R26-mTmG^fl/+; Sox9-creER⁺ (mGFP⁺ HFSCs to make naive HFSC to expose to corpses), Rxra^+/+; Sox9-creER⁺; R26-YFP^fl/fl (Rxra wild type YFP⁺ HFSCs), and Rxra^fl/fl; Sox9-creER⁺; R26-YFP^fl/fl (Rxra cKO YFP⁺ HFSCs). All HFSCs were FACS-isolated (described later) and cultured as described. To generate Rxra wild type and cKO HFSC lines, cells were FACS isolated (described later), and cultures were established prior to activation of Sox9-creER by 4-hydroxytamoxifen (4-OHT). At passage 2, Sox9-creER was activated in culture by 4-OHT in solution (Sigma Aldrich); to do so it was used at a final concentration of 1 μM in E300-Y medium to treat HFSCs. Medium plus 4-OHT was refreshed each day for three consecutive days, and then replaced by E300-Y. HFSCs were allowed to grow for 4 further days prior to FACS isolation of YFP⁺ cells. All HFSC lines, as well as the 3T3/J2 fibroblast line, were functionally and morphologically validated as HFSC or fibroblast lines respectively.

To generate HFSCs carrying the RARE-driven RFP with pGK-driven H2B-GFP (‘RARE reporter’), the RARE-reporter virus was transduced⁵⁷ into wild-type primary HFSCs derived from a second telogen mouse, as described⁵⁶. Stably integrated RARE-reporter HFSCs were FACS sorted on the basis of H2B-GFP, and RARE-driven RFP induction within 4-6 hrs was confirmed by addition of 100 nM 9cRA or 100 nM ATRA ±1 μM AGN 193109.

In pilot experiments nuclear accumulation (by immunofluorescence) of RXRα or RARγ peaked at 30 min post corpse exposure, and so that time point was used to assess immediate effects of corpse-derived signals or recombinant molecules in subsequent experiments. Similarly, the number of corpse-containing HFSCs plateaued at 4–6 h after corpse addition, and so transcriptional activation of the phagocytic programme, surface expression of phagocytic receptors and corpse engulfment (latter two by FACS) were routinely assessed at that time point (Supplemental Fig. 7).

To prepare corpses or secondarily necrotic debris, fully confluent mTomato⁺ HFSCs were treated with 200 μM cisplatin (in 0.9% saline) for 18 h (apoptotic corpses) or 48 h (necrotic debris). Dead and dying cells were collected from the supernatant by pelleting at 700g for 5 min, washed once with E300 medium, and returned to the plate in fresh E300 medium. Prior to returning the floating corpses, dying adherent cells were rinsed with PBS to remove residual cisplatin and a minimal amount of fresh E300 medium was added. Corpses or necrotic debris were allowed to condition the medium for a further 3–4 h. Apoptotic cell corpses were collected by tapping the side of the plate and pipetting their medium over them to detach dying cells. Floating and detached corpses were collected and pelleted by centrifugation as before. Corpse-conditioned medium was carefully removed to a separate tube, before resuspending corpses in a minimal volume of fresh E300 medium. To label any corpses/debris derived from 3T3/J2 fibroblasts, corpses were next incubated with DiI-CM (Invitrogen) for 5 min at 37 °C before pelleting and washing with PBS as before. DiI/mTomato⁺ corpses were resuspended in their corpse-conditioned medium, counted and aliquoted (in their conditioned medium) directly onto experimental plates (medium removed prior) at a ratio of roughly 10 corpses:1 HFSC. For corpse-conditioned medium experiments, corpses were prepared as described and then spun out of the medium at 1,200g for 10 min. Corpse-conditioned medium was further strained through a 0.45-μm syringe filter, prior to use. Corpses, necrotic debris and/or conditioned medium were always prepared immediately prior to their use.

Manipulation of corpse-derived signals was achieved by adding small molecule inhibitors to the corpses after the removal of cisplatin (16 μM BEL (Sigma Aldrich); 100 nM MPA08 (Tocris Bioscience)) or by incubating the corpses with recombinant molecules for 15-20 min prior to adding them to naive HFSCs (1 U Apyrase; 1 ng ml^−1 annexin V (both Tocris Bioscience)). Vehicle treated control corpses were incubated with 1% DMSO. Similar preparations were made for corpse-conditioned medium. To manipulate corpse-sensing mechanisms on HFSCs, naive mGFP⁺ or YFP⁺ or RARE-reporter HFSCs were pretreated with the indicated antagonist (1 nM UVI3003; 1 μM HX531; 1 μM JTE013; 100 nM BMS 777607; 1 μM AGN 193109) (first four: Tocris Bioscience; last one: R&D Systems) for 30 min prior to corpse or corpse-conditioned medium addition. When adding corpses or conditioned medium, the concentration of antagonist was maintained by adding an additional amount of the appropriate compound to the corpses and/or conditioned medium. Experiments were performed in biological duplicate or triplicate and repeated at least twice on separate days. For data visualization, all replicates across independent experiments are represented, unless indicated otherwise in figure legends. Experiments manipulating corpse-derived signals by small molecule inhibitors were performed in parallel, such that a core set of control medium and corpses + Veh experimental replicates exist. For presentation purposes, small molecule inhibitors were grouped according to the step of phagocytosis they affect in separate figure panels, and so the core set of control medium and corpses + Veh is repeated for each (Fig. 4b and Extended Data Figs. 3g, 6d and 7a).

To test the ability of recombinant molecules to recapitulate corpse secreted signals, mGFP⁺ or YFP⁺ or RARE-reporter HFSCs were cultured in E300 medium plus the indicated concentrations of recombinant molecules (see figures and figure legends). Molecules were prepared and stored as stock solutions according to manufacturer’s instructions. In brief, 9cRA and ATRA (both R&D Systems) were each dissolved in 100% DMSO, protected from light, and stored long term at −80 °C with working solutions kept at −20 °C. Recombinant LPC, S1P and AA (all Tocris Bioscience) were each dissolved in 100% ethanol and stored like the retinoids. Free nucleotides were purchased as 100 mM stocks of dATP or dUTP as sodium salts in ultrapure water (NEB) and stored at −20 °C. Stock solutions were diluted individually or in combinations with E300 medium for experiments. Experiments were performed in biological duplicate or triplicate and repeated at least twice on separate days. For data visualization, all replicates across independent experiments are represented, unless indicated otherwise in figure legends. Experiments testing different concentrations and combinations of recombinant molecules for induction of RXRα⁺ or RARγ⁺ nuclear accumulation in cultured HFSCs were performed in parallel, such that a core set of control medium experimental replicates exist. For presentation purposes, sets of concentrations were split across separate figure panels and so the core set of medium control experiments are repeated in each (Fig. 4d and Extended Data Fig. 7e).

To examine the effects of repeated exposure to necrotic damage in vitro, conditioned medium was prepared from live HFSCs, apoptotic corpses or necrotic debris as described above. For proliferation studies, naive mGFP⁺ HFSCs were exposed to new conditioned medium daily three times in 96-well plates on fibroblast feeders and cells were counted by GFP fluorescence in a BioTek Cytation 5 cell imaging multimode reader. Experiments were set up in three biological replicates, and cell numbers were counted daily. Experiment was performed independently three times with one representative experiment shown. To assess the effect of repeated exposure to necrotic debris on colony forming efficiency and size, naive mGFP⁺ HFSCs were plated in 12-well plates on fibroblast feeders and exposed to conditioned medium as described for proliferative studies. Twenty-four after the third exposure to conditioned medium, HFSCs were collected by trypsinization, washed with PBS and replated on fibroblast feeders in E300-Y medium at a density of 10,000 HFSCs per replicate. Medium was replaced 4 days later and changed every second day after that for a total of 14 days of growth. Colony number and size were quantified by GFP fluorescence in a BioTek Cytation 5 cell imaging multimode reader. Experiment was performed twice in experimental triplicates and all six replicates shown.

For colony forming assays on HX531- or annexin V-injected mice, primary HFSCs were FACS-isolated from either uninjected wild type or contralateral skin sites of intradermally injected vehicle versus inhibitor mice as described below. Isolated HFSCs were stained with DiI-CM as described for corpse preparation prior to plating in technical triplicates of 2,000 HFSCs each on 3T3/J2 feeder fibroblasts in 6-well plates in E300-Y medium and allowed to grow for 4 days before medium was changed. At one week post-seeding, colony number and size were counted under an upright fluorescent microscope using bright field and RFP (DiI-CM) fluorescence. Technical triplicates were averaged, and data presented per mouse.

Intradermal injections

Adult mice in early second catagen (CatII) were anaesthetized using isoflurane prior to intradermal injections. Isofluorane anaesthetization was maintained throughout the procedure using a nose cone for delivery. Back skin was shaved using electric clippers and the surface sterilized by wiping with ethanol wipes. Vehicle or small molecule containing solutions were prepared by diluting appropriate chemical in sterile PBS plus 1% FluoSpheres Carboxylate-modified microspheres (1.0 μm, Thermo Fisher Scientific F8816) to assess placement of the injection on tissue sections. To find the injection site for repeated intradermal injections on subsequent days, a small dot was made with permanent marker which the needle was inserted through. One-millilitre insulin syringes with the needle bent to approximately 45° were used to shallowly inject through the epidermis to the dermal space approximately 3–5 mm from the injection site. An injection volume of 25 μl was delivered per injection site, with an average of 4 injection sites per mouse: two anterior and two posterior. Vehicle injections (10% DMSO) were randomly designated to either the left or right side, with the contralateral skin receiving the indicated small molecules. Following the injections, mice were placed in their home cage on a heating pad to recover. Compounds were prepared at 100× the cell culture working solution, from the same stock solutions. To inhibit the phagocytic programme across the course of catagen, intradermal injections were performed three times, separated by 20–24 h. Mice were euthanized by lethal CO₂ administration, 4 h after the final injection (at CatVII–VIII), or allowed to progress into telogen for further analysis. For colony-forming assays, injected mice were euthanized 2 days after the end of catagen, back skin was manually dissected into 10 mm² around sites of injections and HFSCs were FACS isolated as described below. Isolated HFSCs were plated on feeders and analysed as described for colony-forming assays (above). To analyse hair cycling defects upon transient inhibition of corpse engulfment during catagen, intradermally injected mice were followed throughout the course of second telogen as described above.

Tissue collection and sectioning

For immunofluorescence analysis of tissue sections, mice were shaved following lethal CO₂ administration and their back skin dissected. Back skin was stretched onto Whatman paper for stability, and immediately prefixed in 1% or 4% paraformaldehyde (PFA) for 1 h at 4 °C or 30 min at 25 °C, respectively. After fixing, tissue was washed twice with PBS for 10 min at 4 °C, before incubating in 30% sucrose in PBS at 4 °C overnight. Tissue was embedded in OCT medium (VWR) and frozen on dry ice blocks before storage at −80 °C. Alternatively, fresh frozen tissue was prepared without prefixation by directly embedding the skin in OCT after it was placed on Whatman paper. Frozen tissue blocks were sectioned at 20 um on a Leica cryostat and mounted on SuperFrost Plus slides (Thermo Fisher). When necessary, sections were stored at −20 °C prior to use.

Immunofluorescence

Skin sections

Following sectioning, tissue was allowed to dry on the slide for 1 h in a partially closed slide box. Fresh frozen tissue was post-fixed with 4% PFA for 5 min, followed by washing in phosphate-buffered saline (PBS) three times for 5 min each. Pre-fixed tissue sections started with the PBS wash step to remove attached Whatman paper. Following washes, samples were permeabilized and blocked in blocking buffer (5% donkey serum, 2.5% fish gelatin, 1% BSA, 0.3% Triton in PBS) for 1 h at room temperature. Primary antibodies were incubated overnight at 4 °C, samples were washed for 5 min in PBS (three times) at room temperature, and secondary antibodies were incubated together with DAPI (to label nuclei) for 1 h at room temperature. Following three final PBS washes of 5 min each, samples were mounted in Prolong Diamond Antifade Mountant (Invitrogen) for imaging. For TUNEL labelling, the Cell Death Detection Kit (TMR red or FITC; Roche) was used according to manufacturer’s instructions, with application of secondary antibodies. A modification was made to halve the concentration of the substrate labelling component to reduce background fluorescence in the skin. For phospho-STAT3 staining, tissue was incubated in ice-cold methanol for 20 min at −20 °C, followed by three times PBS washes prior to blocking. Antibodies were used as follows: rabbit anti-cleaved-caspase-3 (Cell Signaling, 9661, 1:250), rat anti-RFP (Chromotek, 5F8, 1:1,000), rabbit anti-RFP (MBL, PM005, 1:1,000), chicken anti-GFP/YFP (Abcam, ab13970, 1:1,000), goat anti-P-cadherin (R&D, AF761, 1:250), rabbit anti-keratin14 (Fuchs laboratory, 1:200), rabbit anti-keratin24 (Fuchs laboratory, 1:200), sheep anti-Ki67 (Novus Biologicals, AF7649, 1:200), rabbit anti-MYC epitope (71D10) (Cell Signaling, 2278, 1:250), rat biotinylated anti-CD45 (Biolegend, 5530, 1:200), rabbit anti-RXRα (D6H10) (Cell Signaling, 3085, 1:250), rabbit anti-RARγ (D3A4) (Cell Signaling, 8965, 1:250), rabbit anti-MFGE8 (Invitrogen, PA5-109955, 1:200), rat AlexaFluor647-conjugated anti-F4/80 (BM8) (Biolegend, 123121, 1:200), rat biotinylated anti-ITGA6 (also known as CD49f) (GoH3) (Biolegend, 313603, 1:500), rabbit anti-cJun (60A8) (Cell Signaling, 9165, 1:250), rabbit anti-FosB (5G4) (Cell Signaling, 2251, 1:250), and rabbit anti-phospho-STAT3 (Tyr705)(D3A7) (Cell Signaling, 9145, 1:250). All secondary antibodies used were raised in a donkey host, and conjugated to AlexaFluor488, Rhodamine, or AlexaFluor647 (Jackson ImmunoResearch Laboratory; 1:500). Catalogue numbers (given in order of: AlexaFluor488, Rhodamine, and AlexaFluor647 conjugates) for donkey anti-rabbit antibodies (711-545-152; 711-295-152; 711-605-152), for donkey anti-rat antibodies (712-545-150; 712-295-150; 712-605-150), for donkey anti-chicken antibodies (703-545-155; 703-295-155; 703-605-155), for donkey anti-goat antibodies (705-545-003; 705-295-003; 705-605-003), and for donkey anti-sheep AlexaFluor647 (713-605-003). 4′,6-diamidino-2-phenylindole (DAPI) was used to label nuclei (1:10,000). To co-stain RARγ and RXRα, the rabbit primary antibodies were individually directly conjugated to one of AlexaFluor350, AlexaFluor488, AlexaFluor568, or AlexaFluor647 using the rabbit specific Zenon Antibody Labelling Kit (Thermo Fisher Scientific) and following manufacturer’s instructions.

Iterative bleaching extends multiplexity

Tissue sections were processed, fixed and sectioned as for normal immunofluorescence, with one modification. Tissue cryosections (25-μm) were placed in glass bottom slide wells coated with chrome gelatin alum to securely adhere the tissue to the glass coverslip. The IBEX protocol was followed as described⁶¹, with the following modifications. Following blocking with our blocking buffer (above), tissue was incubated with primary antibodies directly conjugated to fluorophores for 3 h at room temperature, followed by PBS washes and imaging using a spinning disk confocal microscope. DAPI was used as described before, as a fiducial stain to align images from iterative cycles. To bleach fluorophores between iterative cycles of staining and imaging, we exposed tissue to 1 mg ml^−1 of lithium borohydride for 15 min at room temperature, followed by three 1× PBS washes. Antibodies used were as follows: (panel 1) rat anti-Foxp3-AlexFluor488 (FJK-16s) (ThermoFisher, 53-5773-82, 1:100), In situ cell death detection kit, TMR red (Roche), rat anti-CD8-AlexaFluor647 (BioLegend, 100724, 1:150); (panel 2) rat anti-CD206-AlexaFluor488 (MMR) (BioLegend, 141710, 1:500) and rat anti-CD68-AlexaFluor647 (BioLegend, 137004, 1:500); (panel 3) rat anti CD11c-AlexaFluor488 (N418) (BioLegend, 117311, 1:100) and rat anti-Ly6g-AlexaFluor647 (1A8) (BioLegend, 127610, 1:150); (panel 4) rat anti-ITGA6-AlexaFluor488 (BioLegend, 313608, 1:150) and rat anti-Langerin–AlexaFluor647 (929F3.01) (Novus Biologicals, DDX0362A647-100; 1:100); (panel 5) rat anti F4/80-AlexaFluor488 (BioLegend, 123122, 1:150) and rat anti-CD172a (Sirpα)-AlexaFluor647 (BioLegend, 144028, 1:150); (panel 6) hamster anti-TCRgd-AlexaFluor488 (BioLegend, 118128, 1:100) and rat anti-Tim4-AlexaFluor647 (RMT4-54) (BioLegend, 130008, 1:150); (panel 7) rat anti-CD4-AlexaFluor488 (RM4-5) (BioLegend, 100529, 1:100) and rat anti-CD3-AlexaFluor647(17A2) (BioLegend, 100209, 1:100); (panel 8) Avidin–FITC (ThermoFisher Scientific, A821 1:1,000) and rat anti-I-A/I-E (MHCII)–AlexaFluor647 (M5/114.15.2) (BioLegend, 107618, 1:150); (Panel 9) rat anti-CD45-AlexaFluor488 (BioLegend, 103122, 1:150) and rat anti-P-cadherin-AlexaFluor647 (R&D Systems, FAB761R-100UG, 1:200).

Cell culture

For immunofluorescence experiments, feeders were split onto poly-l-lysine coated glass coverslips, seeded in 12-well plates 24 h prior to the addition of HFSCs. HFSCs were grown to confluency before feeders were detached by repeated PBS washes, and corpse or corpse-conditioned medium experiments were performed. At the end of the experiment, cells were washed twice with PBS and prefixed with 4% PFA for 3 min at 25 °C. Cells were washed three times with PBS and stained as for tissue sections.

Microscopy

Images of Sox9-creER; Rosa26^Brainbow2.1 tissue was acquired using a Zen-software driven Zeiss LSM 780 inverted laser scanning confocal microscope and 20× air objective (NA = 0.8), a 40× water immersion objective (NA = 1.2), or a 63× oil immersion objective (NA = 1.4). To separate CFP, YFP, GFP, RFP and AlexaFluor647 fluorophores, excitation with specific laser lines (405, 440, 488, 514, 561, 594, and 633) and narrow wavelength emission cut-offs on 4 detectors were set up as follows: CFP (excitation 440 nm, emission 450 nm–490 nm), GFP (excitation 488 nm, emission 500 nm–515 nm), YFP (excitation 514 nm, emission 525 nm–570 nm), RFP (excitation 561 nm, emission 595 nm–620 nm), and AlexaFluor647 (excitation 633 nm, emission 650 nm–690 nm). Due to their well-separated excitation and emission spectra, GFP and AlexaFluor647 were acquired simultaneously on the same detector. Stacks with a 1-μm step were acquired. Confocal microscopy was performed in The Rockefeller University’s Bio-Imaging Resource Center, RRID: SCR_017791.

Images of Rxra^fl; Sox9-creER; Rosa26^YFP tissue stained using IBEX methodology was acquired on an inverted Dragonfly 202 spinning disk confocal system (Andor Technology Inc.) using the 40× oil immersion objective, a 40-μm pinhole and a Zyla camera. Four laser lines (405, 488, 561 and 625 nm) were used for near simultaneous excitation of DAPI, Alexa-448, RRX and Alexa-647 fluorophores. Tiled images with a 1-μm stack step were acquired using the Andor Fusion software (v 2.3). Images were stitched and aligned using DAPI as a fiducial with Imaris and the SimpleITK Image Registration Pipeline plug-in for Imaris.

Images of other cryosections were acquired using a Zen-software-driven Zeiss Axio Observer.Z1 epifluorescence/brightfield microscope with a Hamamatsu ORCA-ER camera, Axiocam350, and an ApoTome.2 slider (to reduce light scatter in z). Stacks with a 1-μm step were acquired. Apotome acquired images were processed via ‘Apotome Raw Convert’ function, and stitched (if necessary), in Zen software (v 3.1). Subsequent image processing was conducted in ImageJ (v. 2.9.0) and Imaris (v. 10.1) (Oxford Instruments) software. For presentation purposes, images were cropped and assembled in Adobe Illustrator.

Phenotyping corpse engulfment via microscopy

To quantify apoptotic cell clearance in tissue sections, confocal images were acquired with a 40× or 63× oil immersion objective, with a 1-μm step size, of tissue stained with either P-cadherin or KRT14 to mark cell boundaries, DAPI to mark nuclei, and TUNEL to label late-stages of cell death. On single z-plane images, dying cells were scored as engulfed when a small TUNEL+ apoptotic body was nestled inside the cell boundary of a cell with a healthy nucleus, and could be visualized as such across the consecutive z-stack. Apoptotic bodies were generally round and found against health nuclei, in accordance with electron microscopy images of wild-type hair follicle ORS cells. Unengulfed apoptotic cells were either large TUNEL+ signal that completely overlapped a condensed nucleus and occupied roughly 50–75% the area of a healthy cell, or were visualized as small, irregularly shaped TUNEL+ debris pushed to the cell boundary edges. TUNEL+ debris was slightly more prevalent at the dermal–epithelial junction or the epithelia directly adjacent to the hair shaft, but was also visible throughout the ORS at cell boundaries.

Electron microscopy

Dissected back skin was placed on thin paper towel for stability, and fixed in 2% glutaraldehyde, 4% PFA, and 2 mM CaCl₂ in 0.1 M sodium cacodylate buffer (pH 7.2) for 2 h at room temperature, postfixed in 1% osmium tetraoxide and processed for Epon embedding. Ultrathin sections of 60-65 nm were counterstained with uranyl acetate and lead citrate, before images were taken with a transmission electron microscope (Tecnai G2-12;FEI) equipped with a digital camera (AMT BioSprint29). Samples were processed and imaged at The Rockefeller Electron Microscopy Resource Center. The number of engulfed apoptotic corpses per stem cell was quantified via transmission electron micrographs.

Flow cytometry

To obtain single-cell suspensions for fluorescence activated cell sorting (FACS) at all stages of the hair cycle, back skin was excised, and the dermal side scraped with a dull scalpel to remove excess fat prior to incubation with 0.25% collagenase (Sigma-Aldrich) in warm PBS, dermal side down for 45–60 min at 37 °C with gentle rotation in a plastic petri dish. The dermal side was scraped gently with a dull scalpel to mechanically dissociate cells in the lower ORS and hair bulb (“dermal fraction”). The dermal fraction was only kept for late anagen and early-to-mid catagen samples, and was processed separately from the epidermal fraction. To collect the epidermal fraction, the skin was placed dermal side down in 0.25% trypsin-EDTA (Gibco) for 20–25 min at 37 °C with gentle rotation. The hairy side of the skin was scraped against the direction of hair growth with a dull scalpel to release cells in the upper hair follicle (including the hair follicle bulge stem and hair germ progenitor cells). For both dermal and epidermal fractions, the resulting cell suspensions were pipetted up and down with a 5 ml serological pipette for 5 min, before being quenched with FACS buffer (5% fetal bovine serum, FBS, in PBS). Plastic petri dishes were rinsed with 5 ml of FACS buffer 2–3 times, which was collected and added to the appropriate cell suspension. Suspensions were filtered through sequential 70-μm and 40-μm nylon filters (VWR), before being pelleted at 350g for 15 min at 4 °C. Cell pellets were resuspended in ice cold FACS buffer, re-filtered into FACS tubes, and incubated with primary antibodies for 20 min on ice. Secondary antibodies and LysoTracker DeepRed (Invitrogen, 1:4,000) were added directly to FACS tubes, and incubation continued for 10 min on ice. Samples were further diluted with FACS buffer plus DNase (Roche) to minimize cell clumping prior to sorting or analysis. For analysis of RXRα levels by FACS, cells were stained with cell-surface specific primary and secondary antibodies, before being fixed and processed using the BD Cytofix/Cytoperm kit following manufacturer’s instructions. Primary antibodies were used as follows: rat biotinylated anti-CD45 (30-F11) (eBioscience, 13-0451-82, 1:200), rat biotinylated anti-CD117 (2B8) (eBioscience, 13-1171-82, 1:200), rat biotinylated anti-CD140a (APA5) (eBioscience, 13-1401-82, 1:200), rat biotinylated anti-CD31 (390) (eBioscience, 13-0311-82, 1:200), rat anti CD34-FITC (RAM34) (eBioscience, 11-0341-82, 1:200), rat anti CD34–eFluor660 (RAM34) (eBioscience, 50-0341-82,1:200), rat anti ITGA6–PercpCy5.5 (GoH3) (BioLegend, 313617, 1:250), rat anti-Ly6A/E-APC-Cy7(BioLegend, 108125, 1:1,000), rabbit anti-RXRα (D6H10) (CST, 3085, 1:250), rat anti-Tyro3/Dtk-AlexaFluor700 (R&D Systems, FAB759N, 1:200), rat anti-Mertk-AlexaFluor700 (R&D Systems, FAB5912N, 1:200), and rat anti-Axl-AlexaFluor700 (R&D Systems, FAB8541N, 1:200). Secondary antibodies were used as follows: Strepavidin-PE-Cy7 (1:3,000) and donkey AlexaFluor 488 or AlexaFluor568 (1:500). Annexin V-AlexaFluor568 (Invitrogen, A13202, 1:100) and/or DAPI was used to identify apoptotic and dying cells, respectively. For FACS using annexin V, primary and secondary antibody staining was performed in annexin V Binding Buffer (10 mM HEPES, 140 mM NaCl, 2.5 mM CaCl₂, pH 7.4). For AldeFluor activity assay, manufacturer’s instructions were followed (ALDEFLUOR Kit 01700, StemCell Technologies) with the addition of a more specific ALDH1A inhibitor at 1 μM (673-A, R&D Systems 6934).

For FACS analysis, live cell suspensions from back skin were collected and analysed as described above. Alternatively, cultured HFSCs were trypsinized for 7–10 min (as for passaging the cell lines), and pelleted at 300g before resuspension, filtering and incubating with primary antibodies. A minimum of 20,000 HFSCs were analysed per sample using either a BD LSRII Flow Cytometer or a BD Fortessa Flow Cytometer (BD Bioscience). Representative sort schemes pertaining to analysis of TAM-family receptors, LysoTracker expression, RXRa expression and the AldeFluor assay can be found in Supplementary Figs. 3–5 in combination with Extended Data Figs. 4 and 8. For analysis of the Brainbow2.1 HFSCs, a LSRII specially equipped with a 445 nm laser was used to excite CFP, separately from YFP/GFP and RFP. Phagocytic HFSCs were scored as double positive (containing a corpse of one fluor inside a cell of another fluor) following stringent gating against doublets, and separately confirmed as engulfment events via immunofluorescence. A representative sort scheme is shown in Supplementary Fig. 1.

For HFSC isolation for single cell RNA-sequencing cells were sorted according to the scheme shown in Supplemental Fig. 2, with an 85-μm nozzle into 96-well PCR plates (Bio-Rad) containing 2 μl of lysis buffer (0.2% Triton X-100, 2 U μl RNaseOUT (Thermo Scientific), 0.25 μM oligo-dT30VN primer, 1:2 × 10⁶ diluted ERCC spike-in RNAs (Ambion)). For HFSC isolation for in vitro culture and bulk ATAC-sequencing, cells were sorted using a 70-μm nozzle into E300-Y medium and FACS buffer (Supplementary Fig. 3). Representative sort schemes pertaining to Sox9-creER;R26^DTA ectopic corpse response and RXRα overexpression and knockout in HFSCs, are available in Supplementary Figs. 4 and 5, respectively. To isolate primary HFSCs for culture from the Sox9-creER:mTmG^fl/+ treated or not with tamoxifen as described previously, we gated on mTomato⁺ or mGFP⁺ in combination with CD34⁺ITGA6⁺ (Supplementary Fig. 6), using a 70-μm nozzle to sort into FACS buffer. Similar methodology was used to isolate Rxra HFSC lines for culture as described in Supplemental Fig. 3, after which Sox9-creER was activated with 4-OH-tamoxifen in culture. For bulk RNA-sequencing, cells were sorted using a 70-μm nozzle directly into Trizol and FACS buffer (Supplementary Fig. 7). Sorting was performed on a BD FACSAriaII equipped with Diva software (v. 8.0) (BD Biosciences).

Flow cytometry was performed at The Rockefeller University’s Flow Cytometry Resource Center (RRID: SCR_017694). Flow cytometry plots were generated using FlowJo to illustrate the strategies used for cell isolation, and manually compensated for presentation.

scRNA-sequencing libraries

Single-cell RNA-sequencing libraries were prepared from FACS-isolated hair follicle epithelial cells in AnaVI, CatVI, and CatVII, using a slightly modified Smart-Seq2 protocol as previously described^62,63. For each hair cycle stage, cells from 3–6 mice were pooled prior to FACS isolation. In brief, cells were sorted into hypotonic lysis buffer, snap frozen in liquid nitrogen and stored at −80 °C until all samples were collected. Cells were lysed by heating at 72 °C for 3 min, followed by reverse transcription of mRNA using dT30 oligonucleotides, template switching oligonucleotides and Maxima H- reverse transcriptase. The whole transcriptome was amplified (15 cycles) by KAPA HiFi DNA polymerase (Roche), and then size-selected using 0.6× AmpPure XP beads (Beckman Coulter). To exclude cells with poor amplification, and wells containing multiple cells, quantitative PCR (qPCR) for Gapdh was performed. Illumina sequencing libraries were indexed with unique 5′ and 3′ barcode combinations (up to 384 cells) using the Nextera XT DNA library preparation kit (Illumina). Libraries were pooled and size-selected with 0.9× AmpPure XP beads. Prior to sequencing on Illumina NextSeq500 using a 75 bp paired-end read mid-output setting, library quality was assessed by TapeStation (Agilent).

ATAC-sequencing libraries

ATAC-seq was performed on 20,000–75,000 (in vivo samples) or 50,000 (culture samples) FACS-sorted HFSCs, as previously described^64,65,66. In brief, cells were lysed in ATAC lysis buffer for 1 min on ice, washed and nuclei resuspended in transposase buffer. Genomic DNA was transposed using Tn5 transposase (Illumina) for 30 min at 37 °C, at which point the reaction was halted. Samples were uniquely barcoded in batches of 10-12 (in vivo samples) or one batch of 27 (cultured HFSCs), using Buenrostro⁶⁴ or Nextera XT index kit v2 indices. Sequencing libraries were prepared according to manufacturer’s instructions (Illumina). Libraries were sequenced to a depth of 50–100 million sequences, using paired-end runs on an Illumina Novaseq 6000 (at The Rockefeller University Genomics Resource Center).

CUT&RUN sequencing libraries

CUT&RUN sequencing was performed on 500,000 cultured HFSCs, as previously described with minor modifications^67,68. Unless otherwise indicated, steps were performed at room temperature. In addition to biological samples, antibody validation and specificity was verified using (1) a rabbit IgG control antibody; and (2) rabbit anti-RXRα in the Rxra-cKO HFSCs. In brief, cells were trypsinized as for re-plating, then washed with PBS and resuspended in crosslinking buffer (10 mM HEPES–NaOH pH 7.5, 100 mM NaCl, 1 mM EGTA, 1 mM EDTA and 1% formaldehyde) with rotation for 10 min. Crosslinked cells were quenched with 0.125 M (final concentration) glycine for 5 min, then washed with ice cold 1× PBS and resuspended in NE1 buffer (20 mM HEPES–KOH pH 7.9, 10 mM KCl, 1 mM MgCl₂, 1 mM dithiothreitol, 0.1% Triton X-100 supplemented with Roche complete protease inhibitor EDTA-free) and rotated for 10 min at 4 °C. Nuclei were washed twice with CNR wash buffer (20 mM HEPES pH 7.5, 150 mM NaCl, 0.5% bovine serum albumin and 0.5 mM spermidine supplemented with protease inhibitor) and incubated with concanavalin-A (ConA) beads washed with CNR binding buffer (20 mM HEPES–KOH pH 7.9, 10 mM KCl, 1 mM CaCl₂ and 1 mM MnCl₂) for 10 min at 4 °C. ConA-bead bound nuclei were incubated overnight at 4 °C in CNR antibody buffer (CNR wash buffer supplemented with 0.1% Triton X-100 and 2 mM EDTA) and 1:50 RXRα antibody (Cell Signaling Technologies, clone D6H10, 3085). ConA-bead bound nuclei were washed with CNR Triton wash buffer (CUT&RUN wash buffer supplemented with 0.1% Triton X-100) then resuspended and incubated at 4 °C for 60 min in CUT&RUN antibody buffer and 2.5 μl pAG-MNase (EpiCypher). Following this, ConA-bead bound nuclei were washed twice with CUT&RUN Triton wash buffer, resuspended in 100 μl of Triton wash buffer and incubated on ice for 5 min before 2 μl of 100 mM CaCl₂ was added per sample. Samples were incubated on ice for 30 min and the reaction was then stopped by adding 100 μl of 2× stop buffer (340 mM NaCl, 20 mM EDTA, 4 mM egtazic acid, 0.1% Triton X-100 and 50 μg ml^−1 RNaseA) and incubated at 37 °C for 10 min. ConA-bound nuclei were captured on a magnet, and supernatant containing Cut-and-Run DNA fragments was collected. Supernatant was incubated at 70 °C for 4 h with 2 μl 10% sodium dodecyl sulfate and 2.5 μl 20 mg ml^−1 proteinase K, prior to DNA purification using PCI reagent (phenol:chloroform:isoamyl alcohol, Millipore). DNA fragments were precipitated overnight with ethanol and glycogen at −20 °C before resuspension in elution buffer (1 mM Tris–HCl pH 8.0 and 0.1 mM EDTA).

CNR sequencing libraries were generated using NEBNext Ultra II DNA Library Prep Kit for Illumina and NEBNext Multiplex Oligos for Illumina. PCR-amplified libraries were purified using 1× ratio of SPRI beads (Beckman) and eluted in 15 μl EB buffer (Qiagen). All CNR libraries were sequenced on Illumina NextSeq using 40 bp paired-end reads.

RNA isolation

Total RNA was isolated from FACS-isolated HFSCs using the Direct-zol RNA MicroPrep kit (Zymo Research) following manufacturer’s instructions. The optional DNase I treatment was included in all sample preps, and RNA was eluted in DNase/RNase-free water. Quality and concentration of RNA samples were determined using an Agilent 2100 Bioanalyzer. All samples for sequencing had RNA integrity (RIN) numbers >8.5. RNA samples were used for qPCR with reverse transcription (RT–qPCR) or bulk RNA sequencing, as described.

Bulk RNA-sequencing libraries

Comparable amounts of RNA per sample were used to prepare bulk RNA-sequencing libraries using Illumina Trueseq standard mRNA library kit (non-stranded, poly-A selection) following manufacturer’s guidelines. Libraries were then uniquely barcoded, pooled and sequenced on an Illumina Novaseq 6000 using single-end runs (at Weill Cornell Medical College’s Genomic Core Facility).

RT–qPCR

Equivalent amounts of RNA were reverse transcribed using SuperScript III Reverse Transcriptase (Thermo Fisher Scientific) following manufacturer’s instructions. To normalize cDNA amount across samples, primers for B2m were used. cDNAs were mixed with gene specific primers (Supplementary Table 6) and SYBR green PCR MasterMix (Sigma Aldrich) and run on an Applied Biosystems 7900HT Fast Real-Time PCR system.

Single cell and bulk RNA-sequencing analysis

Trimmed FASTQ files were obtained from the Rockefeller University’s Genome Resource Center (scRNA-sequencing, this study), or from the Gene Expression Omnibus (GSE90848 and GSE130850 for previously published telogen and AnaI-II HFSC scRNA-sequencing datasets), or from the Genomic Core Facility (Weill Cornell Medical College; bulk RNA-sequencing), and raw sequencing reads were aligned to the mouse reference genome (UCSC release mm39) using STAR (v2.6)⁶⁹. The expression values of each gene were quantified as both raw counts and transcripts per million (TPM) using Salmon (v.1.4.0)⁷⁰, and compiled in R (v.3.6.1) using RStudio (v.3.4.2) by Tximport (v.1.12.3)⁷¹.

Bulk RNA sequencing

For differential gene expression analysis in R, low detection genes (minimum average read count <10) were filtered before DESeq2 analysis (v.1.16.1)⁷². Differential expression modelling used a negative binomial distribution and Wald test. Genes were differentially expressed for log₂[fold-change]>|1| and adjusted P < 0.05. Heat maps and bar graphs illustrating differential gene expression were constructed in a Python environment (detailed in next paragraph).

scRNA-sequencing

Analysis and visualization of the data were conducted in a Python environment built on Pandas (v.2.0.1), NumPy (v.1.24.2)⁷³, SciPy (v.1.10.1)⁷⁴, scikit-learn (v.1.2.0), SCANPY (v1.9.3)⁷⁵, AnnData (v.0.9.1)⁷⁵, matplotlib (v.3.7.1)⁷⁶ and seaborn (v.0.13.1)⁷⁷ packages. Raw count and metadata matrices for 1,489 single ORS cells across the hair cycle were loaded in SCANPY as an AnnData object. Single cell data was preprocessed to remove lowly detected genes (expressed in <75 cells) and cells with low complexity libraries (<2,000 genes detected). SCANPY was used to normalize counts per cell, and highly variable genes were detected. Prior to dimensionality reduction by principal component analysis (PCA), data were centred and scaled. PCA was performed on highly variable genes, with 100 components and the svd_solver using ‘arpack’ (SCANPY default setting). To construct a k-nearest neighbours graph on Euclidean distance, 41 principal components were used (which captured 25% of the variance in the data). Data was visualized using UMAP in SCANPY, and clustering was done using the Leiden algorithm (with a resolution setting of 0.5). Cluster resolution was chosen after iterating through resolution parameters from 0.1 to 0.75, as best capturing both hair follicle cycle stages and anatomic location (upper bulge region/upper ORS versus hair germ/upper-middle ORS versus lower ORS). Marker gene expression based on the literature^7,11,63,78, together with the FACS markers each population was sorted on, was used to identify clusters. SCANPY was used to visualize selected marker genes in dot plots, or as normalized counts visualized on UMAPs.

Differential gene expression based on cluster identity was used in DESeq2 to identify genes that varied as cells transitioned from late anagen growth phase to catagen. Differential expression was performed as described for bulk RNA-sequencing, with the modification of a threshold of 0.75 to construct Wald tests of significance. Gene set enrichment analysis (GSEA) on differentially expressed genes was performed using GSEA software (v.4.3.2)^79,80, and run with the MSigDB 2022 mouse database. Gene set terms with false discovery rate < 0.1 and showing high normalized enrichment scores in catagen cells were considered interesting. To construct gene set scores based on the GSEA identified terms the corresponding Mus musculus gene lists were obtained by Amigo2 through the Gene Ontology consortium. The SCANPY tl.score_genes function was used to compute the average expression of each gene set across single cells, and normalized to a randomly sampled reference set of genes^81,82. The resulting gene set scores were colour coded on corresponding UMAP visualizations of the data.

ATAC-seq analysis

Trimmed FASTQ files were obtained from the Rockefeller University’s Genome Resource Center and aligned to the mouse reference genome (UCSC release mm39) using Burrows-Wheeler Aligner (BWA, v.0.7.18), using BWA-MEM with default parameters. The output.sam files were name-sorted and duplicate reads were marked and removed using SAMtools (v.1.17)⁸³. Peaks were called on each replicate using MACS3 (v.3.0.0) using the callpeak command, BAMPE, and a mappable genome estimate of 1.87 × 10⁹ (from the ENCODE pipeline). The fraction of reads in peaks was calculated using bedTools (v. 2.31.0)⁸⁴ and used to scale bigwig files equivalently in deepTools (v.2.0.0)⁸⁵. Bigwig files were created from deduplicated, pooled replicate bam files using deepTools, and normalized as reads per genome coverage. Pooled replicate bigwig files were also used to calculate peak coverage matrices to plot heatmaps of centred differential peaks, extended by 1 kb upstream and downstream. Differential peak analysis was done in DESeq2, using read count matrices across each individual replicate from concatenated, merged union peak sets from each replicate. These union peak sets were created separately for in vivo samples and in vitro samples. Differential analysis used negative binomial modelling, and Wald’s test for significance. To assign peaks to nearest expressed gene, part of the Inferelator-prior (v.0.3.8)⁸⁶ package was used. Peaks were assigned to genes if they fell 50 kb upstream or 5 kb downstream of the gene body and were curated for expression using either scRNA-seq (in vivo samples) or bulk RNA-seq (in vitro samples). To make sure that all potential enhancers for genes related to the apoptotic cell clearance programme were identified, any unassigned intergenic peaks within approximately 200 kb of phagocytosis-related genes were manually curated. If no genes were expressed transcriptionally in the interval between phagocytic gene and unassigned intergenic peak, the intergenic peak was considered a potential enhancer for said gene. Peaks of interest were visualized using the integrated genome viewer (IGV) software (v.2.13.2), together with.bed files of differential peaks.

Motif enrichment analysis for in vivo samples was performed in two ways: First, the MEME suite (v. 5.5.2) package XSTREME⁸⁷ in web browser format was used to search for motifs enriched in differential peaks, using as background the union set of all peaks detected, and the JASPAR 2022 vertebrate CORE transcription factor motif database, with lengths of 6–18 bp specified. Both known and de novo enriched motifs were collapsed to clusters based on similarity and ranked based on adjusted P value. Second, the transcription factor occupancy prediction by investigation of ATAC-seq signal (TOBIAS, v.0.14.0)⁸⁸ framework was used to perform chromatin footprinting analysis. In brief, replicate-pooled bam files read coverage across the genome was calculated and corrected for Tn5 transposase cutting bias before footprint scores were calculated within the union set of called peaks. TOBIAS footprint scores were used to compute differential binding between anagen and catagen pooled replicates, or between Rxra wild-type and cKO pooled replicates. RXR-family catagen bound footprints were visualized in IGV by pooling each individual RXR-family member’s bed footprint file.

CUT&RUN sequencing analysis

Trimmed FASTQ files were obtained from the Rockefeller University’s Genome Resource Center and aligned to the mouse reference genome (UCSC release mm39) using Burrows-Wheeler Aligner (BWA), using BWA-MEM with default parameters. The output.sam files were name-sorted and duplicate reads were marked and removed using SAMtools (v.1.17)⁸³ Reads were filtered to less than 121 bp using SAMtools (v.1.3.1). BAM files for each replicate were combined using Samtools. Bigwig files were generated using Deeptools (v.3.1.2) with reads per kilobase of transcript per million mapped reads (RPKM) normalization and presented with Integrative Genomics Viewer software. CNR peaks were called using SEACR (v.1.3)⁸⁹ from bedGraph files generated from RPKM-normalized Bigwig files (bigWigToBedGraph, UCSC Tools) using stringent setting and a numeric threshold of 0.01.

Statistics and reproducibility

All data from every experiment were included for analysis unless an error was detected via failed positive or negative controls; in that case the entire experiment was excluded from analysis. Measurements were taken from independent distinct samples, unless stated otherwise. Statistical methods were not used to predetermine sample size. Experiments were not randomized or blinded, given the lack of ambiguity in phenotypes observed and internal controls used.

Statistical and graphical analyses were performed in Jupyter Notebooks, running a custom Python environment built as described in the single cell sequencing analysis section. Sample sizes, replicates and statistical tests used are indicated in each figure legend. Unless otherwise stated, unpaired two-tailed Student’s t-tests with a 95% confidence interval were performed to test for pair-wise differences among the means. Data are visualized as box-and-whisker plots, with the box representing the first to third quartiles of the data set, the median line inside the box, and the whiskers extending a maximum of 1.5 times the inter-quartile range. Observations that fall outside this range are plotted independently. For clarity, each observation in a data set is also visualized as a point overlaid on the box plot. Whenever representative plots or images are shown, data sets with similar results were generated from additional n > 3 independent biological replicates, from separate litters of mice or two independent cell culture experiments from separate days. All attempts at replication in this study were successful. In general, experiments were not randomized or performed in a blinded manner, due to the complex genetic models and obvious phenotypic differences in samples.

Reporting summary

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