Decoding Activation of ILC2 using Time-Dependent Cell-State Selection

Cell fate determination, a fundamental process of life, is controlled through dynamic intracellular molecular networks. However, the low population of cells at the switching period of fate determination has made it technically dicult to analyze the transcriptome of the stage. Here we developed the Time-Dependent Cell-State Selection (TDCSS) technique, which detects an index of the switching period by live-cell imaging of secretion activity followed by simultaneous recoveries of the indexed cells for subsequent transcriptome analysis. Specically, we used the TDCSS technique to study the switching period of group2 innate lymphoid cells (ILC2s) activation indexed by interleukin (IL)-13 secretion onset. TDCSS newly classied time-dependent genes, including transiently induced genes (TIGs). The nding of IL4 and MIR155HG as TIGs demonstrated their regulatory function of ILC2s activation.

group 2 innate lymphoid cells (mILC2s) 18,19 that select their fate through inducing a rapid type 2 immune response, which produces a tremendous amount of type 2 cytokines. We further demonstrate the power of the TDCSS technique by applying it to the rare human ILC2s (hILC2s) [20][21][22] and obtained the signature of expression in the switching period of the fate determination process.

Development of the TDCSS Technique
The key technical challenge of our method is integrating the destructive single-cell transcriptome analysis and the non-invasive LCI without compromising the bene ts of either one (Fig. 1g). To overcome this challenge, we designed and developed the TDCSS technique to create an index of a switching period using real-time image analysis of LCI and to synchronously harvest the indexed cells for subsequent transcriptome analysis. The TDCSS technique consists of four steps ( Fig. 1g): (1) measurements of an individual cell activity using LCI, (2) real-time analyses of the activity trace, (3) decisions of recovery, and (4) recoveries of the target cells. In the measurement step, an image of a cell is acquired using the LCI technique. The acquired image is then analysed in real-time in the analysis step and indexed for a targeted activity of the cell, which results in the third step that determines whether to recover the cell or to abandon that cell and revert to the step 1 for a different cell. In the recovery step, the agged cell is harvested for transcriptome analysis. Evaluating the effectiveness of the TDCSS As a proof of concept and to evaluate the effectiveness of the TDCSS technique, we selected the cytokine secretion response of mILC2 because of the following reasons. First, mILC2 rapidly produces a large amount of type 2 cytokine, such as interleukin (IL)-5 and IL-13, in the activation process 23 . Therefore, the initiation of the switching period is easily detectable by the LCI technique. Second, mILC2 is activated by a humoral factor, such as IL-33 18 . Therefore, the activation response of mILC2 is easily mimicked on microscopy. Third, mILC2s are non-adherent cells, the agged mILC2s thus can be easily harvested using a glass capillary.
We performed the TDCSS technique as follows. We observed the IL-13 secretion activities of 187 individual ILC2s from mouse fat tissues stimulated by recombinant IL-33 using the real-time single-cell imaging of protein secretion 13 (LCI of secretion activity, LCI-S) (Fig. 2a). The traces of the secretion signal showed that cells started secretion after various latencies over several hours (Fig. 2b). We found that a delay time (⊿t) over 0.5 hours from the onset of secretion (index) is needed to reliably ag the target cells (Extended data Fig. 1a).
To validate the TDCSS technique, we compared mRNA levels of Il13, which encodes IL-13, between "timedependent" recovered cells using the TDCSS technique at ⊿t = 0.5 or 1.0 hours (Extended data Fig. 1b) and "snap-shot" recovered cells only from secretion-positive mILC2s (Fig. 2c) using single-cell quantitative RT-PCR (qRT-PCR) (Fig. 2d). The variation of Il13 expression of "time-dependent" recovered mILC2s at ⊿t = 0.5 hours was signi cantly reduced compared to that of "snap-shot" recovered cells, although the variation in cells at ⊿t = 1.0 hours was larger than the variation in cells at ⊿t = 0.5 hours. These results indicate homogeneity in the Il13 mRNA level in the switching period of mILC2's activation process among cells, which diversi ed within 1.0 hours. In conclusion, we found that using the onset of IL-13 secretion as an index was effective in obtaining the intracellular molecular information of mILC2 in the switching period using the TDCSS technique. The transcriptome of mILC2 in the switching period In order to explore the activation process of the mILC2 comprehensively, mILC2s in the switching period (⊿t = 0.5, 1.0, and 1.5 hours) and in others (pre-stimulation, post-activation, and silent; see Methods) were subjected to single-cell RNA-seq (scRNA-seq). Based on the principal component analysis (PCA) or the trajectory inference (TI) results, the switching period cells was found in a cluster/branch different from that of pre-stimulation cells or post-activation cells (Extended data Fig. 2a, b). This indicates that mILC2 has a characteristic gene expression pattern in its switching period of the activation process.
To classify genes expressed in the switching period, we clustered genes with differential expression (DEG) between states (one-way ANOVA p < 0.05). In the end, we classi ed 1,511 DEGs based on the transition modes in their expression levels into six groups: Early induced, Early reduced, Late induced, Late reduced, Transiently reduced, and Transiently induced (Fig. 3a). The results revealed that a large number of genes were dynamically regulated in the switching period of the mILC2's activation process.
Most importantly, the transiently induced genes (TIGs) could only be identi ed using the TDCSS technique because their expression levels were indistinguishable between pre-stimulation and late stages (Fig. 3b, arrows, the number of DEG was nine between Pre-Stim and Post-activation and was three between Pre-Stim and Silent among 212 genes). The TIGs were composed of genes that were distinct from those induced in the late stage of activation (Fig. 3c). Remarkably, Mir155hg, which is the host gene of miR-155 known to be a critical regulator of mILC2 immune response 24,25 , was among TIGs (Supplementary table 1, 2). To study TIGs' functions, we performed gene ontology (GO) analysis and found the relationships with immune activation related terms such as response to peptide and in ammatory response (Fig. 3d). The result that TIG had a uniquely enriched ontology cluster indicates that the TDCSS technique succeeded in revealing the context of the switching period of mILC2's activation process previously unavailable.
In order to investigate whether the supervision of cellular function related to TIGs was affected by choice of the index, we reassigned the index to the onset of IL-5 secretion and re-identi ed the TIGs (Extended data Fig. 4a, b). Although some GO terms were shared between IL-5 indexed TIGs and IL-13 indexed TIGs, these two groups of TIGs were not the same (Extended data Fig. 4c, d). These results demonstrate that the TDCSS technique can be used to comprehensively characterise the gene expression in the switching period of the mILC2 activation process.

