Single-Cell Analysis Identifies Thymic Maturation Delay in Growth-Restricted Neonatal Mice

Fetal growth restriction (FGR) causes a wide variety of defects in the neonate which can lead to increased risk of heart disease, diabetes, anxiety and other disorders later in life. However, the effect of FGR on the immune system, is poorly understood. We used a well-characterized mouse model of FGR in which placental Igf-2 production is lost due to deletion of the placental specific Igf-2 P0 promotor. The thymi in such animals were reduced in mass with a ~70% reduction in cellularity. We used single cell RNA sequencing (Drop-Seq) to analyze 7,264 thymus cells collected at postnatal day 6. We identified considerable heterogeneity among the Cd8/Cd4 double positive cells with one subcluster showing marked upregulation of transcripts encoding a sub-set of proteins that contribute to the surface of the ribosome. The cells from the FGR animals were underrepresented in this cluster. Furthermore, the distribution of cells from the FGR animals was skewed with a higher proportion of immature double negative cells and fewer mature T-cells. Cell cycle regulator transcripts also varied across clusters. The T-cell deficit in FGR mice persisted into adulthood, even when body and organ weights approached normal levels due to catch-up growth. This finding complements the altered immunity found in growth restricted human infants. This reduction in T-cellularity may have implications for adult immunity, adding to the list of adult conditions in which the in utero environment is a contributory factor.


