extfieldsofstudy: Medicine Biology

text:
Reference-based comparison of adaptive immune receptor repertoires

Summary B and T cell receptor (immune) repertoires can represent an individual’s immune history. While current repertoire analysis methods aim to discriminate between health and disease states, they are typically based on only a limited number of parameters. Here, we introduce immuneREF: a quantitative multidimensional measure of adaptive immune repertoire (and transcriptome) similarity that allows interpretation of immune repertoire variation by relying on both repertoire features and cross-referencing of simulated and experimental datasets. To quantify immune repertoire similarity landscapes across health and disease, we applied immuneREF to >2,400 datasets from individuals with varying immune states (healthy, [autoimmune] disease, and infection). We discovered, in contrast to the current paradigm, that blood-derived immune repertoires of healthy and diseased individuals are highly similar for certain immune states, suggesting that repertoire changes to immune perturbations are less pronounced than previously thought. In conclusion, immuneREF enables the population-wide study of adaptive immune response similarity across immune states.

INTRODUCTION B and T cell receptor (BCR, TCR) repertoires (also called adaptive immune receptor repertoires, AIRR) are continually shaped throughout the lifetime of an individual in response to environ-mental and pathogenic exposure. As of yet, however, there exists only a limited quantitative conception of how immune receptor repertoires differ across individuals and cell populations (Brown et al., 2019;Miho et al., 2018;Raybould et al., 2021). This is primarily because a method for measuring inter-individual MOTIVATION B and T cell repertoires record past and current immune states. Therefore, the majority of immune repertoire studies aim to measure the impact of the immune state on the immune repertoire because it is widely assumed that repertoires change as a function of the immune state. So far, a method to measure and immunologically interpret differences between immune repertoires has remained unavailable. We have addressed the methodological challenge of immune repertoire comparison by implementing a referencebased multidimensional repertoire similarity measure based on in silico and experimental immunologically interpretable ground truth.
(inter-repertoire) similarity is lacking, thus greatly impeding the understanding of how health and disease shape immune repertoires and how disease contributes to the deviation of an individual's baseline repertoire (Cobey et al., 2015). Although it is generally thought that infection or disease induces measurable repertoire changes (even on the antigen-specific agnostic level), this belief remains unproven and, in fact, is counter to current evidence finding, using statistical learning, that even in systemic infections such as cytomegalovirus (CMV) only a comparatively very small number of TCRs are infection associated (DeWitt et al., 2018;Emerson et al., 2017;Pavlovi c et al., 2021). As opposed to machine learning approaches that aim to detect the most differentiating factors (i.e., subsets of a repertoire) between, for example, two different immune states Pavlovi c et al., 2021;Pertseva et al., 2021;Shemesh et al., 2021;Widrich et al., 2020aWidrich et al., , 2020b, we investigate here a method for quantitatively comparing any two repertoires in an unsupervised fashion. We thus seek to understand to what extent individuals differ with respect to their entire repertoire and not just class-associated subsets.
The need for comparing immune repertoires using a quantitative measure has recently been addressed by approaches based on single sequence-dependent and sequence-independent features, which vary in statistical dependency (mutual information) and immunological interpretability (Chiffelle et al., 2020;Miho et al., 2018;Olson et al., 2019). Sequence-dependent approaches range from the measurement of clonal overlap (Bolen et al., 2017;Greiff et al., 2015a;Miho et al., 2018;Rognes et al., 2022;Yaari and Kleinstein, 2015) to more sophisticated algorithms that identify disease-specific enrichment of sequence clusters by testing against VDJ recombination models  or similarity networks of control datasets Shugay et al., 2015). Sequence-independent approaches are mainly represented by entropy-based diversity indices (Alon et al., 2021;Greiff et al., 2015a;Kaplinsky and Arnaout, 2016;Strauli and Hernandez, 2016), which have lately been augmented with a correction for sequence similarity (Arora et al., 2018;Vujovi c et al., 2021). None of the currently available comparative methods, which are based on single repertoire features, however, represent an integrated multi-feature measure of immune repertoire similarity that takes into account the complexity of information encoded in the ensemble of the existing immune repertoire features (Gupta et al., 2015;Nazarov et al., 2020;Shugay et al., 2015). Such an integrated measure, encoding per-feature similarity in one common mathematical structure, is needed to enable a representation of repertoire similarity.
Here, we introduce immuneREF: a measure for quantifying immune repertoire similarity across multiple immune repertoire features. Our framework, implemented in an R package, measures immune repertoire similarity using a combination of features that are immunologically interpretable (clonal expansion, sequence composition, repertoire architecture, and clonal overlap) and that cover largely distinct dimensions of the immune repertoire spaces. Specifically, to interpret immune repertoire similarity scores, immuneREF establishes a self-augmenting dictionary of simulated and experimental datasets where each new dataset analyzed may be used as a comparative reference for scoring and biologically interpreting inter-individual variation (and thus the deviation) of immune repertoire features (Figure 1). We applied immuneREF to >2,400 immune repertoires from humans with varying immune states (healthy, virus infection, autoimmune disease) and found that the similarity of blood-derived immune repertoires is not consistently a function of the immune state.
Overall, immuneREF enables the quantification of repertoire similarity at population scale while still providing single-individual resolution, and it enables answering fundamental questions such as to what extent immune repertoires are robust to perturbations introduced by immune events.

