Introduction

T-cell acute lymphoblastic leukemia (T-ALL) is a malignancy characterized by the proliferation of immature T cells, known as thymocytes, and is subclassified according to the stage of thymic maturation (Terwilliger and Abdul-Hay 2017; De Keersmaecker et al. 2013). Representing approximately 10-15% of pediatric cases and 25% of adult cases of all acute lymphoblastic leukemia (Vadillo et al. 2018), T-ALL originates in the thymus. This organ is the primary site where early T cell progenitors (ETPs) originating from the hematopoietic tissue differentiate into mature, immunocompetent T cells. Thymocytes follow a well-orchestrated migratory path within distinct thymic niches, receiving essential extrinsic signals for their proliferation, differentiation, and selection (Takahama 2006). Thymic epithelial cells (TECs) play a central role in this process, acting as the main cellular component of these niches, and influencing thymocyte migration and differentiation by expressing crucial growth factors, ligands (such as DLL4), or releasing chemokines (e.g. CCL25) and cytokines such as interleukin-7 (IL-7) into the thymic niche (Yui and Rothenberg 2014; Koch et al. 2008; Zamisch et al. 2005).

The etiology of T-ALL is a topic of extensive research. Isolation and genome sequencing of malignant thymocytes from patients have uncovered mutations in over 100 putative driver genes (Liu et al. 2017; Girardi et al. 2017). Notably, over 50% of T-ALL patients harbor gain-of-function mutations in the NOTCH1 gene, which encodes the receptor for DLL4, while around 10% of patients exhibit similar mutations in the IL-7 receptor (IL7R) gene (Liu et al. 2017). The leukemogenic potential of constitutive activation of these receptors has been demonstrated across various vertebrate models (González-García et al. 2009; Silva et al. 2021; Shochat et al. 2011; Zenatti et al. 2011; Oliveira et al. 2022; Chen et al. 2007; Blackburn et al. 2012), emphasizing the evolutionary conservation of thymic developmental mechanisms that become misrouted in T-ALL (Bajoghli et al. 2019). These studies have contributed valuable insights into the impact of genetic alterations and oncogenic pathways on T-ALL development. However, the gene-centric view provided only limited understanding of the initiation and progression of the disease. In recent years, the growing evidence of the role of an abnormal tumor microenvironment in carcinogenesis strongly supports a more integrative view that considers the convergence of cellular genetics and the surrounding malignant niche as crucial elements for disease progression (Vadillo et al. 2018). In the context of T-ALL, however, the impact of the thymic niche or the role of TECs has not been studied.

Computational and mathematical models are powerful tools for integrating complex biological processes and can be used to test new hypotheses in health and disease (Ji et al. 2017; Metzcar et al. 2019; King et al. 2021). Several such models have been developed to study population dynamics during T cell development (Robert et al. 2021; Aghaallaei et al. 2021; Efroni, Harel, and Cohen 2007; Thomas-Vaslin et al. 2008; Vibert and Thomas-Vaslin 2017; Souza-e-Silva et al. 2009). Our recently developed cell-based computational model, for example, simulates individual thymocytes using agents within a spatially-resolved “virtual thymus” to explore cell-level behaviors and their effects on thymocyte population dynamics. This model was originally developed to investigate the impact of thymic niche signals and intrathymic cell localization on the αβ/γδ T cell sublineage outcomes (Aghaallaei et al. 2021).

In this study, we developed new computational tools to enhance the virtual thymus model, enabling us to identify both cell-autonomous and non-autonomous factors in thymocytes and TECs that could drive the proliferative expansion of a lineage derived from a single progenitor cell (hereafter referred to as a clone) within the thymus. We conducted an unbiased systematic analysis of various parameters, including TEC shape and architecture, IL-7 signaling, cell cycle, duration of the proliferative phase, and cell migration. This analysis facilitated a detailed comparison of proliferative expansion between wildtype (WT) and lesioned clones with alterations in the IL7R and NOTCH1 receptors. Through simulations of over 1500 scenarios, we identified factors that promote substantial expansion of malignant clones. In particular, modifications in the TEC network and niche emerged as a previously uncharacterized factor that could synergistically accelerate clonal expansion. Subsequent in vivo experiments using medaka fish (Oryzias latipes) confirmed the outcomes of our in silico analysis, solidifying the predictive accuracy of our virtual model, and revealing the impact of TECs on the initiation and progression of T-ALL.

Results and discussion

A virtual thymus to model the medaka embryonic T cell development

The details of the virtual thymus model implementation and parameter estimation have been extensively described previously (see Supplementary Material of Aghaallaei et al. 2021), and in the following we will only delineate its main features. The model considers both spatial and temporal aspects of early thymus development at multiple scales: Cells are physically represented as one or multiple spheres that mechanically interact with adhesive and repulsive forces to capture cell crowding; each cell has an independent internal state that changes dynamically based on subcellular signaling modeled with ordinary differential equations or with phenomenological rules; and at the tissue level a partial differential equation is used to solve the diffusion of extracellular molecules such as cytokines. To reduce computational cost, the computational representation focuses on a 5 µm deep slice of the lower half of the radially symmetric organ, approximately 1/10 of the total volume of a medaka embryonic thymus (Figure 1 - Supplement 1A, B). Consistent with observations in vertebrates, including in medaka (Bajoghli et al. 2009; 2015; Aghaallaei et al. 2021; 2022), ETPs within our virtual thymus mechanically and biochemically interact with TECs, proliferate and differentiate in response to signals from the thymic niche, and subsequently diverge into two distinct T cell sublineages before selection and exit from the thymus (Figure 1A). In this work, we define cells of a lineage derived from a single founder ETP as a “clone”, and we do not consider differences in immunological clonotypes that emerge after recombination of the T cell receptor.

Figure 1

In the virtual thymus model, factors governing thymocyte motility, including cell speed (Figure 1B) and directionality (Figure 1C) were fine-tuned based on quantitative noninvasive imaging of the thymus using multiple medaka transgenic reporter lines (Bajoghli et al. 2009; 2015; Aghaallaei et al. 2021). Based on these experimental observations, we modeled entry into the organ by creating new cells at the bottom of the simulated domain (ventral to the organ), these cells move into the thymus and, some time after differentiation, either leave the organ by exiting from the ventral side as a naïve T cell or undergo cell death by gradually shrinking and eventually being removed from the simulation (Figure 1 – Supplement Movie 1-2). The model accounts for the time taken for a single progenitor to enter and exit the thymus, which is estimated at approximately 3 days in zebrafish (Hess and Boehm 2012) and medaka (Bajoghli et al. 2015). For time calibration, each simulation step was set to represent 15 seconds. With a total of 48,000 steps per simulation, the entire simulation spans 720,000 seconds, equivalent to about 8.3 days (Figure 1D). This time interval was chosen to recapitulate thymic development in medaka from 2.5 to 11 days post-fertilization (dpf), because at 2.5 dpf the first ETP enters the embryonic thymus (Bajoghli et al. 2019). Correspondingly, each simulation starts with entering of the first ETP into the virtual thymus and ends with a homeostatic cell population of roughly 100 cells (mean ± standard deviation: 98±7; Figure 1 – Supplement Movie 1).