INTRODUCTION
Fetal growth restriction (FGR) is a serious condition affecting between 3 and 7% of infants in developed nations (1) but can be higher than 40% in low-income nations (2). The growth restricted infants suffer from an increased risk of a wide variety of acute complications, including necrotizing enterocolitis, perinatal asphyxia, hypothermia, hypoglycemia, immunodeficiency, and death (3). FGR is commonly attributed to insufficient fetal nutrients. This may be because either the mother is malnourished, as in the prototypical case of the Dutch Hunger Winter (4), or because the placenta is dysfunctional, leading to poor nutrient transport to the fetus (5). Severely growth restricted infants exhibit asymmetric growth wherein resources are reallocated to the fetal brain at the cost of the rest of the body (6), so called "brain sparing." The infants undergo "catch-up growth" for years, although some may remain small for life (7). Unfortunately, this rapid growth often fails to overcome the complications of intrauterine restriction, as a wide variety of morbidities including diabetes (8), polycystic ovary syndrome (9), or perturbed mental health (10) have been linked to FGR or catch-up growth (11). Together, these support the "Developmental Origins of Adult Disease" paradigm, wherein adverse factors during critical periods of fetal growth irreversibly impact development (11,12). To overcome these effects, a clear understanding of the mechanism for each of these sequelae is necessary (13).
Several animal models have been developed for the study of adult diseases originating in the uterus. Most of these involve restricting certain nutritional components, such as protein (14) or oxygen (15) in pregnant dams. However, this does not distinguish between maternal placental effects. The Igf-2 P0 model overcomes this through knock-out of placenta-specific insulinlike growth factor 2. There are four Igf-2 isoforms which result from the use of different promoters although they ultimately generate the same protein. The P 0 promoter is specific to the placenta. This gene is paternally imprinted, allowing for generation of both wildtype and affected offspring within the same litter. Importantly, all offspring develop in a wildtype dam, preventing maternal variables from affecting development. This targeted knock-out reduces placental growth and therefore the nutrient transport to the fetus, resulting in a brain-sparing phenotype reminiscent of human FGR (16). Early hypocalcemia in the fetuses of these mice (17) mimics the hypocalcemia found in human neonates (18). These mice have been shown to develop anxiety later in life (19), which recapitulates known long-term effects of FGR on mental health (20).
While long-term effects are a subject of much interest, most acute FGR complications are simply attributed to a lack of tissue mass and developmental delay: for example, a smaller and less mature kidney (21), pancreas (22), or bowel (23) will simply not function as well. Adaptive immunity is mediated by T-cells which develop in the thymus. However, while the thymus is a short-lived organ which involutes shortly after birth it continues to function well into adult life (24). Deleterious effects on this transient organ could, therefore, have a significant and irreversible impact on immunity in adult life. Initially, infants with FGR have acutely smaller thymi and altered CD4/CD8 ratios of peripheral T cells (25). Later in life, FGR is associated with abnormal responses to vaccines and higher rates of death due to infection (26). For example, indirect evidence comes from a study showing that young adults born in the annual "hungry season" in rural Gambia-and therefore likely to be born with FGR-have a 10fold higher risk of premature death, largely due to infection (27).
At a cellular level, broad definitions classify cells based on discrete cell-surface markers. In T-cell development, lymphoid progenitors travel from the bone marrow through the bloodstream to arrive at the thymus where NOTCH signaling directs them toward the T-cell lineage (28). These cells divide and differentiate through four stages of DN (double negative, referring to lack of either CD4 or CD8 T-cell surface markers) whilst undergoing Tcr β rearrangements in the thymic cortex. With the correct set of signals, a minority of the cells become γδ T cells. The rest gain CD8 and then CD4 surface makers (DP, double-positive) at which point Tcr α rearrangements take place. The rearrangements provide the diverse T-cell receptor repertoire necessary for antigen recognition. The cell is then activated by either MHC Class I or Class II molecules presenting self-antigens. Strong affinity for either MHC Class I or Class II leads to apoptosis to prevent auto-immune disease. Weak affinity to Class I or II leads to differentiation into either CD8+ cytotoxic T cells or CD4+ T helper cells, respectively. Intermediate affinity to Class II molecules leads to differentiation into a T-regulatory cell. Errors anywhere along this pathway of T-cell differentiation, from insufficient Tcr rearrangements to underrepresentation of any of the T cell lineages, can lead to impaired immune function (29,30).
Single-cell RNA sequencing, for example Drop-Seq (31), In-Drop, or the commercial 10X Genomics and Dolomite platforms allow the analysis of the transcriptomes of thousands of singlecells (32). These analyses have been invaluable for identifying immune-cell subtypes within populations traditionally classified by discrete cell surface markers (33) and revealed new regulatory pathways (34).
Here, we used a previously established murine model of FGR in order to assess the effect of an adverse in utero environment on neonatal and adult immunity. The Igf-2 P0 (P 0 ) transcript is imprinted and paternally expressed leading to specific loss of placental Igf2 and growth restriction in the fetuses carrying a P 0 transcript deletion (16). We then used Drop-Seq to profile the transcriptomes of 7,264 cells from neonatal thymi in order to characterize possible immune perturbations at a cellular level.

Mouse Tissue Preparation
Mice were maintained at Central Biomedical Services in accordance with the UK Home Office, Animals (Scientific Procedures) Act 1986 which mandates ethical review. C57BL/6 dams were mated with Igf-2 P0 heterozygous males. Mice were genotyped as previously described (16). For neonatal analysis, mice were weighed 5 days after birth prior to organ collection. For adult mice, body and organ weights were assessed at 8 weeks of age. Spleens were digested with collagenase D (11088866001, Sigma, Gillingham, UK) at 1mg/ml in HBSS (24020091, ThermoFisher, Paisley, UK) for 30 min at 37 • C and washed with PBS (D8537, Sigma). Single cell thymic suspensions were prepared by maceration of the thymi through a 40 µm strainer. Red blood cell lysis was performed with PharmLyse (555899, BD Biosciences, Wokingham, UK) for 3 min at RT for both spleens and thymi. Resultant suspensions were passed through a 40 µm strainer and live cells counted using trypan blue (T8154-20ML, Sigma).

