Cytokine injection
A list of the cytokines studied and their alternative names are shown in Fig. 1d and Supplementary Table 1. We selected 86 cytokines representing most members of major cytokine families, including the following families: IL-1 (IL-1α, IL-1β, IL-1Ra, IL-18, IL-33, IL-36α and IL-36Ra); common γ-chain/IL-13/TSLP (IL-2, IL-4, IL-13, IL-7, TSLP, IL-9, IL-15 and IL-21); common β-chain (GM-CSF, IL-3 and IL-5); IL-6/IL-12 (IL-6, IL-11, IL-27, IL-30, IL-31, LIF, OSM, CT-1, NP, IL-12, IL-23 and IL-Y); IL-10 (IL-10, IL-19, IL-20, IL-22 and IL-24); IL-17 (IL-17A, IL-17B, IL-17C, IL-17D, IL-17E, IL-17F); interferon (type I: IFNα1, IFNβ, IFNε and IFNκ; type II: IFNγ; and type III: IFNλ2); TNF (LTα1/β2, LTα2/β1, TNF, OX40L, CD40L, FasL, CD27L, CD30L, 4-1BBL, TRAIL, RANKL, TWEAK, APRIL, BAFF, LIGHT, TL1A and GITRL); complement (C3a and C5a); growth factor (FLT3L, IL-34, M-CSF, G-CSF, SCF, EGF, VEGF, FGFβ, HGF and IGF-1). Representative cytokines (TGFβ1, GDNF, persephin (PSPN), prolactin (PRL), leptin, adiponectin (AdipoQ), resistin (ADSF), noggin, decorin and thrombopoietin (TPO)) from other protein families were also included.
Every recombinant mouse cytokine was obtained from at least two separate orders. The endotoxin level was <0.1 ng µg–1 of protein for every cytokine per the information from the vendors (Peprotech and R&D). To preserve cytokine activities, carrier-free cytokines were freshly reconstituted according to the manufacturer’s instructions, stored at 4 °C in sterile conditions and used within 28 h after reconstitution. For each cytokine, 5 μg in 100 μl sterile PBS was injected into each animal. Wild-type female C57BL/6 mice were purchased from the Jackson Laboratory and used in studies as 11–15-week-old young adults after resting for at least 1 week in the facility. Mice were maintained on a 12-h light–dark cycle at room temperature (21 ± 2 °C) and 40 ± 10% humidity. Cytokines were injected under the skin (50% subcutaneous, 50% intradermal) bilaterally in the abdominal flank of each mouse. Bilateral skin-draining inguinal lymph nodes were collected 4 h after injection at 6:00–8:00 and pooled for downstream processing. For each of the 86 cytokines, replicate experiments were performed in three independent C57BL/6 mice to ensure reproducibility. As a control, PBS alone was injected into mice for each experimental batch, totalling 14 PBS-injected mice. All experiments were reviewed and approved by the Broad Institute’s Institutional Animal Care and Use Committee.
Data generation quality assurance
All samples were processed using an optimized experimental pipeline to ensure quality. In particular, batch effects that arise from experiments performed on different days are known to be a major source of artefact in transcriptomic studies. Therefore, batch-to-batch consistency was strictly experimentally ensured and then computationally verified. Specifically, the mice were ordered from the same batch and housed in the same environment. Animals were randomly allocated to the experimental groups. Lymph nodes were collected at 6:00–8:00 in all experiments to exclude the impact of circadian clocks on transcriptomic profiles. Samples were processed fresh in every experiment and were kept on ice during processing whenever possible. The same researchers performed the same steps of the sample processing and sequencing pipeline following the same, highly optimized procedures. The investigators performing animal experiments and RNA sequencing were blinded from each other during data collection. The number of batches was minimized whenever possible. The three replicated mice for each cytokine were processed in different batches to ensure that batch effects, if any, would not influence biological interpretations. All samples were sequenced on two sequencing runs, with the first sequencing run containing the first set of replicates and the second containing the second and third set of replicates. PBS controls were included in every batch to ensure comparability, and transcriptomic profiles of PBS samples from different batches were computationally compared to verify batch-to-batch consistency (Extended Data Fig. 2a). In brief, Euclidean distances were calculated for each pair of PBS-treated cells of the same cell type based on the entire transcriptome to ensure that the within-batch distances and between-batch distances were comparable.
Lymph node dissociation and cell sorting for scRNA-seq experiments
An optimized pipeline for viable cell recovery and more balanced cell-type representation was used to process lymph nodes for scRNA-seq. Lymph nodes were enzymatically digested using a protocol that maximizes the recovery of myeloid and stromal cells while maintaining high viability40. In brief, lymph nodes were placed in RPMI with collagenase IV, dispase and Dnase I at 37 °C, and cells were collected once they were detached. The cells were then immediately placed on ice and washed with PBS supplemented with 2 mM EDTA and 0.5% biotin-free BSA, then filtered through a 70 µm cell strainer. Cells were incubated with Fc blocking antibodies 4 °C, then with a biotinylated anti-CD3 and anti-CD19 antibody cocktail. Antibodies were used at a dilution of 1:100. Streptavidin microbeads were then added and the cells were magnetically sorted using MACS MS columns according to the manufacturer’s protocol (Miltenyi Biotec). After cell sorting, a small fraction of the CD3+ or CD19+ cells was pooled with CD3–CD19– cells for more balanced representation of all cell types and proceeded immediately to scRNA-seq.
scRNA-seq
Cell hashing was used to combine multiple samples into the same single-cell emulsion channel41. The mouse cells obtained from different stimulation conditions were stained with TotalSeq antibodies (BioLegend anti-mouse hashtags 1–8; used at 1:100 dilution), washed 5 times at 4 °C and pooled in PBS with 0.04% BSA according to the manufacturer’s protocol. Next, 55,000 cells were loaded onto a 10x Genomics Chromium instrument (10x Genomics) according to the manufacturer’s instructions. The scRNA-seq libraries were processed using a Chromium Single Cell 3′ Library & Gel Bead v3 kit (10x Genomics) with modifications for generating hashtag libraries41. Quality control for amplified cDNA libraries and final sequencing libraries was performed using a Bioanalyzer High Sensitivity DNA kit (Agilent). scRNA-seq and hashing libraries were normalized to 4 nM concentration and pooled. The pooled libraries were paired-end sequenced on a NovaSeq S4 platform targeting an average sequencing depth of 20,000 reads per cell for gene expression libraries, and on a NovaSeq S4 or SP platform targeting 5,000 reads per cell for hashtag libraries.
scRNA-seq data pre-processing
The raw bcl sequencing data were processed using the CellRanger (v.3.0) Gene Expression pipeline (10x Genomics), including demultiplexing and alignment. Sequencing reads were aligned to the mm10 mouse reference genome, and transcriptomic count matrices were assembled. Hashtag library FASTQ files were processed using the CITE-seq-Count tool (v.1.4.3; github.com/Hoohm/CITE-seq-Count). Gene expression and hashtag were matched using the MULTIseqDemux function of the Seurat R package (v.4.1)42. Cells with multiple hashtags were considered multiplets (for example, doublets or triplets) and were excluded from further analysis. The Seurat R pipeline was used to perform quality control to include only cells with >500 genes, >1,000 unique molecular identifiers and <10% mitochondrial gene content. The expression matrix was globally scaled by normalizing the gene expression measurements by the total expression per cell. The resulting values were multiplied by a scale factor of 10,000 and natural log-transformed.
For the initial global analysis of all cells, the top 3,000 variable genes were selected for dimensionality reduction analysis. Principal component analysis (PCA) was then used to denoise and to find a lower-dimensional representation of the data. The top 75 principal components (PCs) were used for global clustering and for visualization using a t-SNE map43. Clusters were identified using the Louvain clustering algorithm. This step resulted in a total of 61 non-singleton global (level 1) clusters (Extended Data Fig. 1d). We removed potential multiplets by removing the cells with the top 2% gene counts in each cluster. As different cell types have variations in the numbers of genes detected on average, this step was done at the cluster level rather than for all data. For each level 1 cluster of cells, we then performed another round of clustering (level 2) to further verify the identity of each cluster and to remove potential doublets. This step resulted in a total of 183 global level 2 clusters. The cell-type identity of each level 2 cluster was assigned on the basis of the expression of 115 known marker genes (Supplementary Table 2). Clusters enriched for marker genes of multiple cell types were considered multiplets and removed. The top DEGs between each cell type and others are listed in Supplementary Table 2.
Quantitative measures of reproducibility across animal replicates
A gene expression vector for each biological replicate was created for each cytokine stimulation condition in a given cell type by taking the difference between the average expression vector of cytokine-treated cells and the average expression vector of PBS-treated cells. Genes were included if they were significantly differentially expressed (FDR < 0.05 and absolute log2(FC) > 0.25) compared with PBS controls and were expressed in >10% of cytokine-treated cells for upregulated genes or >10% of PBS-treated cells for downregulated genes. Rps, Rpl, mitochondrial genes and unlabelled genes were excluded. Pairwise Pearson correlation coefficients were then calculated for these vectors (Extended Data Fig. 2b).
Assessment of access to injected cytokines by each lymph node cell type
To determine whether the injected cytokines can be accessed by each cell type in the draining lymph nodes, we examined ISG expression levels in cells treated with IFNα1 or IFNβ (IFNα1/IFNβ) or PBS for each cell type. Type I interferon was chosen for this analysis given its strong induction of antiviral programmes in a wide range of cell types. A maximum of 100 cells were sampled from each condition. ISG scores were obtained by summing the normalized expressions of ISGs in each cell. These scores were then used to predict whether each cell was treated with IFNα1/IFNβ or PBS, and the accuracy of the prediction was represented as receiver operating characteristics curves (Extended Data Fig. 2c). The ISGs were obtained from the MSigDB hallmark gene set.
Differential gene expression analysis
Analyses of DEGs were performed to identify marker genes for cell clusters or cytokine-responsive genes. Analyses of DEGs were performed between two groups of cells using the two-sided Wilcoxon rank-sum test on normalized gene expression values. The P values obtained from the tests were then adjusted (Bonferroni or FDR) to address multiple testing. Genes were considered DEGs in the cytokine signature if they had FDR < 0.05 and absolute log2(FC) > 0.25 between cytokine-treated and PBS-treated cells of the same cell type, were expressed in >10% of cytokine-treated cells for upregulated genes or >10% of PBS-treated cells for downregulated genes, and satisfied the FC threshold for at least two out of the three mice (to mitigate the influence of potential single-mouse outliers). Rps, Rpl, mitochondrial genes and unlabelled genes were excluded from the signatures. Cell-type-specific cytokine signatures are listed for all major cell types in Supplementary Table 3.
Quantification of overall magnitude of transcriptomic responses to each cytokine in each cell type
We constructed a global reference map that quantified the overall gene expression changes induced by each cytokine in each cell type (Fig. 1d). This maps takes into account two metrics: the number of DEGs and the magnitude of change across the entire transcriptome. The number of DEGs is the total number of genes in each cytokine signature. The overall magnitude of cytokine-induced differential expression was computed as the Euclidean distance between the centroid vectors of cytokine-treated cells and PBS-treated cells. This value was normalized to a scale ranging from 0 (low) to 100 (high). To reduce the impact of outliers in the normalization process, winsorization was applied such that values above the 95th percentile were replaced with values at the 95th percentile before normalization. A maximum of 100 cells from each cytokine treatment condition were sampled for each cell type for the magnitude calculation. A distinct colour ramp was used for each cell type to emphasize that cell types have different properties (for example, different number of genes expressed on average) and were independently analysed. Cytokine–cell type combinations with five or more cells sampled were included in this analysis.
Identification of GPs per cytokine treatment and per cell type
GP analysis was used to identify co-regulated genes for each cytokine treatment across cell types (Fig. 2c and Supplementary Fig. 1) or for each cell type across all treatment conditions (Extended Data Figs. 5i,j, 6i,j, 7i,j and 8i,j and Supplementary Figs. 2i,j, 3i,j, 4i,j, 5i,j, 6i,j, 7i,j, 8i,j, 9i,j, 10i,j and 11i,j). GPs were constructed using the non-negative matrix factorization (NMF) algorithm44 using the R package NMFN (v.2.0). We removed genes associated with tissue dissociation45 and the cell cycle, as well as mitochondrial genes, Rps and Rpl, unlabelled genes, globally overabundant genes and those expressed in fewer than ten cells.
In the cytokine-centric analysis, NMF was used to identify cell-type-specific GPs in response to each cytokine. Cells treated with the specified cytokines or PBS were used for NMF, with a maximum of 100 cells per cell type per condition. NMF was run separately for each cytokine in this analysis, except that IFNα1 and IFNβ were processed together and IL-1α and IL-1β were processed together in Fig. 2c and Supplementary Table 5. We identified 40 GPs per cytokine treatment, with some predominantly corresponding to cell-type identity and others predominantly to cellular responses to cytokine stimulation. To quantitatively identify GPs predominantly corresponding to cellular responses to cytokine stimulation, a two-sided Wilcoxon rank-sum test was used to identify GPs with weights that were significantly different between the cytokine-treated cells and PBS-treated cells. GPs showing significant upregulation in any cell types were displayed. The top 30 genes with the highest weights for each GP were used to identify enriched biological processes using clusterProfiler (v.4.2.1)46 on the Hallmark gene sets from the MSigDB database47. Genes highlighted in the text for biological significance satisfied two criteria: (1) they were among the top 30 highest-weighted genes for a significantly upregulated GP in response to cytokines; and (2) the genes individually also showed significant upregulation in response to cytokines in the DEG analysis of cytokine signatures.
In the cell-type-centric analysis, NMF was used to analyse GPs of each cell type across all cytokine treatment conditions to identify cytokines that induce similar cellular processes. We identified ten GPs for each cell type and visualized the relationship between GPs and subclusters for each cell type using heatmaps in Extended Data Fig. 5–8 and Supplementary Figs. 2–11. The top genes with the highest weights for each GP are shown in Extended Data Figs. 5j, 6j, 7j and 8j, Supplementary Figs. 2j, 3j, 4j, 5j, 6j, 7j, 8j, 9j, 10j and 11j and Supplementary Table 8.
Analysis of secondary responses to induced IFNγ
The IFNγ-induced gene expression signature was used to infer the level of cellular responses to cytokine-induced IFNγ. The IFNγ signature score for each cell type was constructed by summing the expressions of significantly overexpressed genes in IFNγ-treated cells relative to PBS-treated cells (FDR-adjusted P < 0.001 and log2(FC) > 1) in the corresponding cell type. The log2(FC) of the signature scores and FDR-adjusted P value relative to PBS-treated cells were calculated for every cytokine treatment and shown in Extended Data Fig. 4e.
Identification of cellular polarization states
To identify cellular polarization states induced by cytokines, subclustering was performed for cells of each cell type. For heterogeneous cell types (for example, macrophage, MigDC and γδ T cell), the most abundant homogeneous subset was analysed to identify cytokine-induced states instead of re-deriving cell subsets in the polarization state analysis. We used PCs to subcluster on the basis of discriminating genes, defined as genes with a large absolute log2(FC) (between 0.75 and 1.5 depending on cell type) in any cytokine-treated cells compared with PBS-treated cells. We removed genes associated with tissue dissociation45 and the cell cycle, as well as mitochondrial genes and Rps and Rpl. We then performed PCA and visualized the cells using UMAP48. The proportion of cells falling into each cluster was calculated for each cytokine or PBS control. Major polarization states were identified on basis of two criteria: (1) cell clusters with significantly (FDR-adjusted P < 0.01) more than the expected number of cytokine-stimulated cells using a hypergeometric test; and (2) manual verification of biological relevance of the highly expressed genes or GPs in the subcluster and cytokines inducing the changes. To find discriminating markers and biological functions of each state, we analysed DEGs and co-regulated GPs per state relative to all other cells for the cell type. DEGs were identified using the two-sided Wilcoxon rank-sum test between each polarization state and other cells of the same cell type. The significantly overexpressed genes with the largest log2(FC) are shown. The most strongly polarized states are summarized in Fig. 3. The complete landscape, including less-strongly polarized states, in each cell type can be found in Extended Data Figs. 5–8 and Supplementary Figs. 2–11. We compared the polarization states by calculating the pairwise Pearson correlation coefficients between the gene expression profiles of each polarization state after subtracting the profiles of PBS-treated cells of the same cell type to remove cell-type-specific gene expression. These results are displayed in Extended Data Figs. 5b, 6b, 7b and 8b and Supplementary Figs. 2b, 3b, 4b, 5b, 6b, 7b, 8b, 9b, 10b and 11b.
We assigned a unique identifier to each polarization state using the following convention: ‘<cell type abbreviation>-<lower case letters>’. When applicable, the letters a–d were reserved for type I interferon, type II interferon, IL-1α and IL-1β, and TNF, respectively, which are cytokines that induce polarization states across a large number of cell types.
Comparative global view of polarization states across immune cell types
To gain a global view of the 66 polarization states across immune cell types defined in Fig. 3, we used Jaccard similarity index to evaluate similarity between each pair of cell states (Extended Data Fig. 9). The gene expression profile of each polarization state was compared with PBS-treated cells of the same cell type to remove cell-type-specific gene expression. The genes with an absolute log2(FC) > 0.5 compared with PBS-treated cells were used to compute the Jaccard similarity score. Upregulated and downregulated genes were separately calculated. The rows and columns were hierarchically clustered using the average-linkage method on the Euclidean distances to identify groups of similar polarization states. To visualize unique polarization states with low similarity to other states, the same results were illustrated using a force-directed network, with a higher circle size indicating a more unique state, which was calculated on the basis of the inverse of mean Jaccard similarity value with other states.
Pathway enrichment analysis for NK cells
To identify biological processes enriched for the IL-18-treated NK cells, a pre-ranked gene list was computed by subtracting average gene expression values of PBS-treated NK cells from those in IL-18-treated NK cells. Gene set enrichment analysis was performed using clusterProfiler (v.4.2.1) on the gene ontology biological processes gene sets. Gene sets with a FDR-adjusted P < 0.1 are shown. As a comparison, representative cytokines from other NK cell polarization states were analysed using the same method.
Cytokine and cytokine–receptor gene expression maps
A map of cell-type-specific production of cytokines was derived from our dataset. Cytokine genes expressed in at least 50 cells were included in the cytokine expression heatmap. The cells were obtained from all conditions (PBS or cytokine treated) to provide a map of cytokine expression under all unstimulated or cytokine stimulation conditions (that is, to account for induced expression). The gene expression level was then normalized relative to the cell type with the maximum expression level (whereby the maximum level is capped at 1 expression unit before normalization) to account for the variation in the number of transcripts produced or detected for each cytokine. A cytokine was considered expressed in a cell type if more than 0.1 normalized expression units were detected.
The cytokine–receptor expression map was constructed using the same approach. This included signalling receptors, decoy receptors and receptors that form complexes with cytokines. A list of genes encoding known functional receptors for the 86 studied cytokines are listed in Supplementary Table 1. The cytokine expression map and the cytokine–receptor expression map are shown in Fig. 4a, Extended Data Fig. 10c and Supplementary Table 9.
Cell–cell interactome network construction
A cell–cell interactome network was constructed to chart available cytokine-mediated cell–cell communication channels. The network was constructed such that the source and sink nodes are cell types and intermediate nodes are cytokines. The paths between source cell-type nodes and sink cell-type nodes through cytokine nodes were established on the basis of the detectability of the cytokine mRNA in the cell population (normalized expression > 0.1) and the responsiveness of the cell type to the cytokine (more than ten DEGs in the corresponding cytokine signature). For heteromeric cytokines or cytokine complexes composed of two subunits (IL-12, IL-23, IL-27, LTα1/β2 and LTα2/β1), the cytokine is shown as expressed and is annotated with an asterisk if the genes encoding at least one subunit are expressed as there is evidence of extracellular assembly of some components into functional cytokines under healthy or pathological conditions49. The network was plotted separately for each source node for ease of interpretability.
To construct the ligand–receptor interactome, we identified functional cognate receptors for each cytokine from the literature, which is listed in Supplementary Table 1. For a receptor to be considered expressed in a cell type, the normalized expression value of the receptor gene needed to be greater than a cut-off threshold (default of 0.1 expression unit). For heteromeric receptors, all components needed to be expressed for the receptor to be considered expressed. For cytokines with more than one functional receptor, the receptor was considered expressed if any functional receptors are expressed. We then connected the cytokines with the cell types expressing the cognate receptors. The cytokine production portion of the interactome is the same as the one in the ligand–response interactome. The ligand–response and ligand–receptor networks were then compared to generate the cell–cell communication paths that are common or different between these two approaches.
IREA for cytokine response analysis
We offer two types of IREA analysis options to assess cytokine responses in a user’s data depending on the input, which can be a gene set or a gene expression matrix. The cell type in the user data is specified by the user. User data are then compared with the transcriptional cytokine responses of the same cell type from the Immune Dictionary using the following methods:
-
1.
For the gene set input, we first find gene set scores by summing the normalized expression value of all genes in the gene set in each of the cytokine-treated cells or PBS-treated cells. Statistical significance is assessed using a two-sided Wilcoxon rank-sum test between gene set scores on cytokine-treated cells and gene set scores on PBS-treated cells, and an FDR correction is applied to all cytokine calculations. Enrichment can also be calculated using the hypergeometric test on significant DEGs (FDR < 0.01 between cytokine-treated cells and PBS-treated cells), a method commonly used in pathway analyses.
-
2.
For the gene expression matrix input, the expression matrices are first normalized such that the total expression per cell sums to 10,000 units; the expression is then log-transformed. Genes giving significant contribution to the enrichment score, with the default being those having an average of more than 0.25 expression values, were included. Next, the projection score is calculated by finding the cosine similarity score between user input and cytokine-treated or PBS-treated cells. Statistical significance is assessed using a two-sided Wilcoxon rank-sum test between projection scores on cytokine-treated cells and projection scores on PBS-treated cells, and an FDR correction is applied to all cytokine calculations. The effect size is the mean difference between projection scores on cytokine-treated cells and on PBS-treated cells. The effect size and FDR-adjusted P value for each of the 86 cytokines can then be visualized using a compass plot shown in Fig. 5d. Conceptually, this method takes into consideration the direction and magnitude of expression of each gene. That is, a strongly upregulated gene in both the user dataset and Immune Dictionary reference dataset is given a high weight that increases the overall likelihood of enrichment; a strongly upregulated gene in one dataset but not the other is given lower weight; and a gene that is upregulated in one dataset but downregulated in the other is given negative weight that decreases the overall likelihood of enrichment. The genes contributing the highest weights to the enrichment can be visualized using a diverging bar plot shown in Extended Data Fig. 12b.
IREA for cellular polarization analysis
IREA polarization analysis implements the same statistical test as the IREA cytokine response analysis. In IREA polarization, user data are compared with the polarization state gene expression profiles. A polarized radar plot is shown if at least one cellular polarization state is significantly enriched (FDR-adjusted P < 0.05). If no state is significantly enriched, the radar plot shows an enrichment score of 0 for every state, which signifies that the input cells are unpolarized. The enrichment score is normalized to be between 0 and 1 on the radar plot.
IREA for cell–cell communication network construction
We constructed models of cell–cell communication networks by taking into account cytokine production and cytokine response. Cytokine production was obtained by examining the transcripts mapped to each of the 86 cytokines. The cytokine response was assessed using IREA, and cytokines with IREA output of FDR < 0.01 were included. For heteromeric cytokines or cytokine complexes composed of two subunits, the cytokine was displayed as expressed if at least one subunit is expressed as there is evidence of extracellular assembly of some components into functional cytokines49 (same method as for the cell–cell interactome). Cytokine networks can be visualized as shown in Fig. 5e and Extended Data Fig. 12c.
IREA analysis of mouse tumour scRNA-seq data
The scRNA-seq data were downloaded as 10x Genomics data files36. Data were processed using the same approach as described above but with a minor modification, whereby 40 PCs were used for downstream analysis. Cell types were annotated as shown in the publication36. The IREA analysis was done between anti-PD-1 treatment and controls for each cell type using the transcriptome-wide approach with default parameters. A receptor expression threshold of 0.05 was applied to produce the data in the receptor ring in the cytokine enrichment plot.
IREA analysis of human COVID-19 blood scRNA-seq data
The scRNA-seq data were downloaded as a Seurat object from the human COVID-19 blood study39. Cluster annotations were used as defined in the Seurat object. IREA analysis was performed using data from ventilated patients with COVID-19 and compared with healthy individuals for each cell type using the transcriptome-wide approach with default parameters and species specified as human. IREA implements mouse and human homologue gene conversion using the most recent release of the National Center for Biotechnology Information HomoloGene database (release 68). A receptor expression threshold of 0.05 was applied to produce the data in the receptor ring in the compass plot.
Statistical analysis
The statistical tests used are described for each analysis in the corresponding text. Two-sided statistical tests were used unless otherwise specified. FDR or Bonferroni adjustments were made for the analyses for which multiple hypothesis testing applies.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Source link