At 11 dpf, the medaka thymus contains 877±146 (N=7) thymocytes. Considering that the virtual thymus represents 1/10 of the medaka embryonic thymus, the observed cell population in the virtual thymus model reproduces the temporal dynamics in medaka within the margin of error.

Each cell’s internal state and decision to proliferate or differentiate depends on intrinsic and extrinsic factors that are integrated by signaling pathways. The virtual thymus includes Delta-like 4 (Dll4), the ligand of Notch1 receptor, and interleukin-7 (IL-7) cytokine as two cell-extrinsic factors provided by the TECs (Figure 1E). Consistent with medaka WT embryonic thymus (Aghaallaei et al. 2021; 2022), all virtual TECs in this model express the Dll4 (Figure 1E, top-right panel), with a subset spatially releasing IL-7 into the environment (Figure 1E, top-left panel), creating a short-ranged IL-7 cytokine gradient in the thymic cortex niche (Aghaallaei et al. 2021). In terms of cell-intrinsic factors, all ETPs uniformly express the Notch1 receptor (Figure 1E, middle-right panel), and each cell lineage expresses the IL-7 receptor (IL7R) at a constant level chosen randomly between 0 and 1 (Figure 1E, middle-left panel). In each simulation, new ETPs continually enter the thymic niche at regular intervals. Following engagement with TECs and receiving Notch1 and IL-7 signals (Figure 1E, bottom panels), they proliferate and differentiate (step 3 in Figure 1A). The cell cycle was modeled with an average duration of 7 hours and was subdivided into G1, S, G2, and M phases. In line with in vivo observations (Ruijtenberg and van den Heuvel 2016), virtual cells require external pro-proliferative signals during the G1 phase to commit to S through M. Once differentiated, cells can no longer enter the cell cycle. The dynamics of pro-proliferative signals and the proliferation stop after differentiation allow clones in silico to undergo up to four rounds of cell division (Aghaallaei et al. 2021).

In this work, we implemented a code that enables us to meticulously trace the unique identity of each cell upon entering the virtual thymus and subsequently track the identity of its descendants (Figure 1F, left panel; Figure 1 – Supplement 1 C). This tool identified various waves of clones because of the constant influx of new ETPs and negative selection and efflux of fully differentiated T cells (Figure 1F, arrows in the right panel, Figure 1 – Supplement Movie 1). We note a slight variability in the clone sizes, which could arise from three processes. First, random variations in the expression level of IL7R impact IL-7 signal transduction and, consequently, proliferation (Aghaallaei et al. 2021). Second, due to random cell motility and cell crowding effects, certain cells could coincidentally remain in close contact with IL-7-secreting TECs, maximizing their exposure to pro-proliferative signals such as IL-7 and Notch ligand, and thus increasing their chances of entering the cell cycle. Third, each cell determines the duration of its next cell cycle by drawing an Erlang-distributed random variable with a mean of 7 hours and a standard deviation of ≅0.99 hours.

After normalizing all clones based on their time of entry into the thymus, the pattern of developmental phases became clear across all clones (Figure 1G). In this visualization, the now near-synchronous rounds of cell division across clones can be seen as bumps in the graph indicating clonal expansion (Figure 1G, gray arrowheads). During the commitment phase, which has a minimum duration of 24h but can be extended due to insufficient Notch signaling, thymocytes were competent for proliferation, leading to an increase in the total cell number per clone. In the differentiation phase, which has a fixed duration of 24h, clones reached a maximum size. Note that differentiation prevents further entry into the cell cycle but permits cells that passed G1 and therefore committed to S through M phases to complete their division. Finally, clones underwent selection and exit from the virtual thymus, thus leaving the simulation. Together, the enhanced virtual thymus model facilitated the investigation of the heterogeneity and temporal dynamics of individual clones originating from a founder progenitor cell.

Modeling the impact of TEC architecture on thymocyte population size

In vivo, individual TECs exhibit a distinctive star-shaped morphology, and protrusions of neighboring TECs contact each other forming a three-dimensional network (Figure 2A, left panel). Using confocal imaging, we estimated the relative area occupied by thymocytes (53±3%; N=4) and TECs (47±3%; N=4). A similar morphology and tissue density was replicated in our in silico model (Figure 2A, right panel), with thymocytes taking up on average 56±2% and TECs 44±2% of the volume in homeostasis. We wondered to what extent TEC morphology and density could impact thymocyte population size. This aspect is difficult to study experimentally and is, therefore, an ideal use case for the virtual thymus model. In a series of simulations, we combinatorically varied (i) the size of TECs, (ii) the number of their protrusions, and (iii) the TEC cell density (Figure 2B). The average homeostatic thymocyte population size in each simulation was then used as a readout. The results from 26 different tested conditions predicted that a higher density of larger TECs with more protrusions led to an almost two-fold increase in the number of thymocytes (200±13, N=3; scenario 26 in Figure 2C, Figure 2 – Supplement 1A bottom) compared to the reference condition (98±7, N=19; Reference in Figure 2C). Conversely, scenarios where these parameters were reduced, particularly a reduced TEC size, also diminished the homeostatic thymocyte population size (scenarios 1-14 in Figure 2C). The condition with the least amount of cells, scenario 1, had a reduced TEC size, but an unchanged number of protrusions and density. The condition diametrically opposite to scenario 26 was scenario 6, with a reduced TEC size, reduced number of protrusions, and reduced TEC density (Figure 2 – Supplement 1A top). In both scenario 1 (45±1, N=3) and scenario 6 (51±2, N=3), the thymic population size was only about half of that in the reference condition (98∓7, N=19). Statistical analysis confirmed that all three parameters and their interactions were significant predictors of total thymocyte population size (Figure 2 – Supplement 1B-C). The statistically significant interaction can be explained by the fact that higher TEC density increases the number of TECs, and thereby amplifies the effect of an increased TEC radius or number of protrusions. Thus, these three alterations of the thymic niche morphology act synergistically to modulate thymocyte population size.

Figure 2

There are two possibilities that could explain the variations in total thymocyte numbers in simulations. Either over- or under-proliferation of a few thymocyte lineages, or the cumulative impact of subtle changes across several thymocyte clones. To distinguish between these possibilities, we used our clonal analysis tool for scenarios 6 and 26 (Figure 2D). This analysis revealed that variations in total thymocyte numbers were a consequence of changes in proliferation across all clones to a similar degree. These results predict that the TEC architecture has a direct impact on thymocyte proliferation rate, which fits well with the fact that TECs act as the main source of ligands, growth factors, cytokines, and chemokines (Gameiro, Nagib, and Verinaud 2010). Any changes in the TEC architecture could therefore influence the amount of cytokine production. In our virtual thymus model, this affects the extracellular spatial distribution of IL-7 cytokine (Figure 2E). To further explore this, we tested two additional scenarios. In scenario 27, TECs maintained their reference architecture and density but all expressed IL-7, leading to a more uniform cytokine distribution in the thymic environment (Figure 2E). The extracellular level of IL-7 (Figure 2F) and thymocyte population size (155±11, N=5; Figure 2G) in scenario 27 was similarly high as in scenario 26, where the size, protrusion number, and density of TECs were increased. Conversely, when TECs did not express IL-7 (scenario 28), the thymocyte population size (Figure 2G) was low (42±1, N=5), akin to scenario 6, where the size, protrusion, and density of TECs were reduced.