RESULTS
Reference-based comparison of immune repertoires based on immunological features: Constructing a similarity atlas of immune repertoires To derive a similarity measure for immune repertoires, we devised a framework that calculates a repertoire similarity score based on six features that reflect immune repertoire biology ( Figure 1). These features are (1) germline gene diversity (Greiff et al., 2015a;Yaari and Kleinstein, 2015), (2) clonal diversity (Greiff et al., 2015a;Stern et al., 2014), (3) clonal overlap (Greiff et al., 2015a;Yaari and Kleinstein, 2015), (4) positional amino acid frequencies (Mason et al., 2019), (5) repertoire similarity architecture (Bashford-Rogers et al., 2013;Ben-Hamo and Efroni, 2011;Miho et al., 2019), and (6) k-mer occurrence (Greiff et al., 2017b;Thomas et al., 2014) (see the STAR Methods section for a detailed immunological and mathematical description of these features). A similarity score is calculated for each pair of repertoires and each feature (six n x n symmetric matrices, n = number of repertoires), creating a similarity matrix for each feature. This matrix may be viewed as a weighted network, in which the nodes correspond to repertoires and the edges connecting the nodes are the similarity scores. The resulting six single-feature similarity networks enable insight into per-feature similarity. Finally, a composite network of the six feature similarity networks represents an interpretable multidimensional picture of the repertoire landscape. Briefly, the single features are condensed into a multi-feature composite network by taking the mean of all single-feature similarity values resulting in a single repertoire similarity value (for alternative approaches to computing composite networks, see STAR Methods section). By virtue of representing a similarity matrix as a weighted network repertoire, similarity may be computed on selected levels such as one (repertoire) to many (repertoires), many to one, and many to many ( Figure 1). Interpretability stems from all repertoire features being transformed into a similarity measure on a 0-1 scale allowing for direct quantification of their individual contribution to multidimensional immune repertoire similarity.
immuneREF measures immune repertoire similarity with high sensitivity We sought to quantify the sensitivity by which immuneREF can detect differences between immune repertoires with respect to the six repertoire features. The simulated repertoires, varying in a controlled manner, represent a ground truth reference map that enables a more precise assessment of immuneREF sensitivity. For example, simulated repertoires may be used to guide the evaluation of variation between experimental repertoires with respect to each repertoire feature as well as multi-feature combinations. Simulations were performed using the immune-SIM repertoire simulation suite , which was used to create native-like repertoires that were varied across eight parameters. Native-likeness was demonstrated in Weber et al. (2019). The parameters that were varied across simulated repertoires included clone count distribution, V-, (D-), J-gene frequency noise, insertion, and deletion likelihoods, species (human and mouse), and receptor type (IgH, TRB). We constructed additional simulated repertoires with spiked-in motifs (mimicking antigen-binding motifs; Akbar et al., 2019), excluded hub sequences in the sequence similarity network (simulating network architecture variation; Miho et al., 2019), and replaced nucleotide codons with synonymous codons (simulating biases in the k-mer occurrence that are relevant in detectable immunogenomic patterns of public clones; Greiff et al., 2017b) (see STAR Methods; Table S1 lists the parameter variations used for the simulations and how each of the parameters is expected to influence the six immuneREF features). The parameter combinations were chosen so each simulated repertoire varied only along one parameter dimension at a time, allowing us to determine the sensitivity of each feature to each parameter change.
The mathematical structure of the single-feature similarity matrices enables their merging into a composite network that provides the opportunity for a condensed single-score representation of inter-sample repertoire similarity. The composite immu-neREF network (which combines all six repertoire features) recovers major variation in the repertoires including noise The complexity of AIRRs spans the frequency, motif, and feature space to each, of which distinct repertoire features may be attributed: the immune information stored in AIRRs is multidimensional. A longstanding question in the AIRR field is how to quantitatively measure inter-sample (sample, e.g., individual, immune cell population) AIRR similarity by accounting for AIRR feature multidimensionality in the effort to understand the distribution of inter-sample AIRR similarity across different immune events or immune cell populations. (B) We set out to develop an AIRR similarity measure that is sensitive, captures maximal immune information, and is sufficiently flexible to allow future integration of additional repertoire features (extensibility). (C) Each AIRR is represented as a node in a similarity network. The edges connecting the nodes represent the similarity score between the AIRR based on the six repertoire features. The immuneREF approach establishes interpretability on different levels: (1) from a single-feature perspective, the application of spider plots allows for an interpretable comparative analysis between repertoires, enabling the user to interpret the result observed in the condensed network on a per feature basis.
(2) From the condensed feature network perspective, a major novelty introduced by the immuneREF workflow is the ability to combine established repertoire features into a common coordinate system. This transformation allows the combination of trends across features into a single condensed network that represents pairwise-cross-feature similarities. These pairwise similarities allow for the identification of subsets of more similar or aberrant repertoires. Interpretability on both features means allowing comparison to other repertoires and to simulated ones (of which we know the repertoire structure as ground truth), thus creating similarity equivalence classes. Equivalence classes create sets of reference repertoires, which enable interpreting the repertoire structures of other repertoires solely based on the immuneREF similarity score.  (Figures 2A and 2B). immu-neREF also clearly distinguishes repertoires from different receptors and species based on strongly distinguishing features such as V-, (D-), and J-gene usage while allowing the identification of commonalities in amino acid usage, clonal diversity, and architecture across immune receptors and species. This sensitivity analysis also underlines a major advantage of immuneREF, namely its flexibility to accommodate both BCR and TCR repertoires from different species in one single analysis workflow.
We quantified the sensitivity of immuneREF by detecting significant changes in similarity scores corresponding to the variation in simulation parameters across both the single feature (Figures S1-S3) and composite network ( Figure 2A) and found that each feature had a unique sensitivity profile to changes in the simulation parameters, underscoring the value of perfeature similarity evaluation. For example, a change in the alpha parameter of the Hill function (controlling clone count distribution) solely impacted the immuneREF diversity feature. As the immuneSIM parameter controlling the distribution of clone counts only affects the clone count simulation without impacting simulated sequences, the fact that only the feature targeted by the parameter change is impacted shows that immuneREF is robust to random noise in the simulation that is not introduced through parameter changes. An increase in the V(D)J noise parameter, which modifies the frequencies of the germline genes used in the simulation, led to detectable and significant changes in similarities of the germline gene usage and k-mer occurrence features. Modification of the insertion/deletion patterns (dropout of deletions and or insertions) led to a consistent impact in the amino acid frequency feature and, more importantly, the architecture feature, where a lower diversity due to restricted insertions and deletions led to significant changes in network architecture. Implanting motifs at various frequencies led to a significant similarity change in the k-mer occurrence feature. The deletion of hub sequences led to an impact in the architecture feature and also changed the repertoire overlap similarity, C  Article ll OPEN ACCESS thus underlining the importance of public clones in the network architecture as reported previously (Miho et al., 2019). Finally, we modified the repertoires by introducing synonymous codons at various percentages and found that the k-mer occurrence feature was the only one impacted. Therefore, we conclude that immuneREF features largely react as hypothesized to variation in simulation parameters (Table S1). Taken together, we demonstrated that the immuneREF framework is sensitive to even comparatively small repertoire variations.
Mutual information analysis demonstrates no interdependence to limited inter-dependence of immuneREF features While the examined features were initially chosen based on immunological criteria, we also wished to verify whether each feature provides a sufficiently different measurement of the immune repertoire information space ( Figure 2C). Specifically, having integrated all features into a common coordinate system, we were able to compute cross-feature mutual information and found that features show no dependence to limited dependence (range = 0.01-0.57; Figure 2C) indicating largely non-overlapping and distinct spaces of immune information captured. The highest mutual information was found between the positional and sequential sequence-derived features (i.e., positional amino acid frequency and gapped k-mer occurrence, respectively), whereas the lowest mutual information value was found between the diversity and convergence features ( Figure 2C).
Complementarily, we sought to quantify to what extent the addition of new repertoire features leads to diminishing returns (sufficiency analysis). To this end, we computed the mean change in repertoire similarity values when increasing the number of features from one through six. Thereby, we could show that each additional feature added increasingly less information, as shown by the diminishing change of the mean similarity value with each added feature. The saturation of the mean similarity change curve indicated information saturation independent of the order in which features were arranged (Figures 2D and S3G-S3J). As discussed below, mutual information values behaved similarly for experimental repertoire data. Thus, we demonstrated that the immuneREF framework creates information-laden similarity networks, whose topologies capture the immunological similarity landscape of immune repertoires.
The similarity landscape of simulated repertoires defines reference repertoires By calculating the similarity matrix for each of the six immune repertoire features, we embedded the six different immunological features into a common coordinate system, i.e., a network structure. This network (with nodes representing repertoires and weighted edges representing pairwise similarity) situates each repertoire within a similarity landscape allowing quantification of many-to-many repertoire similarity.
A more fine-grained image of the similarity landscape may be gained by examining the similarity from the perspective of every single repertoire ( Figures 3C and 3D). We define the local similarity of a repertoire to its neighboring repertoires as a scaled node strength (see STAR Methods). This local similarity represents the position of the repertoire with respect to its direct neighbors in its cohort (defined by an application-dependent label, e.g., same species and disease) and allows us to distinguish between well embedded and aberrant repertoires. The local similarity measure further acts as a magnifying glass by elucidating finer differences between repertoires, which are diluted by population averages when examining repertoire similarity across the full similarity network. Using this perspective, repertoires that are most (locally) similar to other repertoires in their cohort can be identified, allowing the extraction of repertoires most representative for a given immune state. Such detailed one-to-one feature comparisons highlight, in the most simple case, which features of the simulated repertoires are receptor specific (amino acid frequency, k-mer occurrence, VDJ usage, and convergence) and which are more general to immune repertoire data showing higher similarity across different species and receptors (diversity, architecture) ( Figures 3E and S4).
Having evaluated the similarity of simulated datasets, these may serve as a reference to interpreting similarity score variation of experimental repertoires ( Figure 3C), thus enabling the creation of equivalence classes of immune repertoires not only as previously performed based on clonal expansion (Greiff et al., 2015b) but based on six repertoire features. Furthermore, any evaluated repertoire, be it of experimental or simulation origin, will become a new node in the similarity network and may serve as a valid reference point (just as any other node in the network). This network of self-augmenting repertoire similarity reference points is another source of interpretability as it allows the linking of the repertoire similarity of any number of repertoires with their underlying features. In the next section, we provide such a repertoire similarity network on experimental datasets.
Validation of immuneREF on experimental data: Detection of differences between cell populations in mouse immunization and human COVID-19 datasets To validate immuneREF sensitivity on experimental data, we used antibody repertoire datasets generated from a mouse antigen immunization study, where differences in the similarity between antigen immunization cohorts are expected (Greiff et al., 2017a(Greiff et al., , 2017bMiho et al., 2019). Notably, we were able to recover clear differences between isotypes and cell populations (both with higher within-cohort and lower across-cohort similarity); additionally, we found that the antigen immunization cohorts have more distinct similarity profiles in the plasma cell populations (IgG) compared with the antigen-inexperienced cell populations ( Figure S5). The overall high similarity scores across the full immunological feature range are in agreement with our previous studies where we observed high similarity between these repertoires on a single feature basis (Greiff et al., 2017a).
Similarly, applying immuneREF to TCR repertoires of patients recovered from mild cases of COVID-19  revealed clusters of increased similarity within patients and cell populations ( Figure S5).
Application of immuneREF to >1,500 experimental blood immune repertoires indicates only small similarity-based differences between health and autoimmune disease Having established the sensitivity of our approach in detecting a wide range of differences between simulated repertoires  Figure S5) with respect to immunologically relevant and interpretable repertoire features, we set out to determine the similarity landscape of large-scale experimental TCR repertoire datasets. We evaluated 1,522 human TCR repertoires derived from peripheral blood mononuclear cells (PBMCs) of patients with varying and diverse immune states (PanImmune Repertoire Database (PIRD) dataset containing samples from healthy, rheumatoid arthritis (RA), and systemic lupus erythematosus (SLE) patients; Table S2). We found an even similarity landscape of overall high similarity scores ( Figure 4A). Similarity score distribution was also even in single features, which despite feature-specific differences, show overall high similarity scores between repertoires. We examined networks at three different similarity cutoffs (an edge is drawn between two repertoire nodes if their similarity is in 25%, 50%, and 75% top weights, respectively), and we found that in all three cases, no immune state-specific grouping could be observed ( Figure 4B).
The range of general and local similarities across all samples as well as within each disease cohort was evaluated using an analogous approach to that used for the simulated datasets ( Figures 4C and 4D). While the similarity scores ranged between $0.5 and 0.8 overall, the within-disease cohort spread varied, with the healthy and RA cohorts showing a more restricted range of similarity scores compared with a broader range for SLE ( Figures 4C and 4D).
To quantify per feature similarity and dissimilarity with respect to a reference dataset, we compared the repertoires identified as the ones best connected (highest local similarity) within their cohort to an immuneSIM reference repertoire (human, TRB, standard parameters; see STAR Methods) ( Figure 4E). The similarity scores of all tested immune states largely overlap with respect to the healthy reference repertoire, with convergence being the feature dimension with the largest dissimilarity, meaning there is almost no convergence between the RA or SLE samples and the reference.
Following our observations of high repertoire similarity within the PIRD dataset, we ran immuneREF on another large publicly available dataset (human, TCR) (Emerson et al., 2017)  Similarity score D Figure 3. The similarity landscape of simulated repertoires defines reference repertoires (A) Baseline similarity between replicates for repertoires simulated using default immuneSIM parameters (see Table S1) is R0.96 for five of six features, with the convergence feature being the exception by definition at %0.09. Bar graphs show mean SEM across replicates. (B) Repertoire similarity distribution in a condensed network across the various evaluated parameter range. Across cohorts, similarity scores have a broad range, whereas within cohorts the range is more restricted.
(C) Workflow to determine representative repertoires per cohort going from many-to-many to a one-to-one comparison.
(D) Local similarity distribution per species/receptor combination enables situating each repertoire based on its connectivity with respect to neighbors in the same cohort.
(E) Comparing repertoires with maximal local similarity in their cohort visualizes the commonalities between receptor types; here the Murine IgH repertoire with maximal local similarity serves as a reference repertoire. The plot visualizes the similarities of each non-reference repertoire to the Murine IgH reference. tients, 351 are from CMV-negative patients, and 26 are from patients with unknown CMV status. This dataset has previously been used to showcase immune state classification with high accuracy via the identification of CMV-associated public TCR sequences (sequences shared between individuals). In a similar fashion, immune state-associated public sequences were used to successfully classify RA and SLE samples from the PIRD dataset (Liu et al., 2019). As with the PIRD dataset, we observed high within and across immune state repertoire similarity ( Figure S5). This is in line with the findings of Emerson and colleagues as they found that only a small subset of clones (CMV-associated ones in Emerson et al., 2017) significantly differed in abundance between immune states (CMV+, CMV-) and that that shared antigen exposure to CMV led to a reduced number of shared TCRb clones, even after controlling for individual human leukocyte antigen (HLA) type, indicating a largely private response to a major viral antigenic exposure (Johnson et al., 2021).
In summary, the results of our analysis of human TCR repertoires strongly support the argument that the signal-to-noise ratios, where signal means repertoire features associated with disease status, are unfavorably tilted toward noise, where noise is defined as technological and immunological information, which cannot yet be linked to a given disease state.

