Reverse transcription was performed with M-MLV reverse transcriptase (Promega, USA) according to the manufacturer’s protocol. Real-time PCR was carried out on DNA Engine analyzer (BioRad) using RT-PCR kit (Syntol, Russia). As internal control for efficiency of RNA reverse transcription of Poliovirus strain Sabin I was added in the same concentration to each sample before extraction. Reverse transcription was performed with primers Pow-TBE-3’: 5’- AGCGGGTGTTTTTCCGAGTC-3’ for TBEV and PVR1: 5’-CGAACGTGATCCTGAGTGTT-3’ for Poliovirus. Polymerase chain reaction was performed with primer for TBEV R-TBE: 5’-ACACATCACCTCCTTGTCAGACT -3’, F-TBE: 5’-GGGCGGTTCTTGTTCTCC-3’ and probe TBE: 5’-(FAM*)-TGAGCCACCATCACCCAGACACA-(BHQ1*)-3’ for 3’-end of TBEV genome and by using standards with a known concentration of TBEV RNA described earlier (Gmyl et al., 2003). For internal control PVR1, PVL1: 5’-GGCAGACGAGAAATACCCAT -3’ and probe PVP: 5’- (R6G*)-TTGATTCATGAATTTCCTTCATTGGCA-(BHQ1)-3’were used for qRT-PCR (32). The concentration was expressed as the logarithm of the viral genome copy number (GCN) per ml (log10GCNmL). Experiments were repeated three times for each sample to calculate an average value.

Virus was quantified using a plaque assay conducted in 24-well plates with PEK cells under a 1.26% methylcellulose (Sigma, USA) overlay as described earlier (33). Plaques were visualized by staining with 0.4% gentian violet dye (Sigma, USA) and subsequently counted. The number of plaques was expressed as log₁₀PFU per mL. Two parallel replicates were conducted for each sample.

The neuroinvasiveness index for female BALB/c mice was estimated as the ratio between the number of infectious viral particles determined by plaque assay in PEK cells (PFU) and the 50% lethal dose (LD₅₀) of the virus after subcutaneous (s/c) inoculation, calculated according to the Kerber method (41) Table 1. Virulence is presented as a ratio between the genome copy number (GCN) and the 50% lethal dose (LD₅₀) of the virus after subcutaneous (s/c) inoculation as described previously (42) Table 1.

As a murine model, we have selected the BALB/c strain of mice due to its increased susceptibility to TBEV (38, 43). In experiments, female BALB/c mice (8 weeks old) (FSBI “Scientific Center of Biomedical Technologies” of Russian Academy of Sciences, “Stolbovaya” branch) were divided into groups of 10 and injected subcutaneously with 100 μl of virus suspended in Earle’s balanced salt solution (199 medium) (Chumakov FSC R&D IBP RAS, Russia) at an infectious dose of approximately 6 log¹⁰PFU per mouse. The dose of the immature virus sample was equalized by the number of GCN relative to mature EK-328.