Interestingly, although the average extracellular IL-7 level is at its highest in scenario 27 with ubiquitous expression, the spatially restricted elevated IL-7 concentration in scenario 26 is more effective at driving thymic population growth (Figure 2E-G). Thymocytes first enter the organ at its periphery from below (at the ventrolateral site), where IL-7 levels are highest, and tend to migrate to the IL-7-depleted center of the thymus as they become non-proliferative (Aghaallaei et al. 2021; see also Figure 1 – Supplement Movie 1). Thus, in scenario 26, thymocytes in their proliferative phase immediately encounter elevated IL-7 at the organ periphery, stimulating cell proliferation. In contrast, in scenario 27, thymocytes that enter the organ periphery encounter comparable levels of IL-7 to the reference scenario, and only after entering deeper into the organ can the ubiquitous expression of IL-7 unfold its effects. This delay contingent to random cell migration (which is further exacerbated by cell crowding posing an obstacle), reduces the impact of elevated IL-7 in scenario 27. These results highlight how the spatial structure of the thymic niche can impact on thymocyte population dynamics.

Our simulations predict that the availability of IL-7 provided by TECs has a direct impact on thymocyte population size. To validate this prediction in vivo, we conducted two functional analyses using the medaka model organism (Figure 2H). In one experiment, employing the CRISPR-Cas9 technique, we knocked out the il7 gene (Figure 2-Supplement 2A, B) in a medaka transgenic line, where thymocytes express green fluorescent protein (GFP). Consistent with our in silico outcomes, il7 crispant embryos displayed fewer thymocytes and a smaller thymus size (Figure 2I, Figure 2 - Supplement 2C, D). To estimate the extent of cell proliferation, we counted mitotic cells throughout the entire thymus using the M-phase marker phospho-histone 3 (pH3) and normalized values to the mean of the control. This analysis further supported that a reduction in IL-7 availability in the thymus reduces thymocyte proliferation (Figure 2I, J). In a second experimental setup (Figure 2H), we artificially increased il7 levels in the thymic niche by injecting the ccl25a:il7 construct into embryos (hereafter called TEC:IL-7^HI), resulting in a strong upregulation of the levels of il7 produced by the TECs (Aghaallaei et al. 2021). Compared to WT, the thymus of TEC:IL-7^HI embryos appeared visibly larger, and the number of pH3+ cells was significantly increased (Figure 2I, J). Therefore, the in vivo results confirm the critical role of IL-7 for thymocyte proliferation in medaka, a mechanism that is evolutionarily conserved among vertebrates (Iwanami et al. 2011). Additionally, the alterations in thymic size suggest a regulatory crosstalk between thymocyte proliferation and TEC architecture.

Modeling the impact of IL-7 signaling activity on thymocyte population size

Given the important contribution of extracellular IL-7 on thymocyte population size, we next evaluated the impact of (i) extracellular IL-7 depletion (e.g., via ligand internalization), (ii) the rate of signal transduction activation upon IL7R and IL-7 binding, and (iii) the rate of decay of signaling activity (Figure 3A, B). In our virtual thymus model, the IL-7 signaling activity (σ_IL7) of a thymocyte over time (t) is modeled phenomenologically using a function involving the IL7R concentration ([IL7R]), the average extracellular IL-7 concentration (〈[IL-7_ex]〉), a parameter (a_IL-7) that scales the strength of signal transduction activation, and another parameter (d_IL-7) that scales the rate at which the signaling activity diminishes over time (Figure 3B). As before, we combinatorically perturbed these three factors, simulating 17 different scenarios, to assay their impact on thymocyte population size (Figure 3C). As expected, we observed a positive correlation between the level of signal transduction activation a_IL-7 and thymocyte numbers. Conversely, variations in signaling decay rate d_IL-7 showed the opposite effect. The rate of signal transduction deactivation models cellular short-term memory, thus a lower signaling deactivation rate indicates that cells retain their IL-7 stimulus for a longer duration. Because IL-7 signaling activity σ_IL7 promotes cell cycle entry, we expect that parameter changes that increase stimulus duration should promote cell division. Indeed, the scenario with the highest number of thymocytes occurred when signaling deactivation was low and activation was high (274±11, N=4; scenario 45 in Figure 3C). The diametrically opposite combination had the least amount of cells (42±1, N=4; scenario 29 in Figure 3C). Statistical analysis confirms that signal transduction activation a_IL-7, and the signal transduction deactivation rate d_IL-7 were both statistically significant predictors of thymocyte population size (Figure 3 - Supplement 1). There was no statistically significant effect of parameter interactions, likely because these parameters exert their effect independently. There was only a limited reduction in population size with IL-7 depletion (91±11, N=4; scenario 36 in Figure 3C) compared to no depletion in our reference setting (98±7, N=19). Overall, IL-7 depletion did not have a significant effect on thymocyte population (Figure 3 - Supplement 1). However, we noted that population size difference was more pronounced at higher cell numbers, e.g., comparing between scenario 44 (237±7, N=4) and scenario 45 (274±11, N=4). Similarly, the mean extracellular IL-7 content in the virtual thymus barely changed when comparing the reference to scenario 36 but was noticeably reduced in scenario 44 compared to scenario 45 (Figure 3D). Therefore, thymocyte proliferation could self-inhibit via depletion of extracellular IL-7, but we expect this effect to be small unless cell numbers are massively increased.

Figure 3

Together, the outcomes of our simulations reveal that molecular and cellular changes, specifically those capable of increasing the extrinsic factor IL-7 in the niche — such as a higher density of TECs or an elevated expression level of IL-7 by individual TECs — directly contribute to an increase in thymocyte population size. Tuning thymocytes’ sensitivity to IL-7 by manipulating IL7R signaling can further amplify proliferation. Finally, we expect that cytokine internalization by thymocytes has only a mild inhibitory effect on excessive population growth by limiting the availability of the pro-proliferative IL-7 in the extracellular space.

A systematic approach identified synergistic factors promoting clonal expansion of IL7R-lesioned clones

Our 45 tested scenarios thus far reached a maximum three-fold induction of thymocyte population size when parameter changes promoted increased IL-7 signaling. While the activation of IL-7 signaling is linked to T-ALL development, this pathology is associated with proliferation of clones carrying somatic mutations (Oliveira et al. 2019; 2022; Silva et al. 2021). We therefore modified the code of our virtual thymus model to enable the introduction of parameter changes exclusively in a single clone and all of its progeny (hereafter referred to as the lesioned clone). Inspired by observations that 10% of T-ALL patients exhibit dominant active mutations in the IL7R gene (Liu et al. 2017; Shochat et al. 2011; Zenatti et al. 2011), and a significant subset of T-ALL patients with active IL7R signaling display elevated IL7R expression levels (Silva et al. 2021), two types of lesions were introduced: (1) a dominant active IL7R lesion (hereafter called IL7R^DA lesion), (2) an IL7R overexpression lesion (hereafter called IL7R^HI lesion). For the IL7R^DA lesion, we set the IL7R level of the lesioned clone to the maximum WT level of 1 and allowed it to activate IL-7 signaling regardless of extracellular IL-7 cytokine, mimicking the dominant activation mutations in the IL7R gene (Figure 4A, middle panel). For the IL7R^HI lesion, we mimicked overexpression of the IL7R gene by assigning a value of 10 for the IL7R level in the lesioned clone, i.e. 10-fold higher than the maximum attainable in the WT population (Figure 4A, right panel). In both scenarios, non-lesioned thymocytes retained reference IL7R levels randomly ranging between 0 to 1 (Aghaallaei et al. 2021). In each simulation, a single lesioned clone entered the thymic niche after the 60th hour of a simulation, that is shortly after the establishment of a homeostatic population size (Figure 4B, Figure 4 - Supplement Movie 1).