Flow Cytometry
Isolated cells were incubated with antibodies listed in Table S1 on ice for 30 min in PBS prior to washing. Cell types were identified as in Table S2. Samples were processed on a BD LSR Fortessa and the data was analyzed using FlowJo 10. Neonatal FIGURE 1 | FGR leads to Brain-Sparing at the Expense of Thymus. WT dams were mated with P 0 males to produce P 0 and WT offspring. The number of WT and P 0 animals or organs analyzed for each measurement is given in parentheses. (A) Neonates were analyzed 5 days after birth for absolute body weight (n = 61, 19). Relative body and organ weights from P 0 neonates were compared with littermate controls. Purple dashed line was defined by the average P 0 body weight, to allow relative comparison of organ weights to both P 0 and WT. (Relative Body Weight n = 41, 18; Brain n = 7, 5; Spleen n = 32, 11; Thymus n = 28, 12). Organ/Body ratios were compared between P 0 and WT littermates. (B) Absolute cell counts within the spleen and thymus were compared with littermate controls (Spleen Cells n = 16, 8; Thymus Cells n = 7, 5). Cell counts as a ratio to mouse weights were compared with littermate controls. (C) Spleen cells were analyzed using flow cytometry and the proportions of the indicated immune cell populations were compared with littermate controls (Spleen Flow Cytometry n = 10, 5, with 1 outlier removed using the ROUT method). Datasets were tested for normality using the Shapiro-Wilk normality test. Normally distributed data were compared using unpaired two-tailed t-tests while non-parametric data were compared using Mann-Whitney Rank comparison tests. For all tests, * P < 0.05. and adult samples were analyzed using the same gating strategy (Figures S1-S3).

Statistical Analysis
For fold change calculations, Igf-2 P0 mice were compared to litter mate controls with each metric being divided by the mean value obtained from WT litter mates. Outliers were identified within each data-set using the ROUT method in GraphPad Prism where Q = 1. This identified only one outlier (Splenic immune subsets, Figure 1, P 0 ) which was subsequently removed. Datasets were tested for normality using the Shapiro-Wilk normality test. Normally distributed data were compared using unpaired twotailed Student's t-test while non-parametric data were compared using the Mann-Whitney Rank comparison test. Cell cycle data was analyzed using a two-way ANOVA with Tukey's multiple comparisons test across cell types or multiple t-tests using Twostage step-up method of Benjamini, Krieger, and Yekutiele across genotypes.

