Selection of target genes and sgRNA design
Five hundred sgRNAs were selected as initial targeting library. Fifty out of the 500 sgRNAs were non-targeting control guide RNAs. Three sgRNAs per gene targeted the top 150 most frequently mutated or CNV genes in human HNSCC and sSCC. We selected genes that have a mutation frequency of ≥6% in HNSCC, ≥30% in sSCC and ≥15% of copy number alterations in both. Mutation frequency was assessed using the cBioPortal for Cancer Genomics51,52 with the datasets for sSCC53,54 and HNSCC55,56,57,58,59. The exact selection of genes and their mutation frequency can be found in Extended Data Fig. 1. The overlap and differences between HNSCC and sSCC genes is graphically represented in Fig. 1b (Venn diagram). Individual guides were designed with the Broad Institute sgRNA design tool (sgRNA Designer: CRISPRko (https://portals.broadinstitute.org/gpp/public/analysis-tools/sgrna-design)). The three top-ranking guides per target were chosen. Fifty non-targeting control sgRNAs were randomly picked from the Brie library60 (Addgene #73633). All 500 sgRNA sequences are listed in Supplementary Table 9.
CROP-mCherry vector and library cloning
The puromycin resistance cassette in the original CROP-seq-Guide-Puro vector (Addgene #86708) was replaced with a mCherry sequence, amplified from pAAVS1-NDi-CRISPRi (Gen1) (Addgene #73497) and cloned in via Pfl23II and MluI restriction sites for FACS. The 500 sgRNA sequences were ordered as an oligonucleotide pool from IDT, cloned in batch via Gibson assembly or via BsmbI restriction sites for sgRNAs as previously described5. For all single-cell experiments, three different library batches were cloned and individually sequenced to confirm homogeneous sgRNA representation.
High-titre lentivirus production
Production of vesicular stomatitis virus (VSV-G) pseudotyped lentivirus was performed by calcium phosphate transfection of Lenti-X 293 T cells (TaKaRa Clontech, 632180) with CROP-mCherry and helper plasmids pMD2.G and psPAX2 (Addgene plasmids 12259 and 12260). Sixteen hours after transfection, medium was changed to viral production medium (DMEM (Gibco 11965092), 1% penicillin-streptomycin-glutamine (Gibco 10378016), 1% 100 mM sodium pyruvate (Gibco 11360070), 1% sodium bicarbonate 7.5% solution (Gibco 25080094), 5 mM sodium butyrate (Sigma-Aldrich B5887)) and incubated at 37 °C/7.5% CO2. Viral supernatant was collected 46 h after transfection and filtered through a 0.45-μm filter (Millipore Stericup Quick Release Durapore S2HVU02RE). For in vivo lentiviral transduction, the viral supernatant was concentrated ~2,000-fold using a 100 kDa MW cut-off Millipore Centricon 70 Plus (Merck Millipore; UFC710008). The virus was further concentrated by ultracentrifugation using the SW 55 Ti rotor for the Beckman Coulter Optima L-90 Ultracentrifuge at 45,000 rpm and 4 °C. Final viral particles were resuspended in viral resuspension buffer (20 mM Tris pH 8.0, 250 mM NaCl, 10 mM MgCl2, 5% sorbitol) and stored at −80 °C until used for titration or injection. For low-titre virus, production of lentivirus was performed similarly as described above, and viral supernatant was collected and filtered using a 0.45-μm syringe filter (Sarstedt 83.1826) and stored at −80 °C until used.
In vivo experiments
All animal experiments were conducted in strict accordance with the Swiss Animal Protection law and requirements of the Swiss Federal Office of Food Safety and Animal Welfare (BLV). The Animal Welfare Committee of the Canton of Zurich approved all animal protocols and experiments performed in this study (animal permits ZH074/2019, ZH196/2022).
Genetically modified mice of the strain Tg(B6J.129(Cg)-Gt(ROSA)26Sortm1.1(CAG-cas9*,-EGFP)Fezh/J (denoted as B6.Cas9) were purchased from the Jackson Laboratory (strain 026179). B6.129 S Tnfrsf1atm1ImxTnfrsf1btm1Imx (denoted as Tnfrsf1a-KO, originally from Jackson Laboratory 003243) were acquired through the Swiss Immunology Mouse Repository (SwImMR). Wild-type mice of the strain CD1-IGS (denoted as CD1) were purchased from Charles River. B6.Cas9 were crossbred either with Tnfrsf1a-knockout or CD1 mice to provide embryos heterozygous for the Cas9 allele, suitable for lentivirus injection at embryonic development stage E9.5.
Ultrasound-guided in utero injections were conducted as previously described61. In brief, females at E9.5 of gestation were anaesthetized with isoflurane. Each embryo was injected with 0.5 μl of high-titre lentivirus, and up to 8 embryos were injected per litter. Surgical procedures were limited to 30 min to ensure fast recovery.
All mice were housed at the Laboratory Animal Services Center (LASC) of the University of Zurich in individually-ventilated cages in a humidity- and light-controlled environment (22 °C, 45–50%, 12 h light/dark cycle), and had access to food and water ad libitum. B6.Cas9 males used for crossbreeding and pregnant CD1 females were housed individually. All other mice were group-housed. Successful infection of mice was controlled by eye through excitation with a dual fluorescent protein flashlight (NIGHTSEA, DFP-1). Positively infected mice were then euthanized by decapitation (P4) or with CO2 (P60), shaved and treated with hair removal cream if necessary.
P4 back skin was processed as previously described8, with minor modifications. After being surgically removed and scraped for fat, the back skin was washed once in cold PBS and then placed in dispase (Corning 354235) for 35 min at 37 °C, dermis-side facing down, on an orbital shaker. Next, epidermis was separated from dermis with fine forceps, torn into smaller pieces and placed in 4 ml 0.25% Trypsin-EDTA (1×, Gibco; 25200056) for 15 min at 37 °C with orbital shaking. Epidermis was washed with cold PBS and pipetted vigorously to achieve a single-cell suspension. Suspension was then afterwards filtered through a 70-µm and a 40-µm strainer (Corning; 431750, 431751) consecutively, then centrifuged for 10 min at 400g and resuspended in FACS buffer (PBS + 2% chelexed FBS (FBS(-)) + DAPI).
P60 back skin was processed as previously described62. In brief, back skin was harvested with surgical scissors and placed dermis-facing up on a styrofoam tray with pins. Fat and subcutis was scraped off with a scalpel. Fat-free skin was then washed, dermis-side facing down in 1× PBS. Skin was then placed in the same orientation in 0.5% Trypsin-EDTA (10×, Gibco; 15400054) and incubated at 37 °C on an orbital shaker for 25 (females) or 50 min (males). Using a glass microscopy slide, skin was secured to the bottom of the dish and scraped again with a scalpel until it started to break. Excess cold PBS with 2% FBS(-) was added to neutralize trypsin. Cell suspension was strained on ice first through a 70-µm strainer, washed with an additional 15 ml of cold 1× PBS with 2% FBS(-) and then filtered through a 40-µm strainer. Cell suspension was spun down at 400g for 10 min and resuspended in FACS buffer.
P4: a total of 58 injected P4 mice were collected and assessed in 7 independent single-cell RNA-seq runs. P60: a total of 10 injected mice were collected and assessed in 8 independent single-cell RNA-seq runs. Because the P4 skin contains less cells, more mice had to be pooled to sort a sufficient number of mCherry-positive cells.
mCherry-positive cells were sorted on a BD FACSAria III using a 70-µm nozzle and processed for single-cell capture with a BD Rhapsody Single-Cell Analysis System using BD Rhapsody Cartridge Kit (633733). For P4, 7 separate cartridges with an average of 50,037 cells were loaded. For P60, 8 cartridges with an average of 64,625 cells were loaded. Reverse transcription from beads and sequencing library production was carried out according to the manufacturer’s instructions (BD Biosciences, doc ID: 210967 rev. 1.0).
Dial-out nested PCRs for sgRNA amplification
The sgRNA-containing region was amplified from BD Rhapsody beads in a separate nested PCR reaction. In PCR1, we used 100 µl KAPA HiFi HotStart ReadyMix (Roche 07958935001), 6 µl forward primer (5′-ACACGACGCTCTTCCGATCT-3′, 10 μM), 6 µl reverse primer (5′-TCTTGTGGAAAGGACGA-3′, 10 μM), 12 μl Bead RT/PCR Enhancer reagent from BD Biosciences and 72 μl nuclease-free water for a total volume of 200 μl. Rhapsody beads from each separate cartridge preparation were resuspended in the PCR mix, aliquoted in 4 × 50 μl reactions and quickly moved into a pre-heated PCR machine without allowing beads to settle. The following conditions were used: Initial denaturation at 95 °C for 5 min, followed by 25 cycles of denaturation 95 °C for 30 s, annealing 53 °C for 30 s, extension 72 °C for 20 s and a final extension of 10 min at 72 °C.
PCR1 product was pooled, beads were removed magnetically and amplicons were cleaned up with Agencourt AMPure XP beads (Beckman Coulter A63881) according to the manufacturer’s instructions before the second PCR. For PCR 2, 3 μl of cleaned up PCR1 was used as template along with 2 μl forward primer (5′-ACACGACGCTCTTCCGATCT-3′, 10 μM), 2 μl reverse primer (CAGACGTGTGCTCTTCCGATCTCTTGTGGAAAGGACGAAACA*C*C*G-3′, 10 μM), 18 μl nuclease-free water, along with 25 μl of KAPA HiFi HotStart ReadyMix for a total volume of 50 μl. PCR2 was performed under the following conditions: Initial denaturation at 95 °C for 3 min, followed by 10 cycles of denaturation 95 °C for 30 s, annealing 60 °C for 3 min, extension 72 °C for 60 s, and a final extension of 5 min at 72 °C. PCR2 product was cleaned up with Agencourt AMPure XP beads according to the manufacturer’s instructions before the indexing PCR. Indexing PCR was performed according to BD Biosciences ‘mRNA Targeted Library Preparation’ protocol (doc ID: 210968 rev. 3.0).
sgRNA enrichment and depletion calculation
To calculate enrichment and depletion of sgRNAs, we used sgRNA cell counts for P4 and P60. For T0, read counts from the sequenced pre-injection libraries A, B and C originating from separate batch plasmid preparations were used. As each Px (P60 or P4) sample received one of the three T0 library batches, we normalized each sample individually as follows. Counts at Px and T0 were first transformed to proportions over the total counts, either of each entire sample or library batch. Then, each Px sample guide proportions were normalized by dividing against its T0 proportion (that is, its initial representation in the matching pre-injection library). This yields a fold change (of ratios). Fold changes from triplets of guides targeting the same gene were then averaged together (and indicated with the symbol of the targeted gene), and so were the 50 control guides (grouped for a total of 17 pseudo-triplet controls: ctrl_1 for control sgRNA from 1 to 3, ctrl_2 for control sgRNA from 4 to 6, etc, and ctrl_17 which comprised only control sgRNA 49 and 50). This resulted in 167 fold change values in each Px sample: 150 for the target genes and 17 for the pseudo-triplet controls. Ranks were defined for each Px sample based on these fold changes in decreasing order, therefore with the highest fold change ranking first. Ties (typically, fold changes with value zero) were assigned to the same lowest rank (see min_rank function from R package dplyr). These ranks were used to compute the correlation of the top 20 and bottom 20 perturbations (guide triplets) among the eight P60 samples (biological replicates), as shown in Extended Data Fig. 1h.
Finally, the fold changes obtained were also averaged across the same-time point samples (either P4 or P60). For this average, a weighted mean was used, with the weights being the total number of guide-positive cells for each sample. These results are reported in log2 scale in Fig. 2a and in Extended Data Figs. 1i and 3. For the P60 versus P4 guide representation (Extended Data Fig. 3c), the fold changes shown in Extended Data Fig. 3a were divided by those in Extended Data Fig. 3b (by matched perturbation identity) and re-ranked. Epidermal stem cell-specific sgRNA representations (as shown in Extended Data Fig. 3d) were computed with the same approach, after subsetting the P4 or P60 datasets to include only epidermal stem cells (cluster 0, in Fig. 1e,h).
DMBA/TPA chemical carcinogenesis
Two-stage cutaneous chemical carcinogenesis was performed as previously described63. In brief, the back skin of 6- to 8-week-old female mice were shaved with electric hair clippers and treated once with 400 nmol DMBA (Sigma-Aldrich D3254) per 100 µl (dissolved in acetone) as an initiator mutagen. After a 2-week rest period, 40 nmol TPA (Sigma-Aldrich 79346) dissolved in 100 µl 100% ethanol was applied twice a week. Tumour formation was monitored regularly by visual inspection and mice were euthanized after 12 weeks of treatment. The maximum cumulative tumour size of 2 cm, as permitted by the Animal Welfare Committee of the Canton of Zurich (animal permits ZH074/2019, ZH196/2022), was not exceeded.
Tumour preparation
Tumours and 1–2 mm of surrounding tissue from DMBA/TPA treated mice were excised and finely minced with a scalpel. Tumour pieces were submerged in 5 ml pre-warmed DMEM with 0.25% Trypsin/EDTA and 3 U ml−1 DNAse I at 37 °C on an orbital shaker with rigorous resuspension every 5 min for a total of 30 min. The reaction was quenched with 1 ml FBS and strained through a 70-µM mesh, pelleted at 400g for 10 min with excess PBS and resuspended in 0.5 ml FACS buffer before sorting. Single-Cell capture was performed with BD Rhapsody as described above using Rhapsody Enhanced Cartridge Kit (BD Biosciences 664887). Reverse transcription from beads and sequencing library production was carried out according to manufacturer’s instructions (BD Biosciences, doc ID: 210967 rev. 1.0).
Tumour preparation for sgRNA amplification
DMBA/TPA-induced tumours that were not used for either scRNA-seq or spatial analysis were minced. Genomic DNA was extracted using the Qiagen DNeasy Blood and Tissue Kit. sgRNA region was amplified in a nested PCR. In PCR1, we used 12.5 µl KAPA HiFi HotStart ReadyMix (Roche 07958935001), 0.75 µl forward primer (5′-CTTGTGGAAAGGACGAAACACCG-3′, 10 μM), 0.75 µl reverse primer (5′-GTGTCTCAAGATCTAGTTACGCCAAGC-3′, 10 μM), 1 µl of extracted genomic tumour DNA (with a concentration between 50–100 ng μl−1) and 10 μl nuclease-free water for a total volume of 25 μl per tumour. PCR1 was performed using the following conditions: Initial denaturation at 98 °C for 2:30 min, followed by 22 cycles of denaturation 98 °C for 20 s, annealing 62 °C for 30 s, extension 72 °C for 30 s. PCR1 products for each tumour were individually column-purified with a GenElute PCR Clean-Up Kit (Sigma-Aldrich NA1020-1KT). For PCR2, 1 μl of cleaned up PCR1 was used as template along with 0.75 μl forward barcoded i5 primer (see Supplementary Table 1, 10 μM), 0.75 μl reverse barcoded i7 primer (see Supplementary Table 1, 10 μM), 10 μl nuclease-free water, along with 12.5 μl of KAPA HiFi HotStart ReadyMix for a total volume of 25 μl. Eight different i5 and i7 primers allowed for a total of 64 different combinations. PCR2 was performed under following conditions: initial denaturation at 95 °C for 3 min, followed by 32 cycles of denaturation 98 °C for 30 s, annealing 62 °C for 30 s, extension 72 °C for 30 s. PCR2 products for a maximum of 64 individually barcoded tumour amplifications were pooled and cleaned up with Agencourt AMPure XP beads (Beckman Coulter A63881) according to manufacturer’s instructions before sequencing.
Single-cell RNA sequencing
With the exception for whole-transcriptome amplification analysis 7 (sent to Novogene and sequenced on a S4 flow cell of an Illumina Novaseq), all single-cell or amplicon sequencing was prepared as ready-made libraries and sequenced either on Illumina NovaSeq or Nextseq 500 instruments at the Functional Genomics Center Zurich. Sequencing libraries were checked for peak size and concentration on a Bioanalyzer 2100 using DNA High Sensitivity Chip. Concentrations were additionally measured using a Qubit Fluorometer. Whole-transcriptome amplification analyses from Rhapsody Cartridges were sequenced in full SP flow cells at 1.8 nM concentration with 20% PhiX and the following read configuration: Read 1 60, I1 8, Read 2 62 cycles. Tumours were sequenced with 10% PhiX, dual indexing i5, 8 bp; i7, 8 bp; and Read 1, 64 cycles.
Bulk RNA sequencing and processing
mCherry-positive or negative keratinocytes were FAC-sorted on a BD FACSAria III directly into TriZol LS and column-purified with a Direct-zol RNA Miniprep kit (Zymo Research). Sequencing libraries from 1,000 ng of RNA per sample were prepared with the TrusSeq mRNA library kit and sequenced on a NovaSeq X 2 × 150 bp at the Functional Genomics Center Zurich. Quality control of the paired-end bulk RNA-seq data was performed using the FastQC program. Adapter trimming of the Illumina Universal Adapter (AGATCGGAAGAG) was carried out by Cutadapt. Transcript level quantification was performed by Salmon64 using the Gencode Mouse M25 reference transcriptome. DESeq2 was utilized to identify the differentially expressed genes between different conditions65. Initial quality control of the count data was carried out by PCA and hierarchical clustering using top 500 variable genes after variance stabilizing normalization. The gene level count data were normalized by the median of ratios method and dispersion was estimated. Next a generalized linear model was fitted to identify differentially expressed genes. The P values were obtained by the Wald test and corrected for FDR using the Benjamini–Hochberg procedure. FDR cut-off <0.05 was used to detect significant differentially expressed genes for subsequent analysis.
Immunofluorescence
Back skin from P4 and P60 mice was scraped to remove fat, placed on Whatman paper and cut into strips. Sections were embedded in OCT (Tissue-Tek; 4583), frozen by placing them in a liquid nitrogen-cooled isopentane filled metal beaker, stored at −80 °C and sectioned at 10–12 μm thickness on a Leica CM1900. Sections were immobilized on Superfrost glass slides, fixed for 10 min at room temperature with 4% paraformaldehyde (PFA), washed twice in 1× PBS before blocking in blocking solution (1% BSA, 1% gelatin, 2.5% normal goat serum, 2.5% normal donkey serum, 0.30% Triton X-100, 1× PBS). Primary antibodies (TNF [D2D4] XP[R] Rabbit anti-mouse monoclonal antibody, CST 11948; 1:300) and purified rat anti-mouse CD104 (also known as β4-integrin) (BD Pharmingen 553745, 1:300) were incubated overnight at 4 °C. After washing twice with 1× PBS, slides were incubated with secondary antibodies (Alexafluor488, Cy3, or AlexaFluor647, Jackson ImmunoResearch Laboratory; 1:500–1:1,000) and 0.5 µg ml−1 4′,6-diamidino-2-phenylindole (DAPI) at room temperature for 1 h. Sections were then washed again with 1× PBS, dried, covered with home-made mounting medium and sealed with nail polish before image acquisition. Pictures were acquired with a 20× objective on a ZEISS Axio Observer microscope controlled by the ZEN microscopy software (version 3.1).
Whole-mount immunofluorescence and antibodies
Dissected back skin from P4 and P60 mice injected with CROP-mCherry library were fixed in 4% PFA for 1 h at room temperature. Following fixation, samples were permeabilized in 0.8% PBS-Triton overnight. All steps during staining were carried out at room temperature. Primary antibodies were diluted into blocking buffer (5% donkey serum, 2.5% fish gelatin, 1% BSA, 0.8% Triton in PBS) and were incubated for at least 16–20 h at room temperature. Samples were then washed for 3–4 h in 0.8% PBS-Triton, and incubated with secondary antibodies (in blocking buffer) together with DAPI (to label nuclei, 0.25 mg ml−1) for at least 16–20 h. After staining, samples were extensively washed with 0.8% PBS-Triton every hour for 4–6 h. Samples were then dehydrated in increasing concentrations of Ethanol: 30%, 50% and 70% in doubled-distilled water (with the pH adjusted to 9.0 with NaOH/HCL) for 1 h each, and finally in 100% ethanol for 1 h. For tissue clearing, samples were transferred to Eppendorf tubes with 500 μl ethyl cinnamate (ECi, Sigma 112372) and shaken overnight at room temperature under dark conditions. Fresh ECi was used for mounting for imaging.
For fixed sections, back skin was cut into strips, embedded and frozen in OCT (Leica), and sectioned with a Leica cryostat (10-μm sections). Sections were permeabilized in 0.3% PBS-Triton for 15 min at room temperature, blocked in blocking buffer (5% donkey serum, 2.5% fish gelatin, 1% BSA, 0.3% Triton in PBS) for 15 min at room temperature. Primary antibodies in blocking buffer were incubated on slides overnight at 4 °C, washed off with 0.3% PBS-Trition before adding secondary antibodies, with DAPI for 1 h at room temperature. After washing off the secondary, samples were mounted for imaging using ProLong Diamond Antifade Mountant (P36961, Invitrogen).
Antibodies used were as follows: rat anti-RFP (Chromotek, 5F8; 1:200), rabbit anti-RFP (MBL, PM005; 1:200), chicken anti-GFP (Abcam, ab13970; 1:200), rat anti-CD45-biotin (Biolegend, 103104; 1:200), goat anti-TNFR1 (R&D, AF-425-PB; 1:200). All secondary antibodies used were raised in a donkey host and were conjugated to Alexa Fluor 488, Cy3 or Alexa Fluor 647 (Jackson ImmunoResearch Laboratory; 1:500- 1:1000). 4′,6-diamidino-2-phenylindole (DAPI) was used to label nuclei (0.25 mg ml−1).
Whole-mount microscopy
Whole-mount images were acquired using an LSM980/Airyscan module combined with Tiling (where indicated) 40× multi-immersion objective (LD LCI Plan-Apochromat 40×/1.2 Imm Corr DIC (silicone oil, glycerol or water immersion), WD 0.41 mm) with glycerine as an immersion medium for cleared samples, or a 40× oil immersion objective (Plan-Apochromat 40×/1.4 Oil DIC, WD 0.13 mm) for sections.
Whole-mount image processing and analysis
All processing was done in Zen Blue, with a custom batch processing macro (Thomas Peterbauer). All images depicted are maximum intensity projections.
In vivo EdU labelling
One hour before they were euthanized, mice were injected intraperitoneally (50 µg per g body weight) with 5-ethynyl-2′-deoxyuridine (EdU, Targetmol, TMO-T17341-500mg) dissolved in PBS at a concentration of 10 mg ml−1. After cryo-sectioning, click-chemistry was performed for staining with the In Vivo EdU Click Kit 488 (Merck, BCK488-IV-IM-S) according to manufacturer’s instructions. In brief, 10–12 µm thin frozen tissue sections were fixed with 4% PFA in PBS for 15 min at room temperature. Cells were washed twice with 1× PBS before permeabilization with 0.3% Triton X-100 in PBS for 15 min at room temperature. Five hundred microlitres of EdU reaction cocktail was added to each section and incubated for 30 min at room temperature, protected from light. Afterwards, sections were washed in PBS and antibody staining was carried out as described above.
Keratinocyte culture
Newborn, primary mouse epidermal keratinocytes derived from B6-LSL-Cas9-eGFP mice were isolated as previously described66. In brief, isolated epidermal keratinocytes were cultured on 3T3-S2 feeder layer previously treated with Mitomycin-C in 0.05 mM Ca2+ E-media supplemented with 15% serum. After 3 passages on 3T3-S2 feeder layer, cells were cultured in 0.05 mM Ca2+ E-media, made in house as previously described14. The B6-LSL-Cas9-eGFP keratinocytes were transiently transfected with 1 µg of a Cre-expressing plasmid in 6 wells using Lipofectamine 2000 (Invitrogen 11668-027) to activate Cas9–GFP. Cells were then sorted for GFP-positive signal on a BD Aria III and put back in culture using 0.05 mM Ca2+ E-media. For this project, all further modified cell lines derived from the Cas9–GFP activated cell line and were grown at standard conditions, 37 °C and 5% CO2. The SCC cell line was previously isolated from DMBA/TPA chemically induced skin tumours and brought into culture. Cell lines were tested for mycoplasma using the Mycoplasma PCR detection kit (Sigma; D9307).
In vitro lentivirus infections
For lentiviral infections in culture, B6.Cas9 keratinocytes (see above) were plated in 6-well plate (Thermo Scientific Nunclon TM Delta Surface 140675) at 1.5 × 105 cells per well and infected with 100–300 µl of low-titre virus in the presence of infection mix (1 in 10 dilution of polybrene (10 mg ml−1 Sigma 107689-100MG in PBS) in FBS[-]), by centrifuging plates at 1,100g for 30 min at 37 °C in a Thermo Heraeus Megafuge 40 R centrifuge. Infected cells were sorted by FACS for mCherry on a BD FACSAria III. Viruses used for generation of knockout keratinocytes carried the following sgRNAs: Trp53_sgRNA.3, Fat1_sgRNA.3, Notch1_sgRNA.2, ctrl_44. Additionally, we used constructs carrying sgRNAs and Tnf cDNA: Trp53_sgRNA.3 + Tnf, Notch1_sgRNA.2 + Tnf, ctrl44-Tnf.
Enzyme-linked immunosorbent assay
Mouse TNF uncoated ELISA Kit from Invitrogen (88-7324-22) was used according to the manufacturer’s instructions to quantify TNF concentration in supernatants from confluent 6 wells of keratinocyte cultures after 24 h of incubation. Measurements of 96-well plates in technical and biological triplicates from separate wells were taken on a Tecan Infinite M1000Pro plate reader. Concentrations in supernatant were calculated according to included recombinant TNF standard curve.
RT–qPCR
In vitro-infected keratinocytes and SCC cells were lysed in 6 wells with 1 ml TRIzol Reagent and the RNA was purified using the Direct-zol RNA Miniprep kit (Zymo Research). The procedure was performed according to the manufacturer protocol except for an additional 1 min centrifugation after the last washing step to completely remove residual ethanol. cDNA was synthesized using the Promega GoScript Reverse Transcription Mix, Oligo Protocol. In this procedure, 500–1,000 ng of RNA were converted into oligo(dT)-primed first-strand cDNA. iTaq Universal SYBR Green Supermix was used according to the manufacturer’s protocol for RT–qPCR reaction in a Quant Studio 7 Flex (Applied Biosystem by Life Technologies). RT–qPCR primers can be found in Supplementary Table 9. The delta-Ct method in Quant Studio Real-Time PCR software (v1.3) was used for analysis and to calculate fold changes based on Ct values.
Western blot
Keratinocytes were cultured as described above, washed with PBS and lysed in protein sample buffer (100 mM Tris-HCl pH 6.8, 4% SDS, 20% glycerol, and 0.2 M DTT). The lysate was heated for 5–10 min to 95 °C and vortexed briefly to shear genomic DNA. Isolated proteins were stored at −20 °C for later use or loaded directly (10–20 µg). Proteins were separated by 4–12% Bis-Tris SDS–PAGE electrophoresis (at 180 V, 300 W, 55 min) and transferred onto a nitrocellulose membrane (GE Healthcare). Membranes were blocked in 5% BSA (Sigma A3059-100G) in TBS-Tween (0.1%) for 1 h at room temperature, incubated with primary antibodies in 3% BSA at 4 °C overnight and washed with TBS-Tween. The secondary antibodies were added at room temperature in 3% BSA for 2 h. Western blots were developed with freshly mixed ECL solutions (GE Healthcare). The following antibodies were used: anti-TNF (Cell Signaling Technology 11948), anti-vinculin (Abcam, ab129002). Secondary Goat Anti-Rabbit IgG HRP Linked Antibody (1:1,000 dilution, Cell Signaling Technology 7074S).
Scratch assay
Cells were seeded in an ibidi Culture-Insert 2-well (81176) according to the manufacturer’s recommendations. This insert creates a clean cell-free gap in a 6- or 12-well plate (Thermo Scientific, 150239). The assay was performed according to the manufacturer’s protocol (ibidi). In brief, once confluency was reached, the insert was removed, and detached cells were washed away with PBS. New medium (with or without 10 ng ml−1 recombinant TNF (R&D Systems 410-MT) was added, and the wells were imaged every hour for 22 h in an incubation chamber at 37 °C and 5% CO2. The images were taken with a Zeiss Axio Observer controlled via the ZEN software. The data analysis was performed with an Image J plugin67 (version 1.54 h) and default settings.
TNF and TNF module overexpression constructs
Full-length mouse Tnf was amplified by PCR from GFP-TNF-alpha plasmid (Addgene #28089) via Phusion PCR according to the manufacturer’s protocol with the following primer: forward: 5′-agcctgctgaagcaggccggcgacgtggaggagaaccccggccccatgagcacagaaagcatgatccgcgacgt-3′, reverse: 5′-ATATaacgcgtCTAcagagcaatgactccaaag’. Primers created a 5′-P2A-site overhang with the opposing site added to the 3′ end of mCherry. mCherry-P2A-TNF was then ligated into the CROP-seq vector via Pfl23 and MluI restriction sites. Expression was tested with quantitative PCR, western blot and immunofluorescence (see Fig. 5).
Full-length cDNA for Jun, Ccn1 and Fos was Phusion PCR amplified from Addgene plasmids #47443, #12519 and #193089 with the following primers: Fos fw (5′-ATGGATCCGGGgccaccaacttcagcctgctgaagcaggccggcgacgtggaggagaaccccggccccATGATGTTCTCGGGTT-3′), Fos rv (5′-taacgcgtTCACAGGGCCAGCAGCGTG-3′); Jun fw (GGATCCGGGgccaccaacttcagcctgctgaagcaggccggcgacgtggaggagaaccccggccccatgactgcaaagatggaaac-3′), Jun rv (5′-ataacgcgttcaaaacgtttgcaact-3′); Ccn1 fw (5′-ataGGATCCGGGgccaccaacttcagcctgctgaagcaggccggcgacgtggaggagaaccccggccccATGAGCTCCAGCACCTTC -3′), Ccn1 rv (5′-taaacgcgtTTAGTCCCTGAACTTGTGGATGTC-3′).
PCR amplicons were ligated via BamHI and MluI restriction site overhangs into the same EF1-mCherry-2A expression vector as TNF.
MTT assay
Proliferation rates for different treatments and overexpression constructs were tested using the MTT Cell Growth Assay Kit (Merck CT01). In short, 10,000 keratinocytes were seeded per 96-well of a flat-bottom tissue culture plate and incubated in 100 μl regular E-low medium for 24 h at 37 °C, 5% CO2 before adding 10 μl of freshly prepared MTT solution for four hours. The resulting formazan crystals were dissolved by adding 100 μl isopropanol with 0.04 N HCl to each well with vigorous pipetting. Absorbance at 570 nm with a reference wavelength of 630 nm was measured directly on a Tecan Infinite M1000Pro plate reader.
CRISPR sgRNA efficiency
Eleven sgRNAs representing top, middle and bottom enrichment cohorts were selected and sgRNA efficiency was assessed by infecting B6.Cas9 keratinocytes. Lentivirus was derived previously with low-titre lentiviral preparations of the CROP-Seq-Puro vector and single guides. Cells were harvested after 10 days of puromycin selection (1 µg ml−1, Gibco A11138-03). Genomic DNA was extracted using the DNeasy Blood and Tissue Kit (QIAGEN 69504), followed by Taq (NEB M0273) PCR following the manufacturer’s instructions. Amplicons were 800–1,000-bp-long and amplified by primers, designed with primer 3 listed in Supplementary Table 1. T7 endonuclease assay was performed as described in the ALT-R genomic editing detection kit (IDT). Samples were quantified on a Bioanalyzer 2100 with a DNA high-sensitivity chip. Primers for amplification of genomic target regions of guides can be found in Supplementary Table 1.
Visium spatial gene expression
The workflow of Mirzazadeh et al.68 for RNA-rescue spatial transcriptomics was followed to overcome low RNA integrity number score. A microarray of fresh-frozen DMBA/TPA-induced squamous cell carcinoma samples was cryo-sectioned with 10 µm thickness (Leica, CM3050S) and placed on the capture area of a spatial gene-expression slide (10X Genomics, 1000188). Samples were stored in −80 °C before processing. The samples were fixed with 4% PFA and H&E staining was performed. The spatial libraries were then generated from the probe hybridization step (10X Genomics, 1000365) according to Visium Spatial Gene Expression Reagent Kits for FFPE (10X Genomics, User Guide CG000407 rev D, 1000361). The resulting libraries were sequenced by the Genomics Facility Basel. The sequencing was performed in a paired-end manner with dual indexing (10X Genomics, 1000251) on a NovaSeq 6000 (Illumina) with a S4 or SP (PE 51) flow cell. S4 libraries were then sequenced with the following cycle settings: Read 1 101 cycles, i7 10 cycles, i5 10 cycles, Read 2 101 cycles.
Visium spatial expression data processing
Sequenced libraries were processed using Space Ranger software (version 2.0.1, 10X Genomics). Reads were aligned with the ‘spaceranger count’ pipeline to the pre-built mouse reference genome provided by 10X Genomics version 2020-A (comprising the STAR-indexed mm10 mouse genome and the GTF gene annotation from GENCODE version M23 (Ensembl 98)), and the probe set reference CSV file (Visium Mouse Transcriptome Probe Set v1.0) was provided for filtering of valid genes. The resulting filtered count matrix of features (that is, genes) per spot barcodes were employed for the downstream analysis. Proportions of cell-type mixtures in each Visium spot were deconvoluted with the robust cell-type decomposition (RCTD)69 method from the spacexr package (version 2.2.1; https://github.com/dmcable/spacexr), leveraging as single-cell reference the union of mCherry-positive (Fig. 4c) and mCherry-negative (Extended Data Fig. 6c) tumour cells (merged dataset in Extended Data Fig. 7e). Seurat R package (version 4.9.9.9041) was used to plot these cell-type deconvolution results. Spatial co-occurrence of cell types was analysed with the ISCHIA package40 (version 1.0.0.0), and resulting probabilities for positive (Pgt) or negative (Plt) co-occurrence were plotted in −log10 scale.
Bioinformatic processing of the scRNA-seq data
The raw sequencing data comprising the BCL files were demultiplexed using the Illumina bcl2fastq program v2.20 with default options, allowing one mismatch of the sample barcode sequences. The resulting fastq files were processed by the BD Rhapsody Pipeline version 1.9.1 hosted at the Seven Bridge’s cloud platform (https://www.sevenbridges.com/). A custom mouse reference genome based on the Gencode version 25 was utilized and combined with our 500 sgRNA sequences. Since the BD pipeline uses STAR aligner70 in the backend, the custom STAR reference genome was generated by the genomeGenerate command in STAR. The Refined Putative Cell Calling option was disabled while running the pipeline. The output unique molecular identifier (UMI) counts data corrected by the BD Genomics RSEC (Recursive Substitution Error Correction) method were used for further downstream analysis.
Processing of the PCR sgRNA dial-out data
The PCR dial-out data was processed by a custom in-house Python script. In brief, the script extracted three 9-nucleotide (nt) long cell barcodes out of possible 97 barcodes at specific locations from the Read 1, allowing one mismatch per barcode (Hamming distance of 1). It also obtained 8-nt UMI sequences from the Read 1, as provided in the BD manual. For valid barcodes, it counted the presence of sgRNAs in Read 2 by an exact match of the 20-nt sgRNA sequence. Next, the UMI counts per cell were deduplicated and sgRNA UMI counts per cell were obtained. If multiple sgRNAs were detected in one cell, the cell was only assigned to a specific sgRNA if the sgRNA UMI count was higher than the quantile 0.99.
Downstream processing of the scRNA-seq data
To remove the doublets, the Scrublet71 Python package was utilized. Doublets were removed from each dataset separately with initial filtering of genes expressed in minimum 5 cells and expected doublet rate of 20% as detected by the BD Rhapsody. Next, 7 and 8 scRNA-seq datasets for P4 and P60, respectively, were merged in Seurat 472. Cells with detected sgRNAs and that were not doublets were selected for further analysis. Cells with UMI counts >500, UMI counts lower than quantile 0.99 of total UMI counts (to remove outliers) and cells <20% mitochondrial genes were filtered. Afterwards, the data were processed by the standard Seurat pipeline of normalizing the data with NormalizeData function, which scales the counts to 10,000 per cell followed by natural log transform using log1p. Two thousand variable genes were detected by the FindVariableGenes function. The normalized data was scaled by the ScaleData function followed by PCA analysis with the RunPCA function. The contribution of each principal component in explaining the variance of the data with the ElbowPlot function and selected top principal components for each dataset of P4, P60 and tumours. The selected principal components were used for finding nearest neighbours with the FindNeighbors function. The clustering of the data with FindClusters function with modularity optimization by Louvain algorithm. UMAP dimensional reduction technique was used to visualize the data in two dimensions. No batch effect in the P4 and P60 datasets was observed in the UMAP projections. Hence no batch effect correction was performed. The marker genes for each cluster were obtained by the FindAllMarkers function and the clusters were annotated with cell types based on the expression of known marker genes. The differential expression analysis of the scRNA-seq data for each perturbation was carried out using the MAST package73 and Wilcoxon rank-sum test, which yielded comparable results. The kernel density estimates for gene expression were inferred by the Nebulosa package.
Analysis of average expression and perturbation-perturbation matrix
To obtain an overview of the transcriptional phenotype of each perturbation, the average gene expression of cells for each perturbation in basal cells and all cells of P60 was calculated. Next the top 2,000 most variable genes in those perturbations were selected and heat maps of correlation matrices were computed using the pheatmap package with ward.D2 method for clustering. Moreover, we utilized the preserve_neighbors function from PyMDE package (pymde.org) for minimum-distortion embedding of the data in two dimensions.
sgRNA counts in tumour data
One-hundred and sixty eight individual tumours were sequenced to quantify sgRNA representation in each tumour. The raw sequencing data were demultiplexed with bcl2fastq function with corresponding i5 and i7 indexing barcodes. The number of sgRNAs in the sample was counted by the count command from the MAGECK74 package. The resulting sgRNA representations in tumours were used for total representation and selection scores in tumours.
Immune cell depletion
CD4+ or CD8+ T cells were depleted by intraperitoneal injection of a cocktail of monoclonal rat anti-mouse CD4 (YTS 191.1, Hoelzel Diagnostika LEIN-C3210) and rat anti-mouse CD8 (YTS169.4 Hoelzel Diagnostika LEIN-C2850) neutralizing antibodies at 10 mg per kg body weight. Antibodies were injected every two weeks from P4 to P60. Macrophages were depleted by intraperitoneally injecting rat anti-mouse CD115 (CSF-1R, Hoelzel Diagnostika LEIN-C2268) at an initial dose of 0.5 mg followed by 25 mg per kg body weight intraperitoneally twice per week from P4 to P60. Control mice received rat IgG isotype control (Hoelzel Diagnostika, LEIN-I-1177). We used 6 mice for T cell depletion, 7 mice for macrophage depletion and 7 control IgG mice. At P60, the epidermis was prepared as described above and split into two parts (anterior back skin and posterior back skin) for separate amplicon sequencing runs. Following genomic DNA isolation, PCR amplification of sgRNA cassettes was carried out with the same protocol as the sgRNA amplification protocol from tumours. Amplicons were sequenced on a Novaseq 6000 SP flow cell with 10% PhiX, dual indexing i5, 8 bp; i7, 8 bp; and R1 for 64 cycles. The number of sgRNAs in the fastq files was counted by MAGECK. DESeq2 was utilized for differential enrichment analysis of the sgRNAs at gene level65. The count data for the samples were normalized using the median of ratios method using the 50 control sgRNAs and the dispersion across samples was estimated. A negative binomial generalized linear model was fitted for each guide to detect differentially enriched or depleted guides. The P values obtained by the Wald test were corrected by the Benjamini–Hochberg multiple testing procedure.
WGCNA analysis
For the WGCNA24 analysis, the dataset was downsampled to a maximum of 500 cells per perturbation using the subset function with downsample = 500 option in Seurat. Next, the data were normalized and 2,000 most variable genes were selected with NormalizeData and FindVariableFeatures functions. The power calculation was performed by the pickSoftThreshold function from the WGCNA package for powers 1 to 30. Next, we ran WGCNA analysis with blockwiseModules function with signed network, minimum module size of 10 and mergeCutHeight value 0.15 and correlation function bicor for robust correlations using the power estimate from the previous step. The cells with control sgRNA were not included for WGCNA analysis unless specified otherwise. WGCNA analysis was carried out for each cluster separately. To calculate the perturbation score, moduleEigengenes function was run to obtain the module eigengenes (first principal component) of the given module. Next, the values of the first principal component were normalized using the scale function in R. Next, linear modelling was performed in R to obtain perturbation score for each perturbation over control with the formula: normalized gene score ~ perturbation. The P values, standard error and effect sizes were extracted from the linear model for plotting. This approach was similar to the one described in Jin et al.25. The P values of differentially expressed module genes over control by the Fisher’s method were combined using the metap R package with the average of log fold change to compute the average perturbation effect.
Cell–cell signalling analysis
In order to decipher cell-to-cell signalling between each cluster, the Cellchat33 package with truncated mean approach with minimum ten cells for inferring signalling was used.
Integration of scRNA-seq data
The integration of P60 cluster 0 and tumour mCherry-positive clusters 0 and 7 was performed using Seurat IntegrateData function. First, each dataset was normalized and 2,000 variable genes were obtained. Next SelectIntegrationFeatures was used to select features that are repeatedly variable across datasets for integration. Next, the datasets were scaled and PCA was performed using the common features. After that, Seurat FindIntegrationAnchors function was run with the RPCA algorithm utilizing 30 principal components followed by IntegrateData to integrate the datasets. Results were compared with the Harmony75 algorithm, which yielded similar integration.
SEACells analysis
To identify meta cell states in the scRNA-seq data, the SEACells approach was employed, which identifies a group of cells defined as metacells to uncover cell states. First, the P60 EpSC dataset was downsampled to a maximum of 200 cells per sgRNA to obtain a similar number of metacells for each perturbation. Around 50 cells were chosen to define one meta cell (SEACell) and used the PCA to build the kernel. The SEACells function was run with the option n_waypoint_eigs= no of SeaCell +1 and convergence_epsilon = 1e-5 for each perturbation separately resulting in 1–4 metacells per perturbation. The SEACells model was fitted with a maximum of 100 iterations. Furthermore, the average normalized expression was calculated for each meta cell by the AverageExpression function in Seurat for plotting the gene-expression heat maps.
General data analysis
The analysis was carried out using in-house Python 3.9 and R 4.1 scripts. Data wrangling was done with the pandas library in Python and tidyverse library in R. Heat maps were created using the ComplexHeatmap and pheatmap packages. Other plots were made using the ggplot2 library in R and seaborn library in Python.
Pathway enrichment
The overrepresentation analysis for the pathway enrichment was performed by EnrichR31 (https://maayanlab.cloud/Enrichr), which is based on Fisher’s exact test, and with custom R scripts utilizing the enrichR library on MSigDB Hallmark 2020 and other gene sets. A cut-off of FDR < 0.05 was set to select differentially expressed genes. Pathway enrichment analysis was carried out using the GSEA31 pre-ranked method using the GSEApy76 (https://github.com/zqfang/GSEApy), which enables the analysis of up- and downregulated genes simultaneously. 10,000 permutations were performed at a FDR cut-off value of 0.05 to determine enriched gene sets.
Survival analysis
TCGA data was obtained from UCSC Xena Toil77 (http://xena.ucsc.edu). Survival analysis was conducted for HNSCC and other TCGA tumour types in R utilizing the survival and survminer library. The Cox Proportional Hazards regression model was used to fit the gene-expression data to survival to obtain the Hazard Ratio. Kaplan–Meier analysis was performed on samples with lower and upper tertiles (top 1/3 versus bottom 1/3 mRNA) of gene expression (in transcripts per million (TPM)) and the P values were computed by log-rank test. To obtain the association of survival with the TSK signature genes, the average expression levels of MMP9, MMP10, PTHLH, FEZ1, IL24, KCNMA1, INHBA, MAGEA4, NT5E, LAMC2 and SLITRK6 were computed. Kaplan–Meier analysis was performed with lower and upper tertiles of average expression of these genes.
Statistics and reproducibility
Whole-mount stainings shown in Figs. 1c and 3e are representative images from three fields of view from two mice each. Immunofluorescence as shown in Fig. 4 are representative images from six fields of view from four different tumours. Figure 5j–n are representative images from five fields of view from two mice each.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.