Figure 4

First, we evaluated the impact of IL7R^HI and IL7R^DA lesions in the reference condition, meaning no additional modifications were made to our model. The clonal analysis of simulations revealed that introducing a lesioned clone did not markedly alter the division behavior of non-lesioned clones in the same simulation (scenario with no lesion: 1.04±1.26 rounds of cell division; scenario with IL7R^DA: 0.98±1.22 rounds of cell division; scenario with IL7R^HI: 1.04±1.23 rounds of cell division; Figure 4B-C). The large variance in non-lesioned clones is consistent with our previous work (Aghaallaei et al. 2021), which showed that clones expressing endogenously high levels of IL7R undergo more rounds of cell division in the virtual thymus, whereas clones with lower levels of IL7R proliferate very little or not at all. In contrast, lesioned clones averaged just above 4 rounds of cell division (4.20±0.34 for IL7R^DA and 4.10±0.13 for IL7R^HI). Consequently, the size of IL7R^DA and IL7R^HI clones was nearly 8-fold larger than an average non-lesioned clone (Figure 4C). This result suggests that lesioned clones acquired a distinct proliferative advantage over non-lesioned clones in our virtual thymus model. Nevertheless, since the expansion of lesioned clones amounted to, at best, only one additional round of cell division compared to the upper range of the non-lesioned clone distribution (Figure 4C), our virtual thymus model predicts that a single lesion in the IL7R gene will not be clinically impactful.

Therefore, we next explored whether additional modifications in both cell-autonomous and non-autonomous factors might substantially enhance the clonal expansion. To identify these factors, we decided to undertake a systematic approach and tested combinations of scenarios affecting the TEC architecture and IL-7 signaling parameters, as tested in Figures 2 and 3, together with changing other parameters affecting proliferation and differentiation of all clones (Figure 4D). Furthermore, we considered scenarios in which only lesioned clones exhibited additional modifications, such as reduced cell motility, slower differentiation time, or autocrine IL-7 production. The inclusion of the latter was inspired by a study showing that certain malignant T cells derived from T-ALL patients possess the ability to ectopically express IL-7 cytokine (Buffière et al. 2019). In theory, there are 104976 possible permutations of these scenarios. To reduce the space of possibilities and thus the computational cost, we prioritized alterations expected to increase proliferative potential based on outcomes shown in Figures 2 and 3. We also decided to focus on simulating the extremes; e.g., low and high levels shown in Figure 4D. In addition, we simulated a small sample of scenarios expected to decrease proliferative potential, such as a reduction in TEC density. In the end, 1580 scenarios were simulated. The clonal analysis tool was used to compare clone size and the number of cell division rounds between lesioned and non-lesioned clones under the same conditions. Overall, we observed a similar proliferation pattern among IL7R^HI and IL7R^DA clones in most of the tested scenarios. The lesioned clones consistently showed a higher proliferation rate than their counterpart clones (Figure 4E; Figure 4 – Source Data 1). This trend persisted even when parameters known to generally increase cell division for all cells, such as mean cycle duration and proliferation duration, were modified (Figure 4 - Supplement 1, Figure 4 – Supplement Movie 2). In fact, lesioned clones disproportionally gained up to 6 rounds of cell division by modulating proliferative parameters, while normal clones only gained 2 rounds of cell division in the same scenarios. This difference results from cell crowding creating a disproportional disadvantage to non-lesioned clones: In the scenario with increased proliferation, the first clones to colonize the organ can massively proliferate and completely surround the TECs, preventing later arriving clones from accessing TEC-derived pro-proliferative signals (Figure 4 – Supplement Movie 2, right column). Thus, these few “lucky” early colonizers prevent late cells from attaining their full proliferative potential. Lesioned clones carry an intrinsic pro-proliferative advantage due to their lesions in the IL7R, enabling proliferation despite a lack of direct contact with TECs. In contrast, modifying parameters that affected IL-7 activity only enhanced the proliferation rate of non-lesioned clones but showed no impact on either IL7R^HI or IL7R^DA clones (Figure 4 - Supplement 2). While this outcome is expected for IL7R^DA clones, which are insensitive to IL-7, for IL7R^HI clones it indicates that 10-fold receptor overexpression is enough to activate downstream pro-proliferative effects of the pathway regardless of the parameter permutations that we tested. In contrast, proliferation of normal clones is strongly affected in these conditions, gaining up to 3 rounds of cell division between least and most proliferative conditions.

In our systematic approach, we did not detect an added effect on clonal expansion from combining reduced cell motility or autocrine IL-7 production with IL7R^HI or IL7R^DA lesions (Figure 4 - Supplement 3). Added IL-7 secretion by lesioned clones with the autocrine lesion had only a small effect on proliferation of non-lesioned clones in the same simulation. The most critical lesion in conferring a proliferative advantage to lesioned clones was a delay in differentiation. This delay essentially doubled the duration of the proliferative phase, and indeed we observed up to 8 rounds of cell division, double the level of clones having only a single lesion in the IL7R. This result fits well with in vivo data suggesting that perturbations in differentiation have been implicated in contributing to the initiation and progression of cancer (Ruijtenberg and van den Heuvel 2016).

Surprisingly, the expansion of lesioned clones was amplified when the TEC niche was at its sparsest – a result opposite to the non-lesioned clones which instead profit from a dense TEC niche (Figure 4 - Supplement 4). While non-lesioned clones gained 2 rounds of division with an increase in TEC density, lesioned clones lost 2 rounds of division. Indeed, parameter combinations with reduced density and size of TECs were overrepresented in the scenarios that most increased proliferation (Figure 4 - Supplement 5A-C). For example, in the most extreme scenario observed in our systematic approach, IL7R^HI or IL7R^DA clones underwent up to 18 rounds of cell division – far above the 12 rounds predicted from modifying each group of parameters in isolation (+6 from proliferation, +4 from added lesions, +2 from reduced TEC density, assuming an additive effect; Figure 4 – Supplement 1-4). This result suggests that a defective niche amplifies the effect of other pro-proliferative modulations such as slower differentiation and shorter cell cycle, conferring a synergistic advantage to IL7R^HI and IL7R^DA clones compared to their normal counterparts. Note that at 18 rounds of division lesioned clones were so massive that the total lesioned cell volume was over a 1000-fold larger than the volume of the simulated organ slice (Figure 4 - Supplement 5D). In comparison, in the homeostatic condition, the total volume occupied by both thymocytes and TECs amounted to only 0.7-fold of the available tissue volume. Though our simulations by default include physical volume exclusion effects that prevent dense packing, we did not implement density-dependent feedback on proliferation, enabling cells to duplicate their volume at cell division and proliferate uncontrolled despite the lack of available space. Presumably, similar conditions in vivo would lead to thymus hyperplasia. Indeed, our experimental manipulations of IL-7 levels in the thymus had a direct effect on organ size (Figure 2I, Figure 2 - Supplement 2C, D).