Identi cation of TIGs of human ILC2
Because the e cacy of the TDCSS technique was demonstrated by the experiments using mILC2, we next analysed the expression signature of human ILC2 (hILC2) 21,22,26 during the switching period. Besides being a rare cell type in peripheral blood (50-150/mL) 21 , hILC2 starts secretion after various latency over several days (Fig. 5a, Extended data Fig. 5a). The populational rarity and temporal heterogeneity make hILC2 in the switching period hard to be recovered by using a "snap-shot" method (Extended data Fig. 5b), which has been an obstacle for studying the activation mechanism of hILC2.
We applied the TDCSS technique to hILC2 by indexing the onset of IL-13 secretion and harvested the hILC2s in the switching period (⊿t within 3 hours). We compared the transcriptomes among hILC2s in the pre-stimulation, in the post-activation (over three days after stimulation), and in the switching period obtained from three donors. We classi ed 453 differentially expressed genes into six groups in the same manner described in the mILC2 study (Fig. 4b). We obtained 10 TIGs (CARD16, HPGDS, STK17B, C9orf135, RTKN2, FRY, CS, STAT4, MIR155HG and IL4) (Fig. 4c). Remarkably, TIGs of hILC2s also included MIR155HG, same as those of mILC2s, and IL4 which encodes IL-4. IL-4 is a key inducer of type 2 immune response 27 and has an important role in the hILC2 function, including proliferation and maintaining of the CRTh2 expression 28 . Although hILC2 has been thought as a source of IL-4 20 , it was di cult to characterize the IL-4 secretion dynamics of hILC2 due to its low production of IL-4. Therefore, we examined IL-4 secretion dynamics from individual hILC2s by LCI-S and con rmed their transient secretion for several hours (Fig. 4d).

