A disease-associated gene desert directs macrophage inflammation through ETS2

Analysis of existing data relating to chr21q22

IBD GWAS summary statistics³ were used to perform multiple causal variant fine-mapping using susieR⁵¹, with reference minor allele and LD information calculated from 503 European samples from 1000 Genomes phase 3 (ref. ⁵²). All R analyses used v.4.2.1. Palindromic SNPs (A/T or C/G) and any SNPs that did not match by position or alleles were pruned before imputation using the ssimp equations reimplemented in R. This did not affect any candidate SNP at chr21q22. SuSiE fine-mapping results were obtained for ETS2 (identifier ENSG00000157557 or ILMN_1720158) in monocyte/macrophage datasets from the eQTL Catalogue⁵³. Co-localization analyses were performed comparing the chr21q22 IBD association with summary statistics from other chr21q22-associated diseases^3,4,5,6 and monocyte/macrophage eQTLs^{54,55,56,57,58} to determine whether there was a shared genetic basis for these different associations. This was performed using coloc (v.5.2.0)⁵⁹ using a posterior probability of H4 (PP.H4.abf) > 0.5 to call co-localization.

Raw H3K27ac ChIP–seq data from primary human immune cells were downloaded from Gene Expression Omnibus (GEO series GSE18927 and GSE96014) and processed as described previously⁶⁰ (code provided in the ‘Code availability’ section).