Identifying a sparse TEC network as a new factor that could influence the clonal expansion of lesioned clones was an unexpected outcome. This is because, in the WT situation, a sparser TEC network results in decreased extracellular IL-7, leading to reduced clonal size (Figure 2). While being disadvantageous for non-lesioned clones, this compromised thymic niche proved to be optimal for a massive expansion of lesioned clones. Conversely, we observed that a denser TEC network – characterized by a higher density of TECs with more protrusions and larger size – positively influenced the population size of non-lesioned cells but had a negative impact on the IL7R^HI and IL7R^DA lesioned clones (Figure 4 - Supplement 4). Similarly, in conditions of high cell crowding, lack of contact with TECs put non-lesioned clones at a disadvantage, while lesioned clones appeared to benefit (Figure 4 – Supplement 1, Figure 4 – Supplement Movie 2). This outcome might be explained by the signaling dynamics within our virtual model, where the engagement of the NOTCH1 receptor on thymocytes with the DLL4 ligand on TECs, on the one hand, promotes proliferation and, on the other hand, is a prerequisite for the differentiation process until commitment to a non-dividing fate. In our virtual model, if the NOTCH1 receptor fails to engage with DLL4 on TECs due to their low density or cell crowding effects, cells will remain in an undifferentiated state for a longer time. Closer inspection of time-normalized clones in scenario 6 (very sparse TECs) and in scenario 26 (very dense TECs) indeed confirms that the thymocyte population has a longer turnover time when TECs are sparse (almost 80h, Figure 2D middle panel), and a shorter turnover when TECs are dense (less than 60h, Figure 2D right panel), indicative of the effect of Notch signaling. Together, cells in an environment with fewer TECs and thus lower Notch signaling will experience a prolonged proliferative phase but still require pro-proliferative signals to commit to the cell cycle. Owing to the higher levels of IL7R signal intrinsic to the lesion, lesioned clones will be more easily competent to proliferate even without Notch signaling and, therefore, benefit from a lower TEC density.

The interplay between IL7R and NOTCH1 signals in the clonal expansion of thymocytes

In our virtual thymus model, ETPs simultaneously assess whether the combined sum of IL7R and NOTCH1 signals surpasses a threshold required for entry into the cell cycle ((Aghaallaei et al. 2021); Figure 5A). Besides its effect on proliferation, NOTCH1 signaling also promotes ETP differentiation (Aghaallaei et al. 2021; 2022), and we used this fact and the observations in notch1b medaka mutant phenotypes to include an accelerating effect of NOTCH1 on differentiation in our model (Aghaallaei et al. 2021). Thus, NOTCH1 has a dual effect in thymocytes, which we implemented in our model as an incoherent feed-forward loop (Figure 5A): NOTCH1 signaling promotes proliferation and accelerates differentiation, but differentiation inhibits proliferation by preventing further entry into the cell cycle. IL7R signaling promotes proliferation independently of NOTCH1 and has no direct or indirect effect on the differentiation process. To further explore the interplay between these two factors, we decided to compare scenarios where lesioned clones exhibited only constitutive activation of the NOTCH1 receptor (hereafter called NOTCH1^DA) or in combination with either IL7R^DA or IL7R^HI. The simulations predicted that lesioned clones with only NOTCH1^DA modification displayed a very small advantage compared to non-lesioned clones; this increase in proliferation was only marginally higher than expected for a clone expressing the highest endogenous levels of IL7R (hereafter IL7R^WT; Figure 5B). In contrast, lesions that specifically delay differentiation led to higher proliferation. Indeed, the addition of NOTCH1^DA modification to IL7R^DA or IL7R^HI lesioned clones did not show any effect on clonal expansion but combining delayed differentiation with IL7R^DA or IL7R^HI doubled proliferation. This result indicates that modulating the effect of NOTCH1 globally does not affect proliferation in our model, which is likely due to feedback inhibition, but specifically targeting the differentiation process leads to clonal expansion.

Figure 5

To test the outcome of the virtual model, we performed a series of in vivo experiments to induce NOTCH1 and IL7R signaling in thymocytes. To constitutively activate NOTCH1 signaling in thymocytes, we took advantage of the previously cloned medaka notch1b intracellular domain (NICD) construct (Aghaallaei et al. 2021; 2022), which was driven by a thymocyte-specific promoter (Bajoghli et al. 2015). Furthermore, the destabilizing PEST domain from the NICD was removed in this construct (hereafter referred to as NICD^ΔPEST; Figure 5C). This genetic modification was made due to the known impact of the PEST domain on protein stability, as nonsense mutations lacking the PEST domain of the human NOTCH1 gene have been frequently observed in T-ALL patients (Breit et al. 2006; Ferrando 2009; Weng et al. 2004) (Figure 5 - Supplement 1A). To mimic the in silico IL7R^HI clones, we employed a previously developed construct wherein a thymocyte-specific promoter drives the medaka full length il7r cDNA (Aghaallaei et al. 2021). We then performed mutation in the extracellular juxtamembrane-transmembrane region of the medaka il7r cDNA (Figure 5 - Supplement 1D) to develop a dominant active IL7R form, akin to the NPC mutation found in the IL7R gene in some T-ALL patients (Zenatti et al. 2011; Oliveira et al. 2022), hereafter called il7r^NPC. In our experimental setup, the promoter also co-expressed GFP, which allowed us to (i) identify thymocytes expressing the oncogenes, (ii) determine the clonal expansion of cells expressing the oncogene, and (iii) assess thymus hyperplasia and infiltration into other organs, a characteristic feature for T-ALL. DNA constructs were then injected into blastomeres of embryos at one-cell stage and they were observed during their development using live imaging, with a focus on two specific phenotypes: Firstly, we assessed whether the thymus exhibited enlargement beyond the normal size for its developmental stage. In particular, thymus lymphoma was defined as hyperplasia with an increase of more than twice the organ’s typical size. Secondly, we identified T-ALL by observing infiltration of GFP-co-expressing malignant cells in other organs, including the brain, intestine, or heart (Figure 5D). To align our in vivo experiments with the in silico conditions, we limited our observations to 11 days post-fertilization (dpf), concluding the experiments at the freshly hatched yolk-sac larval stage. Since thymopoiesis begins at 3 dpf in medaka embryos (Bajoghli et al. 2015; 2009), this means we monitored thymus growth for a time span of 8 days.

None of the embryos injected with the NICD^ΔPEST construct (N=32) displayed thymus hyperplasia at 11 dpf (Figure 5E), supporting our simulation outcome showing that global activation of NOTCH1 signaling in thymocytes alone is insufficient to result in clonal expansion within a short time frame of 8 days. The MYC oncogene is an endogenous downstream target of NOTCH1 and plays a major role in NOTCH1-induced transformation (Sanchez-Martin and Ferrando 2017). Therefore, to attempt to shift the balance of NOTCH1 action towards proliferation, we decided to introduce the medaka mycn cDNA into our construct. Consequently, we observed that 25% of injected embryos with NICD^ΔPEST and mycn (N=36) displayed thymus hyperplasia at 11 dpf (Figure 5E). Among them, 8% also exhibited a massive infiltration of GFP-expressing cells in other organs such as the brain and gut. Further analysis of sorted GFP-expressing thymocytes revealed ectopic il7 expression, while the expression level of endogenous il7r remained unchanged (Figure 5 - Supplement 1B). Whole-mount in situ hybridization (WISH) analysis further confirmed a robust upregulation of il7 expression in the thymus of these embryos (Figure 5 - Supplement 1C). MYCN is a basic helix-loop-helix transcription factor that is downstream of several pro-proliferative signaling pathways (Ruiz-Pérez, Henley, and Arsenian-Henriksson 2017), including NOTCH1 in the context of thymocytes (Sanchez-Martin and Ferrando 2017)). Together with our experimental results, it is therefore likely that MYCN could enact the proliferative effect of NOTCH1, shifting the balance in the incoherent feed-forward loop towards proliferation. Interestingly, our results also indicate that transformed thymocytes acquire the ability to release IL-7 (Figure 5 - Supplement 1B, C), which may further promote proliferation in an autocrine fashion.

