Comparison of T-cell Receptor Diversity of people with Myalgic Encephalomyelitis versus controls

doi:10.21203/rs.3.rs-3164397/v1

Objective:

Myalgic Encephalomyelitis (ME; sometimes referred to as Chronic Fatigue Syndrome or CFS) is a chronic disease without laboratory test, detailed aetiological understanding or effective therapy. Its symptoms are diverse, but it is distinguished from other fatiguing illnesses by the experience of post-exertional malaise, the worsening of symptoms even after minor physical or mental exertion. Its frequent onset after infection might indicate that it is an autoimmune disease or that it arises from abnormal T-cell activation.

Results:

To test this hypothesis, we sequenced the genomic loci of a/d, b and g T-cell receptors (TCR) from 40 human blood samples from each of four groups: severely affected people with ME/CFS; mildly or moderately affected people with ME/CFS; people diagnosed with Multiple Sclerosis, as disease controls; and, healthy controls. Seeking to automatically classify these individuals’ samples by their TCR repertoires, we applied P-SVM, a machine learning method. However, despite working well on a simulated data set, this approach did not partition samples into the four subgroups, beyond what was expected by chance alone. Our findings do not support the hypothesis that blood samples from people with ME/CFS frequently contain altered T-cell receptor diversity.

The host immune system is able to respond to an extremely diverse array of foreign antigens, in part due to recombination of variable, diversity, and joining (V, D and J) chromosomal segments within T-cell receptor (TCR) a, b, , and loci in thymocytes. Additionally, random nucleotide insertions and deletions occur at the CDR3 (complementarity-determining region 3) where segments are joined, together providing the variation necessary to recognise the diverse range of immunogenic peptides. Rearrangement of human TCR loci yields at least 10¹⁵ different TCR chains [1] dispersed among about 10¹¹ naïve T-cells circulating at any one time [2].

When T-cells are sampled from blood, most TCR sequences are observed only once, although some are found multiply in part due to identical recombination events occurring in the thymus [3] and in part due to clonal expansion of T-cells whose TCR binds to an antigen-bound major histocompatibility complex protein. Clonal expansion of T-cells occurs in disease states, such as Alzheimer’s disease [4], amyotrophic lateral sclerosis-4 [5], IgG4-related disease [6], Kawasaki Disease [7], multiple sclerosis [8], and Ras-associated lymphoproliferative disease [9].

Clonal expansion of T-cells also occurs in response to infection, for example with Epstein-Barr virus in acute infectious mononucleosis [10]. Expansion of a T-cell population occurs after a viral antigen-MHC complex binds its naïve T-cell’s heterodimeric TCRαβ or TCR with sufficient affinity. This increased abundance of T-cells that recognise a specific pathogen elicits an enhanced immune response.

A substantial minority (~10%) of individuals with acute viral or bacterial infection experience a prolonged illness lasting 12 months [11] or longer. This proportion is observed with many infections, including by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [12]. Many people with myalgic encephalomyelitis (ME, also known as chronic fatigue syndrome [CFS]) report an onset of symptoms after infection [13-15]. ME/CFS is a multi-systemic disabling disease that results in a health-related quality of life that is worse than 20 other conditions including chronic inflammatory disorders, cancer and multiple sclerosis [16, 17]. A quarter of people with ME/CFS report being house- or bed-bound [18]. ME/CFS is not rare: it affects, for example, between 836,000 and 2.5 million Americans [16] and an estimated 0.2-0.4% of the UK population [19, 20]. Despite its high prevalence and burden of illness, no reliable biomarker or diagnostic test exists and its aetiology is unknown.

Given that its onset is frequently ascribed to infection and hypotheses that it may be an autoimmune condition [16] we set out to determine whether blood samples from people with ME/CFS tend to show an expansion of TCR clonotypes compared with healthy or disease controls.

Samples and entropy metrics

Forty human peripheral blood mononuclear cell (PBMC) samples were received from each of four groups: (i) Severely affected people with ME/CFS (either house- or bed-bound; MEsa); (ii) Mildly or moderately affected people with ME/CFS (MEmm); (iii) people diagnosed with Multiple Sclerosis (MS; disease controls); and, (iv) Healthy controls (HC). Samples were sourced from the CureME Biobank [21] from female donors aged 40-60 years, chosen to reduce possible age- or sex-dependent confounding effects, although their limited availability among MEsa required us to source samples from younger donors (Table 1).

Table 1. PBMC samples

Age at sample collection	18-29 years	30-39 years	40-49 years	50-60 years
Severely affected people with ME/CFS, MEsa	8	9	6	17
Mildly- or moderately-affected people with ME/CFS, MEmm	0	0	20	20
People with Multiple Sclerosis, MS	0	0	12	28
Healthy controls, HC	0	1	16	23

Some samples were removed from the analysis due to quality or process issues. Six CD8⁺ samples were discarded: one had insufficient PBMCs, and another five libraries had been cross-contaminated. Among CD4⁺ samples, 12 were discarded: two had insufficient PBMCs, four had insufficient DNA (<80ng), three failed the UMI ligation step, two had a barcoding error, and one failed the pooling stage. This left 154 CD8⁺ samples and 148 CD4⁺ samples. Analysis of variance (ANOVA) for the experimental parameters (Table 2) showed that the only association with sample group was CD8⁺ cell number (), although this was not significant after accounting for multiple tests.

Table 2. Mean values of cell count and purity, and numbers of CD8⁺ or CD4⁺ samples.

	PBMCs ()	CD8⁺ sample number	CD8⁺ cells ()	CD8⁺ Purity (%)	CD4⁺ sample number	CD4⁺ cells ()	CD4⁺ Purity (%)
MEsa	5.27	38	0.632	70.7	35	0.665	72.0
MEmm	4.50	38	0.541	70.3	37	0.703	76.0
MS	4.62	39	0.473	68.0	39	0.652	73.5
HC	4.72	39	0.585	69.8	37	0.698	73.6
All samples	4.78	154	0.558	69.7	148	0.679	73.7