Extensibility of immuneREF: Integration of gene expression with immune repertoire data
The mathematical structure of the composite network obtained from immuneREF allows the extensibility of the immuneREF framework to other features. As proof of principle of this immu-neREF capability, we show here an integrative analysis of immune repertoires and gene expression. This integration is of high interest to RNA-seq experiments that include both receptor and global transcript sequences, or even repertoire experiments paired with transcriptomics (Rubio et al., 2022;Song et al., 2021). Integration of immune repertoire with gene expression is challenging due to the multidimensional nature of both kinds of datasets and the discrepancy in their data structure. Previous attempts of integration are still over-simplistic, such as the calculation of correlation between the number of distinct CDR3 amino acid sequences and gene expression of some marker genes such as CD3, CD4, CD8, HLA class I, and class II genes (Brown et al., 2015).
immuneREF includes the option to evaluate similarity based on a gene expression matrix and add it to the composite   Table S1).
Cell Reports Methods 2, 100269, August 22, 2022 7 Article ll OPEN ACCESS network. Briefly, immuneREF first filters all genes with low variation between experimental conditions and then calculates the pairwise correlation between observations to construct a single gene expression feature (similarity matrix). Once the seven features (six from immune repertoires and one for gene expression) are calculated, they may be condensed into a multi-feature network as described above. Our solution for integrating receptors with gene expression confers immuneREF the advantage of overlaying dual biological information ( Figure S6A).
As an example, we analyzed bulk RNA-seq gene expression of pre-B cell line B3 from the published STATegra project (Gomez-Cabrero et al., 2019). This is a time-course experiment that collects samples at six time points using an inducible Ikaros system where B cell progenitors undergo growth arrest and differentiation ( Figure S6B). Principal-component analysis (PCA) showed clear differences at gene expression level when control and Ikaros groups were compared but also within the Ikaros group across time, being t0 the nearest to controls ( Figure S6B). To generate the single-feature similarity matrix of gene expression that better collects these differences, we tested the three available correlation-based methods implemented in immuneREF ( Figures S6C-S6E). All of them perfectly separated control (blue) and Ikaros (red) groups. Additionally, ''Pearson correlation'' and ''PCA scores'' nearly recovered correctly the time series pattern (purple to yellow degradation), while mutual rank matched perfectly.