To test whether constitutive IL7R activation could enhance the development of thymus hyperplasia and T-ALL, we injected various constructs to overexpress either il7r or the dominant active il7r^NPC alone, or in combination with NICD^ΔPEST and mycn. Thymus hyperplasia was found in a slightly higher frequency of 26% (N=49) or 59% (N=29) of embryos, respectively, when only il7r or il7r^NPC was overexpressed in thymocytes. Additionally, we FACS-sorted thymocytes and performed a qPCR for lck, since this gene was shown to be upregulated in zebrafish after IL7R^DA T-ALL development using RNASeq data (Oliveira et al. 2022). Indeed, we detected a significant upregulation of lck upon il7r^NPC overexpression when compared to the WT counterparts. However, upregulation was not significantly different when the WT was compared to il7r^HI condition (Figure 5 - Supplement 1D). Combining il7r^NPC and NICD^ΔPEST and mycn yielded a notably elevated frequency of ∼73% (N=33) for thymus hyperplasia (Figure 5E). Of this group, 33% also exhibited a T-ALL phenotype, suggesting an additive effect of these factors in the T-ALL development (Figure 5 - Supplement Table 1). GFP and pH3 double staining of embryos exhibiting the T-ALL phenotype further confirmed that the combination of NICD^ΔPEST, mycn and il7r overexpression in thymocytes resulted in a higher number of mitotically active cells within the thymus (Figure 5F). Notably, we observed many pH3-stained cells in the inner zone of the medaka thymus (Figure 5G, arrows), an area where WT thymocytes are mitotically quiescent (Bajoghli et al. 2015; Aghaallaei et al. 2021).

The prediction that higher NOTCH1 signaling does not lead to increased proliferation is likely due to the incoherent feed-forward loop downstream of NOTCH1, which shuts down excessive proliferation by triggering differentiation to a non-dividing cell fate (Figure 5B). We observed a similar effect in our in vivo experiments, where activation of NOTCH1 alone did not lead to thymus hyperplasia within the observed time window of 8 days (Figure 5E). In contrast, combining NOTCH1 with its downstream pro-proliferative effector mycn led to thymus hyperplasia. Intriguingly, in vivo, continuous activation of NOTCH1 signaling in thymocytes is one of the main drivers of T cell leukemogenesis and T-ALL development (Weng et al. 2004; Lin et al. 2006; Chen et al. 2007; Liu et al. 2017; Neumann et al. 2014). Several in vivo studies have demonstrated that constitutive activation of NOTCH1 receptor leads to leukemia development after several weeks in zebrafish (Chen et al. 2007; Blackburn et al. 2012) or months in mice (Chiang et al. 2008; Hu et al. 2009; Sharma et al. 2006; Wendorff and Ferrando 2020). Indeed, we did not observe thymus hyperplasia when upregulating NOTCH1 signaling alone, which may stem from the shorter time window of our in vivo experiments (i.e., 8 days). However, in the virtual thymus model, a sole upregulation of NOTCH1 cannot produce thymus hyperplasia regardless of the duration. Different possibilities could reconcile the discrepancy between in vivo and the computational model. Firstly, the balance between pro-proliferative and pro-differentiation effects of NOTCH1 is different from what we modeled. Secondly, there is additional negative feedback between proliferation and differentiation, which our model does not consider. Thirdly, proliferation in vivo can continue even past differentiation if NOTCH1 signaling is sufficiently stimulated.

Overall, our in vivo results support the outcomes of our simulations, showing that overexpression of IL7R alone, but not NOTCH1 alone, is sufficient to induce thymus hyperplasia in a very short time period. However, dual activation of IL7R and the NOTCH1-MYCN axis can promote clonal expansion in vivo, leading to the rapid development of thymus hyperplasia and T-ALL in the medaka model system within eight days.

An abundance of thymic IL-7 milieu could accelerate the T-ALL development

One of the most unexpected outcomes of our simulations was that a defective TEC network provides lesioned clones with a substantial advantage in proliferation. Given the crucial role of IL-7 in proliferation of thymocytes, and considering that TECs, and not thymocytes, serve as the primary source of IL-7 cytokine in the normal thymic milieu, we next asked whether elevated il7 expression in TECs alone is sufficient to induce thymus hyperplasia and to enhance T-ALL development in a thymus where thymocytes express high levels of il7r. This question was also driven by the results of scenario 27, which showed a twofold increase in thymocyte population size when all TECs express IL-7 (Figure 2E-G). In this thymic niche enriched with IL-7, lesioned clones either with IL-7^HI or extended differentiation times showed a significant advantage in clonal expansion compared to their WT counterparts (Figure 6A). Further, the impact of alterations in the thymic niche on T-ALL development has not been previously addressed. Therefore, we monitored the thymus development of the TEC:IL-7^HI embryos until 11 dpf to mimic scenario 27 in vivo. We found that 12% of these embryos (N=24) displayed signs of thymus hyperplasia (Figure 6B), however, we did not observe massive infiltration of GFP-expressing cells into other organs. The frequency of thymus hyperplasia was increased to 41% (N=32) when DNA constructs designed to overexpress il7 in TECs were co-injected with the construct overexpressing il7r in thymocytes (Figure 6B). Comparably, thymus hyperplasia was only observed in 26% of embryos injected with the IL7R overexpression construct in thymocytes alone (N=49). A synergistic effect was observed when il7 was overexpressed in TEC along with NICD and mycn in thymocytes (Figure 6B).

Figure 6

Interestingly, the frequency of embryos with T-ALL development was increased to 25% (N=48) in this group. Taken together, our in vivo results reveal that an excess of IL-7 in the thymic environment, combined with alterations that affect the differentiation status of thymocytes, triggers the T-ALL development in a short period of time.

Overall, several of our experimental observations indicate that there might be additional negative feedback from pro-proliferative signaling to differentiation (Figure 6C). Most strikingly, manipulations in IL7R or in thymic IL-7 levels are sufficient to induce hyperplasia and T-ALL (Figure 5E and Figure 6B) – an outcome that is much more severe than the model predicts (Figure 5B and Figure 6A). In simulations, we could only obtain such massive overproliferation by modifying multiple parameters, including delaying differentiation, accelerating cell divisions, and reducing TEC density (Figure 4 – Supplement 5), as differentiation in our model is a hard stop to proliferation regardless of pro-proliferative stimuli. Similarly, the experimental observation that NOTCH1-MYCN upregulation leads to thymus hyperplasia and T-ALL development also suggests that proliferation can proceed unbridled by differentiation in vivo. A negative feedback loop from proliferation (or cell cycle entry) on differentiation could explain these discrepancies between the virtual thymus model and in vivo observations. We therefore propose that continued thymocyte proliferation can downregulate differentiation pathways and maintain cells in a constant proliferative state.