Human herpesvirus infection is proposed to be a trigger for ME/CFS [22]. Cytomegalovirus (CMV, a member of the Herpesviridae family) also affects T-cell clonal diversity [23, 24]. CMV infection is very common in the UK [25], and thus could affect TCR clonotype diversity analyses. Nevertheless, CMV seropositivity, available for 98 of the 160 samples, was not significantly associated with sample group status (, χ²-test; Table 3).

Table 3. Cytomegalovirus seropositivity/seronegativity results, shown by group. Equivocal samples were those with indeterminate IgG levels.

Seropositivity	MEsa	MEmm	MS	HC
Positive	14	4	11	12
Negative	19	9	14	13
Equivocal	2	0	0	0

We chose to adopt Rényi entropy as our T-cell receptor diversity metric. Many metrics have been proposed for this purpose, most adapted from standard measures used in ecology [26]. Although Shannon entropy [27] and Simpson's diversity [28] are commonly used in the T-cell literature, these metrics are often applied uncritically, without considering what is being measured [29]. To avoid biasing results by choosing a specific diversity metric prior to analysis, we decided to use the more general, and thus more appropriate, Rényi entropy (Methods). For each cell type and a/b/g-chain combination we defined a vector of clonotype counts in a sample. Here, a clonotype is defined as the CDR3 region, plus the full V, D and J gene segments without considering pairing between a-b or g-d chains. Note that -chain data is not included, because recombination of the locus, which occurs first, preferentially removes the TCR-δ locus [30]. Next, we constructed a matrix of Euclidean distances between each vector pair, ensuring that the pair contained the same number of recombinants by randomly down-sampling the more-populous sample 1,000-times (Methods). We needed to down-sample because the TCR recombinant count varied over two orders of magnitude among study samples. Distances were calculated over a pre-set optimised range of a, the order of Rényi entropy. Next, for cell type and a/b/g-chain combinations we partitioned samples by their clonotypes, adapting the machine learning approach described in Greiff et al. [31] (Methods). Once distances were precomputed, investigators were unblinded to the group identity (e.g. MEmm or MS) of each CD8⁺ sample. CD4⁺ data were acquired following unblinding.

Multidimensional scaling (MDS) plots

The distance matrix for all three CD8⁺ chain types (TCR a-, b-, and g-chains) considered together was visualised using a multidimensional scaling (MDS) plot (Figure 1). The four groups (MEsa, MEmm, MS and HC) are not clearly separated in this plot’s two dimensions. Nevertheless, we noted the presence of seven MEsa, two MEmm and three MS disease cases located away from the main cluster at high x-values (). Age, rather than disease status, likely accounts for these apparent outliers. This is because among 21 variables (covering group, age, PBMC count, CD8⁺ counts pre- and post-selection, purity, DNA amounts, fractions of a/b/g-chains and counts of a/b/g-coding junctions; Table 4), only age (binned by decade) was strongly negatively and significantly correlated with x-axis values (). This likely reflects a narrowing of immune repertoire diversity due to reduced thymic activity (thymopoiesis) as we age [32], which explains our finding that age group is a significant predictor () of CD8⁺ TREC count. Indeed, retrieving the ages at donation of the MEsa outliers (Figure 1) showed that most (four of seven) were from individuals then aged 18-29 years.

Table 4: Variables tested as potential confounders with MDS plot x-axis using a generalised linear model fit.

Variable	Estimate	Std. Error	t-value		Significance
(Intercept)	1.25E+11	2.22E+11	0.564	0.574
ME/CFS (mm)	1.89E+02	3.20E+02	0.590	0.556
ME/CFS (sa)	4.51E+01	3.68E+02	0.123	0.903
MS	4.73E+01	3.23E+02	0.147	0.884
Age	-1.36E+03	4.29E+02	-3.16	0.00199	**
PBMC Number	1.35E+02	1.02E+02	1.32	0.190
CD8+ % Input	-4.54E+01	2.87E+01	-1.58	0.116
CD8+ Cells Post-selection	-7.90E+02	7.85E+02	-1.01	0.316
Purity	6.47E+00	8.34E+00	0.776	0.439
DNA Extracted (μg)	-9.18E+01	1.51E+02	-0.608	0.544
DNA Yield (ng)	6.11E-03	3.90E-01	0.016	0.988
Demultiplexed Reads (641)	1.23E-04	5.47E-05	2.25	0.0263	*
PEAR Assembled (641)	1.19E+01	1.69E+01	0.701	0.485
Demultiplexed Reads (782)	-1.56E-05	4.27E-05	-0.366	0.715
PEAR Assembled (782)	-1.39E+01	1.87E+01	-0.741	0.460
a Coding Junctions	5.74E-01	2.38E-01	2.41	0.0174	*
b Coding Junctions	-4.18E-01	3.32E-01	-1.26	0.211
g Coding Junctions	-8.88E-01	3.07E-01	-2.90	0.00452	**
a %	-1.25E+11	2.22E+11	-0.564	0.574
b %	-1.25E+11	2.22E+11	-0.564	0.574
g %	-1.25E+11	2.22E+11	-0.564	0.574

Significance levels for the t-test statistic are shown: for (***), (**), and (*). ME/CFS(sa), ME/CFS(mm) and MS were testing whether being a member of the given subgroup was correlated with the axis, relative to healthy controls (HC). Age (binned by decade) was the linear term of the model t for that covariate. Other covariates: PBMC Number, total number of cells per sample; CD8⁺% Input, percentage of input PBMCs that were CD8⁺ cells; CD8⁺ cells post-selection, number of CD8⁺ cells after the SureSelect process; Purity, corresponding CD8⁺ purity at this stage; DNA Extracted (μg), total amount of DNA initially recovered from each sample; DNA Yield (ng), amount remaining after size selection; Demultiplexed Reads [Library], total number of reads recovered for the two Illumina flowcells of sufficient quality (library numbers 641 and 782); PEAR Assembled [Library], percentages assembled into clonotypes for the same libraries using the Paired-End Read Merger tool [40]; a/b/g Coding Junctions, total numbers of coding junctions for each chain type; a/b/g %, percentages of each chain type, relative to the other two types (such that the sum was 100%).