DISCUSSION
Combining methods from both immune repertoire and network analysis, we have provided a framework for flexible referencebased quantification of immune repertoire similarity. Using ground truth simulations, we show that immuneREF is sensitive to inter-repertoire differences in all immunological features. Taking advantage of information theory, we showed using both simulated and experimental data that the features selected for immuneREF cover a large extent of immune repertoire biology. We introduced the concepts of full-network repertoire similarity and local similarity, which allow complementary quantification of the impact of the differences in the repertoire similarity landscape. Specifically, while the more general repertoire similarity evaluated on the entire network provides insight into the range of similarity within and across conditions, local similarity shows a particular advantage of the network approach, as the embedding of a repertoire in its neighborhood can markedly differ from what can be expected by its pairwise connections.
immuneREF not only provides a framework for measuring immune repertoire similarity but also for interpreting it. Specifically, it enables the creation of equivalence classes of immune repertoires lacking from existing methods. For example, once the similarity observed within a given set of experimentally obtained immune repertoires has been computed, such repertoires may function as reference points that in turn enable the interpretation of relative similarity in other repertoires (Figures 1 and 3C). Of note, the concept of diversity measures creating equivalence classes has been noted previously for Hill diversity measures (Greiff et al., 2015b) and is here extended to include additional repertoire features immuneREF unifies as single and composite features, frequency-dependent, and sequence-dependent similarity measures into one computational framework. Beyond quantifying the repertoire similarity of experimental immune repertoires, immuneREF also enables the comparison of simulated Marcou et al., 2017;Safonova et al., 2015;Weber et al., 2020;Marcou et al., 2017;Safonova et al., 2015;Weber et al., 2019) and in vitro synthetic immune repertoires used for therapeutic antibody discovery (Mason et al., 2018). Furthermore, immuneREF may be used for data curation purposes in immune repertoire databases such as iReceptor (Corrie et al., 2018), VDJserver (Cowell et al., 2015, PIRD (Zhang et al., 2019), and Observed Antibody Space (Kovaltsuk et al., 2018). Specifically, upon the integration of an immune repertoire into a database, the similarity of the repertoire with all other stored repertoires may be computed. Beyond immunological insight, immuneREF may reveal unexpected technological variation, thus motivating follow-up inspection (Barennes et al., 2021). Since immuneREF has been built to work across species, cell populations, and receptor types and experimental or simulated data (all-in-one comparative framework), it enables rapid distinction of cohort-specific and cohort-unspecific features. This is also important for comparative immunological approaches not centered on health versus disease comparison but, for example, the evolution of adaptive immunity (Pancer and Cooper, 2006).
The ease of use of the immuneREF approach opens new possibilities for large-scale comparative studies as shown on the PIRD dataset, which may yield additional insight into the challenges of predicting immune state based on repertoire profiling. Indeed, we found that the population average quantified by im-muneREF may ''conceal'' relevant immunological phenotype signals, despite the fact that the sensitivity of immuneREF was shown to be high in simulated and experimental data (Figures 2  and S5). Given the lack of large-scale (antigen-specific) data, it remains unclear how the information of the immune state is distributed across immunological features. Specifically, our finding-that repertoire similarity does not differ across immune states-is strictly only valid for unsorted PBMC TCR repertoire data as examined in this study. As known from previous studies (Amoriello et al., 2020(Amoriello et al., , 2021Csepregi et al., 2021;Ghraichy et al., 2021;Greiff et al., 2017b;Li et al., 2020;Ota et al., 2022;Riedel et al., 2020;Rosati et al., 2021), different cell populations (in different lymphoid organs) may behave in a highly different manner ( Figure S5). On the other hand, it did not escape our attention that this broad similarity in human blood samples might suggest the maintenance of lymphocyte homeostasis even in the event of chronic disease.
Our results reinforce the notion that while some diseases may introduce abnormalities into the immune repertoire, others result in a comparatively normal one (Bashford-Rogers et al., 2019), a result that suggests the absence of a signature unique to health. If this is true, then blood-based immune repertoire diagnostics will require even more advanced methodologies to be developed (Arnaout et al., 2021;Dahal-Koirala et al., 2022;Widrich et al., 2020a). For example, for simulated repertoires, motif implants in R10% of sequences were required to affect the amino acid frequency and architecture features, suggesting that even in the case of high clonal expansion, the impact on the repertoire might not be sufficient to significantly change major repertoire features. This is reinforced by results showing that the diseasedriving response in multiple autoimmune diseases is only to a small part antigen specific (Christophersen et al., 2019;Dahal-Koirala et al., 2022). More generally, our paper advances the state of the art of the immune repertoire field by changing the null hypothesis. Specifically, currently, the predominant thinking is that any immune state changes measurably the immune repertoire in a systematic fashion. Our paper challenges this view by finding that, a priori, we should not expect to see differences (Figures 4 and S10), and any substantial change must be proven. This change of perspective is highly valuable to the field as it pushes it toward more sensitive and robust approaches to immune repertoire and machine learning analysis (Arnaout et al., 2021;Kanduri et al., 2021;Slabodkin et al., 2021). Specifically, the usefulness of global features for diagnostics is severely limited, and to detect single-sequence-level differences (Emerson et al., 2017;Kanduri et al., 2021;Widrich et al., 2020a), single-sequence-level statistical and machine learning approaches are needed Schattgen et al., 2021).
In the future, ultra-deep (Briney et al., 2019;Soto et al., 2019Soto et al., , 2020 and population-wide, large-scale immune repertoire projects such as Human Vaccines Project (Crowe and Koff, 2015) may benefit from using immuneREF for identifying immune event-driven aberrations from a baseline repertoire similarity. Furthermore, large-scale database initiatives such as the iReceptor gateway (Corrie et al., 2018) may benefit from immu-neREF functionality for on-the-fly computation of inter-dataset similarity.
Limitations of the study Although we consider the usefulness of the six chosen features to be established (Figures 2 and S3G-S3J), we concede that the asymptotic nature of the sufficiency calculation leaves the door open to the introduction of additional features. The proposed set of immuneREF features denotes in this sense a minimally sufficient set for the analysis of immune repertoire datasets. It ensures sufficient coverage of the major variationintroducing aspects. It is for that reason we devised immu-neREF as inherently modular, allowing single-and multi-feature analysis as well as encouraging the addition of new features relevant for particular problems such as transcriptome analysis (Schneider-Hohendorf et al., 2018; Figure S6), HLA typing for TCR studies (DeWitt et al., 2018;Emerson et al., 2017;Francis et al., 2021), single-cell omics information Setliff et al., 2019;Sturm et al., 2020;Yermanos et al., 2021), gene-specific substitution profiles for somatic hypermutation analysis (Sheng et al., 2017), lineage-specific information (Hoehn et al., 2016(Hoehn et al., , 2021, and antigen-specific and antigenassociated motifs identified by sequence clustering and machine learning (Akbar et al., 2019;Dash et al., 2017;Friedensohn et al., 2020;Glanville et al., 2017;Greiff et al., 2017b;Horst et al., 2021;Mason et al., 2019;Mayer-Blackwell et al., 2021;Meysman et al., 2018;Quiniou et al., 2020;Sidhom et al., 2019;Wong et al., 2020;Yohannes et al., 2021). In particular, a future extension of immuneREF may be a feature that reliably identifies antigen-specific sequences, thus increasing the amount of immune information recovered. More generally, adult repertoires are very complex and contain hidden informa-tion of many antigens at different time points that might have been shared by different individuals. For instance, repertoire fingerprints of influenza infection might be present on most studied individuals and could explain the difficulty to distinguish healthy and diseased individuals. New features including (single-cell-based) antigen specificity patterns may help separate shared infection marks on the immune repertoire.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:  Evenness profiles were calculated as described previously (Greiff et al., 2015b) on the CDR3 nucleotide level. Briefly, we calculated the Hill-diversity for alpha values 0-10 in steps of 0.1 with alpha = 1 being defined as the Shannon evenness. Each entry in the profile varies between z0 and 1, where higher values indicate an increasingly uniform clonal frequency distribution. We determined evenness profiles for each repertoire and evaluated cross-repertoire evenness similarity by Pearson correlation of the repertoires' evenness profiles as described previously (Amoriello et al., 2020;Greiff et al., 2015bGreiff et al., , 2017a. immuneREF feature: Positional amino acid frequencies The positional amino acid frequencies were calculated separately for each CDR3 sequence length. To decrease bias by extraordinarily short or long CDR3 sequences, we limited this analysis to a range of the most common lengths (8-20 amino acids) (Greiff et al., 2017a;Raybould et al., 2019). Briefly, per position amino acid frequencies were calculated for each length. Subsequently, the resulting per length frequency vectors of each repertoire were Pearson-correlated by length and the mean correlation was calculated. Unlike in the case of k-mer occurrences, no positions are excluded, making AA frequency more sensitive to VDJ usage perturbations. Relative frequencies were used for all positional amino acid frequency calculations.
immuneREF feature: Sequence similarity network architecture As previously described (Miho et al., 2019), we constructed a sequence similarity network for each immune repertoire: nodes represent amino acid CDR3 sequences connected by similarity edges if they had a Levenshtein Distance of 1 (LD = 1). The igraph R package was used to calculate networks (v.1.2.4.1, Csardi and Nepusz, 2006), that were analyzed with respect to four measures representing different aspects of network architecture: (i) cumulative degree distribution, (ii) mean hub score (Kleinberg hub centrality score), (iii) fraction of unconnected clusters and nodes and (iv) percent of sequences in the largest connected component. An LD = 1 network was constructed for each repertoire and the similarity between the repertoires' resulting network was evaluated with respect to their differences in the cumulative degree distribution, mean hub-score, outlier sequence occurrence, and largest network components; these metrics have been shown to be defining repertoire characteristics that are robust to subsampling (Miho et al., 2019). The similarity of the architecture between two repertoires A and B was calculated as the mean of four components: (i) the cumulative degree distribution (Pearson correlation between repertoires), (ii) mean hub scores ð1 À jMeanHubScore A À jMeanHubScore B Þ, (iii) the fraction of unconnected components, and (iv) the fraction of sequences in the largest component ð1 À jPercLargestComponents A À PercLargestComponents B jÞ. Unlike many of the other features, the network feature combines multiple single measures, which rendered it difficult to perform Pearson correlation analysis involving all four investigated network measures. Therefore, we adopted the network feature comparison approach described above.
immuneREF feature: Repertoire overlap (convergence) The pairwise repertoire clonal overlap (clones defined based on 100% similarity of CDR3 amino acid sequence), was calculated across repertoires, as previously described (Greiff et al., 2017b): This clonal sequence overlap measure represents the similarity value between repertoires with respect to clonal convergence.
immuneREF feature: Germline gene diversity The relative frequency of germline genes (defined by the ImmunoGenetics Database, IMGT) (Giudicelli et al., 2004) across clones in each repertoire was calculated for each repertoire depending on species and immune receptor class (Ig, TR). The germline gene usage allows insight into deviations from a baseline recombinational likelihood and thereby captures the potential impact of disease, vaccine, or other events on the immune state (Avnir et al., 2016;Greiff et al., 2017a). To determine germline gene usage similarities, we examined the V-and J-gene frequencies across clones for each individual. The Pearson correlation coefficient was determined for each of the frequency vectors (V-, D-, J-gene) with entries of all IMGT variants in a pairwise fashion between samples as described previously (Greiff et al., 2017a;Weber et al., 2019). Specifically, the correlations are calculated per germline gene, leading to separate V_cor, D_cor, J_cor values (and additionally VJ_cor for each V_J combination). The resulting correlation values are combined into a single value by calculating a weighted mean of these components. The weight vector used for the results in the manuscript is c(V = 1,D = 1,J = 1,VJ = 0).
immuneREF feature: Gapped k-mer occurrence For a given k-mer size k and maximal gap length m, the nucleotide-based gapped-pair-k-mer occurrences were counted for all gap sizes % m (Palme et al., 2015). The parameters k and m were chosen based on previous research (Greiff et al., 2017b), where defining parameters k = 3, m % 3 was shown to lead to an encoding sufficient for sequence classification. The counts were normalized by the total number of gapped k-mers found across all gap sizes such that short-gap gapped-k-mers were weighted higher than larger gap sizes. While the amino acid frequency distribution contains positional information, the gapped k-mer occurrence represents shortand long-range sequential information encoded in the repertoire. We counted the occurrence of gapped k-mers (k = 3, m % 3) across Article ll OPEN ACCESS all CDR3 sequences of a repertoire and correlated the resulting distributions between repertoire pairs using Pearson correlation as described previously .
immuneREF feature: Transcriptome integration In order to keep the most informative genes from the genes obtained in a transcriptome experiment, immuneREF firstly applies a low variation filter (Hackstadt and Hess, 2009). Specifically, the standard deviation (SD) is calculated per gene across samples, and all genes above a certain threshold (default, SD > 1) are preserved for subsequent analysis.
To construct the gene expression feature similarity matrix, the Pearson correlation was calculated between samples. Additional approaches for the calculation of the gene expression feature similarity matrix implemented in the immuneREF package (mutual rank, PCA) are described in the package documentation.
Calculating repertoire similarities per feature The calculation of the similarity values between a pair of repertoires was performed in a feature-specific manner as described in the methods section of each feature.

Repertoire similarity -Condensing features into a composite network
The single features are condensed into a multi-feature network by taking the mean of all single-feature similarity values resulting in a single repertoire similarity value. The resulting condensed network represents a weighted composite of the single-feature similarity networks. Additional approaches to obtain a composite network (max similarity, min similarity, SNF (Wang et al., 2014)) are implemented in the R-package as described in the package documentation.

Mutual information
Mutual information is a measure that quantifies to what extent one random variable explains another. Mutual information was defined as IðX; YÞ = HðXÞ À HðXjYÞ = X x;y P XY ðx; yÞlog P XY ðx; yÞ P x ðxÞP y ðyÞ Where, H(X) is the marginal entropy, H(X|Y) the conditional entropy P XY the joint probability distribution of X and Y and P X and P Y the respective marginals. Mutual information was calculated using the R packages entropy (v.1.2.1, Hausser, 2014) and infotheo (v.1.2.0, Meyer, 2014). The values were normalized to the range [0,1] by dividing the mutual information by the sum of the entropies H(X)+H(Y). This normalized mutual information, also known as redundancy, is zero when both are independent and maximal when knowledge of one of the variables becomes redundant given the other.

Quantification of mutual information across ensembles of repertoire features
The mutual information between two features was calculated across all values in the similarity matrix, whereas the similarity matrix represents all pairwise similarity values between repertoires for a given feature. For the V(D)J diversity feature, values were set to zero by definition (i.e., the similarity between repertoires of different species/receptors) and were excluded from this calculation. We ensembled immune information captured by the repertoire features ( Figure 2) as the extent to which repertoire features collectively cover immune repertoire complexity. Specifically, we evaluated the change in mutual information between subsequently added features. Features were added one by one (1-feature network / 2-feature network, 2-feature network / 3-feature network, and so forth, where n-feature means n features combined into a composite network), with the next feature to be chosen randomly (500 permutations of feature combinations per ''n-features / n+1-feature'' step).

Local repertoire similarity
To determine a single value measure for how connected a repertoire is within a subgraph (e.g. the repertoires of healthy human IgH repertoires and the similarity values between them), we defined the local similarity measure. It is calculated by dividing the node strength of each repertoire within a subgraph (sum of all edge weights connecting it to the other nodes in the subgraph) by the sum of all node strengths in the subgraph.
Local similarity gives the ratio of node strength that is connected to each repertoire in a subgraph and thus allows the identification of the most and least representative node of any category (the one most and least strongly connected within that category, respectively, see Figure 2C). The local similarity is dependent on the number of nodes within the subgraph and is therefore only used to compare repertoires within the subgraph. To enable comparison of local similarity values across different subgraphs, local similarity can be scaled by dividing by the number of nodes in the subgraph to correct for varying subgraph sizes in cases where the number of repertoires per subgraph differs.


The genome sequence of the square-spot rustic, Xestia xanthographa (Schiffermuller, 1775)

We present a genome assembly from an individual male Xestia xanthographa (the square-spot rustic; Arthropoda; Insecta; Lepidoptera; Noctuidae). The genome sequence is 934 megabases in span. The majority of the assembly (99.94%) is scaffolded into 31 chromosomal pseudomolecules, with the Z sex chromosome assembled.

Background
Xestia xanthographa (square-spot rustic) is a widespread noctuid moth found across much of the Palearctic, Europe, North Africa and North America; its larvae are nocturnal feeders on various grasses. In the UK, adults are abundant in late summer from August to September and the species overwinters as a larva. Xestia xanthographa was a key species in a recent study revealing the detrimental effects of street-lighting on caterpillar abundance in the UK (Boyes et al., 2021). The species has also been recorded as a common prey species for autumn-flying bats (Razgour et al., 2011), and, as an adaptation to facilitate bat avoidance, the auditory sensitivity of X. xanthographa is broadly tuned with an optimal frequency of 30 kHz (Norman & Jones, 2008).
The genome of X. xanthographa, was sequenced as part of the Darwin Tree of Life Project, a collaborative effort to sequence all of the named eukaryotic species in the Atlantic Archipelago of Britain and Ireland. Here we present a chromosomally complete genome sequence for X. xanthographa, based on one male specimen from Wytham Woods, Oxfordshire, UK.

Genome sequence report
The genome was sequenced from a single male X. xanthographa collected from Wytham Woods, Oxfordshire, UK (latitude 51.772, longitude -1.337) (Figure 1). A total of 27-fold coverage in Pacific Biosciences single-molecule long reads (N50 12 kb) and 40-fold coverage in 10X Genomics read clouds were generated. Primary assembly contigs were scaffolded with chromosome conformation Hi-C data. Manual assembly curation corrected 86 missing/misjoins and removed 17 haplotypic duplications, reducing the assembly size by 1.41% and scaffold number by 44.86%, and increasing the scaffold N50 by 1.61%.
The final assembly has a total length of 934 Mb in 59 sequence scaffolds with a scaffold N50 of 31 Mb (Table 1). Of the assembly sequence, 99.94% was assigned to 31 chromosomallevel scaffolds, representing 30 autosomes (numbered by sequence length), and the Z sex chromosome (Figure 2- Figure 5; Table 2). The assembly has a BUSCO (Simão et al., 2015) completeness of 98.7% using the lepidoptera_odb10 reference set. While not fully phased, the assembly deposited is of one haplotype. Contigs corresponding to the second haplotype have also been deposited.

Sample acquisition and nucleic acid extraction
A male X. xanthographa (ilXesXant1) and a second specimen of unknown sex (ilXesXant2) were collected from Wytham Woods, Oxfordshire, UK (latitude 51.772, longitude -1.337) by Douglas Boyes, University of Oxford, using a light trap. The specimens were identified by the same individual and snapfrozen on dry ice. DNA was extracted from whole organism tissue of ilXesXant1 at the Wellcome Sanger Institute (WSI) Scientific Operations core from the whole organism using the Qiagen Figure 1. Image of the Xestia xanthographa specimens taken prior to preservation and processing. Above, ilXesXant, used for genome and Hi-C sequencing; below, ilXesXant2, used for RNA-Seq. and Illumina HiSeq 4000 (RNA-Seq) instruments. Hi-C data were generated from abdomen tissue of ilXesXant1 using the Arima v1.0 kit and sequenced on HiSeq X.

Genome assembly
Assembly was carried out with Hifiasm (Cheng et al., 2021). Haplotypic duplication was identified and removed with purge_dups (Guan et al., 2020). One round of polishing was performed by aligning 10X Genomics read data to the assembly with longranger align, calling variants with freebayes  and Pretext. The mitochondrial genome was assembled using MitoHiFi (Uliano-Silva et al., 2021), which performed annotation using MitoFinder (Allio et al., 2020). The genome was analysed and BUSCO scores generated within the BlobToolKit environment (Challis et al., 2020). Table 3 contains a list of all software tool versions used, where appropriate.

Xin Liu
State Key Laboratory of Agricultural Genomics, BGI (Beijing Genomics Institute)-Shenzhen, Shenzhen, China The data note by Douglas et al. described the dataset of the assembled square-spot rustic genome. The manuscript clearly presented the methods used for sampling, sample preparation, sequencing, and data analysis, with the related statistics presented. The genome assembly is of good quality, with long scaffold N50 and anchored to chromosomes, making it valuable for later genome and other studies. I think the current manuscript is good enough to be published as a data note. Meanwhile, I have the following suggestions for the authors to consider: In the 'Genome sequence report' section on Page 5, the statement of BUSCO results might need to be modified. Interpretation of the BUSCO result should be like "98.7% of the BUSCOs in the lepidoptera_odb10 dataset were found to be complete", as you can also observe the other 1% to be either in multiple copies (duplicated), or fragmented. This reflected the completeness of the genome assembly.

1.
Also, in the part mentioned above, I would suggest providing some statistics (at least the length maybe) of the assembled second haplotype, which would be informative for understanding the heterozygosity in the genome.

2.
I found all the data notes in the same gateway directly used figures from the BlobToolkit Viewer. I would suggest revising the figures in the data note thus they would be a better fit. For example, the axis labels might be better to be revised in cases and the font sizes and the dataset label can also be removed. I don't know what the 'Scale' means in Figure 2, and I think Figure 3 and Figure 4 are quite difficult to understand, at least for me.

3.
Last, I would like to express my condolence on Douglas Boyes' death. After noting this from the manuscript and searched online, I would like to thank his contributions in providing this valuable dataset.

Is the rationale for creating the dataset(s) clearly described? Yes
Are the protocols appropriate and is the work technically sound? Yes

Are sufficient details of methods and materials provided to allow replication by others? Yes
Are the datasets clearly presented in a useable and accessible format? Yes


Higher Daily Air Temperature Is Associated with Shorter Leukocyte Telomere Length: KORA F3 and KORA F4

Higher air temperature is associated with increased age-related morbidity and mortality. To date, short-term effects of air temperature on leukocyte telomere length have not been investigated in an adult population. We aimed to examine the short-term associations between air temperature and leukocyte telomere length in an adult population-based setting, including two independent cohorts. This population-based study involved 5864 participants from the KORA F3 (2004–2005) and F4 (2006–2008) cohort studies conducted in Augsburg, Germany. Leukocyte telomere length was assessed by a quantitative PCR-based method. We estimated air temperature at each participant′s residential address through a highly resolved spatiotemporal model. We conducted cohort-specific generalized additive models to explore the short-term effects of air temperature on leukocyte telomere length at lags 0–1, 2–6, 0–6, and 0–13 days separately and pooled the estimates by fixed-effects meta-analysis. Our study found that between individuals, an interquartile range (IQR) increase in daily air temperature was associated with shorter leukocyte telomere length at lags 0–1, 2–6, 0–6, and 0–13 days (%change: −2.96 [−4.46; −1.43], −2.79 [−4.49; −1.07], −4.18 [−6.08; −2.25], and −6.69 [−9.04; −4.27], respectively). This meta-analysis of two cohort studies showed that between individuals, higher daily air temperature was associated with shorter leukocyte telomere length.

INTRODUCTION
Climate change is a major public health concern that is becoming increasingly important worldwide. The global surface temperature for July 2021 was the highest for July in the 142-year record of the National Centers for Environmental Information of National Oceanic and Atmospheric Administration (NOAA), which dates back to 1880. 1 The climate is warming quickly in the World Health Organization European Region, which is experiencing accelerated rates of temperature increase and increased frequency and intensity of heat waves. 2 Previous studies have shown that the frequency and severity of extreme weather events are increasing as a consequence of climate change. 3 Prominent examples are the heat waves in America and Canada and the floods in Germany and China in 2021. Importantly, accumulating evidence shows that increases in air temperature endanger human health and well-being in numerous ways and are associated with increased morbidity and mortality, especially in age-related vulnerable population subgroups, such as the elderly. 4−10 A recent worldwide study based on data from 750 locations in 43 countries showed that the global heat-related excess death ratio increased by 0.21% between 2000−2003 and 2016−2019. 7 The global population is aging rapidly and the prevalence of age-related diseases is quickly increasing. 11 Telomeres are highly conserved tandem repetitive nucleotide sequences (TTAGGG), which provide a protective cap at the ends of chromosome to maintain genome stability. 12,13 Telomeres are fundamental for cell division and shorten after each round of cell division. 14 Consequently, leukocyte telomere length shortening is evaluated as a potential biomarker for biologic aging, and telomere shortening has been associated with increasing numbers of age-related diseases such as stroke, cardiovascular disease, or cancer. 15−18 A recent large cohort study from the United Kingdom Biobank found that shortened leukocyte telomere length was associated with increased overall cardiovascular, respiratory, digestive, musculoskeletal, and COVID-19 mortality. 19 A further review study showed that telomere attrition was influenced by genetic and environmental factors. 12 Evidence is rapidly growing that telomere shortening can be accelerated by exposure to nonoptimal environmental factors like air pollution, 20,21 and pesticides. 22 Air temperature is an important environmental factor, especially in the context of climate change. A study in the three largest English cities (Greater London, Greater Manchester, and West Midlands) found that heat and cold exposure increased the risk of mortality and years of life lost. 23 In addition, high air temperature or heat stress has been associated with higher levels of oxidative stress and inflammatory biomarkers. 24−27 As oxidative stress and inflammation can speed up telomere attrition, 28,29 these findings suggest that air temperature might affect telomere length. However, only one birth study has reported that prenatal high and low air temperature exposures were associated with shorter cord blood telomere length. 30 So far, the effect of air temperature on leukocyte telomere length has not been investigated among an adult population.
Therefore, we aimed to examine the short-term associations between air temperature and leukocyte telomere length in the region of Augsburg, Germany, within two independent adult cohorts. The study area is located in the city of Augsburg and two surrounding counties in Southern Germany. 31 The study design, standardized sampling method, and data collection have been described in detail elsewhere. 31−33 Of 2974 and 3080 individuals who participated in KORA F3 and KORA F4, the current analysis included 2865 F3 and 2999 F4 participants, with no missing information on leukocyte telomere length, air temperature, or covariates of interest.

MATERIALS AND
The study was approved by the ethics board of the Bavarian Chamber of Physicians (Munich, Germany) in adherence with the declaration of Helsinki. All participants gave written informed consent.
2.2. Exposure Assessment. Countrywide high-resolution (1 km × 1 km) minimum (Tmin), mean (Tmean), and maximum (Tmax) daily air temperature data were estimated using hybrid spatiotemporal regression-based models, following the approach from Kloog et al. 34,35 Several data sets from multiple sources were incorporated in the modeling process, including air temperature measurements from weather stations as well as satellite-derived land surface temperature (LST), elevation, vegetation, urban fabric, arable land, pastures, forests, and inland waters, which were harmonized into 366,536 grid cells of 1 km × 1 km based on the European INSPIRE (Infrastructure for Spatial Information in the European Community) standard using the Lambert Azimuthal Equal-Area projection, EPSG: 3035 (GeoBasis-DE/BKG (2021)) before modeling. We trained three-stage models to achieve air temperature predictions with complete temporal and spatial coverage. In the first stage, a linear mixed effects model with daily random intercepts and slopes for LST and adjusted for the aforementioned spatial predictors was trained in the combinations of days and grid cells, where air temperature measurements and satellite LST were both available. In the second stage, the first stage model was used to predict air temperature for grid cells without air temperature measurements but with available LST data. For the remaining days and grid cells with neither LST data nor air temperature measurements available, we regressed the second-stage air temperature predictions against thin plate spline interpolated air temperature values to obtain fully covered air temperature countrywide. We performed internal and external out-ofsample 10-fold cross-validation to quantify the prediction accuracy of our models. All models achieved excellent performance (0.91 ≤ R 2 ≤ 0.98) and low errors (1°C < root-mean-square error < 2°C). The individual subjectspecific data have been linked to the environmental exposure data with daily resolution over the pseudonymized unique participant ID and over the grid ID of the exposure. The datalinkage processes have been performed using the geocoded participants′ residential addresses. Care was taken to remove all spatial information from the data set. Pseudonymized data were provided to the research team. Finally, we matched phenotype data with exposure data by date, according to different lag days.
The daily concentrations of relative humidity (RH), ozone (O 3 ), nitrogen dioxide (NO 2 ), particulate matter with an aerodynamic diameter < 2.5 μm (PM 2.5 ), and black carbon (BC) were assessed via fixed monitoring sites within Augsburg, Germany. 36−38 O 3 and RH were obtained from an official urban background monitoring site operated by the Bavarian Environment Agency (LfU, Bayerisches Landesamt fur Umwelt), which was located about 5 km south of the city center. NO 2 was measured at an urban background measurement station located approximately 2 km north of the city center [also operated by the Bavarian Environment Agency (LfU, Bayerisches Landesamt fur Umwelt)]. PM 2.5 and BC were measured at a single urban background site located 1 km south of the city center and assessed via a tapered element oscillating microbalance (TEOM model 1400A, Thermo Fisher Scientific) equipped with the filter dynamics measurement system (FDMS, model 8500b; Thermo Fisher Scientific) and an aethalometer (model series 8100; Thermo Fisher Scientific), respectively. Daily 24-h average relative humidity, daily maximum 8-h average O 3 , daily 24-h average NO 2 , daily 24-h average PM 2.5 , and daily 24-h average BC were calculated if at least 75% of the hourly measurements were available.

Measurement of Leukocyte Telomere Length.
Standardized procedures for measuring leukocyte telomere length have been described previously in detail. 39 In brief, blood samples were collected on the same day the questionnaires were completed. Before the KORA study center visit, participants were asked to fast for at least 8 h and avoid exercising and smoking the day before and the morning before blood sampling. Blood samples were taken in a sitting position. Samples for measurement of leukocyte telomere length were stored at −80°C until analysis.
Telomere length in samples of F3 and F4 was assayed at the same time and normalized together. Leukocytes from peripheral blood samples were used to extract genomic DNA, and the telomere length was assessed using a quantitative PCR-based method and expressed as the ratio between the telomere repeat copy number (T) and a singlecopy gene: 36B4 (S) (T/S ratio). 39,40 We ran DNA samples in duplicate in 25 μL reactions on a Rotorgene-Q real-time thermal cycler and CAS-1200 liquid handling system (Qiagen, U.K.). A no-template control and a calibrator sample (genomic Environmental Science & Technology pubs.acs.org/est Article DNA from the K562 cell line) in duplicate were included in each run to standardize the measurement across PCR plates. The PCR output was analyzed using comparative quantification (Qiagen Rotorgene analysis software, Qiagen, U.K.), and quantification was relative to calibrator DNA. For quality control, all samples were checked for concordance between duplicate values and to ensure that the experiment ran within the linear range of the assay established. Furthermore, to ensure the reproducibility of the assay, samples were regularly re-run at random on a different day and/or machine. All samples reproduced well, and the inter-run coefficient of variability (CV) for the T/S ratio was 2.63% in KORA F3 and 3.08% in KORA F4, and the intra-run CV for the T/S ratio was 2.84% in KORA F3 and 3.07% in KORA F4. Correlation coefficients (r 2 ) of the measurements from different days were 0.85 for KORA F3 and 0.93 for KORA F4.

Assessment of Covariates.
Participants completed a computer-assisted personal interview and a self-administered questionnaire with information regarding demographic characteristics, including age (year), sex (male, female), education (years), and body mass index (BMI [kg/m 2 ]); lifestyle characteristics, including smoking status, alcohol consumption (g/day), and physical activity; medical history of the participants, including hypertension (no/yes), angina pectoris (no/yes), myocardial infarction (no/yes), and stroke (no/ yes); and current medications, including antihypertensive (no/ yes), lipid-lowering (no/yes), and antidiabetic medication (no/yes). Body weight and height were obtained during physical examinations. BMI was calculated as weight in kilograms divided by height squared in meters.
Four categories of smokers were defined in the questionnaire: never smokers, former smokers, occasional smokers, and regular smokers. Participants were divided into three categories for the present analysis: never-smokers, former smokers, and current smokers (including occasional and regular smokers). Alcohol consumption was rated by asking participants how many alcoholic drinks (beer, wine, or spirits) they consumed during the previous weekday and the previous weekend. Physical activity levels were estimated by asking subjects how much time each week they spent engaging in physical activity during their leisure time in summer and winter (almost no or no physical activity, irregularly about 1 h per week, regularly about 1 h per week, and regularly more than 2 h per week). For this analysis, participants were categorized into three categories: low (almost no exercise), medium (regularly or irregularly approx. 1 h per week), and high (regularly more than 2 h per week). Seasons at blood draw were defined as spring (March−May), summer (June−August), autumn (September−November), and winter (December−February).
Participants with a history of hypertension were defined by measured blood pressure (over 140/90 mm Hg) or reported use of antihypertensive medications. Participant history of myocardial infarction or stroke was based on self-reported physician diagnoses treated in hospitals. The history of angina pectoris was based on self-reported physician diagnoses.
Data on the use of medications were collected using the standardized software IDOM (an instrument for databasesupported online medication registration). 41 The following classes were considered as antihypertensive medication: antihypertensive medication, diuretics, calcium antagonists, βblockers, ACE inhibitors, and angiotensin antagonists; lipidlowering medication: lipid-lowering, including herbal sub-stances, statins, and fibrates; and antidiabetic medication, including antidiabetics, insulin, and oral antidiabetics.

Statistical Analysis.
Characteristics of the study population were given as mean and standard deviation (SD) for continuous variables and frequency and percentage for categorical variables. We compared the demographic characteristics of two cohort study populations using chi-squared tests for categorical variables and t-tests for continuous variables. The levels of meteorological variables and air pollutants were presented as mean, SD, 5, 25%, median, 75, 95%, and interquartile range (IQR), and Spearman's correlation was used to evaluate their correlations.
We applied cohort-specific generalized additive models (GAMs) to explore the short-term effects of mean air temperature on leukocyte telomere length. To explore the potential cumulative effects of mean temperature, we investigated the moving averages of daily mean temperature at lags 0−1, 2−6, 0−6, and 0−13 days before blood draw. In preliminary analyses, we investigated the shape of the association by including the temperature term as a nonlinear function (spline with three degrees of freedom). Since none of the curves indicated significant deviations from linearity (Supporting Information, Figure S1), we included air temperature as a linear term for the main analyses. Leukocyte telomere length was natural log-transformed to increase normality of residuals. We controlled for age, sex, BMI, education, smoking status, alcohol consumption, physical activity, day of the week, season, time trend (cubic spline with ten degrees of freedom), and relative humidity with the same lag period as the air temperature. The time trend was calculated as the day of the year. We included the time trend in our model to adjust for unmeasured confounding, such as variables varying over time, even in this short period.
We conducted effect modification analysis by including an interaction term between air temperature and the potential effect modifier: sex (male vs female), age (<65 years vs ≥ 65 years), physical activity (low vs medium or high), obesity (BMI <30 kg/m 2 vs ≥30 kg/m 2 ), smoking status (current vs former or never smoker), hypertension (yes vs no), cardiovascular disease (defined as a history of angina pectoris, myocardial infarction, or stroke [yes vs no]), season (warm: April−September vs cold: October−March), and O 3

(low [<median O 3 ] vs high [≥median O 3 ]).
To assess the robustness of the results, we performed several sensitivity analyses. First, we applied four different confounder models. Model 1 only adjusted for age, sex, BMI, day of the week, time trend, and relative humidity. Model 2 included all covariates of the main model and additionally adjusted for total cholesterol, triglyceride, LDL, and HDL. Model 3 extended Model 2 by additionally adjusting for medication intake (antihypertensive, lipid-lowering, or antidiabetic medication). In Model 4, we further included air pollutants (O 3, PM 2.5, NO 2 , and BC) with the same lag period as the air temperature. These four pollutants were additionally included in the model separately to avoid collinearity. Further, we only additionally adjusted for BC only in KORA F4 because there was no BC data available in KORA F3. Second, instead of Tmean, we used Tmin and Tmax as alternative exposure metrics. Third, we excluded outliers in leukocyte telomere length values less than Environmental Science & Technology pubs.acs.org/est Article the first quartile of the data (Q1) − 1.5 × IQR or more than the third quartile of the data (Q3) + 1.5 × IQR to avoid overestimating the effects induced by extreme values. Fourth, to avoid the effects of extreme air temperature and the influence of elevated environmental temperature on the quality of blood samples, we excluded subjects exposed to air temperature greater than 95% of temperature. Finally, we extended the lag days and investigated the effect of 4-week moving air temperature averages (lags 0−27 days) on leukocyte telomere length. All models were first fitted separately for each cohort, and the results from the individual cohorts were then pooled using fixed-effect meta-analysis. Heterogeneity was examined using the I-squared (l 2 ) statistic. A P-value > 0.05 and/or l 2 <50% was considered homogeneous.
To compare the effect estimates across different exposure windows, effect estimates were presented as percent changes of the geometric outcome mean with 95% confidence intervals (CIs) per IQR increase in air temperature. P < 0.05 was considered to be statistically significant for all statistical tests. All statistical analyses were performed using R (Version 4.1.2) with "mgcv" and "metafor" packages.

Study Population and Exposure Data.
The characteristics of the study population and level of leukocyte telomere length are described in Table 1. The overall geometric mean level of leukocyte telomere length was lower in KORA F3 with 1.7 T/S compared to 1.8 T/S in KORA F4. The KORA F3 population was older, had a lower level of education, higher proportions of current smokers and hypertension, and higher levels of triglycerides, HDL, and total cholesterol than the KORA F4 population. The KORA F4 population had a higher level of LDL and a higher proportion of antihypertensive medication intake than the KORA F3 population.
Mean daily mean temperature was 7.2°C in KORA F3 and 7.8°C in KORA F4 (Table 2). Spearman's correlation coefficients (r) between meteorological variables and air pollutants (O 3, PM 2.5, NO 2 , and BC) were similar for both studies (Supporting Information, Figure S2). The correlations between the different temperature variables and PM 2.5, NO 2, and BC were generally high (r ≥ 0.72), whereas the correlations between temperature variables, RH, and air pollutants were only weak to moderate.

Short-Term Effects of Air Temperature on Leukocyte Telomere Length.
Our meta-analyses showed that between individuals, higher air temperature was associated with shorter leukocyte telomere length for all investigated lags ( Figure 1). We found that between individuals, an IQR increase in daily mean air temperature (10.77, 10.11, 10.15, and 9.54°C, respectively) was significantly associated with shorter leukocyte telomere length at lags 0−1, 2−6, 0−6, and , respectively). The cohort-specific associations (KORA F3 and KORA F4) showed similar results with slightly more pronounced estimates for F4 and lags 2−6 days for F3 being the only weaker association (Figure 2). Effect estimates expressed as absolute changes between individuals with 95% CIs per IQR increase in air temperature are shown in Figure S3 (Supporting Information).

Effect Modification.
Our analysis mainly indicated no significant effect modification by the investigated potential effect modifiers (Figure 3). Only participants examined in the cold season showed significantly shorter leukocyte telomere length for air temperature at lags 0−1 days. In contrast, no association was observed for participants examined in the warm season.
Effect estimates were presented as percent changes of the geometric outcome mean with 95% CIs per IQR increase in air In general, the associations between air temperature and leukocyte telomere length between individuals remained robust in the sensitivity analyses (Supporting Information, Figure S5). We observed similar associations when using different sets of confounder adjustment. Also, alternative exposure metrics Tmin and Tmax showed comparable effects as Tmean. Moreover, the exclusion of outliers and the exclusion of subjects exposed to air temperatures greater than 95% of temperature did not affect the results. Finally, we still found a significant effect of daily air   46, P = 0.28. All models were adjusted for age, sex, BMI, education, smoking status, alcohol consumption, physical activity, day of the week, season, time trend (cubic spline with ten degrees of freedom), and relative humidity with the same lag period as the air temperature . Effect estimates were presented as percent changes of the geometric outcome mean with 95% CIs per IQR increase in air temperature. The respective IQR increases were 10.77°C for lags 0−1 days, 10.11°C for lags 2−6 days, 10.15°C for lags 0− 6 days, and 9.54°C for lags 0−13 days. Figure 2. Estimated effects (percent change [95% CI] of geometric mean) of air temperature on leukocyte telomere length between individuals at lags 0−1, 2−6, 0−6, and 0−13 days per IQR increase in air temperature in KORA F3 and KORA F4. All models were adjusted for age, sex, BMI, education, smoking status, alcohol consumption, physical activity, day of the week, season, time trend (cubic spline with ten degrees of freedom), and relative humidity with the same lag period as the air temperature. Effect estimates were presented as percent changes of the geometric outcome mean with 95% CIs per IQR increase in air temperature. KORA F3: the respective IQR increases were 12.42°C for lags 0−1 days, 11.62°C for lags 2−6 days, 11.80°C for lags 0−6 days, and 10.74°C for lags 0−13 days; KORA F4: the respective IQR increases were 9.12°C for lags 0−1 days, 8

DISCUSSION
To our knowledge, this is the first study to explore the shortterm effects of air temperature on leukocyte telomere length within an adult population using data from two cohort studies.
Pooling the results of the two cohort studies showed that between individuals, higher daily air temperature was significantly associated with shorter leukocyte telomere length. Leukocyte telomere length measurement is increasingly recognized as a clinical indicator of the risk of age-related disease. 42 Our study is the first to show significant short-term associations between higher air temperature and shorter leukocyte telomere length between individuals in an adult population. To date, only one birth study on prenatal temperature exposure was able to demonstrate an association, finding that the strongest effect of a 1°C increase in air temperature above the heat threshold (19.5°C) took place at week 36 of gestation, resulting in a significantly associated 3.29% (95% CI: 1.88, 4.67%) T/S shorter cord blood telomere length. However, the association with a 1°C decrease in air temperature below the cold threshold (5.0°C) was strongest at week 10 of gestation, with a 0.72% (95% CI: 0.46, 0.97%) T/S longer cord blood telomere length. 30 Another study explored the association between a normal body temperature range (35.0−37.5°C) and leukocyte telomere length in middle-aged and older adults, reporting a significant negative correlation between baseline body temperature and telomere length. However, there was no association between body temperature and the follow-up leukocyte telomere length and the 6-year longitudinal differences in telomere length. 43 We found both significant immediate and lagged effects of air temperature on leukocyte telomere length at lags 0−1, 2−6, 0−6, and 0−13 days. The results indicate that high air temperatures over the preceding 2 weeks are associated with shorter telomere length. A recent epidemiological study based on 1792 participants with overweight/obesity reported a  . It is conceivable that the immediate effects on telomere length on a day-to-day basis are at least partly compensated. However, the findings of this study suggest that repeated exposure to high temperatures may have compounding effects on telomere length. The reversibility of these effects still needs to be explored, in particular, to determine if there is an exposure level or frequency at which effects become partially or totally irreversible. Except for season with lags 0−1 days, no significantly different effects were found across investigated population subgroups concerning age, sex, smoking status, obesity, physical activity, hypertension, cardiovascular disease, or different O 3 -levels. Furthermore, it is noteworthy that we found between individuals significant short-term effects of higher air temperature on shorter leukocyte telomere length in all of these subgroups. This means that the significant effects of air temperature on leukocyte telomere length are robust across all population subgroups. Moreover, we found similar results for different temperature metrics (Tmean, Tmin, and Tmax), again suggesting the stability of our results.
Global warming is an increasingly critical global challenge. The Intergovernmental Panel on Climate Change (IPCC) Working Group II Sixth Assessment Report 2022 reported that there is at least a greater than 50% probability that global warming will reach or exceed 1.5°C in the near term, even for the very low greenhouse gas emissions scenario. 45 Furthermore, leukocyte telomere length shortening is a clinical gauge for biological aging and age-related disease risk. 15−18 Our results showed that between individuals, an IQR/1°C increased daily air temperature was associated with shorter leukocyte telomere length. The findings of our study support the growing evidence that climate change and subsequent temperature increases can lead to adverse human health effects. In light of these findings, public policies should be implemented to decrease the rate of global warming and prevent heat-induced health risks, which, to a certain extent, may help increase lifespan as well as delay or reduce agerelated diseases.
The exact underlying mechanisms by which air temperature impacts leukocyte telomere length are still far from being understood. Telomere lengths tend to shorten with increasing age, and telomere attrition can be accelerated by oxidative stress and inflammation. 28,29 Because of rich guanine content in the 5-TTAGGG-3 repeat sequence, telomeres are highly sensitive to oxidative stress, which causes telomere DNA to be deficient in the repair of single-strand breaks. 46 A recent study among bakery workers found that heat stress can increase the level of malondialdehyde (oxidative stress biomarker). 24 Heat stress also has been shown to induce the production of reactive oxygen species and the increased expression of the interleukin (IL)-8 and IL-8 receptor genes in human dental pulp cells. 25 Increasing short-term air temperature was significantly associated with hypomethylation TLR-2, which may activate the expression of TLR-2 gene and lead to biological responses activating the C-reactive protein expression. 26,27 Another previous study showed that short-term apparent temperature increase was associated with an increase in the acute inflammation factor high-sensitivity C-reactive protein. 47 A significant decrease in leukocyte telomere length was found as high-sensitivity C-reactive protein levels increased. 48,49 Thus, oxidative stress and inflammation may be among the mechanisms by which exposure to increasing temperature could lead to telomere shortening. Furthermore, telomerase, an enzymatic activity that plays a pivotal role in stabilizing telomeres by which telomeric repeats are added to the ends of chromosomes, 50 has been shown to be temperature-dependent, peaking at 37°C under experimental conditions, followed by a decrease in the human telomerase′s processivity with increasing temperature, primer concentration, and potassium ion (K + ). 51 This study has several strengths. We estimated the air temperature at each participant′s residential address for each calendar day by a sophisticated statistical modeling method using satellite, meteorological, and land-use data, which introduced sufficient spatial and temporal air temperature variability and reduced exposure misclassifications compared to a fixed monitoring site. However, no personal temperature exposure measurements have been considered. Second, this study combined two independent cohort studies with a large sample size (total N = 5,864), which should be considered more representative and generalizable than a single cohort study. Third, the KORA cohort (F3 and F4) provided extensive information regarding subject characteristics, which allowed us to adjust for a large number of confounders and explore several potential effect-modifying factors. Fourth, this study gives insight into the short-term effect of air temperature on shorter leukocyte telomere length among an adult population, something that has not been previously studied. However, our study also has some limitations. First, since this is an observational study, we cannot ignore the influence of potential residual confounding. However, we applied several different confounders to our model, and the results stayed robust. Second, our results are based on single-center cohort studies in Augsburg, Germany, and may not be generalizable to other populations with different ethnic, climatic, or demographic conditions. Furthermore, we only used ambient area level air temperature, but we do not know the extent to which individuals were exposed to outdoor temperatures, which increases the likelihood of misclassification of exposures. Finally, due to the lack of data in the present study, the association between air temperature and leukocyte telomere length was not investigated for different leukocyte subtypes, which may be characterized by different telomere lengths.
In conclusion, we found that between individuals, higher daily air temperature was associated with shorter leukocyte telomere length. Our findings add to the burgeoning evidence regarding how increased air temperature can adversely impact human health. ■ ASSOCIATED CONTENT