Drop-Seq
Thymi were dissected from postnatal day (PND) 6 neonates and homogenized in PBS before passing through a 40 µm filter. Cells were counted and diluted to 200,000 cells/mL in PBS/0.1% BSA (A7906-100G, Sigma) for processing according to Drop-Seq protocol (31). Single cells were encapsulated in droplets with primer barcoded beads (MACOSKO-2011-10, ChemGenes, Wilmington, USA). After RNA binding to the oligo-dT portion of the primer, the droplets were broken to pool beads for reverse transcription (EP0743, ThermoFisher). Exonuclease (EN0581, ThermoFisher) treatment removed unused primers from the beads prior to PCR amplification at up to 8,000 beads/tube (∼400 cells/tube). Samples were pooled with up to 4,000 cells/sample and the DNA purified before determining concentration. DNA was fragmented and tagged with Illumina adapter and index sequences using Nextera Sequencing kit (FC-131-109, Illumina, Little Chesterford, UK) and 600 pg DNA/sample. Samples were diluted 1:1 in ddH 2 0 water before double-sided purification with Agencourt Ampure beads (A63881, Beckman Coulter, High Wycombe, UK), using 0.6× beads to remove large fragments and 1× beads twice to remove contaminating primers. Samples were analyzed for size distribution using an Agilent HS DNA kit (5067-4626, Agilent, Stockport UK) and concentration and sequenced on an Illumina HiSeq4000 with an Index read of 8bp, Read 1 of 25 bp and Read 2 of 100 bp. Samples were collected from two separate litters (two WT and one or two P 0 animals per litter). Single-cell barcoded cDNA was generated for each sample with littermates being processed on the same day. Library preparation was done for all samples on the same day using Nextera indices to distinguish between each neonate. All samples were sequenced together in one lane of an Illumina HiSeq4000 flow cell.
Raw Fastq files are demultiplexed (dropseq_demultiplex.sh) using the Nextera indices and then converted to uBAM using PicardTools:FastqToSam (v2.9.0). Quality control, alignment (STAR v020201) gene quantification and final matrix generation were performed using DropSeqTools (v1.12, http://mccarrolllab. com/dropseq). Alignments were performed against the mouse reference genome (mm10 available from http://mccarrolllab. com/dropseq). Initial thresholds of a minimum 200 genes per cell and genes must be present in at least 3 cells were applied. The resulting digital expression matrix (DEM) was imported into Seurat (35) (v2.3.0) for downstream analysis and is split into two scripts. The first, dropseq_seurat_splitDEMs.R, performs the more computationally intensive tasks intended to be run on high performance computers, the Seurat object is saved in Robj format to be imported in to second script, dropseq_seurat_splitDEMs_Plots.R, for differential transcript analysis and figure creation.
Two separate DEMs were calculated for the WT and WT+P 0 samples. The WT only samples were used to calculate variable genes (FindVariableGenes), which were then used as input to generate the PCA (RunPCA), find clusters (FindClusters) and produce a tSNE (t-distributed stochastic neighbor embedding) visualization (RunTSNE) from the combined WT and P0 sample DEM. FindClusters is run across multiple resolutions (0.2, 0.4, 0.6, 0.8, and 1.0), each stored on the Seurat Object. Normalization (NormalizeData), UMI and MT regression (FilterCells) were performed using Seurat including a more stringent threshold of a minimum 300 genes per cell and genes must be present in at least 3 cells was applied. Cell cycle assignments were performed using SCRAN (36) (v1.6.9) on the combined WT+P0 DEM, using an intermediate SingleCellExperiment (v1.0.0) data structure, and then added back to the Seurat Object. Cell cycle genes were regressed out using a subtraction of G2M from S cell cycle scores per cell. The resulting Seurat data object is saved as an RObj for input into the plotting and differential analysis part of the pipeline.
The RObj generated from first dropseq_seurat_splitDEMs.R script are imported into dropseq_seurat_splitDEMs_Plots.R to extract (e.g., with GetCellEmbeddings) the required data for each of the plots in the figures. Custom tSNE plots were generated using ggplot2. Transcript abundance dotplots were generated from AverageExpression extracted from the Seurat object and ggplot2. Cluster trees were generated using clustree (37). Differential transcript analysis was performed by comparing each cluster in turn to all others (FindAllMarkers) and using a log fold change threshold of > 0.7 and adjusted p-value < 0.01. The heatmap (pHeatmap), used the same thresholds, and the top 20 genes for each cluster selected.
A custom tool was created to classify whether ribosomal proteins are localized on the surface or are internal to the ribosome structure (https://github.com/darogan/Ribosomal-Protein). The bioinformatics tool scores each ribosomal subunit by its distance from the subunit core. We used this to identify the subunits present on the surface of the ribosome. The tool also provides features to display and color the subunits in Pymol, and was used to confirm they are accessible on the surface. A further

Placental Insufficiency Leads to Brain Sparing at Expense of Thymus
When crossing wildtype (WT) C57BL/6 females with heterozygous mutant Igf2 P 0 males, both WT and P 0 offspring develop in the WT mother. The P 0 neonates exhibit poor placental and fetal growth. Animals were assessed 5 days after birth for body and immune organ weights, with littermates used to determine fold changes ( Figure 1A). This time-point was chosen as it is early within the murine neonatal period but provided sufficient tissue. Consistent with previous data (16), the overall body weight was reduced while brain weight was only marginally affected. Spleens of FGR P 0 neonates were smaller than WT but this change was proportionate to body weight. However, overall cell number was reduced ( Figure 1B). Further investigation using flow cytometry revealed a number of perturbed lineages, where neutrophils and macrophages were the most affected by FGR (Figure 1C). The thymus suffered markedly more growth restriction than the spleen with disproportionately reduced size and a nearly 70% reduction in cells compared to littermate controls. Given the overall severity in thymocyte reduction, we used Drop-Seq to investigate the effects of FGR at the single-T-cell level. Within the Igf-2 locus, neither WT nor P 0 thymocytes cells contained reads aligned to the P 0 -specific upstream exon U2 (Figure S4), confirming that placental-specific P 0 was not expressed in T cells. This suggests, therefore, that the observed phenotype was a not result of a direct effect of loss of local Igf-2 but rather, insufficient placental supply of nutrients. The Igf2 locus contains a number of transcript variants from the placental P 0 through to P 4 (38). We also examined reads mapping to the other isoforms and did not notice any distinct patterns except that in the P 0 mutant there were a number of reads mapping to an untranscribed region upstream of P 1 .

Drop-Seq Identifies Developing Thymocytes
To investigate the effect of FGR on T-cell development, we performed Drop-Seq analysis of 4,679 WT and 2,585 P 0 thymic cells (n = 4 and 3 animals respectively) with an average read depth of 32,399 reads/cell resulting in an average of 994 transcripts per cell (sequencing statistics are in Table S1). We performed a post-hoc power estimate based on our sequencing metrics (39). Our 7,264 cells have a mean reads (per gene per cell) of 5.8 which is well in excess of the minimum threshold of 0.02. This also holds for the WT (6.6) and P 0 (4.0) experiments if considered separately.
The transcript signatures of each WT cell were used to identify clusters and the genes driving that clustering. We regressed out the signals from cell cycle regulator genes when determining the cluster map, which they otherwise affected due to the intrinsic link between cell cycle and T-cell differentiation (40) (Figure S5). Using the remaining cluster determinants, we mapped both WT and P 0 cells into 8 clusters of 7,264 single cells (Figure 2A). Initially, we used known differentiation markers to define these clusters (Figures 2B,C). We identified the known T-cell maturation pathway, beginning with cells containing high levels of Il2ra (Double Negative, DN) through to Cd4/8 containing cells (Double Positive, DP) and ending with high levels of Ccr7 and Itm2a (T Mature ). In this map dominated by DP cells, we were unfortunately not able to distinguish CD4 and CD8 single positive mature T cells. Bias in identifying clusters was checked using various degrees of resolution, from 0 (identifying only 1 cluster) to 1 (identifying 11; Figure 2D and Figure S6). This shows distinct DN1 separation even at the lowest resolution but heterogeneous DP subclustering. Regardless of resolution, cells clustered into distinct T-cell populations based upon differentiation status (Figure 2E). Of note, we found that the DP cells were highly heterogeneous and could be divided into multiple clusters, suggesting that the traditional definition using two markers is insufficient to fully characterize these cells. Cells in Cluster 7 were identified as contaminating red blood cells due to high hemoglobin mRNA levels, while Cluster 8 contained macrophages, identified by abundant Aif1, Lyz2, Complement, and H2 complex transcripts. As both of these clusters contained fewer than 100 cells, we focused our attention on the larger T-cell clusters.

Drop-Seq Identifies Transcript Signatures in T-Cell Maturation
To further characterize these cell clusters, we identified the top marker genes in each cluster by performing a differential analysis of each cluster compared to all others in an unsupervised manner (adjusted P-value <0.01; Table S2). These genes included both known T cell specific markers such as Cd8 and Ccr9 and cytoskeletal components (Figure 3). We prevented cell cycle genes driving the clustering (by regression), but in the differential expression analysis the unregressed data were used. A number of cell cycle regulators were found to differ between mapped clusters, demonstrating a clear connection between cellular phenotype and the stage of the cell cycle. Cluster 1, containing the DN cells, was characterized by early T-cell markers Il2ra and Hes1. Tcr -C1, and C2 were also found within these cells, likely because only DN cells contain precursor T δ cells. Cluster 2, the early DP cluster, contained T-cell genes Granzyme A and Ctage5, as well as Ranbp1, Nolc1, and Hspd1, which are associated with nuclear pore transport, rRNA processing and protein folding, respectively. Plac8 was also found, which was previously implicated in T-cell development (41) but here specifically found within DN cells. Clusters 3 to 5, which we have identified as increasingly mature DP cells, contained high levels of Tarpp (also known as Arpp21), Rag1 and Rag2, whose proteins are present during TCR rearrangements. During T-cell positive selection, Tarpp levels are reduced (42), which is shown here as the T Mat cluster contains few Tarpp transcripts. These cells also showed increasing levels of Ccr9, the CCL25 receptor expressed on late DN and DP T cells. These known T-cell markers were found alongside Btg2, a cell cycle inhibitor, and Rmnd5a, which is necessary for mitotic division. Clusters 2 & 3 both contained Rrm2, Tyms, and Esco2, which synthesize nucleotides and bind sister chromatids, respectively, and are required for S phase of the cell cycle. The variations of Clusters 2 through 5 then are likely in part due to progressive rounds of proliferation as the T cells mature.
The mature lymphocyte markers Itm2a and Cd52 are specific to Cluster 6, alongside a number of cytoskeletal (Cd2, Tmsb4x, and Tmsb10). Negative growth regulators Slfn2, Gimap4, Ms4a4b, and Shisa5 are also present in high levels in the T Mat cluster, likely involved with positive and/or negative selection. Interestingly, the DN Cluster 1, DP Cluster 4, and T Mat Cluster 6 contained heterogeneous mixtures of cells with high and low levels of transcripts coding for ribosomal proteins.

Cell Cycle and Ribosomal Protein Structure Drive Thymocyte Maturation
We further investigated the trends in cell cycle and ribosomal protein genes. Although we prevented known cell cycle genes from determining the cluster map, we did not remove these genes from the transcript profile of each cell. When comparing transcript signatures between clusters, we detected differences in cell cycle stage, likely due to a link between cell type and proliferation state (see Figure S5 for cluster mapping with and without cell cycle regression). We scored each cell for its G 1 , S, or G 2 /M phase transcript signature ( Figure 4A) and compared cycle phase between cell types (Figure 4B). DN and T Mat stages had lower levels of proliferation than the DP stage. There were no differences in cell cycle overall between WT and P 0 genotypes ( Figure 4C).
A differential transcript analysis for each cluster against all others revealed that within the DP stage, ribosomal protein transcripts were varied, most notably between Clusters 3 and 4. We made the same observation of a differential ribosomal protein transcript population in the unsupervised clustering in Figure 3. The ribosomal protein mRNAs that significantly changed in Cluster 3 were compared against those found in Cluster 4 ( Figure 5A). In each case, cells in Cluster 3 were notable for low ribosomal protein transcript abundance while in Cluster 4 such transcripts were increased. The majority of the transcripts (9 out of 13, 69%) found by this comparison between clusters coded for proteins that contributed to the ribosome surface. Surface proteins are more likely to play a functional or regulatory role than internal proteins. The protein subunits of a ribosome are known to play a key role in directing translation of specific mRNA pools (43). Looking at the overall cluster mapping, we found that the ribosomal transcripts appeared in waves dominating Cluster 4 and dividing Cluster 6 (Figures 5B,C Movie S1).

FGR Prevents Thymocyte Maturation
Having identified each cluster and compared the top gene signatures, we then compared the WT and P 0 thymocytes within each cluster (Figure 6A). Even after accounting for the higher number of WT cells present, we still found a higher proportion of DN cells in the P 0 thymus than WT (Figures 6B,C). The opposite held true with the late stage T Mat cluster, as there were proportionally more T Mat cells in the WT thymus than the P 0 . Overall, there was a general skew of P 0 cells toward earlier, more immature clusters in the DP pathway, while the WT thymus favored the later stage clusters. Using flow cytometry for surface markers, we detected discrete increases in DN cells at the expense of mature CD4+ cells in the P 0 (Figure 6D). Within the smaller DN cell population, there were proportionally more DN3 cells in the P 0 over the WT, which is noteworthy given that rearrangement of V and D chains as well as β selection occur at this stage in T-cell differentiation.

FGR Transcript Signatures
Given the apparent waves of ribosomal transcripts appearing in Figure 3, we looked more deeply within clusters to determine whether FGR changed transcript signatures of seemingly equivalent cell types (Figure 7, Supplemental Tables 7A-E). We first identified genes that distinguished major cell type (DN vs. DP or DP vs. T Mat , Figures 7A,B, individual cluster comparisons in Figure S7, Supplemental Tables 7A1-15). As expected, Cd8 and Cd4 genes dominated for the DP to DN comparison, whilst Cd52 and Itm2a distinguished T Mat from DP in this unsupervised comparison. We then looked within clusters for changes due to FGR (Figures 7C-E, individual cluster comparisons in Figure S8, Supplemental Tables 8A1-6). Within each cluster, ribosomal protein transcripts were found at significantly higher levels in the WT over the P 0 . Xist was also found, which indicates the presence of a female in the P 0 group ( Figure S9 shows the consistent overlap of the female P 0 with the male P 0 neonatal T-cells indicating that sex has little effect on these cells). Interestingly, Igf2 was present in lower levels in the WT DN and T Mat clusters, note that this refers to total Igf2 rather than the placental specific isoform. We then investigated overarching trends across clusters by comparing genes which differ between WT and P 0 cells, within each cluster to another cluster (Figures 7F,G). The ribosomal protein transcripts aligned clearly, demonstrating an overall trend as opposed to a clusterspecific phenomenon.  . Body (n = 7, 4), spleen (n = 7, 4), and thymus (n = 5, 3) weights were compared with littermate controls. Spleen/body and thymus/body weight ratios were also compared relative to littermate controls to account for lower overall body weight in the P 0 mice. (B) Total cells within the spleen (n = 4, 4) and thymus (n = 5, 3) were compared directly with littermate controls as well as relative to overall mouse weight. (C) Surface markers were used to classify T cells within the thymus and the proportions were compared to littermate controls (n = 8, 7). (D) Splenic immune subset proportions were determined from flow cytometry and compared with littermate controls (n = 5, 7). Datasets were tested for normality using the Shapiro-Wilk normality test. Normal data were compared using unpaired two-tailed t-tests while nonparametric data were compared using Mann-Whitney Rank comparison tests. For all tests, *P < 0.05.

Immune Defect Lasts into Adulthood
Thymus size is correlated with thymus function (44). Humans with DiGeorge syndrome have little or no thymus, resulting in a significant immunodeficiency (45,46). We wondered whether the maturation delay in the thymus of the 6-day old neonate would affect the immune cell population in the adult. Infants with FGR exhibit increased growth and "catch up" relative to their unrestricted counterparts (47). Adult male mice displayed signs of catch-up growth, although they were still found to be significantly lighter than their littermate controls ( Figure 8A). However, the spleen and thymus weights were not different from littermate controls. Despite this, the cellularity of the thymus was still significantly reduced (58% ± 19%, P = 0.019, Figure 8B). However, as both DN and CD4+ lineages recovered to WT proportions the long-term effects of poor placentation on cell distributions in the adult thymus are less clear ( Figure 8C). While we did not evaluate the repertoire diversity or clonality, we did observe that DN subtypes were statistically distinct from WT ( Figure 8D).

DISCUSSION
We present the first single-cell analysis of the neonatal thymus, in the context of a disease model for FGR. Recently, interest in fetal immunity has led to the generation of a cell atlas for the developing thymus (48). Here, we add to this atlas with an insight into thymic development after birth. We found a number of genes in common with these early developmental snapshots, such as H19 and Plac8, but were further able to resolve multiple developing T-cell populations due to the seemingly increased thymic activity after birth.
The "fetal programming of adult disease" theory has largely been associated with non-communicable diseases such as obesity or heart disease. Within this wider context, we present the first genetic model suggesting that an adverse in utero environment perturbs immune cell development. This may lead to poor immune function later in life and consequently a higher risk of communicable diseases.
Our study has identified a novel sequela of FGR in the mouse: a lasting deficit in T cells, resulting from a maturation delay found in the growth restricted neonate. At the cellular level, there were more DN cells and a skew toward early DP cells. The high heterogeneity of DP cells found here suggests that current definitions of DP are not specific enough given the large proportion of the thymic T cells in this category. Our 26 supplemental tables from the single-cell dataset provide numerous candidates for future investigation of subcluster markers and some of these may be functionally relevant to the differentiation process. Of particular note is the high number of ribosomal protein transcripts driving clustering within this group. Ribosome biogenesis has recently been shown to play a key role in stem cell differentiation (49) and translation of specific pools of mRNAs depends on the ribosomal protein composition (43). This is particularly noticeable with defects in specific ribosomal protein genes, which can lead to discrete effects such as an erythropoiesis block in hematopoietic stem cells (50) or a loss of body hair (51), as opposed to an expected generalized protein deficit. Clusters 1, 4, and 6-comprised of DN, late DP and maturing T cells-seem to subcluster based on ribosomal protein transcript levels. The ribosomal proteins encoded by these transcripts, located largely on the outside of the ribosome, are greatly reduced in the P 0 mice. It is interesting to speculate that within this T-cell differentiation pathway, ribosome biogenesis, and hence function, plays a role in directing maturation.
The clusters were defined by transcript signatures after excluding cell cycle regulators. Nevertheless, distinct cell cycle signatures were present within the clusters, confirming that proliferation and T-cell maturation are inexorably linked, which is consistent with a similar finding in specifically CD4+ T-cell differentiation (40).
How this early phenotype leads to low thymic cellularity in the adult remains unclear. Infants with FGR typically undergo "catch-up growth" to reach normal or near normal weight. This was the case here for all parameters except thymus cellularity. T-cell maturation, is linked at a molecular level to cell cycling because RAG2 is degraded during G2, S, and M phases (52). Perhaps thymic "catch-up growth" and T cell diversification and maturation are incompatible because proliferation inhibits T-cell maturation. Indeed, the similar cell-cycle profiles of the P 0 and WT thymic cells is indicative of a lack of catch-up growth in the P 0 thymus. It is also possible that catch-up growth cannot occur because the thymus is programmed to involute after birth, rather than expand. In either case, growth restriction early in life would be irreversible: the number of T-cell clones generated in the neonatal thymus restricts the number of T-cell clones present in the adult.
In both mice and humans, thymus output declines with age, with the memory T cells dominating the peripheral pool in old age (53). The diversity and clonality of naive T cell also change with age (54). Immune function declines with age and there is an increased risk of infection and malignancy late in life (44). Genetically manipulated mice with enlarged thymi have increased thymic activity and enhanced recovery from γ irradiation (55). In contrast, patients with DiGeorge syndromes where the thymus is small display immune vulnerability early in life (45,46). Thus, a thymic deficit early in life could have dramatic effects on adult immunity which would only become exacerbated by normal thymic involution.
Overall, we have identified key immune phenotypes caused by placental insufficiency. Future studies will ultimately determine the role of FGR in both neonatal and adult immunity. Nevertheless, we have shown a distinct and negative impact of FGR on immune cells in the mouse.

DATA AVAILABILITY STATEMENT
The single-cell datasets generated for this study can be found in the ArrayExpress database at EMBL-EBI under accession number E-MTAB-6945, https://www.ebi.ac.uk/arrayexpress/ experiments/E-MTAB-6945. All other raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
FC, JK, DC-J, and WB contributed conception and design of the study. WB, JK, DH, and AK performed the research. RH and WB analyzed the data. RH contributed novel bioinformatic tools. ZY and CA contributed microfluidic equipment and expertise. WB wrote the first draft of the manuscript. RH and DC-J wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.