Table 5: Variables tested as potential confounders with MDS plot y-axis using a generalised linear model fit.

Variable	Estimate	Std. Error	t-value		Significance
(Intercept)	4.93E+10	8.91E+10	0.553	0.581
ME/CFS (mm)	-1.48E+01	1.29E+02	-0.115	0.908
ME/CFS (sa)	1.77E+02	1.48E+02	1.20	0.235
MS	1.50E+02	1.30E+02	1.16	0.249
Age	3.49E+01	1.72E+02	0.202	0.840
PBMC Number	-2.28E+01	4.11E+01	-0.554	0.580
CD8+ % Input	2.70E+00	1.15E+01	0.234	0.815
CD8+ Cells Post-selection	3.02E+02	3.15E+02	0.959	0.340
Purity	-4.29E+00	3.35E+00	-1.28	0.203
DNA Extracted (μg)	2.42E+00	6.07E+01	0.04	0.968
DNA Yield (ng)	1.88E-01	1.57E-01	1.20	0.232
Demultiplexed Reads (641)	5.19E-05	2.20E-05	2.36	0.0199	*
PEAR Assembled (641)	-5.77E+00	6.79E+00	-0.85	0.397
Demultiplexed Reads (782)	-6.02E-06	1.72E-05	-0.351	0.726
PEAR Assembled (782)	-1.04E+00	7.51E+00	-0.138	0.890
a Coding Junctions	1.56E-01	9.56E-02	1.63	0.105
b Coding Junctions	-4.50E-01	1.33E-01	-3.38	0.000998	***
g Coding Junctions	-2.47E-01	1.23E-01	-2.00	0.0474	*
a %	-4.93E+10	8.91E+10	-0.553	0.581
b %	-4.93E+10	8.91E+10	-0.553	0.581
g %	-4.93E+10	8.91E+10	-0.553	0.581

Significance levels are shown for (***), (**), and (*). Parameters are as for Table 4.

MDS y-axis values were negatively correlated with b (TRB) locus coding junction count (; Table 5). This observation motivated us to perform similar analyses that separated CD8+ recombinants by each of the three loci: TCR a-, b- or g-loci. Again, no clear separation was visible in two dimensions (Figure 2A,B,C), or when a third dimension was added, or when different a-ranges and step size were used (Figure 2D,E). Finally, we undertook this analysis for CD4⁺ (helper T) cells, separating by TCR a-, b-, or g-chain. Although some structure was apparent, again we observed no clear separation between groups in two dimensions (Figure 3).

Potential Support Vector Machine (P-SVM)

Whether clonotypes of one group (e.g. ME/CFS cases) clustered separately from another (e.g. healthy controls) was then assessed using the Potential Support Vector Machine (P-SVM) approach of Hochreiter et al. [33] using leave-one-out cross-validation. Classification boundaries between groups were defined by the P-SVM and classifications were tested against random expectation using permutation testing (up to label shuffles; Methods; Figure 4).

Our three primary aims were to test for a difference in the expansion of T-cell clonotypes inferred from TCR sequencing between: (i) ME/CFS cases (MEsa+MEmm) and healthy controls, or (ii) ME/CFS cases and MS controls, or (iii) MEsa and MEmm cases. Recognition that each of the three hypotheses was tested six times (both CD4⁺ and CD8⁺ cells; and each of a-, b- or g-chains) resulted in a Bonferroni multiple testing correction [34] of for the 18 primary tests. Our six secondary aims were to test for T-cell clonotype expansion differences between: healthy controls and (a) ME/CFS or (b) MEmm or (c) MEsa cases; or between MS cases and (d) healthy controls or (e) MEmm or (f) MEse cases for each of the two cell types and three TCR chains, i.e. tests. Accounting for this total of 54 tests, our Bonferroni-corrected threshold for testing secondary aims was .

Comparison of group classifications against random permutations yielded single-test p-values, each representing the likelihood that two groups had been separated by the P-SVM algorithm by chance alone. Four comparisons (1 of 18 testing a primary hypothesis, and 3 of 36 testing secondary hypotheses) achieved significance at (Table 6). Nevertheless, once the number of tests was accounted for (as above) no comparison remained significant after Bonferroni correction.

We conclude that none of the null hypotheses that we investigated should be rejected.

		CD8⁺			CD4⁺
		a	b	g	a	b	g
Primary Hypotheses	ME vs HC	0.224	1.00	1.00	0.137	0.030	0.053
	ME vs MS	0.236	1.00	1.00	0.663	0.351	0.293
	MEsa vs MEmm	0.395	0.224	0.433	0.396	0.875	0.735
Secondary Hypotheses	Cases vs HC	1.00†	1.00†	1.00†	0.301†	0.104†	0.024†
	MS vs HC	0.998	0.715	0.709	0.379	0.827	0.140
	MEmm vs HC	0.547	0.207	0.301	0.454	0.385	0.231
	MEsa vs HC	0.020	0.597	0.133	0.064	0.193	0.318
	MEmm vs MS	0.193	0.451	0.553	0.135	0.491	0.514
	MEsa vs MS	0.129	0.597	0.267	0.662	0.020	0.588

Table 6. Uncorrected p-values from permutation testing for CD4⁺ and CD8⁺ cells, by a/b/g chain type. p-values are indicated in bold. The † symbol indicates that fewer than 1,000 permutations were used to generate the p-value due to computational constraints (<25 days computation time).