Conclusion

This study shows the powerful synergy between computational modeling and in vivo experimental validation in elucidating the complex cellular and molecular interactions within an organ that contribute to the initiation and progression of diseases. Here, we highlight a previously underexplored area – the role of TECs and the thymic niche in the development of T-ALL – bridging a significant gap in our understanding of the initiation phase of the disease. By simulating over 1500 scenarios, we identified potential drivers such as alterations in the shape, density, and network of TECs that significantly influence the clonal expansion of lesioned clones, aspects which are challenging to manipulate in experimental setups. Consistent with the outcome of our simulations, our in vivo results reveal that an enriched thymic IL-7 milieu can accelerate disease progression, particularly when thymocytes display enhanced IL7R levels or continuous NOTCH1 activation and MYCN overexpression. This challenges traditional gene-centric views of the disease and highlights the importance of spatial structure of the thymic niche for understanding the mechanisms of T-ALL initiation and progression. It is important to note the limitations of our cell-based computational model, calibrated to mimic the embryonic thymus of medaka (Bajoghli et al. 2015; 2009; Aghaallaei et al. 2021), which might not scale accurately to different species. Therefore, further experimental studies, particularly using murine models, are essential to comprehensively investigate the role of TECs in T-ALL development. Despite simplifying the biological complexity inherent in modeling T-cell development, our integrated approach significantly enhances our understanding of environmental influences on disease progression. It also could open new avenues for exploring how the environment contributes to chemoresistance and disease relapse. This could lead to innovative therapeutic strategies targeting the leukemic-supportive niche, offering hope for improved treatment outcomes.

Materials and methods

In silico model

To develop, implement, and simulate the virtual thymus model, we used the modeling and simulation software EPISIM (Sütterlin et al. 2017; 2013). This multiscale simulation software uses an agent-based paradigm to represent cells as spheres or ellipsoids in three-dimensional space, enabling each of the cells to individually perform internal processes based on flow diagrams and logic rules (e.g. if/else statements) or differential equations, and also implements a partial differential equation solver to simulate diffusion of chemicals in the extracellular space. The virtual thymus model has been comprehensively described in the supplementary material to our previous study (Aghaallaei et al. 2021). In this work, we used this model as a baseline to introduce a small number of additions, explained in the following.

Rate of IL-7 pathway signal transduction

The equation describing the IL-7 signaling rate in each cell was rescaled to introduce the parameter a_IL-7 for the IL-7 signal transduction activation rate. The default value of a_IL-7 was chosen such that the model behavior did not change. The rescaled equation for IL-7 signal transduction reads:

where σ_IL-7 is the IL-7 pathway signal transduction activity, [IL7R] is the IL7R concentration in the given cell,〈[IL-7_ex]〉is the mean extracellular IL-7 concentration in the cell’s microenvironment, and d_IL-7 is the signal transduction deactivation rate.

In the scenario where we introduced a lesioned clone with dominant-active IL7R, IL-7 signal transduction activity was calculated as

Thus, in IL7R^DA clones, neither extracellular IL-7 nor the signal transduction activation rate a_IL-7 had an impact on cells’ IL-7 signaling pathway activity.

Depletion of extracellular IL-7

In our original model (Aghaallaei et al. 2021), IL-7 ligand is released back to the extracellular space after binding to IL7R. Thus, thymocytes do not affect the extracellular IL-7 gradient, which is given by the partial differential equation

where D_IL-7 is the IL-7 diffusion constant, k_IL-7 is a baseline level of extracellular IL-7 degradation, and s_i is a source term accounting for IL-7 secretion, e.g. by a subset of TECs. In the modeling platform we use, secretion occurs throughout a cell’s volume; more specifically, for solving the partial differential equation, space is discretized into voxels, and all voxels i that intersect with the ellipsoids that are used to represent a cell contribute to the source term.

In this work, we introduce the option, controlled via a Boolean flag DEPL that is set by the user before the simulation starts, to simulate a scenario where thymocytes internalize IL-7 ligand after it is bound to IL7R. In this scenario, the extracellular IL-7 concentration is given by:

Note that the additional sink term in Equation 4 represents extracellular IL-7 depletion and is simply the IL-7 signal transduction activity from Equation 1 summed over all voxels j that intersect with thymocytes. We assume that all internalized IL-7 is permanently removed from the extracellular pool.

For lesioned IL7R^DA clones, we used the alternative sink term

which is summed over all voxels k that intersect with thymocytes of the lesioned clone. This term corresponds to the IL-7 signal transduction activity as shown in Equation 2. Reduction of extracellular IL-7 below zero was prevented by a computational check.

Determining clonal lineages

The EPISIM software we used assigns to each cell a unique cell identifier. We made use of this unique cell identifier to generate unique clonal identifiers as follows. For each thymocyte that entered the simulation via migration from outside the thymus, we used the value of their unique cell identifier to assign the value of a new clonal identifier variable. This cell-specific value of the clonal identifier variable was passed on to any daughter cells produced via cell division. Thus, using the clonal identifier, we essentially generated an in silico clonal label which we used to determine clonally related cells.

Lesioned clone

To generate a lesioned clone in a given simulation, we created a computational label to mark a single cell selected randomly among new thymic immigrants once the thymic cell population reached its homeostatic size (t >= 60h). This label marked this cell and its progeny as a “lesioned clone”, such that we could specifically modulate parameters and introduce new rules only for this single clonal lineage. We tested the effect of several alterations specific to the lesioned clone (Figure 4 - Supplement Table 1).

Parameter variations

Please refer to the supplementary material of (Aghaallaei et al. 2021) for a comprehensive explanation of all parameters and the choice of the reference values. Altered parameters used in this work are listed in Figure 1 - Supplement Table 1.

Calculation of simulated clone size

Clone size was defined as the number of terminal leaves in a cell lineage tree (see scheme in Figure 4C); terminal leaves were counted regardless of the ultimate fate of those cells (i.e. positive or negative selection). All clones of the same genotype (wildtype or lesioned) belonging to the same simulated scenario were averaged across all replicate simulations of that scenario. This calculation gives the average number of cells per clone in a given scenario. The number of cells was log2-transformed to represent rounds of cell division.

In vivo model

Medaka (Oryzias latipes) husbandry was performed in accordance with the German animal welfare standards (Tierschutzgesetz §11, Abs. 1, Nr. 1, husbandry permit no. 35/9185.46/Uni TÜ.) The transgenic line (tg) lck:gfp was described previously (Bajoghli et al. 2015). All experiments conducted in medaka embryos were performed prior to the legal onset of animal life stages under protection, utilizing both males and females of the wildtype line Cab (Loosli et al. 2000).

Cloning of DNA constructs