Discussion
We succeeded in a holistic exploration of the potential regulators of the fate determination by developing and using the TDCSS technique to harvest individual cells at the switching period even though such potential regulators existed in a small population of cells and only for a short duration. Discovering a uniquely enriched gene ontology cluster of the transiently induced genes (Fig. 3d) indicates that the TDCSS technique is able to exploit the IL-13 secretion onset to identify the characteristic state during the activation process of the ILC2.
The process of how activated ILC2s determine their fate to produce tremendous amount of cytokine 18,29 was unknown. IL4 and MIR155HG, whose respective products IL-4 and miR-155 are two positive regulators of the cytokine production and proliferation of the ILC2 24,28 , are upregulated in the switching period of the activation process, suggesting that they are key regulators of the following activation process. MiR-155 is known to target and repress IL-33 signalling repressors, including INPP5D 29-31 and DUSP10 32,33 . Although these genes were not expressed immediately after the secretion onset of the IL-13, they were expressed at the late stage of the activation process under the IL-33 stimulation (Extended data Fig. 5g). MIR155HG might be prepared to silence these negative regulators from the switching period of the activation process to maintain ILC2s' proliferation and their cytokine production activity. Furthermore, miR-155 also targets SOCS1 34,35 and DUSP4 36,37 , both encoding inhibitors of IL-2 signalling that support ILC2's activation 38,39 , and IFNGR1 40 , encoding a receptor of IFN-g that suppresses ILC2's activation 41 .
This indicates that miR-155 contributes to the fate determination in several ways. IL-4 28 , the other positive regulator we focused on, was con rmed to be transiently produced at the protein level. Our data show a constant expression of IL4R, which encodes an IL-4 receptor. This implies that IL-4 produced by ILC2 early in its activation contributes to the positive regulation of activation in an autocrine manner though further studies are required to support this. Genes expressed in the switching period of the activation process but currently with unknown function in the ILC2 should be studied in the future in order to better understand the fate determination mechanism of the ILC2.
While we focused on the activation process of ILC2 in this study, many LCI techniques have been developed to trace various fate determination processes, including immune response 16 , cell differentiation 42,43 , cell death 2 , and stimulus response 10,44 . Using these LCI techniques, the TDCSS technique will help elucidating contexts of these fate determination processes as well. The gene expression analysis using the TDCSS technique has the advantage of being able to robustly investigate the fate determination mechanism at a targeted stage but it is not suitable for understanding the continuous aspect of the entire fate determination process. On the other hand, the computational inference methods allow us to obtain a global perspective of the fate determination process, though it suffers shortcomings from an algorithm-dependency 45 . A combined approach of the TDCSS technique and computational inference methods will robustly de ne the fate determination process, leading to a comprehensive understanding of the whole picture of the fate determination processes.

Isolation of human ILC2s
Peripheral blood was obtained from healthy volunteers at Keio University School of Medicine. This study was approved by the Institutional Review Board of Keio University School of Medicine. All subjects provided their written informed consent.