We have applied a method that determines statistical significance of the difference in TCR clonotype repertoires between cases and controls. Using simulated data, the approach perfectly distinguished sample groups at a 1% p-value threshold (Methods). Nevertheless, on experimentally-determined data the method did not detect statistically significant differences between TCR clonotype diversities of people with ME/CFS and others with MS or healthy control individuals, and similarly between the TCR clonotypes of severely affected versus mildly- or moderately-affected people with ME/CFS (MMsa and MMmm). Despite this study comparing a relatively large number of samples from people with ME/CFS ($n=80$) with those from healthy controls ($n=40$) or from disease controls (MS; $n=40$), and despite a large number of recombinants ($\sim{10}^{3}-{10}^{4}$) being sequenced from each sample, no differences between groups were detected. There are four possible explanations for this.

First, given its diverse initiating onsets and symptoms, ME/CFS may yet be found to have heterogeneous biomarkers, each predictive of only a subset of patients. Without such biomarkers, research studies of this size may lack adequate power to identify molecular signatures, such as TCR repertoires, that are predictive of only a minority of patients. Discovery of biomarkers or clinical tests that aid ME/CFS diagnosis remains an urgent priority for future research. Second, a larger study restricted to individuals with particular HLA alleles might have had greater predictive power, because HLA risk alleles are known to modulate autoimmunity risk by increasing the frequency of autoreactive TCRs [35].

Third, even if ME/CFS is aetiologically homogeneous, our negative results could reflect a lack of statistical power to identify true differences. If sampling noise in our experiment’s clonotype counts obscured a clear signal, recombinant numbers may need to be increased by several orders of magnitude, from ${10}^{4}- {10}^{5}$ sequences per donor investigated here, towards the $\sim{10}^{11}$ naïve T-cells in human repertoires [2]. Additional predictive power could be gained by investigating the pairing of TCR chains (e.g. α−β) in single cells, although this would be at the inevitable expense of greater experimental cost.

Finally, this study’s results could also reflect an absence of clonotype diversity differences between groups. If so, then the TCR clonotype sequences themselves, rather than their diversity, could be predictive of disease status, or else TCR repertoire differences are manifest not in blood, but in other more disease-relevant tissues as in MS [36]. Our negative results could also reflect causal mechanisms of ME/CFS that do not result in T-cell repertoire change.

Future work could consider refining the analytical pipeline using samples from a disease with a more clearly established incidence of clonal expansion. A good candidate may be T-cell leukaemia/lymphoma, where research has shown that more than 50% of sample repertoires can consist of a single clonotype [37]. Refining the approach in this way could then allow for the method to detect the potentially subtler differences in disease cases here. However, it is likely that all current TCR repertoire studies are severely underpowered at present, due to current technology greatly undersampling the true T-cell repertoire. This issue is most relevant for investigating specific clonotypes [38], but has relevance to diversity analysis.

Donor sample cell enrichment

160 samples of peripheral blood mononuclear cells (PBMC) were acquired from the CureME Biobank [21] from the same number of anonymised donors. A unique pseudo-anonymised identifier was assigned to each sample upon receipt to allow the study to be blinded. Samples were stored at -180^oC. CureME provides samples from people diagnosed with ME/CFS according to either or both of the Canadian Consensus and Fukuda criteria [21]. Four groups of 40 samples were provided: (i) Severely affected people with ME/CFS (either house- or bed-bound; MEsa); (ii) Mildly or moderately affected people with ME/CFS (MEmm); (iii) people diagnosed with Multiple Sclerosis (MS; disease controls); and, (iv) Healthy controls (HC). As much as possible, these groups were age-matched, with the majority of donors being between 40 and 60 years old at the time of sample collection; all were from female donors (Table 1). Samples’ CMV seropositivity or negativity status was provided by the CureME Biobank based on the level of CMV immunoglobulin G (IgG) measured in their plasma [21]. Data generation and analysis protocols were developed and optimised using two CD8⁺ samples from healthy controls (American Type Culture Collection, Manassas, VA, USA).

PBMCs were sorted by Magnetic Activated Cell Sorting [39] using MACS Micro beads (Miltenyi Biotec) first to positively select CD8⁺ (cytotoxic) cells (CD8 MicroBeads, human 130-045-201), then CD4⁺ (helper/regulatory) cells (CD4 MicroBeads, human 130-045-101), and finally γδ T-cells (TCRγ/δ⁺ T-cell Isolation Kit human 130-092-892). After each selection step, material was retained for verification by Flow Cytometry (FACS) staining using fluorescent conjugated antibodies specific to each of the sorted CD3⁺, CD4⁺, CD8⁺ and TCRγδ⁺ populations (Miltenyi Biotec). Sort results were visualised with plots of co-receptor expression using manually set gates to define discrete cell populations. A minor batch effect was noted, but variability was more pronounced between individuals than between batches, so was considered to have minimal impact.

Sequencing Library Synthesis

After sorting, cells were lysed and DNA extracted using a QIAmp DNA Micro Kit (Qiagen). DNA was fragmented by sonication with a BioRuptor (Diagenode) using 30 rounds at 20 second intervals, size-separated using gel electrophoresis and fragments between 250–350 base pairs excised and recovered using a Clean Gel DNA Recovery kit (Zymo). Yield was measured using a Quant-IT Picogreen DNA assay kit (Thermo Fisher Scientific). Fragments were polished to produce blunt ends using a NextNEB End Repair kit (New England Biolabs).