To overexpress il7 in TECs, we used the construct ccl25a:il7, which was described previously (Aghaallaei et al. 2021). To overexpress il7r in thymocytes, the full-length medaka il7r cDNA (Aghaallaei et al. 2021) was cloned into vectors containing a medaka thymocyte-specific promoter (Bajoghli et al. 2015) that drives a fluorescent protein (sfGFP or TagRFP). To generate a dominant active form of il7r, we introduced three amino acids, Asparagine (N), Proline (P), and Cysteine (C) (known as NPC mutation, as shown in (Zenatti et al. 2011) and (Oliveira et al. 2022)), in the medaka il7r extracellular juxtamembrane-transmembrane interface region after position 266 using site-directed mutagenesis (Figure 5 - Supplement Figure 1D), this corresponds to the human NPC mutation in the position 242 of IL7R gene. To activate Notch signaling in thymocytes, we utilized the medaka notch1b intracellular domain (NICD), as described previously (Aghaallaei et al. 2021), and removed its PEST domain to enhance protein stability, mimicking the situation reported in most T-ALL patients harboring gain-of-function mutations in the NOTCH1 gene (Weng et al. 2004). For overexpressing the mycn oncogene, we first isolated the full-length medaka mycn cDNA (accession number: ENSORLG00000022362) and introduced a mutation at position 44 from Proline to Leucine (P44L) using site-directed mutagenesis, as identified in T-ALL patients (Liu et al. 2017). The DNA fragments were then cloned into vectors containing a medaka thymocyte-specific promoter, either alone or in combination (Bajoghli et al. 2015).

DNA micro-injection

Plasmids at concentration 10-25 ng/µl together with I-SceI meganuclease and NEB buffer (NewEngland BioLabs), were co-injected into the blastomere at one-cell stage embryos. Fluorescent signals, based on the constructs either sfGFP, tagRFP or mTurquoise, were used to select positive embryos.

Generation and genotyping of medaka il7 mutant

The CRISPR-Cas9 approach was employed to generate medaka il7 crispant. Crispr RNA (crRNA AGTAGACTGATGCAAAGAAG) was designed using the CCTop website (Stemmer et al. 2015; 2017) and ordered from IDT. RNP complexes were prepared as described before (Hoshijima et al. 2019) using AltR tracrRNA (IDT) and Alt-R S.p. Cas9 nuclease (IDT). The injection mixture contained a total of 25 µM crRNA:tracrRNA duplex and 25 µM Cas9. Injection was performed into the blastomere at the one-cell stage transgenic embryos carrying the lck:gfp reporter. Injected embryos were raised until 8 dpf and then fixed in 4% PFA/2xPBS + 0.1% Tween20 for whole-mount immunostaining. To correlate the phenotype with the genotype, each embryo was genotyped using PCR and Sanger sequencing

Phenotype assessment

The development of injected embryos was monitored using a NIKON SMZ18 stereo-fluorescent microscope. Only embryos with fluorescent signals in the thymus region were selected for further analysis. Thymuses of embryos before hatching or freshly hatched yolk sac larvae were examined. Thymuses with at least a two-fold increase in size compared to transgenic embryos carrying the lck:gfp reporter construct (Bajoghli et al. 2015) were considered as having thymic hyperplasia. The identification of cells expressing fluorescent proteins and oncogenes outside the thymus, such as in the brain, gut, and heart, was considered indicative of the T-ALL phenotype.

Whole-mount immunostaining

Immunostaining procedures followed established protocols (Inoue and Wittbrodt 2011; Aghaallaei et al. 2021). Briefly, mitotically active cells were identified using a rabbit anti-phosphohistone-3 antibody (Ser10, Millipore 06–570, 1:500 dilution), with a Cy3-donkey anti-rabbit immunoglobulin G secondary antibody (the Jackson laboratories, 711-165-152; 1:500 dilution). GFP expression was detected using a goat Anti-GFP antibody (Abcam, ab5450, 1:500 dilution), with Alexa488-donkey anti-goat IgG (Abcam, ab15129, 1:600 dilution) as the secondary antibody. To quantify the number of pH3-expressing cells, the entire thymus region was imaged using a LSM710 (Zeiss) confocal microscope with z-stacks (z=1 µm).

Whole-mount in situ hybridization (WISH)

Whole-mount in situ hybridization in medaka embryos was carried out using digoxigenin-labeled RNA probes, as previously described (Aghaallaei et al. 2005). Probes used in this study for il7r and il7 were previously described (Aghaallaei et al. 2021).

Cell sorting and quantitative RT-PCR

Injected embryos were first smashed through a 40 µM strainer (Greiner) and cells were collected in 0.9x PBS with 500 U.I. Heparinum Natricum (Liquemin). Sorting was then performed using the Cell Sorter MA900 (Sony Biotechnology). WT embryos were used as negative sorting control. Only GFP positive cells were collected and used for RNA isolation. RNA of sorted cells was isolated by using NucleoSpin RNA XS (Macherey-Nagel) and 20 ng carrier RNA, following the manufacturer’s protocol. RNA was treated with rDNase and eluted in 14 µl dAqua. The first-strand cDNA synthesis was carried out with random hexamer primers and SuperScript III Reverse Transcriptase (Thermo Fisher Scientific) by following the manufacturer’s protocol. SYBR Green Kit (Applied Biosystems) was used for Quantitative PCR on the LightCycler 480 (Roche). The data was evaluated in Microsoft Excel using the ΔCt method and normalized to the housekeeping gene ef1a. The primers used for RT-PCR have been described previously (Aghaallaei et al. 2021).

Statistical analysis

Statistical analysis was performed using RStudio version 2024.4.2.764 (Posit team 2024) using R version 4.4.0 (R Core Team 2024), the R library jtools (Long 2022) and GraphPad Prism version 8.0.2. The two-tailed Fisher’s exact t-test was used to compare the thymus hyperplasia phenotype with normal thymus. Welch’s t-test was used for thymus volume size comparison in WT and il7 crispants, multiple t-test was used for comparison of normal thymus with thymus hyperplasia and T-ALL phenotype. Mann-Whitney test was used to evaluate the qPCRs. p-values <0.05 were considered as statistically significant. The number of biological samples for each experiment (N) is indicated in the figures, and bar graphs present the absolute numbers with mean ± SD, ± SEM or the percentage as indicated in the figure legends.

Acknowledgements

The authors thank Anna-Sophia Hellmuth for assistance with il7 gRNA injection; the institute of Medical Virology and Microbiology for continuous support in confocal microscopy; Larissa Doll for support and encouragement. Initial computational work presented here was performed using the ALICE computer resources provided by Leiden University, later work was performed using the computer lab in the Theoretical Biology and Bioinformatics department of Utrecht University. This work was supported by the Madeleine Schickedanz Kinderkrebsstiftung (grant number D.30.28666), and Deutsche Forschungsgemeinschaft (BA 5766/5-1). ET is funded by the Dutch Research Council (NWO) in the NWO Talent Programme with project number VI.Veni.222.323. We thank members of the Theoretical Biology and Bioinformatics division for their constructive feedback on the manuscript.

Author contributions

E.T. designed and developed the computational model. A.M.D. performed the in vivo experiments and analyzed their results. B.B. designed and supervised the study. All authors jointly interpreted the data and wrote the manuscript.

Competing interests

The authors declare that they have no competing interests.

Data and materials availability

All data needed to evaluate the conclusion in the paper are present in the paper and/or the Supplementary Materials. The computational model is available on Zenodo: doi:10.5281/zenodo.11656320 (Tsingos 2024). Additional materials related to this paper such as plasmids or DNA constructs may be requested from the authors.