Scan intervals were 15 minutes for hILC2s experiments and 4 minutes for mILC2s experiments. Detection of the secretion onset
For detection of increased secretion signals from each cell, acquired images were sequentially analysed using Python 3 software developed in-house (https://github.com/TanaYumi/TDCSS). The program calculated the mean intensity of each well and compared the intensity to the rst image intensity. The resulting value was used as the secretion signal. When the last three secretion signals were higher than the value estimated using linear regression of the past signal by one standard deviation, the program noti ed the operator of potential secretion onset.
After the experiments, we reanalysed the secretion signals to determine the accurate time between secretion onset and cell recovery (⊿t). The secretion onset time was determined by tting the values to the following formula:
Pre-stimulation cell was recovered before the stimulus addition from the LCI-S platform. In mouse ILC2 experiments, post-activation and silent cell was randomly recovered from cells with positive or negative secretion signal 8 hours after starting the stimulation. In human ILC2 experiments, post-activation and silent cell was randomly recovered 60 hours (2 specimens) or 120 hours (1 specimen) after starting the stimulation.

qRT-PCR
To increase the accuracy of single-cell mRNA quanti cation, we added the same amount of ERCC Spike-In RNA (622,819 copies per cell) to each sample and used it as a reference. cDNA was synthesised using CellsDirect One-Step qRT-PCR Kits (Thermo Fisher Scienti c) using TaqMan probes targeting four genes purchased from Thermo Fisher Scienti c. The PCR product was diluted ten-fold in RNase-free water, and 1 μL of the PCR product was analysed using qPCR using TaqMan Fast Universal PCR Master Mix (Thermo Fisher Scienti c). Relative gene expression was calculated using the ΔΔCt method using ID136 as a reference. Samples, where Gapdh or Rplp0 could not be detected, were omitted from the analysis.

RNA sequencing
We synthesised cDNA libraries using SMART-Seq v4 3'-DE Kits (Takara Bio Inc., Shiga, Japan) according to the manufacturer's instructions, with some modi cations, adding the synthetic oligo RNA/DNA to suppress undesirable concatemers. A total of 3,892 copies of ERCC spike-in RNA were added to each sample. The cDNA was puri ed using Agencourt AMPure XP magnetic beads (Beckman Coulter). Library quality check was performed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and Agilent High Sensitivity DNA Kits (5067-4626). Degraded or low-yield samples were removed. Qubit High Sensitivity assays (Thermo Fisher Scienti c) were performed to quantify the cDNA in each library, and 5 or 6 libraries with different indexed primers were evenly pooled. Pooled cDNA (400 pg) was tagmented using Nextera XT DNA Library Prep Kits (Illumina, San Diego, CA, USA), as described in the protocol. Library size and cDNA amount were quanti ed using an Agilent High Sensitivity DNA Kit and Qubit High Sensitivity assays, respectively. Pooled libraries were sequenced using 91-bp paired-end sequencing on a MiSeq instrument (Illumina). After library demultiplexing and adaptor trimming, we aligned read 1 to reference sequences (human: GRCh38.90, mouse: GRCm38) using TopHat 2 , and the read counts were calculated using HTSeq 3 .
The normalisation of count data Acquired read count data were converted to counts per million (CPM) and log2-transformed. Then, log2(CPM+1) values were normalised using the trimmed mean of M values (TMM) method. Genes whose mean log2(CPM+1) value was over ten were used for normalisation.
Principal component analysis PCA was performed using the Python 3 module "Scikit-learn 4 ". We used all genes for the PCA, except for genes with log2(CPM+1) > 1 under 9 cells.

Hierarchical clustering analysis
The hierarchical clustering was performed using the Python 3 module "Seaborn".

Pseudo-time analyses
The pseudo-time analysis was performed using Monocle 6 with default parameters. K-Means clustering K-Means clustering was performed with the Python 3 module "Scikit-learn 4 ". In the human ILC2s analysis, mean log2-transformed expression values were calculated among cells in the same state (Fig. 4b).

Enrichment analysis
We used Metascape34 for enrichment analysis using default settings.

Statistical analysis
Variety of the expression levels acquired by qRT-PCR experiment was tested by F-test.
To identify the differentially expressed genes among groups, one-way ANOVA was performed. Genes with P < 0.01 (human ILC2) or < 0.05 (mouse ILC2) were classi ed as highly variable. We used these highly variable genes for subsequent analyses.   Whiskers represent min/max. ND: not detected. **P < 0.01, ***P < 0.001, F-test.