A custom designed UMI adapter was ligated to each end of the polished gDNA fragments (Integrated DNA Technologies). dsDNA adapted design incorporated, in order from the ligatable blunt end: (i) a 4 base-pair (bp) validation barcode, (ii) a 6 bp random Unique Molecular Identifier (UMI), (iii) a 4 bp library ID barcode, (iv) a 2 bp AcuI target cleavage site, (v) a 14 bp filler sequence, (vi) a 6 bp AcuI binding site, (vii) a 5 thymine, 2 uridine, 5 thymine single base run forming a closed hairpin cap preventing adapter concatemerisation. Following adapter ligation (New England Biolabs), Axyprep Solid Phase Reversible Immobilisation (SPRI) beads were used to size select and clean-up the DNA (Corning). The hair-pinned adapter ends were then opened with a NEB USER kit (NEB) to enable five rounds of enrichment PCR with a Kappa Hifi Polymerase $2\times$ Mastermix kit (Roche). Amplicon yields were measured with PicoGreen (Thermo Fisher Scientific) and re-size selected by gel electrophoresis.

Target Enrichment and Sequencing

A target enrichment strategy was used to select for sequencing library fragments containing signal and coding recombinations of the human T-cell receptor loci on chromosomes 14q11 (TRAD, α/δ locus), 7q34 (TRAB, β locus) and 7p14 (TRAG, γ locus). Biotin-conjugated 120mer ssDNA SureSelect enrichment baits were custom designed (Agilent) to give $2\times$ tiling depth across a 500 bp region centred on each predicted recombination signal sequence (RSS) breakpoint motif (IMGT, https://www.imgt.org). Baits for V regions (1,680) and J regions (996) were synthesised as separate libraries. Each library also contained 17 common baits designed across a 1kb section of the 14q11 α/δ locus TRAC constant region with $2\times$ tiling depth for process quality control and assessment of input library coverage. The bait capture library selecting for fragments with J region homology was employed initially to enrich samples batched into 1.5µg pools of five previously UMI indexed sample libraries, 300ng per donor, using manufacturer’s standard protocols. The output of each J region selection was amplified with five rounds of PCR using Kappa Hifi Polymerase $2\times$ Mastermix (Roche). To enrich for recombinant library fragments containing both V and J elements this product was then used as input for a second round of selection with the V region bait library using the manufacturer’s standard protocols. Final capture products were amplified with ten rounds of PCR using Kappa Hifi Polymerase $2\times$ Mastermix (Roche).

Enriched libraries were prepared for Illumina NovaSeq sequencing. AcuI (NEB) was first used to digest the post-selection amplification products to remove the majority of the custom adapter whilst retaining the 4 bp validation sequence, 6 bp UMI, 4 bp library ID barcode and a predictable, 2 bp “sticky" overhang. Custom designed barcoded Illumina compatible adapters with complementary 2 bp overhangs were then ligated to each library pool (IDT). Following ‘adapter switching’ final library pools were enriched using five cycles of PCR, size selected by gel electrophoresis and sequenced on an Illumina NovaSeq platform using a $2\times$ 250 bp PE sequencing protocol (Arizona Genetics Core).

Data Analysis

A custom suite of Java programmes was used for library demultiplexing and assigning clonotypes. Read pairs were quality and length filtered (all bases with Phred $>20$, >250 bp) then demultiplexed using both Illumina pool indexes and library specific barcodes to identify high quality donor specific subsets of reads. Successfully validated 6 bp UMI sequences were identified from each read in a pair, concatenated to form a fragment specific 12 bp UMI sequence, recorded in the ID field of each read along with barcoding information. These sequences were then trimmed from the fragment ends and reads collapsed to make a library of $1\times$ 250 bp SE QC filtered, UMI annotated and trimmed reads for each donor.

A bioinformatics library of 20 base unique sequence tags was curated to unambiguously identify each V, D or J element of the human T-cell receptor loci on chromosomes 14q11 (TRAD, α/δ locus), 7q34 (TRAB, β locus) and 7p14 (TRAG, γ locus). Tags were positioned 15 bp from the predicted element RSS site and reference sequences recorded (Additional File 1). Where local homology levels at an element prevented unique tag production, non-unique tags (rtags) shared by multiple elements were generated. Secondary discriminating tags (alt-tags) located at the nearest point of divergence were then used to differentiate between rtag elements during processing. Separate libraries of tags were produced for ‘coding’ and ‘signal’ flanks of each RSS. Libraries of donor specific $1\times$ 250 bp SE reads were scanned bioinformatically to identify reads with homology to tags in pairwise conformations compatible with VDJ signal or coding recombination. Parallel comparison of VDJ junction identification was conducted using PEAR [40] pre-assembly of overlapping 250 bp PE reads to reconstitute, where possible, the full-length sequenced fragment.

CDR3 regions were defined for such reads as nucleotide sequences located between the predicted position of identified V, D or J element RSS cleavage sites or where homology to the expected reference approaching the RSS site diverged. Individual V-CDR3-J cassettes were determined to be productive or non-productive by scanning for a single ORF containing the known coding frames of each V or J element identified. Predicted in-frame stop codons were also noted and only productive recombinants were analysed further. Each unique V-CDR3-J sequence was assumed to represent a discrete T-cell clonotype and sequencing reads with identical V-CDR3-J recombination were clustered within each donor library. PCR derived duplicates were identified and removed from these clusters using degenerate matching of fragment 12 bp UMI sequences, allowing for up to two base mismatches. Clonotype duplicates within a donor library with unique UMI sequences were assumed to represent T-cell clonal expansions. The number of unique clonotypes and extent of their expansion was then scored for each donor library to give an index of T-cell diversity.

T-cell Receptor Diversity Estimation

Many metrics used for estimating the diversity of the human T-cell repertoire are specific instances of a more general measure, the Rényi entropy [26].

The Rényi Entropy of order α is:

$${H}_{\alpha }=\frac{1}{1-\alpha }\text{log}\sum _{i=1}^{s}{p}_{i}^{\alpha }$$

where $\alpha$ is the weighting parameter and ${p}_{i}$ is the frequency of the ${i}^{th}$^h category:

$${p}_{i}=\frac{{x}_{i}}{{\sum }_{i}{x}_{i}}$$

and the input is a vector of counts for a sample $x$: $x=({x}_{1}, {x}_{2}, {x}_{3}, \dots , {x}_{s})$. Taking the exponent of the Rényi entropy for input vector $x$ results in a diversity index profile ${D}_{\alpha }={e}^{{H}_{\alpha }}$ for a given array of α-values. Choice of $\alpha$ yields different TCR diversity profiles: higher values of $\alpha$ increase the weight assigned to the most frequently occurring clonotypes, whereas when $\alpha =0$ all clonotypes are weighted equally, irrespective of how many instances of each are found.

To define a generalised distance between each pair of TCR diversity vectors we used their Euclidean distance calculated over a range of $\alpha$ (${\alpha }_{min}\le \alpha \le {\alpha }_{max}$, step 0.1).

The range of $\alpha$-values was optimised before unblinding using a zero-point calibration. For this we first selected a sample $i$, and drew a set of 9 random sub-samples of recombinants, at 10% size intervals, from the original set (90%, 80%, etc) which, together with the original data sample, made a set of 10 (sub-)samples. For each sample pair $(i, j)$, we calculated the Euclidean distance (as above) between pairs of their 10 (sub-)samples, summing over different $\alpha$-interval ranges $[X,Y]$, where $X<Y$, with step size 0.5. The $\alpha$-range that yielded the minimum non-zero median distance between the two samples $(i, j)$ was then considered to be optimal. Based on these results $\alpha$-values were chosen for CD8⁺ cells ($30\le \alpha \le 50$) reflecting the mean values of the lower and upper bounds of 27.1 and 47.1, respectively. An $\alpha$-value step-size of 0.1 was chosen, providing twice the resolution used by Greiff et al. [31]. For consistency, the same parameters were later used for the CD4⁺ cell data.

Using sequenced data from two CD8⁺ healthy control samples, we found that Rényi entropy estimates are not robust to subsampling, in contrast to previous claims [31]. To account for the wide variation in samples’ recombinant counts, we adopted a subsampling approach. For each sample pair $[a,b]$ (with sizes $a>b$) we randomly subsampled $a$ to the size of $b$, repeatedly on $n$ occasions, before calculating the Euclidean distance between $b$ and each $a$ subsample and then taking the median over the $n$ values to define their diversity profiles ${D}_{\alpha }$. Results from the two CD8⁺ healthy control samples led us to choose the number of resamples, $n=\text{1,000}$.

Classification using the Potential Support Vector Machine (P-SVM)

We considered various types of non-parametric supervised algorithms to separate cases and controls into groups by the distances calculated between diversity profiles [41]. We chose to apply a SVM approach as it is a distance-based algorithm that draws a classification boundary between two labelled sets of data, and as it accepts a kernel function to allow non-linear boundaries [33, 41–43]. Further, we chose to implement the P-SVM approach of Hochreiter et al. [33], a modified SVM which was formulated to work on pairwise relationship data, such as distance matrices [33]. This has been shown to typically outperform similar SVM approaches [33], and has the added advantage of having been tested previously on T-cell clonotype counts data [31].

We used the P-SVM to attempt to accurately classify the clonotype repertoires into the four groups: severe ME/CFS, mild/moderate ME/CFS, multiple sclerosis, and healthy controls. The P-SVM handles multiple classes iteratively, finding one classification boundary at a time, using leave-one-out cross-validation (rather than training and testing datasets), with optimised hyperparameters [33, 44]. For this step, a single data point is removed, the model is trained on the remaining data, and then the model predicts the status of the removed data point. This process is repeated for each data point in turn. Another leave-one-out loop is used to select the optimal model hyperparameters, $\epsilon$ and $C$, where $\epsilon$ is the level of tolerance for errors in the training data (if too high the model will fail to capture nuance in the training set, and if it too small it will overfit) and $C$ controls the maximum weighting of the support vectors. If left at the default value ($C=1$) then there are few support vectors, leaving the model more prone to outliers. Setting $C$ too low normally provides a less accurate classification boundary [44]. Optimisation of these two hyperparameters was achieved by minimising the generalization error, here the mean squared error [44].

Significance testing was achieved by randomly shuffling data labels (up to 1,000 times) and then training the machine on these shuffled datasets. This generated a distribution of accuracy scores (sum of true positive and true negative classifications). Statistical significance is achieved at the 5% level when fewer than 5% of scores from the shuffled datasets exceed the original model’s accuracy score. An advantage of this implementation of the P-SVM is that the hyperparameters are set per individual permutation. Whilst this dramatically increases the computation time, it is more robust to false positives, because applying the hyperparameters from the initial data for all permutations will lead to a better fit on the actual data than on the permutations. Applying the hyperparameter values acquired from the optimal fit to real data to permutations was found to yield an erroneous significant result.

Simulation

To demonstrate that the P-SVM could differentiate between clonotype repertoires with or without signs of expansion, we created a simulated data set of 160 samples, each with 10⁴ recombinants. Half were ‘severe cases’, whose recombinant counts were drawn from power law distributions with exponents − 1.5 (40 samples) or -1.4 (40 samples), 40 were ‘mild/moderate cases’ using exponent − 0.2, and 40 were ‘healthy controls’ whose counts were drawn from a uniform distribution. Applying the P-SVM workflow described above (using 100 resamples, $30\le \alpha \le 50$ and $\alpha$ step-size of 1) resulted in perfect discrimination of sample labels at a 1% significance threshold.

Generalized linear models (GLMs)

To investigate MDS axes we implemented a generalized linear model using glm() in R. Similarly, glm() was used to model total CD8⁺ T-cell receptor excision circles (TREC) counts per sample against age group.

Ethics approval and consent to participate. Samples were obtained from the UK ME/CFS Biobank, and ethical approval was given by the University College London Biobank Ethical Review Committee (RFL B-ERC) (ref. EC.2018.006). The study was sponsored by the University of Edinburgh.

Consent for publication. Not applicable.

Availability of data and materials. The CD4⁺ and CD8⁺ sequencing read datasets generated during and/or analysed during the current study are available in the Sequence Read Archive, https://www.ncbi.nlm.nih.gov/sra: SAMN33477947-SAMN33478250 (304 files).

Competing interests. The authors declare that they have no competing interests

Funding. JJD was funded by Action for ME and the Chief Scientist Office, Scotland [AME/CSO/18/01]. BF, MF, MR, NP and SM were funded by the Fischer Family Foundation. CPP was supported by the Medical Research Council University Unit award to the MRC Human Genetics Unit, University of Edinburgh, grant number MC_UU_00007/15.

Authors' contributions

The study was devised by M.D.F., N.J.P. and C.P.P. Experiments were performed by B.F., M.R., M.D.F. and N.J.P. Sequence data analysis was undertaken by S.M., and N.J.P. Machine learning aspects were performed by J.J.D. The manuscript was written by J.J.D., N.J.P. and C.P.P. All authors reviewed the manuscript.

Acknowledgements.

This study was carried out using samples from the UK ME/CFS Biobank, which were processed and stored at the UCL-RFH Biobank. We thank Prof Georg Holländer and Simon J McGrath for valuable discussions during the planning of this project. We also thank Mark M. Davis (Stanford University School of Medicine) whose unpublished TCR studies initially inspired this work.

Nikolich-Zugich J, Slifka MK, Messaoudi I. The many important facets of T-cell repertoire diversity. Nat Rev Immunol. 2004;4:123–32.
Jenkins MK, Chu HH, McLachlan JB, Moon JJ. On the composition of the preimmune repertoire of T cells specific for Peptide-major histocompatibility complex ligands. Annu Rev Immunol. 2010;28:275–94.
de Greef PC, Oakes T, Gerritsen B, Ismail M, Heather JM, Hermsen R, Chain B, de Boer RJ. The naive T-cell receptor repertoire has an extremely broad distribution of clone sizes. Elife 2020, 9.
Gate D, Saligrama N, Leventhal O, Yang AC, Unger MS, Middeldorp J, Chen K, Lehallier B, Channappa D, De Los Santos MB, et al. Clonally expanded CD8 T cells patrol the cerebrospinal fluid in Alzheimer's disease. Nature. 2020;577:399–404.
Campisi L, Chizari S, Ho JSY, Gromova A, Arnold FJ, Mosca L, Mei X, Fstkchyan Y, Torre D, Beharry C, et al. Clonally expanded CD8 T cells characterize amyotrophic lateral sclerosis-4. Nature. 2022;606:945–52.
Mattoo H, Mahajan VS, Maehara T, Deshpande V, Della-Torre E, Wallace ZS, Kulikova M, Drijvers JM, Daccache J, Carruthers MN, et al. Clonal expansion of CD4(+) cytotoxic T lymphocytes in patients with IgG4-related disease. J Allergy Clin Immunol. 2016;138:825–38.
Choi IH, Chwae YJ, Shim WS, Kim DS, Kwon DH, Kim JD, Kim SJ. Clonal expansion of CD8 + T cells in Kawasaki disease. J Immunol. 1997;159:481–6.
Babbe H, Roers A, Waisman A, Lassmann H, Goebels N, Hohlfeld R, Friese M, Schroder R, Deckert M, Schmidt S, et al. Clonal expansions of CD8(+) T cells dominate the T cell infiltrate in active multiple sclerosis lesions as shown by micromanipulation and single cell polymerase chain reaction. J Exp Med. 2000;192:393–404.
Levy-Mendelovich S, Lev A, Rechavi E, Barel O, Golan H, Bielorai B, Neumann Y, Simon AJ, Somech R. T and B cell clonal expansion in Ras-associated lymphoproliferative disease (RALD) as revealed by next-generation sequencing. Clin Exp Immunol. 2017;189:310–7.
Callan MF, Steven N, Krausa P, Wilson JD, Moss PA, Gillespie GM, Bell JI, Rickinson AB, McMichael AJ. Large clonal expansions of CD8 + T cells in acute infectious mononucleosis. Nat Med. 1996;2:906–11.
Hickie I, Davenport T, Wakefield D, Vollmer-Conna U, Cameron B, Vernon SD, Reeves WC, Lloyd A. Dubbo Infection Outcomes Study G: Post-infective and chronic fatigue syndromes precipitated by viral and non-viral pathogens: prospective cohort study. BMJ. 2006;333:575.
Tran VT, Porcher R, Pane I, Ravaud P. Course of post COVID-19 disease symptoms over time in the ComPaRe long COVID prospective e-cohort. Nat Commun. 2022;13:1812.
Jason LA, Yoo S, Bhatia S. Patient perceptions of infectious illnesses preceding Myalgic Encephalomyelitis/Chronic Fatigue Syndrome. Chronic Illn. 2022;18:901–10.
Ghali A, Richa P, Lacout C, Gury A, Beucher AB, Homedan C, Lavigne C, Urbanski G. Epidemiological and clinical factors associated with post-exertional malaise severity in patients with myalgic encephalomyelitis/chronic fatigue syndrome. J Transl Med. 2020;18:246.
Underhill RA. Myalgic encephalomyelitis, chronic fatigue syndrome: An infectious disease. Med Hypotheses. 2015;85:765–73.
IOM (Institute of Medicine). Beyond myalgic encephalomyelitis/chronic fatigue syndrome: Redefining an illness. Washington, DC: The National Academies Press; 2015.
Falk Hvidberg M, Brinth LS, Olesen AV, Petersen KD, Ehlers L. The Health-Related Quality of Life for Patients with Myalgic Encephalomyelitis / Chronic Fatigue Syndrome (ME/CFS). PLoS ONE. 2015;10:e0132421.
Pendergrast T, Brown A, Sunnquist M, Jantke R, Newton JL, Strand EB, Jason LA. Housebound versus nonhousebound patients with myalgic encephalomyelitis and chronic fatigue syndrome. Chronic Illn. 2016;12:292–307.
Baker R, Shaw EJ. Diagnosis and management of chronic fatigue syndrome or myalgic encephalomyelitis (or encephalopathy): summary of NICE guidance. BMJ. 2007;335:446–8.
Nacul LC, Lacerda EM, Pheby D, Campion P, Molokhia M, Fayyaz S, Leite JC, Poland F, Howe A, Drachler ML. Prevalence of myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) in three regions of England: a repeated cross-sectional study in primary care. BMC Med. 2011;9:91.
Lacerda EM, Mudie K, Kingdon CC, Butterworth JD, O'Boyle S, Nacul L. The UK ME/CFS Biobank: A Disease-Specific Biobank for Advancing Clinical Research Into Myalgic Encephalomyelitis/Chronic Fatigue Syndrome. Front Neurol. 2018;9:1026.
Rasa S, Nora-Krukle Z, Henning N, Eliassen E, Shikova E, Harrer T, Scheibenbogen C, Murovska M, Prusty BK. European Network on MC: Chronic viral infections in myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS). J Transl Med. 2018;16:268.
Nguyen TH, Rowntree LC, Pellicci DG, Bird NL, Handel A, Kjer-Nielsen L, Kedzierska K, Kotsimbos TC, Mifsud NA. Recognition of distinct cross-reactive virus-specific CD8 + T cells reveals a unique TCR signature in a clinical setting. J Immunol. 2014;192:5039–49.
Wang GC, Dash P, McCullers JA, Doherty PC, Thomas PG. T cell receptor alphabeta diversity inversely correlates with pathogen-specific antibody levels in human cytomegalovirus infection. Sci Transl Med. 2012;4:128ra142.
Vyse AJ, Hesketh LM, Pebody RG. The burden of infection with cytomegalovirus in England and Wales: how many women are infected in pregnancy? Epidemiol Infect. 2009;137:526–33.
Renyi A. On Measures of Entropy and Information. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability 1961, 1:547–561.
Shannon CE. A Mathematical Theory of Communication. Bell Syst Tech J. 1948;27:379–423.
Simpson EH. Measurement of Diversity. Nature. 1949;163:688.
Alatalo RV. Problems in the Measurement of Evenness in Ecology Oikos 1981, 37:199–204.
Verschuren MC, Wolvers-Tettero IL, Breit TM, Noordzij J, van Wering ER, van Dongen JJ. Preferential rearrangements of the T cell receptor-delta-deleting elements in human T cells. J Immunol. 1997;158:1208–16.
Greiff V, Bhat P, Cook SC, Menzel U, Kang W, Reddy ST. A bioinformatic framework for immune repertoire diversity profiling enables detection of immunological status. Genome Med. 2015;7:49.
Blackman MA, Woodland DL. The narrowing of the CD8 T cell repertoire in old age. Curr Opin Immunol. 2011;23:537–42.
Hochreiter S, Obermayer K. Support vector machines for dyadic data. Neural Comput. 2006;18:1472–510.
Noble WS. How does multiple testing correction work? Nat Biotechnol. 2009;27:1135–7.
Ishigaki K, Lagattuta KA, Luo Y, James EA, Buckner JH, Raychaudhuri S. HLA autoimmune risk alleles restrict the hypervariable region of T cell receptors. Nat Genet. 2022;54:393–402.
Amoriello R, Chernigovskaya M, Greiff V, Carnasciali A, Massacesi L, Barilaro A, Repice AM, Biagioli T, Aldinucci A, Muraro PA, et al. TCR repertoire diversity in Multiple Sclerosis: High-dimensional bioinformatics analysis of sequences from brain, cerebrospinal fluid and peripheral blood. EBioMedicine. 2021;68:103429.
Farmanbar A, Kneller R, Firouzi S. RNA sequencing identifies clonal structure of T-cell repertoires in patients with adult T-cell leukemia/lymphoma. NPJ Genom Med. 2019;4:10.
Warren RL, Freeman JD, Zeng T, Choe G, Munro S, Moore R, Webb JR, Holt RA. Exhaustive T-cell repertoire sequencing of human peripheral blood samples reveals signatures of antigen selection and a directly measured repertoire size of at least 1 million clonotypes. Genome Res. 2011;21:790–7.
Schmitz B, Radbruch A, Kummel T, Wickenhauser C, Korb H, Hansmann ML, Thiele J, Fischer R. Magnetic activated cell sorting (MACS)--a new immunomagnetic method for megakaryocytic cell isolation: comparison of different separation techniques. Eur J Haematol. 1994;52:267–75.
Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics. 2014;30:614–20.
Rashidi HH, Tran NK, Betts EV, Howell LP, Green R. Artificial Intelligence and Machine Learning in Pathology: The Present Landscape of Supervised Methods. Acad Pathol. 2019;6:2374289519873088.
Bzdok D, Krzywinski M, Altman N. Machine learning: supervised methods. Nat Methods. 2018;15:5–6.
Noble WS. What is a support vector machine? Nat Biotechnol. 2006;24:1565–7.
Potential - Support Vector Machine (P-SVM). website https://www.tu.berlin/en/ni/research/software-and-data/psvm Accessed 10 July 2023. Tilman Knebel and Sepp Hochreiter.

No competing interests reported.

AdditionalFile1.xlsx
Additional File 1. (Additional File.xlsx) Human T-cell Receptor loci recombination signal sequence (RSS) and Tag sequences used to unambiguously identify each V, D or J element.

Comparison of T-cell Receptor Diversity of people with Myalgic Encephalomyelitis versus controls

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Results

Discussion

Materials and Methods