This study was carried out according to the international guidelines for the handling of laboratory animals (CIOMS Recommendations, 1985, Directive 2010/63/EC and Annex A of the European Convention ETS No. 123). The study protocol was approved by the Ethics Committee of the FSBSI “Chumakov FSC R&D IBP RAS (protocol #29012018 from 29.01.2018).

The blood from the mice was collected from the submandibular vein into Eppendorf tubes containing 10% (v/v) 0.5 M EDTA (Sigma, USA) at -7 days and 7 days post-virus inoculation. Each blood sample was divided into two aliquots of 100–125 μl. Subsequently, 900 μl of TRIzol (Invitrogen, USA) was added to each aliquot using active pipetting. Subsequently, the samples were frozen at -70°C until further analysis.

Each mouse was anesthetized using 100 μl of Zoletil 100 (Verbac S.A., France) diluted to a 1:10 ratio. After anesthesia, the mouse was secured in a supine position. Perfusion was performed with 5 mL of 0.9% NaCl solution by introducing the solution into the posterior aspect of the left ventricle. After puncturing the left ventricle with a needle, a small incision was made in the right atrium to facilitate perfusion (44). The quality of perfusion was assessed by observing the lightening of the color of the brain and liver.

Spleen was collected at 7 days after virus inoculation to the Eppendorf with RPMI 1640 solution (Chumakov FSC R&D IBP RAS, Russia) and 2% fetal bovine serum (FBS) (Gibco, Invitrogen, USA), 2mM L-Gln (Chumakov FSC R&D IBP RAS, Russia), and 0,1% gentamicin (Complete medium). Each spleen was washed several times with a complete medium, eliminated fat and gently homogenized without loss of splenocyte functionality. Splenocytes were precipitated by centrifugation 10’ at 400g on 4 °C in Eppendorf 5804R. Erythrocytes removal was carried out using ACK Lysing buffer (Gibco, Invitrogen, USA) according to manufacturer’s instructions. The cells were counted in a hemocytometer and diluted to 10⁶/mL. CD4⁺ and CD4^- populations were isolated from splenocytes using magnetic beads (Mouse Dynabeads CD4⁺, Invitrogen, USA) according to the manufacturer’s protocol. 500 ml of Trizol was added directly onto CD4⁺ and CD4^- cells populations(Invitrogen, USA). The samples were frozen at -70°C until further investigation.

To detect the virus in the brain tissue of 8-week BALB/c mice 7 days after infection 10% brain suspension underwent a single passage in PEK cell culture to amplify the virus. The presence of virus in the culture fluid was subsequently determined using (ELISA) commercially available kits (Vector-Best, Novosibirsk, Russia).

Blood samples were collected from mice by decapitation on day seven after infection with the viruses under study. The blood was left at +4 °C for 1 hour, after that it was centrifuged at 300 g for 10 minutes. The serum was aliquoted and stored at -20 °C until further analysis. The plaque reduction neutralization test (PRNT50) was performed on PEK cell monolayers in 24-well following the protocol described previously (38). The proceeding steps were identical to those used for the plaque assay. Each experiment included negative and positive murine sera controls with established antibody titers. Neutralizing antibody (NAb) titers were calculated according to the modified Reed and Muench method (45).

Total RNA was isolated from all samples using TRIzol reagent. Complementary DNA (cDNA) synthesis was performed utilizing 5’-RACE template switch technology, which included double-barcoding. Each cDNA molecule was labeled with a sample barcode and Unique Molecular Identifiers (UMIs). This approach helped not only to introduce universal primer binding sites, but to avoid contamination among samples and to count the amount of start mRNA molecules (46). Mouse TCR profiling kit (MiLaboratory LLC) was used for cDNA TCRβ libraries preparation. Bulk TCRβ sequencing was performed on Illumina HiSeq platform with 2X150 bp sequencing length. Detailed information regarding the procedure can be found in SI Appendix.

Raw data was processed during two stages. At the first stage PCR-errors were eliminated and row reads were clustered by UMIs (unique molecular identifiers) using open-source software tool MiGEC (47). At the second stage the MiGEC output were processed by MiXCR (48) to correct PCR and sequencing errors and to extract CDR3 TCRβ clonotypes with marked V, D, and J segments.

We used EdgeR package (49) to identify the significantly expanded clonotypes after infection with TBEV as previously described (12). For each timepoint we collected and analyzed two biological replicates of whole blood samples. Following the edgeR manual TMM-normalization and trended dispersion estimates was performed. To identify clones significantly expanded between pairs of timepoints we used an exact test based on the quantile-adjusted conditional maximum likelihood (qCML). We assigned clonotype to the group of significantly expanded, if its log₂ fold change estimate log₂FC>7 and its p-value after multiple testing correction was lower than 0.05.

To assign phenotype to expanded clones we used repertoires from two fractions of splenocytes (CD4⁺ fraction and CD4^- fraction) for each mouse individually. The circulating T-cells are either CD4⁺ or CD8⁺, thus here we use CD4^- fraction as a proxy for CD8⁺ fraction. For each clone we calculated its concentration in CD4^- fraction to concentration in CD4⁺ fraction ratio. We assigned clone a CD8 phenotype if this ratio were more than 3, otherwise the clone was identified as CD4⁺. The expanded clones are abundant in both blood and spleen samples after the TBEV infection, therefore there were no expanded clones with undetermined phenotype. Detailed information regarding the procedure can be found in Supplementary Figure 3.

The similarity network was processed by the Igraph package (50). To create similarity network TBEV-associated TCRβ CDR3 amino acid sequences from each mouse were combined, additionally V segment and phenotype (CD4⁺/CD8⁺) were added. TCRβ CDR3 amino acid sequences were used as nodes. Edges are aimed to connect nodes with similar CDR3 sequences. An undirected graph object was constructed from a TBEV-associated TCRβ CDR3 amino acid sequences file. The Hamming distance is used to calculate similarities. Edges were created only for TCRβ CDR3 amino acid sequences with Hamming distance 2 or less. Nodes additional information: phenotype (CD4⁺/CD8⁺), TBEV-variant, V-gene were added as attributes. Visualization of the network was conducted using Gephi (51). The ForceAtlas2 algorithm layout was used. Nodes were colored by their TBEV-variants specificity or phenotype, and sized proportionally to their degree centrality.

Structural motifs of TCRβ CDR3 amino acid sequences from detected clusters were processed by ggseqlogo R package (52). Each letter is specific to amino acids. The probability for each position was calculated using ggseqlogo R package, conserved amino acids correspond to bigger letters with probability equal to 1 and the variable positions are shown as columns of letters, letter size corresponds to probability.

P_gen (probability of generation) is the probability of a certain TCR chain to be generated during V(D)J recombination. P_gen integrates the probability of selecting a particular V, D, and J gene, deletions and insertions during V(D)J recombination. The P_gen for expanded TCRβ CDR3 amino acid sequences was calculated using OLGA software with murine TCRβ settings (53).

In addition to identified TBEV-associated TCRβ clonotypes we analyzed published data of TCRβ CDR3 amino acid sequences associated with Far Eastern TBEV strain (54). Firstly, the data previously published by Fujii et al. were preprocessed to achieve uniformity in V-segments nomenclature. Afterwards the similarities between this data and TBEV-associated CDR3 amino acid sequences found by us were measured using Hamming distance. As no similarities were found, an additional graph with TCRβ CDR3 amino acid sequences previously published by Fujii et al. were created and combined with existing graph of TBEV-associated clonotypes for improved visualization.

Statistical analysis was conducted using the R programming language. The Mann–Whitney U test was employed to compare the quantities of expanded clonotypes. To account for multiple comparisons, resulting p-values were adjusted using the Benjamini–Hochberg (BH) method. The significance of shifts towards CD8⁺ T cells was evaluated using the Wilcoxon Signed-Rank Test.

To identify significant changes in V-segment usage, two-sided t-test p-values were calculated for differences in segment frequencies between TBEV-associated clonotypes and the baseline TCR repertoire at day -7. Given that hypotheses related to specific segments can be considered independent, we applied the Benjamini-Hochberg (BH) procedure with two-sided t-tests to control for multiple testing errors. The false discovery rate during normality testing was maintained at ≤0.05 by setting an upper bound for adjusted p-values at 0.05.

We aimed to compare T cell immune response elicited by three virus variants (991/58, EK, and M) with differing levels of pathogenicity. The 991/58 and EK variants are highly virulent and neuroinvasive, in contrast to the low-virulent, low-neuroinvasive M variant. Additionally we compared the immune response to the mature (EK) and immature (IM) form of the same virus. We studied immune response by identifying specific T cell receptor sequences responding to each virus variant. Each group of mice (n=4) received a subcutaneous injection of one virus variant (Figure 2A). A control group consisting of four animals was injected with PBS. The characteristics of the virus dose used to infect animals are given in the table (Table 2). Blood samples were collected from all mice seven days prior to infection (day -7), and each sample was divided into two replicates for parallel processing. At day +7 after virus inoculation (day 0), all mice were bled again (two replicates) and sacrificed to harvest the spleen. Splenocytes were separated into CD4⁺ and CD8⁺ T cell subpopulations. Total RNA from whole blood and splenocyte CD4⁺ and CD8⁺ fractions was isolated and used for TCRβ library preparation.

Although the blood draw samples for a mouse are very small (typically around 100 μl) we were able to identify on average 48,366 (median = 43526, range 19,726–105,848, n = 40) TCRβ clonotypes per sample replicate taken at day -7 and 43,250 (median = 42400, range 78–101,342, n = 40) TCRβ clonotypes at day +7. CD4+ T cells from the spleen contain on average 79,796 (median = 71508, range 44,482–217,039, n = 20) clonotypes while CD8+ fraction 77,416 (median = 73970, range 38,717–139,784, n = 20) TCRβ clonotypes. Unfortunately, reconstructed day +7 repertoire of two mice infected with 991/58 variant contained insufficient number of clones to detect any alterations (Supplementary Figure 2). These mice were excluded from further analysis. To enable comparative analysis of TCR repertoire architecture, sequencing reads were normalized across samples. Repertoire diversity was quantified using the Gini-Simpson index, and clonality was assessed via the DE50 metric (Figure 2B). The results show that the values for both metrics are close to 1 indicating highly diverse repertoires with balanced clonotype frequencies, reflecting high population evenness. With the exception of the EK variant, no significant changes in diversity or clonality were observed. The infection with EK variant exhibited a significant reduction in both diversity and evenness, indicative of an antigen-driven immune response characterized by the clonal expansion of specific T cells. To investigate potential minor changes in repertoire architectures that were not captured by these broad metrics, further in-depth analysis was performed.

The identification of virus-responding clonotypes for each mouse was performed using a previously developed approach (12). TCRβ repertoires from blood replicates of the same day showed strong convergence, demonstrating high reproducibility. Data for single mouse infected with IM is shown (Figure 2C, left and right panel). This reproducible data allowed for a comparison of TCRβ repertoires from whole blood between day -7 and day +7 using EdgeR software. The normalized frequency of each clonotype was calculated at 7 days before and at 7 days after infection, a statistical Generalized Linear Model was employed to assess the probability of clone expansion in response to viral infection. We determined the TCRβ clonotype was expanded only if its concentration significantly changed after infection (Materials and Methods). Every mouse underwent individual analysis according to this approach (Supplementary Figure 5). Some expanded clones were undetectable at day -7, while the abundance of the others drastically increased after TBEV-infection (Figure 2C, middle panel).

To analyze the magnitude of the T cell response, we compared the number of the detected TCRβ clonotypes between groups of mice infected with different TBEV variants. The number of such clonotypes varied from 24 (IM variant) to 213 (991/58 variant) for infected animals, and from 2 to 9 for the control animals. Infected animals exhibited a greater number of expanded clonotypes compared to the control group (Mann–Whitney test, p = 0.003) (Figure 2D, first panel). For the selected comparison group we combined data from mice by virulence and neuroinvasiveness. A statistically significant difference was observed in the number of expanded TCRβ clonotypes between animals infected with the low-virulent M-variant and those infected with the high-virulent TBEV variants (991/58 and EK) as determined by Mann–Whitney tests (p = 0.038) (Figure 3D, second panel). Analysis of each variant independently shows the same trend (Figure 2D, third panel). No statistically significant differences in the magnitude of the T cell response were recorded between the mature and immature forms of the viruses. However, a trend toward an increased number of responded T cell clonotypes to the mature form of the EK virus was detected (Figure 2D, fourth panel).

Thus, we show that TBEV infection induces proliferation of a specific group of TBEV-associated T cells, and infection caused by high-virulent and high-neuroinvasive TBEV variants induce stronger T cell response than infection caused by low-virulent and low-neuroinvasive TBEV variants.

In order to identify CD4⁺/CD8⁺ phenotype of the TBEV-associated clonotypes we used magnetic separation of splenocytes at day +7. Based on the presence of the responding clonotype sequence in CD4⁺ or CD4^- TCRβ repertoire generated from spleen T cells, we assigned each expanded clone to CD4⁺ or CD8⁺ T cell subset (Supplementary Figure 6). In examining the shift towards a cytotoxic T cell response, we calculated the average count of CD8⁺ cells by comparing the non-TBEV-specific TCRβ repertoire of equivalent size to that of the TBEV-associated repertoire for each mouse individually (Figure 3A). This approach allowed us to assess whether there was an enrichment of CD8⁺ cells in response to TBEV infection relative to background levels, thereby providing insights into the role of cytotoxic T cells in this context. Our results reveal an increased presence of CD8⁺ cells within the TBEV-associated repertoires. Statistical significance was observed exclusively for TCRβ clonotypes associated with the EK and IM variants, as determined by Wilcoxon tests (p = 0.03 for both). This finding suggests a trend towards preferential recruitment and expansion of cytotoxic CD8⁺ T cells towards virus-infected cells.

We observed significant alterations in the repertoire structure within the most abundant TCRβ clonotypes following TBEV infection (Figure 3B). Within the top 100 most abundant clones at day -7 none were identified as TBEV-associated. By the day +7 these clones exhibited little changes in their abundance levels and largely remained among the top of the TCRβ repertoire. These clonotypes represent the overall stability of TCR repertoire. In contrast, by day +7 post-infection, the top repertoire of infected mice contained between 21 and 72 novel expanded clones, whereas only between 2 and 9 expanded clones were observed in the control group (Mann–Whitney, p = 0.03). Notably, some of the most frequent clones at day +7 had low concentrations or were undetectable at day -7. In conclusion, these findings suggest that strong immunogenic stimuli like TBEV infection can profoundly impact the range of T cell clones within a given repertoire.

The selection of particular V genes in TCR sequence has been observed in response to immune challenges (8, 14, 56–59). To further evaluate the structural characteristics of virus-responding clones, we performed an analysis of V-segment frequency among these clonotypes within groups of animals infected with different TBEV variants. For each group of mice, we computed the frequency of each V-segment among the TBEV-associated clones and among total T cell clone repertoire (Figure 3C, Supplementary Spreadsheet 2). Our findings indicate a statistically significant increase in the utilization of TRBV13 and TRBV16 for mice infected with the 991/58 virus variant, TRBV3 and TRBV16 for EK, and TRBV29 for mice infected with M or IM variants. We consistently observed a pattern in responding TCRβ clonotypes to neuroinvasive 991/58 and EK TBEV variants utilizing TRBV16. Surprisingly we did not find any similarities in V-segment preferences for mature EK and its immature form IM. These observations may indirectly imply a similarity in TCR sequences that respond to distinct TBEV variants associated. It also suggests that different epitopes might be involved in triggering T cell reactions for these variants, potentially influencing the nature of the immune response and disease outcomes.

In order to identify structural similarities among TBEV responding T cells we constructed TCRβ CDR3 amino acid sequence similarity networks for all clonotypes expanded after TBEV infection (Figure 4A). The largest cluster (#1) includes 10 CDR3 sequences identified in response to each of the investigated TBEV variants. These sequences are characterized by TRBV29 and TRBJ2-1. We also observed two clusters of 7 TCRβ CDR3 sequences formed by clones characterized for 3 TBEV variants. Our analysis revealed that clusters containing 4–5 CDR3 sequences were composed of TCRβ clonotypes that are specific to the 991/58 and EK variants.

The majority of TBEV-specific T cells, whose TCRβ CDR3 sequences form clusters, are CD8⁺ cells. Notably, the largest cluster identified consists entirely of CD8⁺ T cells (Figure 4B). This clustering pattern suggests that these CD8⁺ T cell populations play a significant role in antiviral activity against TBEV. We also identified amino acid consensus motifs within the TCRβ CDR3 of the found clustered TBEV-specific T cells (Figure 4C). We computed the probability of generation (P_gen) for TBEV-associated TCRβ amino acid sequences to evaluate their likelihood during V(D)J recombination. We observed a substantial range in P_gen values across the TBEV-associated TCRβ sequences, moreover TCRβ sequences, which are forming clusters, are characterized by P_gen 10^–5 to 1.3×10^-9. This broad distribution indicates that TBEV-associated TCRβ amino acid sequences from clusters are statistically common products of V(D)J recombination, indicating them as “public” clones (Supplementary Figure 7). Clusters with higher P_gen values contribute to the public TCR repertoire, but the likelihood of random (independent) generation of such receptors across different thymocytes or organisms remains relatively low (53). The appearance of considerable quantities of clones carrying these TCRβ sequences requires selective processes and clonal expansion, which are driven by antigen recognition and immune response dynamics.

To demonstrate that the identified TCRβ CDR3 clusters were formed as a result of TBEV infection rather than being a feature of the TCR repertoire, we randomly extracted TCRβ clonotypes from day -7 samples of each mouse. The number of clones extracted was equivalent to the size of the TBEV-associated TCRβ repertoire. We then calculated the number of edges formed by their amino acid sequences. This process was repeated 1,000 times (Figure 4D). The analysis revealed that the TBEV-associated TCRs formed clusters containing more than twice the average number of clusters observed in the randomized day -7 repertoires, supporting their association with the infection.

We also conducted an analysis of published data on murine TCR amino acid sequences expanded following infection with another TBEV strain (54). Fujii et al. infected C57BL/6 mice with Oshima 5–10 strain. On day 13 post-infection RNA was isolated from the brain and spleen. Nucleotide sequences of TCR were obtained with the use of specific primers annealed to three only V-segments and constant region, afterwards cloned and sequenced by Senger method. Authors observed an increased level of certain V-genes in the brain of infected mice, also identified differences in V-gene families’ frequencies from spleen samples between dying and surviving mice. Amino acid sequences of TCR CDR3 regions of cDNA clones derived from TBEV-infected mouse brains are published by the authors. In order to find similarities between our and published data we looked for exact matches, afterwards matches with Hamming distance 1 and 2. However, no similarities were identified between these sequences and those observed in our study (Supplementary Figure 8).