Processed promoter-capture Hi-C data⁶¹ from 17 primary immune cell types were downloaded from OSF (https://osf.io/u8tzp) and cell type CHiCAGO scores for chr21q22-interacting regions were extracted.

Monocyte-derived macrophage differentiation

Leukocyte cones from healthy donors were obtained from NHS Blood and Transplant (Cambridge Blood Donor Centre, Colindale Blood Centre or Tooting Blood Donor Centre). Peripheral blood mononuclear cells (PBMCs) were isolated by density centrifugation (Histopaque 1077, Sigma-Aldrich) and monocytes were positively selected using CD14 Microbeads (Miltenyi Biotec). Macrophage differentiation was performed either using conditions that model chronic inflammation (TPP)¹⁶: 3 days GM-CSF (50 ng ml⁻¹, Peprotech) followed by 3 days GM-CSF, TNF (50 ng ml⁻¹, Peprotech), PGE₂ (1 μg ml⁻¹, Sigma-Aldrich) and Pam₃CSK4 (1 μg ml⁻¹, Invivogen); or, to produce resting (M0) macrophages: 6 days M-CSF (50 ng ml⁻¹, Peprotech). All cultures were performed at 37 °C under 5% CO₂ in antibiotic-free RPMI1640 medium containing 10% FBS, GlutaMax and MEM non-essential amino acids (all Thermo Fisher Scientific). Cells were detached using Accutase (BioLegend).

Identifying a model of chronic inflammatory macrophages

Human monocyte/macrophage gene expression data files (n = 314) relating to 28 different stimuli with multiple durations of exposure (collectively comprising 67 different activation conditions) were downloaded from the GEO (GSE47189) and quantile normalized. Data from biological replicates were summarized to the median value for every gene. Gene set variation analysis⁶² (using the GSVA package in R) was performed to identify the activation condition that most closely resembled CD14⁺ monocytes/macrophages from active IBD using disease-associated lists of differentially expressed genes⁶³.

CRISPR–Cas9 editing of primary human monocytes

gRNA sequences were designed using CRISPick and synthesized by IDT (Supplementary Table 3). Alt-R CRISPR–Cas9 negative control crRNA 1 (IDT) was used as a non-targeting control. Cas9–gRNA ribonucleoproteins were assembled as described previously⁶⁰ and nucleofected into 5 × 10⁶ monocytes in 100 μl nucleofection buffer (Human Monocyte Nucleofection Kit, Lonza) using a Nucleofector 2b (Lonza, program Y-001). After nucleofection, monocytes were immediately transferred into 5 ml of prewarmed culture medium in a six-well plate, and differentiated into macrophages under TPP conditions. The editing efficiency was quantified by PCR amplification of the target region in extracted DNA. All primer sequences are provided in Supplementary Table 3. The editing efficiency at the chr21q22 locus was measured by quantification of amplified fragments (2100 Bioanalyzer, Agilent) as previously described⁶⁰. The editing efficiency for individual gRNAs was assessed using the Inference of CRISPR Edits tool⁶⁴ (ICE, Synthego).

PrimeFlow RNA assay

RNA abundance was quantified by PrimeFlow (Thermo Fisher Scientific) in chr21q22-edited and unedited (NTC) cells on days 0, 3, 4, 5 and 6 of TPP differentiation. Target probes specific for ETS2 (Alexa Fluor 647), BRWD1 (Alexa Fluor 568) and PSMG1 (Alexa Fluor 568) were used according to the manufacturer’s instructions. Data were collected using FACS Diva software and analysed using FlowJo v10 (BD Biosciences).

MPRA

Overlapping oligonucleotides containing 114 nucleotides of genomic sequence were designed to tile the region containing chr21q22 candidate SNPs (99% credible set) at 50 bp intervals. Six technical replicates were designed for every genomic sequence, each tagged by a unique 11-nucleotide barcode. Additional oligonucleotides were included to test the expression-modulating effect of every candidate SNP in the 99% credible set. Allelic constructs were designed as described previously⁶⁰ and tagged by 30 unique 11-nucleotide barcodes. Positive and negative controls were included as described previously⁶⁰. 170-nucleotide oligonucleotides were synthesized as part of a larger MPRA pool (Twist Biosciences) containing the 16-nucleotide universal primer site ACTGGCCGCTTCACTG, 114-nucleotide variable genomic sequence, KpnI and XbaI restriction sites (TGGACCTCTAGA), an 11-nucleotide barcode and the 17-nucleotide universal primer site AGATCGGAAGAGCGTCG. Cloning into the MPRA vector was performed as described previously⁶⁰. A suitable promoter for the MPRA vector (RSV) was identified by testing promoter activities in TPP macrophages. The MPRA vector library was nucleofected into TPP macrophages (5 µg vector into 5 × 10⁶ cells) in 100 μl nucleofection buffer (Human Macrophage Nucleofection Kit, Lonza) using a Nucleofector 2b (program Y-011). To ensure adequate barcode representation, a minimum of 2 × 10⁷ cells was nucleofected for every donor (n = 8). After 24 h, RNA was extracted and sequencing libraries were made from mRNA or DNA input vector as described previously⁶⁰. Libraries were sequenced on the Illumina HiSeq2500 high-output flow-cell (50 bp, single-end reads). Data were demultiplexed and converted to FASTQ files using bcl2fastq and preprocessed as previously described using FastQC⁶⁰. To identify regions of enhancer activity, a paired t-test was first performed to identify genomic sequences that enhanced transcription and a sliding-window analysis (300 bp window) was then performed using the les package in R. Expression-modulating variants were identified using QuASAR-MPRA⁶⁵, as described previously⁶⁰.

BaalChIP

Publicly available PU.1 ChIP–seq datasets from human macrophages were downloaded from GEO, and BAM files were examined (IGV genome browser) to identify heterozygous samples (that is, files containing both A and G allele reads at chr21:40466570; hg19). Two suitable samples were identified (GSM1681423 and GSM1681429) and used for a Bayesian analysis of allelic imbalances in PU.1 binding (implemented in the BaalChIP package⁶⁶ in R) with correction for biases introduced by overdispersion and biases towards the reference allele.

Allele-specific PU.1 ChIP genotyping

A 100 ml blood sample was taken from five healthy rs2836882 heterozygotes (assessed by Taqman genotyping; Thermo Fisher Scientific). All of the participants provided written informed consent. Ethical approval was provided by the London–Brent Regional Ethics Committee (21/LO/0682). Monocytes were isolated from PBMCs using CD14 Microbeads (Miltenyi Biotec) and differentiated into inflammatory macrophages using TPP conditions¹⁶. After differentiation, macrophages were detached and cross-linked for 10 min in fresh medium containing 1% formaldehyde. Cross-linking was quenched with glycine (final concentration 0.125 M, 5 min). Nucleus preparation and shearing were performed as described previously⁶⁰ with 10 cycles sonication (30 s on/30 s off, Bioruptor Pico, Diagenode). PU.1 was immunoprecipitated overnight at 4 °C using a polyclonal anti-PU.1 antibody (1:25; Cell Signaling) using the SimpleChIP Plus kit (Cell Signaling). The ratio of rs2836882 alleles in the PU.1-bound DNA was quantified in duplicate by TaqMan genotyping (assay C 2601507_20). A standard curve was generated using fixed ratios of geneblocks containing either the risk or non-risk allele (200-nucleotide genomic sequence centred on rs2836882; Genewiz).

PU.1 MPRA ChIP–seq

The MPRA vector library was transfected into TPP macrophages from six healthy donors. Assessment of PU.1 binding to SNP alleles was performed as described previously⁶⁰, with minimal sonication (to remove contaminants without chromatin shearing). Immunoprecipitation was performed as described above. Sequencing libraries were prepared as for MPRA and sequenced on the MiSeq system (50 bp, single-end reads).

ATAC–seq analysis

ATAC–seq in ETS2-edited and unedited TPP macrophages was performed using the Omni-ATAC protocol⁶⁷ with the following modifications: the cell number was increased to 75,000 cells; the cell lysis time was increased to 5 min; the volume of Tn5 transposase in the transposition mixture was doubled; and the duration of the transposition step was extended to 40 min. Amplified libraries were cleaned using AMPure XP beads (Beckman Coulter) and sequenced on the NovaSeq6000 system (100 bp paired-end reads). Data were processed as described previously⁶⁸. Differential ATAC–seq analysis was performed as described previously using edgeR and TMM normalization⁶⁹. Allele-specific ATAC–seq analysis was performed in 16 heterozygous monocyte datasets from healthy controls and patients with ankylosing spondylitis⁷⁰ and in 2 deeply sequenced heterozygous TPP macrophage samples. For these analyses, sequencing reads at rs2836882 were extracted from preprocessed data using splitSNP (https://github.com/astatham/splitSNP) (see the ‘Code availability’ section).

H3K27ac ChIP–seq

H3K27ac ChIP–seq was performed as described previously⁶⁰ using an anti-H3K27ac antibody (1:250, Abcam) or an isotype control (1:500, rabbit IgG, Abcam). Sequencing libraries from TPP macrophages from major and minor allele homozygotes at rs2836882 (identified through the NIHR BioResource, n = 4) were sequenced on the HiSeq4000 system (50 bp, single-end reads). Sequencing libraries from ETS2-edited and unedited TPP macrophages (n = 3) or resting M0 macrophages overexpressing ETS2 or control mRNA (n = 3) were sequenced on the NovaSeq6000 system (100 bp, paired-end reads). Raw data were processed, quality controlled and analysed as described previously using the Burrows-Wheeler Aligner⁶⁰. Unpaired differential ChIP–seq analysis, to compare rs2836882 genotypes, was performed using MEDIPS⁷¹ by dividing the 560 kb region around rs2836882 (chr21:40150000–40710000, hg19) into 5 kb bins. Paired differential ChIP–seq analyses, to assess the effect of perturbing ETS2 expression on enhancer activity, were performed using edgeR with TMM normalization^69,72 (with donor as covariate). Genome-wide analyses used consensus MACS2 peaks. Superenhancer activity was evaluated using Rank-Ordering of Super-Enhancers (ROSE). Chr21q22-based analyses used the enhancer coordinates that exhibited allele-specific activity (chr21:40465000–40470000, hg19). Code is provided for all data analysis (see the ‘Code availability’ section).

Assays of macrophage effector functions

Flow cytometry

Expression of myeloid markers was assessed using flow cytometry (BD LSRFortessa X-20) with the following panel: CD11b PE/Dazzle 594 (BioLegend), CD14 evolve605 (Thermo Fisher Scientific), CD16 PerCP (BioLegend), CD68 FITC (BioLegend), Live/Dead Fixable Aqua Dead Cell Stain (Thermo Fisher Scientific) and Fc Receptor Blocking Reagent (Miltenyi). All antibodies were used at a dilution of 1:40; Live/Dead stained was used at 1:400 dilution. Data were collected using FACS Diva and analysed using FlowJo v.10 (BD Biosciences).

Cytokine quantification

Supernatants were collected on day 6 of TPP macrophage culture and frozen. Cytokine concentrations were quantified in duplicate by electrochemiluminescence using assays (Meso Scale Diagnostics, DISCOVERY WORKBENCH v.4.0).

Phagocytosis

Phagocytosis was assessed using fluorescently labelled Zymosan particles (Green Zymosan, Abcam) according to the manufacturer’s instructions. Cells were seeded at 10⁵ cells per well in 96-well round-bottom plates. Cytochalasin D (10 μg ml⁻¹, Thermo Fisher Scientific) was used as a negative control. Phagocytosis was quantified by flow cytometry, and a phagocytosis index was calculated (the proportion of positive cells multiplied by their mean fluorescence intensity).

Extracellular ROS production

Extracellular ROS production was quantified using the Diogenes Enhanced Superoxide Detection Kit (National Diagnostics) according to the manufacturer’s protocol. Cells were seeded at a density of 10⁵ cells per well and prestimulated with PMA (200 ng ml⁻¹, Sigma-Aldrich).

Western blotting

Western blotting was performed as described previously⁷³ using the following primary antibodies: mouse anti-gp91phox (1:2,000), mouse anti-p22phox (1:500; both Santa Cruz), rabbit anti-C17ORF62/EROS (1:1,000; Atlas), mouse anti-vinculin (Sigma-Aldrich). Loading controls were run on the same gel. Secondary antibodies were as follows: goat anti-rabbit IgG-horseradish or goat anti-mouse IgG-horseradish peroxidase (both 1:10,000; Jackson Immuno). Chemiluminescence was recorded on the ChemiDoc Touch imager (Bio-Rad) after incubation of the membrane with ECL (Thermo Fisher Scientific) or SuperSignal West Pico PLUS (Thermo Fisher Scientific) reagent. Densitometry analysis was performed using ImageJ.

RNA-seq analysis

RNA was isolated from macrophage lysates (AllPrep DNA/RNA Micro Kit, Qiagen) and sequencing libraries were prepared from 10 ng RNA using the SMARTer Stranded Total RNA-Seq Kit v2 Pico Input Mammalian (Takara) according to the manufacturer’s instructions. Libraries were sequenced on either the NextSeq 2000 (50 bp paired-end reads: CRISPR, roxadustat and PD-0325901 experiments) or NovaSeq 6000 (100 bp paired-end reads: overexpression experiments) system and preprocessed using MultiQC. Reads were trimmed using Trim Galore (Phred score 24) and filtered to remove reads <20 bp. Ribosomal reads (mapping to human ribosomal DNA complete repeating unit; GenBank: U13369 .1) were removed using BBSplit (https://sourceforge.net/projects/bbmap/). Reads were aligned to the human genome (hg38) using HISAT2 (ref. ⁷⁴) and converted to BAM files, sorted and indexed using SAMtools⁷⁵. Gene read counts were obtained using the featureCounts program⁷⁶ from Rsubread using the GTF annotation file for GRCh38 (v.102). Differential expression analysis was performed in R using limma⁷⁷ with voom transformation and including donor as a covariate. Differential expression results are shown in Supplementary Tables 1 and 2.

GSEA

GSEA was performed using fGSEA⁷⁸ in R with differentially expressed gene lists ranked by t-statistic. Gene sets were obtained from GO Biological Pathways (MSigDB), experimentally derived based on differential expression analysis or sourced from published literature^{31,42,70,79,80,81,82,83,84,85,86}. Specific details of disease macrophage signatures (Fig. 3f) are provided as source data. GO pathways shown in Figs. 2–5 are as follows: GO:0002274, GO:0042116, GO:0097529, GO:0006909, GO:0071706, GO:0032732, GO:0032755, GO:0032757, GO:2000379, GO:0009060, GO:0006119 and GO:0045649. Statistical significance was calculated using the adaptive multilevel split Monte Carlo method.

IBD BioResource recall-by-genotype study

IBD patients who were rs2836882 major or minor allele homozygotes (n = 11 of each) were identified through the NIHR IBD BioResource. Patients were matched for age, sex, treatment and disease activity, and all provided written informed consent. Ethical approval was provided by the London–Brent Regional Ethics Committee (21/LO/0682). A 50 ml blood sample was taken from all patients and M0 monocyte-derived macrophages were generated as described. After 6 days, cells were collected, lysed and RNA was extracted. Quantitative PCR analysis of a panel of ETS2-regulated genes was performed in triplicate after reverse transcription (SuperScript IV VILO, Thermo Fisher Scientific) using the Quantifast SYBR Green PCR kit (Qiagen) on the Roche LightCycler 480. Primer sequences are provided in Supplementary Table 3 and PPIA and RPLP0 were used as housekeeping genes. Expression values for each gene (\({2}^{\Delta {c}_{T}}\)) were scaled to a minimum 0 and maximum 1 to enable intergene comparison.

In vitro transcription

The cDNA sequence for ETS2 (NCBI Reference Sequence Database NM005239.5) preceded by a Kozak sequence was synthesized and cloned into a TOPO vector. This was linearized and a PCR amplicon generated, adding a T7 promoter and an AG initiation sequence (Phusion, NEB). A reverse complement (control) amplicon was also generated. These amplicons were used as templates for in vitro transcription using the HiScribe T7 mRNA Kit with CleanCap Reagent AG kit (NEB) according to the manufacturer’s instructions, but with substitution of N1-methyl-pseudouridine for uridine and methylcytidine for cytidine (both Stratech) to minimize non-specific cellular activation by the transfected mRNA. mRNA was purified using the MEGAclear Kit (Thermo Fisher Scientific) and polyadenylated using an Escherichia coli poly(A) polymerase (NEB) before further clean-up (MEGAclear), quantification and analysis of the product size (NorthernMax-Gly gel, Thermo Fisher Scientific). For optimizing overexpression conditions, GFP mRNA was produced using the same method. All primer sequences are provided in Supplementary Table 3.

mRNA overexpression

Lipofectamine MessengerMAX (Thermo Fisher Scientific) was diluted in Opti-MEM (1:75 v/v), vortexed and incubated at room temperature for 10 min. IVT mRNA was then diluted in a fixed volume of Opti-MEM (112.5 µl per transfection), mixed with an equal volume of diluted Lipofectamine MessengerMAX and incubated for a further 5 min at room temperature. The transfection mix was then added dropwise to 2.5 × 10⁶ M0 macrophages (precultured for 6 days in a six-well plate in antibiotic-free RPMI1640 macrophage medium containing M-CSF (50 ng ml⁻¹, Peprotech), with medium change on day 3). For GFP overexpression, cells were detached using Accutase 18 h after transfection and GFP expression was measured using flow cytometry. For ETS2/control overexpression, either 250 ng or 500 ng mRNA was transfected and low-dose LPS (0.5 ng ml⁻¹) was added 18 h after transfection, and cells were detached using Accutase 6 h later. Representative ETS2 expression in untransfected macrophages was obtained from previous data (GSE193336). Differential H3K27ac ChIP–seq analysis in ETS2-overexpressing macrophages was performed using 500 ng RNA transfection (see the ‘Code availability’ section).

PRS

Plink1.9 (https://www.cog-genomics.org/plink/1.9/) was used to calculate a polygenic risk score (PRS) for patients in the IBD BioResource using 22 ETS2-regulated IBD-associated SNPs (β coefficients from a previous study³). Linear regression was used to compare PRSs with age at diagnosis, and logistic regression to estimate the effect of PRSs on IBD subphenotypes, including anti-TNF primary non-response (PNR), CD behaviour (B1 versus B2/B3), perianal disease and surgery. For variables with more than two levels (for example, CD location or UC location), ANOVA was used to investigate the relationship with PRS. For analyses of age at diagnosis, anti-TNF response and surgery, IBD diagnosis was included as a covariate.

SNPsea

Pathway analysis of 241 IBD-associated GWAS hits³ was performed using SNPsea v.1.0.4 (ref. ³⁴). In brief, linkage intervals were defined for every lead SNP based on the furthest correlated SNPs (r² > 0.5 in 1000 Genomes, European population) and were extended to the nearest recombination hotspots with recombination rate > 3 cM per Mb. If no genes were present in this region, the linkage interval was extended up- and downstream by 500 kb (as long-range regulatory interactions usually occur within 1 Mb). Genes within linkage intervals were tested for enrichment within 7,660 pathways, comprising 7,658 GO Biological Pathways and two lists of ETS2-regulated genes (either those significantly downregulated after ETS2 disruption with gRNA1 or those significantly upregulated after ETS2 overexpression, based on a consensus list obtained from differential expression analysis including all samples and using donor and mRNA quantity as covariates). The analysis was performed using a single score mode: assuming that only one gene per linkage interval is associated with the pathway. A null distribution of scores for each pathway was performed by sampling identically sized random SNP sets matched on the number of linked genes (5,000,000 iterations). A permutation P value was calculated by comparing the score of the IBD-associated gene list with the null scores. An enrichment statistic was calculated using a standardized effect size for the IBD-associated score compared to the mean and s.e.m. of the null scores. Gene sets relating to the following IBD-associated pathways were extracted for comparison: NOD2 signalling (GO:0032495), integrin signalling (GO:0033627, GO:0033622), TNF signalling (GO:0033209, GO:0034612), intestinal epithelium (GO:0060729, GO:0030277), Th17 cells (GO:0072539, GO:0072538, GO:2000318), T cell activation (GO:0046631, GO:0002827), IL-10 signalling (GO:0032613, GO:0032733) and autophagy (GO:0061919, GO:0010506, GO:0010508, GO:1905037, GO:0010507). SNPs associated with PSC^5,87, ankylosing spondylitis^4,87, Takayasu arteritis^6,88,89 and schizophrenia⁹⁰ (as a negative control) were collated from the indicated studies and tested for enrichment in ETS2-regulated gene lists.

ETS2 co-expression

Genes co-expressed with ETS2 across 67 human monocyte/macrophage activation conditions (normalized data from GSE47189) were identified using the rcorr function in the Hmisc package in R.

¹³C-glucose GC–MS

ETS2-edited or unedited TPP macrophages were generated in triplicate for each donor and on day 6, the medium was removed, cells were washed with PBS, and new medium with labelled glucose was added. Labelled medium was as follows: RPMI1640 medium, no glucose (Thermo Fisher Scientific), 10% FBS (Thermo Fisher Scientific), GlutaMax (Thermo Fisher Scientific), ¹³C-labelled glucose (Cambridge Isotype Laboratories). After 24 h, a timepoint selected from a time-course to establish steady-state conditions, the supernatants were snap-frozen and macrophages were detached by scraping. Macrophages were washed three times with ice-cold PBS, counted, resuspended in 600 µl ice-cold chloroform:methanol (2:1, v/v) and sonicated in a waterbath (3 times for 8 min). All of the extraction steps were performed at 4 °C as previously described⁹¹. The samples were analysed on the Agilent 7890B-7000C GC–MS system. Spitless injection (injection temperature of 270 °C) onto a DB-5MS (Agilent) was used, using helium as the carrier gas, in electron ionization mode. The initial oven temperature was 70 °C (2 min), followed by temperature gradients to 295 °C at 12.5 °C per min and to 320 °C at 25 °C per min (held for 3 min). The scan range was m/z 50–550. Data analysis was performed using in-house software MANIC (v.3.0), based on the software package GAVIN⁹². Label incorporation was calculated by subtracting the natural abundance of stable isotopes from the observed amounts. Total metabolite abundance was normalized to the internal standard (scyllo-inositol⁹¹).

Roxadustat in TPP macrophages

ETS2-edited or unedited TPP macrophages were generated as described previously. On day 5 of culture, cells were detached (Accutase) and replated at a density of 10⁵ cells per well in 96-well round-bottom plates in TPP medium containing roxadustat (FG-4592, 30 μM). After 12 h, cells were collected for functional assays and RNA-seq as described.

CUT&RUN

Precultured TPP macrophages were collected and processed immediately using the CUT&RUN Assay kit (Cell Signaling) according to the manufacturer’s instructions but omitting the use of ConA-coated beads. In brief, 5 × 10⁵ cells per reaction were pelleted, washed and resuspended in antibody binding buffer. Cells were incubated with antibodies: anti-ETS2 (1:100, Thermo Fisher Scientific) or IgG control (1:20, Cell Signaling) for 2 h at 4 °C. After washing in digitonin buffer, cells were incubated with pA/G-MNase for 1 h at 4 °C. Cells were washed twice in digitonin buffer, resuspended in the same buffer and cooled for 5 min on ice. Calcium chloride was added to activate pA/G-MNase digestion (30 min, 4 °C) before the reaction was stopped and cells incubated at 37 °C for 10 min to release cleaved chromatin fragments. DNA was extracted from the supernatants using spin columns (Cell Signaling). Library preparation was performed using the NEBNext Ultra II DNA Library Prep Kit according to a protocol available at protocols.io (https://doi.org/10.17504/protocols.io.bagaibse). Size selection was performed using AMPure XP beads (Beckman Coulter) and the fragment size was assessed using the Agilent 2100 Bioanalyzer (High Sensitivity DNA kit). Indexed libraries were sequenced on the NovaSeq 6000 system (100 bp paired-end reads). Raw data were analysed using guidelines from the Henikoff laboratory⁹³. In brief, paired-end reads were trimmed using Trim Galore and aligned to the human genome (GRCh37/hg19) using Bowtie2. BAM files were sorted, merged (technical and, where indicated, biological replicates), resorted and indexed using SAMtools. Picard was used to mark unmapped reads and SAMtools to remove these, re-sort and re-index. Bigwig files were created using the deepTools bamCoverage function. Processed data were initially analysed using the nf-core CUT&RUN pipeline v.3.0, using CPM normalization and default MACS2 parameters for peak calling. This analysis yielded acceptable quality metrics (including an average FRiP score of 0.23) but there was a high number of peaks with low fold enrichment (<4) over the control. More stringent parameters were therefore applied for peak calling (–qvalue 0.05 -f BAMPE –keep-dup all -B –nomodel) and we applied an irreproducible discovery rate (IDR; cut-off 0.001) to identify consistent peaks between replicates, implemented in the idr package in R (see the ‘Code availability’ section). Enrichment of binding motifs for ETS2 and other transcription factors expressed in TPP macrophages (cpm > 0.5) within consensus IDR peaks was calculated using TFmotifView⁹⁴ using global genomic controls. The overlap between consensus IDR peaks and the core promoter (−250bp to +35 bp from the transcription start site) and/or putative cis-regulatory elements of ETS2-regulated genes was assessed using differentially expressed gene lists after ETS2 disruption (gRNA1) or ETS2 overexpression (based on a consensus across mRNA doses, as described earlier). Putative cis-regulatory elements were defined as shared interactions (CHiCAGO score > 5) in monocyte and M0 and M1 macrophage samples from publicly available promoter-capture Hi-C data⁶¹. Predicted ETS2- and PU.1-binding sites were identified at the rs2836882 locus (chr21:40466150–40467450) using CisBP⁹⁵ (database 2.0, PWMs log odds motif model, default settings).

Intestinal scRNA-seq

Raw count data from colonic immune cells⁴¹ (including healthy controls and Crohn’s disease) were downloaded from the Single Cell Portal (https://singlecell.broadinstitute.org/single_cell). Myeloid cell data were extracted for further analysis using the cell annotation provided. Raw data were preprocessed, normalized and variance-stabilized using Seurat (v.4)⁹⁶. PCA and UMAP clustering was performed and clusters annotated using established markers and/or previous literature. Marker genes were identified using the FindAllMarkers function. Modular expression of ETS2-regulated genes (downregulated after ETS2 editing, gRNA1) was measured using the AddModuleScore function.

Spatial transcriptomics

Formalin-fixed paraffin-embedded sections (thickness, 5 μm) were cut from two PSC liver explants and two controls (healthy liver adjacent to tumour metastases), baked overnight at 60 °C and prepared for CosMx according to manufacturer’s instructions using 15 min target retrieval and 30 min protease digestion. Tissue samples were obtained through Tissue Access for Patient Benefit (TAP-B, part of the UCL-RFH Biobank) under research ethics approval: 16/WA/0289 (Wales Research Ethics Committee 4). One case and one control were included on each slide. The Human Universal Cell Characterization core panel (960 genes) was used, supplemented with 8 additional genes to improve identification of cells of interest: CD1D, EREG, ETS2, FCN1, G0S2, LYVE1, MAP2K1, MT1G. Segmentation was performed using the CosMx Human Universal Cell Segmentation Kit (RNA), Human IO PanCK/CD45 Kit (RNA) and Human CD68 Marker, Ch5 (RNA). Fields of view (FOVs) were tiled across all available regions (221 control, 378 PSC) and cyclic fluorescence in situ hybridization was performed using the CosMx SMI (Nanostring) system. Data were preprocessed on the AtoMx Spatial Informatics Platform, with images segmented to obtain cell boundaries, transcripts assigned to single cells, and a transcript by cell count matrix was obtained⁹⁷. Expression matrices, transcript coordinates, polygon coordinates, FOV coordinates and cell metadata were exported, and quality control, normalization and cell-typing were performed using InSituType⁹⁸—an R package developed to extract all the information available in every cell’s expression profile. A semi-supervised strategy was used to phenotype cells, incorporating the Liver Human Cell Atlas reference matrix. Spatial analysis of macrophage phenotypes was performed according to proximity from cholangiocytes (anchor cell type). Radius and nearest-neighbour analyses were performed using PhenoptR (https://akoyabio.github.io/phenoptr/) with macrophage distribution from cholangiocytes binned in 100 µm increments up to 500 µm. Nearest-neighbour analysis was performed to determine the distance from cholangiocytes to the nearest inflammatory and non-inflammatory macrophage and vice versa.

To generate overlay images, raw transcript and image (morphology 2D) data were exported from AtoMx. Overlays of selected ETS2-target genes (CXCL8, S100A9, CCL2, CCL5) and fluorescent morphology markers were generated using napari (v.0.4.17, https://napari.org/stable/index.html) on representative FOVs: FOV287 (PSC with involved duct), FOV294 (PSC background liver) and FOV55 (healthy liver).

Chr21q22 disease datasets

Publicly available raw RNA-seq data from the affected tissues of chr21q22-associated diseases (and controls from the same experiment) were downloaded from the GEO: IBD macrophages (GSE123141), PSC liver (GSE159676), ankylosing spondylitis synovium (GSE41038). Reads were trimmed, filtered and aligned as described earlier. For each disease dataset, a ranked list of genes was obtained by differential expression analysis between cases and controls using limma with voom transformation. For IBD macrophages, only IBD samples with active disease were included. fGSEA using ETS2-regulated gene lists was performed as described.

LINCS signatures

A total of 31,027 lists of downregulated genes after exposure of a cell line to a small molecule was obtained from the NIH LINCS database⁷ (downloaded in January 2021). These were used as gene sets for fGSEA (as described) with a ranked list of genes obtained by differential expression analysis between ETS2-edited and unedited TPP macrophages (gRNA1) using limma with voom transformation and donor as a covariate. Drug classes for gene sets with FDR-adjusted P < 0.05 were manually assigned on the basis of known mechanisms of action.

MEK inhibition in TPP macrophages

TPP macrophages were generated as described previously. On day 4 of culture, PD-0325901 (0.5 μM, Sigma-Aldrich) or vehicle (DMSO) was added. Cells were collected on day 6 and RNA was extracted and sequenced as described.

Colonic biopsy explant culture

During colonoscopy, intestinal mucosal biopsies (6 per donor) were collected from ten patients with IBD (seven patients with ulcerative colitis, three patients with Crohn’s disease). All had endoscopically active disease and were not receiving immunosuppressive or biologic therapies. All biopsies were collected from a single inflamed site. All patients provided written informed consent. Ethical approval was provided by the London–Brent Regional Ethics Committee (21/LO/0682). Biopsies were collected into Opti-MEM and, within 1 h, were weighed and placed in pairs onto a Transwell insert (Thermo Fisher Scientific), designed to create an air–liquid interface⁹⁹, in a 24-well plate. Each well contained 1 ml medium and was supplemented with either DMSO (vehicle control), PD-0325901 (0.5 μM) or infliximab (10 μg ml⁻¹; MSD). Medium was as follows: Opti-MEM I (Gibco), GlutaMax (Thermo Fisher Scientific), 10% FBS (Thermo Fisher Scientific), MEM non-essential amino acids (Thermo Fisher Scientific), 1% sodium pyruvate (Thermo Fisher Scientific), 1% penicillin–streptomycin (Thermo Fisher Scientific) and 50 μg ml⁻¹ gentamicin (Merck). After 18 h, the supernatants and biopsies were snap-frozen. The supernatant cytokine concentrations were quantified using the LEGENDplex Human Inflammation Panel (BioLegend). RNA was extracted from biopsies and libraries were prepared as described earlier (n = 9, RNA from one donor was too degraded). Sequencing was performed on the NovaSeq 6000 system (100 bp paired-end reads). Data were processed as described earlier and GSVA was performed for ETS2-regulated genes and biopsy-derived signatures of IBD-associated inflammation⁴⁶.

Chr21q22 genotypes in archaic humans

Using publicly available genomes from seven Neanderthal individuals^{100,101,102,103}, one Denisovan individual¹⁰⁴, and one Neanderthal and Denisovan F1 individual¹⁰⁵, genotypes were called at the disease-associated chr21q22 candidate SNPs from the respective BAM files using bcftools mpileup with base and mapping quality options -q 20 -Q 20 -C 50 and using bcftools call -m -C alleles, specifying the two alleles expected at each site in a targets file (-T option). From the resulting .vcf file, the number of reads supporting the reference and alternative alleles was extracted and stored in the ‘DP4’ field.

Inference of Relate genealogy at rs2836882

Genome-wide genealogies, previously inferred for samples of the Simons Genome Diversity Project¹⁰⁶ using Relate^107,108 (https://reichdata.hms.harvard.edu/pub/datasets/sgdp/), were downloaded from https://www.dropbox.com/sh/2gjyxe3kqzh932o/AAAQcipCHnySgEB873t9EQjNa?dl=0. Using the inferred genealogies, the genealogy at rs2836882 (chr21:40466570) was plotted using the TreeView module of Relate.

Data presentation

The following R packages were used to create figures: GenomicRanges¹⁰⁹, EnhancedVolcano¹¹⁰, ggplot2 (ref. ¹¹¹), gplots¹¹², karyoploteR¹¹³.

Statistical methodology

Statistical methods used in MPRA analysis, fGSEA and SNPsea are described above. For other analyses, comparison of continuous variables between two groups was performed using Wilcoxon matched-pairs tests (paired) or Mann–Whitney U-tests (unpaired) for nonparametric data or a t-tests for parametric data. Comparison against a hypothetical value was performed using Wilcoxon signed-rank tests for nonparametric data or one-sample t-tests for parametric data. A Shapiro–Wilk test was used to confirm normality. Two-sided tests were used as standard unless a specific hypothesis was being tested. Sample sizes are provided in the main text and figure captions.

Reporting summary

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