Phenotyping the virulence of SARS-CoV-2 variants in hamsters by digital pathology and machine learning

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has continued to evolve throughout the coronavirus disease-19 (COVID-19) pandemic, giving rise to multiple variants of concern (VOCs) with different biological properties. As the pandemic progresses, it will be essential to test in near real time the potential of any new emerging variant to cause severe disease. BA.1 (Omicron) was shown to be attenuated compared to the previous VOCs like Delta, but it is possible that newly emerging variants may regain a virulent phenotype. Hamsters have been proven to be an exceedingly good model for SARS-CoV-2 pathogenesis. Here, we aimed to develop robust quantitative pipelines to assess the virulence of SARS-CoV-2 variants in hamsters. We used various approaches including RNAseq, RNA in situ hybridization, immunohistochemistry, and digital pathology, including software assisted whole section imaging and downstream automatic analyses enhanced by machine learning, to develop methods to assess and quantify virus-induced pulmonary lesions in an unbiased manner. Initially, we used Delta and Omicron to develop our experimental pipelines. We then assessed the virulence of recent Omicron sub-lineages including BA.5, XBB, BQ.1.18, BA.2, BA.2.75 and EG.5.1. We show that in experimentally infected hamsters, accurate quantification of alveolar epithelial hyperplasia and macrophage infiltrates represent robust markers for assessing the extent of virus-induced pulmonary pathology, and hence virus virulence. In addition, using these pipelines, we could reveal how some Omicron sub-lineages (e.g., BA.2.75 and EG.5.1) have regained virulence compared to the original BA.1. Finally, to maximise the utility of the digital pathology pipelines reported in our study, we developed an online repository containing representative whole organ histopathology sections that can be visualised at variable magnifications (https://covid-atlas.cvr.gla.ac.uk). Overall, this pipeline can provide unbiased and invaluable data for rapidly assessing newly emerging variants and their potential to cause severe disease.


Introduction
As the coronavirus disease-19 (COVID-19) pandemic progressed over the past three years, the virus responsible for the disease, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has continued to evolve giving rise to a number of variants, some of which were defined as "variants of concern" (VOCs) by the World Health Organisation (WHO) [1,2].These VOCs contain mutations, especially but not solely in the spike encoding S-gene, which may confer a selective advantage for example by increasing their transmissibility and/or immune evasion compared to the progenitor virus [1].To date, the WHO has recognised five VOCs: B.1.1.7 (Alpha); B.1.351(Beta); P.1 (Gamma); B.1.617.2 (Delta) and B.1.1.529(Omicron; henceforth referred as BA.1 to differentiate it from other sub-lineages) [1,[3][4][5].Since the emergence of the original BA.1 in November 2021 [5] different sub-lineages have emerged.Soon after the BA.1 emergence, BA.2 became the predominant variant followed by a variety of BA.2 descendants, including BA.5 and BA.2.75, which became predominant in some geographical regions [1,6].Although both BA.5 and BA.2.75 diversified from BA.2, these two Omicron sub-lineages are phylogenetically separated from each other, suggesting that BA.5 and BA.2.75 emerged independently.BQ.1.18is a sub-lineage of BA.5 [7] while the XBB variant is a recombinant between BJ.1 (a BA.2.10 derivative) and BM.1.1.1 (a descendant of BA.2.75), which was first detected in September 2022 in India and spread significantly at the end of 2022 [8].EG.5.1 is a descendent lineage of XBB.1.9.2 but carries two additional spike mutations F456L (defining EG.5 lineage) and Q52H.It was first detected in early March 2023 and was then designated as a variant of interest (VOI) by WHO on 9 August 2023 along with EG.5 and its sub-lineages.EG.5.1 has been detected globally with a rapid increased prevalence to around 50% as of 2nd October 2023 [9].
In general, each new variant spreading globally tends to be more transmissible than the previous dominant variant [1].As population immunity increases, due to either vaccination or continuous virus exposure, it is likely that COVID-19 will adopt an endemic pattern, possibly with seasonal peaks, primed by variants evading pre-existing immunity in the population [10].
Understanding in "real-time" the degree of vaccine escape of any new variant is critical to determine vaccination policies.To this end, in vitro seroneutralisation assays have proven to be a useful surrogate to predict vaccine escape of SARS-CoV-2 variants.Virus virulence is another key phenotypic characteristic of any new variant requiring early assessment.The risk of severe disease and hospitalisation varies, with Alpha, Gamma and Delta VOCs carrying an increased risk of intensive care unit (ICU) admission compared to Beta and BA.1.Hence, the combination of increasing pre-existing immunity in the population together with the intrinsic attenuated characteristics of Omicron, has led to an overall decrease in the incidence of severe disease and mortality associated with COVID-19 [11][12][13][14].
It is however more difficult to predict the trajectory of new variants with respect to virus virulence.Although BA.1 has been shown to be attenuated, and assuming that its transmission potential is maintained, there are no universal evolutionary pressures that may keep this phenotypic trait in newly emerging sub-lineages or new VOCs.
We and others have shown that the observed reduction in virulence of the BA.1 variant correlates to a change in the virus entry pathways in vitro, and importantly in reduced virulence in experimentally infected hamsters [15][16][17][18].Throughout the pandemic, small animal models have been used extensively to assess the virulence of wild type SARS-CoV-2 and emerging VOCs [8,15,[17][18][19][20][21][22][23][24][25][26][27].These studies have provided invaluable data on disease pathogenesis, virus transmission and the efficacy of different anti-viral compounds or vaccines [20,[28][29][30][31][32][33].Importantly, hamsters have been shown to be naturally susceptible to SARS-CoV-2 infection and to be able to transmit the virus to humans [34].In hamsters, BA.1 is unable to infect lung epithelial cells (unlike the original B.1 virus, Delta and other variants) [15,17].In addition, in experimentally infected hamsters it is also possible to recapitulate the increased virulence of the Delta variant compared to B.1 shown in the human population [24,35].
Hence, although no animal model can fully recapitulate a human disease, hamsters represent an excellent model to dissect SARS-CoV-2 pathogenesis and determine the degree of virulence of newly emerging variants.To this end, many studies using experimentally infected hamsters, have focused on measuring in vivo viral replication, on identifying virus infected cells, and on examining pathogenic potential by measuring weight loss and assessing various histopathological criteria in general by qualitative scores [8,[18][19][20][21].Here, we have developed unbiased and automated "digital pathology" methods to assess SARS-CoV-2 virulence.Digital pathology is a broad term that refers to a variety of systems to digitize pathology slides and associated meta-data, their storage, review, analysis, and enabling infrastructure [36].Computational analysis of whole scanned tissue sections provides the opportunity to quantify cells or histological features in wide representative areas of infected organs.We applied these pipelines with recently evolved variants (BA.5, BQ.1.18,BXX, BA.2 and BA.2.75) [1,[6][7][8] and showed that some of them have gained a more virulent phenotype compared to the parent BA.1.This pipeline can contribute to the rapid assessment of newly emerging variants and should prove invaluable as the pandemic enters the next phase.Furthermore, we created an online repository to share with the scientific community high resolution digitized whole organ scanned slides from this study, providing a wider context to histopathology micrographs for experimental models of COVID-19.

Host transcriptional response to SARS-CoV-2 infection
To assess the complex host responses during SARS-CoV-2 infection, we first experimentally infected Golden Syrian hamsters with either the Delta (B.1.617.2) or the BA.1 variants.These variants are on the opposite spectrum of the phenotype associated with the clinical outcome of SARS-CoV-2 infection, both in humans and in experimental models.Hence, Delta and BA.1 can provide the baseline for the development of quantitative pipelines to assess virus virulence.As expected, between 2-and 6-days post-infection (dpi) the Delta-infected animals lost significantly more weight and had higher welfare scores than both the BA.1-and mock-infected hamsters (S1A Fig), confirming the expected phenotype for both variants.We culled animals at either 2 or 6 dpi and collected tissues of both the upper and lower respiratory tract, in addition to peripheral blood, for downstream analyses.
To understand the overall host responses to SARS-CoV-2 infection and identify potential markers of virus virulence, we performed bulk RNAseq on lungs and blood of both infected and mock-infected hamsters.Principal component analysis (PCA) indicated distinct clustering of both the Delta-and BA.1-infected groups at 2 dpi in lungs and blood but showed less distinctive separation of BA.1 and mock-infected animals at 6 dpi (S1B Fig) .Both in the lungs and in peripheral blood, Delta, and BA.1 induced significant differential gene expression compared to the mock-infected samples (Fig 1A).As expected, both variants induced multiple cell responses pathways associated with antiviral mechanisms, immune system activation, cytokine and chemokine responses and interferon signalling as evident by gene ontology (GO) analysis (Fig 1B and 1C).Direct GO analysis between the two variants found a limited number of pathways on day 6 upregulated only in delta-infected animals.These pathways include those associated with remodelling and lipoprotein regulation, suggesting that the differences in disease outcome caused by these two variants may be attributed to lesions and activation of tissue repair pathways in the lungs.Comparison of ~200 interferon stimulated genes showed a general activation of the type-I IFN response in infected animals (as suggested by the GO analysis and previously published studies) both in the lungs and in the blood [37][38][39].Delta-infected animals showed a more robust upregulation of interferon stimulated genes (ISGs) than BA.1-infected hamsters in both lungs and blood (Fig 1D and 1E).Overall, these analyses suggest that markers of type-I IFN response, immune system activation and tissue repair may be useful to characterise the extent of pathology induced by SARS-CoV-2 variants.

Imaging and quantification of SARS-CoV-2 replication in tissues
Next, we characterised the spread of SARS-CoV-2 infection in infected tissues.Throughout this study, for the detection of virus and cellular markers, we aimed to develop unbiased quantitative methods using software assisted whole section imaging and downstream automatic analyses including those enhanced by machine learning approaches (henceforth referred with the broad term of digital pathology) [40].The development of the digital pathology pipeline presented in this study has been evaluated by a board-certified veterinary pathologist (VH).In addition, side-by-side comparisons between standard histopathological scores and the digital pathology pipeline was carried out by three additional board-certified veterinary pathologists (FH, AH, CA).Also, in order to share the histopathology features shown in this study as comprehensively as possible, we have developed an online repository ("CVR Virtual Microscopy"; https://covid-atlas.cvr.gla.ac.uk)where whole organ scanned images can be accessed by users in their entirety and at variable magnification as if they were observing slides under a microscope.

PLOS PATHOGENS
Assessing SARS-CoV-2 virulence by digital pathology First, we compared the detection levels of SARS-CoV-2 nucleocapsid protein in the lungs of infected hamsters at 2 dpi by immunohistochemistry (IHC) with its spike RNA by in situhybridisation.The two methods resulted in essentially identical staining patterns (S2 Fig) .Background staining in mock-infected samples was higher using IHC and therefore we used RNA in situ-hybridisation for the remaining part of the study.Next, we assessed virus replication for Delta and BA.1 along the entire respiratory tract.As expected, tissues collected from animals culled at 2 dpi showed both Delta-and BA.1-infected cells within the respiratory tract (Fig 2A and 2B).We found instead little evidence of virus infected cells in animals culled at 6 dpi.At 2 dpi, we found no significant differences in the number of infected cells in the nose, larynx, and trachea of Delta-and BA.1-infected animals.In our samples, the nose represents the inner mucosa of the small cartilaginous tissue surrounding the nasal cavities of the hamster.However, both nasal turbinates of Delta-infected animals and the lungs showed a significantly higher number of infected cells than the same tissues in BA.1-infected hamsters.Importantly, significant spread of SARS-CoV-2 in the lung parenchyma was evident only in Delta-infected hamsters.Delta infected both epithelial cells in the bronchioles and alveoli forming large foci of infected tissues.Conversely, BA.1 infected only cells in the bronchioles and at most formed small foci of infection in the lung parenchyma in some animal (Fig 2A and 2B).As expected, there was variability between animals within each group with respect to the number of infected cells.This variability was especially evident in the trachea, with 2 of 8 Delta-infected animals showing a number of infected cells more than 10-fold the number in the remaining animals of the group.Except for the two outliers indicated above, the tracheas of the remaining Delta-and BA.1-infected hamsters showed a similar number of infected cells.
Given the data obtained by RNAseq, where the interferon response is a key differentially activated pathway in infected hamsters, we also assessed expression of MX1, as a representative core ISG [41].MX1 expression in the nose and lungs was in general lower at 2 dpi compared to the levels observed at 6 dpi.The larynx showed high MX1 expression on 2 dpi but little on 6 dpi, while the turbinates and trachea showed similar levels at both timepoints (Fig 2C).Bulk RNAseq data suggested a more robust type-I IFN response elicited by the Delta variant compared to BA.1 in lungs from infected hamsters (Fig 1B -1E).We aimed therefore to spatially resolve ISG expression in infected hamsters using in situ RNA hybridisation of serial sections of lungs collected at 2 dpi.We used five sections with the middle section probed for spike, while in the other sections we used RNA probes for the following ISGs: RSAD2, IFIT1, MX1 and OAS1 (Fig 3).The percentage of ISG-positive pixels was directly related to the percentage of spike-positive pixels in the lungs of infected hamsters.There were clear overlapping virusand ISG-positive areas in both Delta-and BA.1-infected lungs.Hence, the more robust type-I IFN response observed by RNAseq in the lungs and blood of infected animals correlates with higher replication levels of Delta in the lungs at 2 dpi.
The differences between Delta-and BA.1-infected hamsters were also present at earlier timepoints.Analysis of tissues collected at 1 dpi of the turbinates, trachea, and lungs (S3A Fig) showed increased levels of spike RNA (S3B Fig) in Delta-infected animals.MX1 expression showed a trend of higher levels of expression in Delta-infected animals but differences were not statistically significant (S3C Fig) .Similarly, nasal washes, throat swabs and whole lung RT-qPCR indicated higher levels of genomic RNA in Delta-infected hamsters, but only in the latter differences were statistically significant (S3D Fig).

Quantifying the extent of pulmonary lesions by digital pathology
We showed above that at 6 dpi, at the peak of clinical signs in experimentally infected hamsters, most of SARS-CoV-2 has been cleared by the host (Fig 2A).Our RNAseq analysis  Hence, we next aimed to image and quantify the extent of virus-induced pulmonary pathology by first evaluating and comparing the immune cells infiltrate in the lungs of Delta-and BA.1-infected hamsters at 6 dpi.We specifically assessed T cells (CD3 + ) and macrophages (IBA1 + ) and found a significant increase in the number of these cells in Delta-infected hamsters compared to those infected with BA.1 (Fig 5A and 5B).By histopathology, both cell types represent most immune cell infiltrates in the lungs of Delta-infected hamsters (Fig 5A and 5B).
We next developed a method to quantify alveolar epithelial hyperplasia.We used the thyroid transcription factor (TTF1) [44,45], a critical factor required for the expression of the surfactant protein in the respiratory epithelia.As expected, we found that TTF1-positive cells included the hyperplastic alveolar epithelial cells but also normal/isolated type 2 pneumocytes in the lungs and epithelial cells in the terminal bronchioles.To quantify only the hyperplastic areas in the lung parenchyma, we used software assisted imaging detection.Using supervised machine learning approaches, we "trained" the HALO software (Indica Labs) to detect clusters of TTF1-positive nuclear areas representing proliferating type-2-pneumocytes while ignoring isolated type 2 pneumocytes or TTF1 + cells in the terminal bronchi.During the training, we made sure that the software excluded hyperplastic bronchial epithelium (S6 Fig) .We found no hyperplastic lesions at 2 dpi in any of the hamster groups while at day 6 we found significantly more hyperplastic type 2 pneumocytes in Delta-infected hamsters compared to those infected with BA.1 (the latter had only values just above background in most animals; Fig 5C and 5D).

Assessing the virulence of Omicron sub-lineages
The pipelines developed so far allow us to provide an automatic, unbiased, and quantitative method to assess the degree of virulence of SARS-CoV-2 in hamsters.Hence, we next used this method to assess virulence of recently emerged omicron sub-lineages.
We first investigated BA.5, as other studies, although carried out with chimeric BA.2/ BA.5 viruses, suggested that this variant had an increased virulence compared to BA.1 [23].Lungs collected from hamsters infected with BA.5 showed in comparison to those infected with BA.1, a consistent trend of increased values for (i) macrophage infiltrate (IBA-1 + cells), (ii) cells expressing MX-1 and (iii) alveolar epithelial hyperplasia.Differences however lacked statistical significance (Fig 6A).BA.5-infected animals also showed a small decrease in body weight at 6 dpi, while animals infected with Delta or BA.1 showed the expected phenotype (slight increase in weight for BA.1-infected animals and weight loss for Delta from 2 dpi; Fig 6B).
Given that the digital pathology pipeline was able to show a possible intermediate virulence phenotype for BA.5, we proceeded to assess the virulence of other Omicron sub-lineages such as BQ.1.18,XBB and BA.2.75 (Fig 6C -6E).None of the variants induced major weight loss in infected hamsters (Fig 6C).Hamsters infected with BA.2.75 showed a modest weight loss between dpi 2 and 5 but neither hamsters infected with BA.2.75 nor BQ.1.18show, unlike animals infected with BA.1, weight increase until 5 dpi (Fig 6C).
We found an increase in macrophages infiltrating the lungs in hamsters infected with the other variants compared to those infected with the original BA.1 (Fig 6D).BA.2.75 showed the highest levels of IBA1 + cells in the lungs, although differences were statistically significant only with BA.1 and XBB, but not with BQ.1.18.Lungs of animals infected with BA.2.75 also displayed a significant increase of alveolar epithelial hyperplasia compared to lungs of BA.1-infected animals (Fig 6E).Overall, the data suggest that virulence of Omicron sub-lineages, and especially BA.5 and BA.2.75 has increased compared to BA.1.We also compared BA.2.75 to its  predecessor BA.2 to determine whether there had been further evolutionary adaptations in BA.2.75 that favoured a more virulent phenotype.No significant differences were observed for weight losses, amount of macrophages and alveolar epithelial proliferation between BA.2.75 and BA.2 (Fig 6F -6H).However, animals infected with BA.2.75 or those infected with BA.2 did not gain weight as steadily as mock-or BA.1-infected animals.We saw however a trend for BA.2.75-infected hamsters to show an increase in the levels of infiltrating macrophages in the lungs compared to BA.2, although differences were not statistically significant ( During the revision of this manuscript, the spread of EG.5.1 around the world substantially increased [9].We compared the virulence of EG.5.1 to BA.2.75, as the latter showed relatively higher virulence in hamsters compared to other Omicron-like viruses.We found the virulence of EG.5.1 to be significantly higher than BA.1 but comparable to BA.2.75 (Fig 6I -6L).
For comparative purposes, semiquantitative scoring [39] of the lung pathology was also performed simultaneously by three board certified pathologists on the lungs of hamsters infected with BA.1, XBB, BQ.1.18and BA.2.75.Overall data showed a similar trend compared to our digital pathology pipelines (S8 Fig).

Discussion
The COVID-19 pandemic has entered a phase characterised by the periodic emergence of immune escape variants, which may differ in their virulence from their predecessors.Assessing the virulence of any newly emerged variant will be therefore one of the key features requiring near real-time monitoring.Animal models have been used throughout the pandemic to unveil many aspects of SARS-CoV-2 pathogenesis [43].So far, there has been an excellent correlation between the virulence of SARS-CoV-2 variants such as Delta and BA.1 in humans and in experimentally infected hamsters [7, 18-21, 26, 27, 42, 46, 47].
Histopathological features in infected animals such as alveolar damage and inflammatory lesions in the respiratory tract are a consequence of virus replication (and host immune responses) and therefore directly reflect virus virulence.In various studies, lung pathology is often characterised by qualitative scores determined by trained pathologists on lesions such as bronchiolitis, haemorrhages, alveolar damage and others (S9 Fig) .Qualitative histopathology scores can vary between individuals, and therefore these types of data are difficult to compare unequivocally between different laboratories.Implementation of digital-pathology pipelines assisted by machine learning not only increase objectivity by eliminating inter-and intraobserver variability, but it also increases the speed of data analysis [49].In addition, assessing lesions in a semi-quantitative fashion with a limited number of options (i.e.scoring with 0, 1, 2, 3) may mask the variability between individuals that can be better appreciated with a more quantitative method used in this study.In this study, we aimed to develop a framework for quantitative unbiased methods to phenotype the relative virulence of SARS-CoV-2 variants by quantifying the extent of pulmonary pathology in experimentally infected hamsters.In line with previous studies [39,[50][51][52], our transcriptomic analysis showed that SARS-CoV-2 infection induces an infiltrate of immune cells in the lungs and immune activation in general.Indeed, here we found that expression of ISGs, and infiltrate of T cells and macrophages are common in lesions induced by the hypervirulent variant Delta but are barely above background in BA.1-infected hamsters.
In addition, our RNAseq data suggest that tissue remodelling was another key parameter distinguishing the transcriptional profile of Delta-and BA.1-infected hamsters.Virulent SARS-CoV-2 variants induce alveolar damage, with necrosis of type 1 and type 2 pneumocytes.Lung damage is subsequently repaired by alveolar epithelial hyperplasia, which is a key feature of lesions induced by SARS-CoV-2 in experimentally infected hamsters [30,42,43], and it is also found in some cases in post-mortem samples of human patients dying as result of COVID-19 [53][54][55].Lung respiratory epithelium repair is due to proliferation of either type 2 pneumocytes following mild injuries [56], or other lung progenitor cells in response to severe injury with abundant loss of type 1 pneumocytes [57][58][59].
Overall, we found that infiltrating macrophages and proliferating epithelial cells constitute most of the cells in the affected areas of the lung parenchyma of hamsters infected with virulent SARS-CoV-2 variants.We consistently observed epithelial hyperplasia and macrophage infiltrates also in a separate study based on histopathology of the lungs of patients who died as result of the first wave of COVID-19 [60].These two features represent therefore exceedingly good markers that by themselves provide an unbiased quantification of the virus-induced pulmonary lesions, and by extension virus virulence.In addition, the use of whole lung sections can also reduce bias by providing a spatial overview of the inflammatory response.For example, our investigation of the spatial distribution of the ISG response and how it directly correlated to the presence of virus in an inflammatory lesion proved to be particularly insightful (Fig 3).Whole scanned imaging of tissue sections and downstream analyses including those based on artificial intelligence approaches ("digital pathology") have been used extensively in diagnostic pathology and other fields including cancer research and infectious diseases [49,61,62].Artificial intelligence-based pipelines to diagnose lung pathology in humans have emerged recently [63].In these approaches, the software is trained to distinguish between different diseases such as tuberculosis and lung cancer as well as fibrotic conditions.
In our study, by performing immunohistochemistry of whole scanned lungs sections to identify IBA1 + cell (as convenient marker for macrophages), we were able to quantify infiltrating macrophages in the lung parenchyma of experimentally infected hamsters.In this manner, the relative number of infiltrating macrophages in the lung parenchyma can be acquired in an unbiased fashion.The same approach to measure alveolar epithelial hyperplasia by simply quantifying TTF1 + cells did not initially provide satisfactory results, due to the expression of this marker in both proliferating and non-proliferating type 2 pneumocytes, as well as bronchial epithelial cells.However, supervised machine learning approaches allowed us to train the software used to detect clusters of TTF1 + cells (hyperplastic type-2 pneumocytes) and exclude isolated type 2 pneumocytes (representing the normal type 2 cells in the lungs) or TTF1 + cells in the bronchiolar epithelium.Our approach allowed us not only to clearly distinguish the pulmonary lesions between those induced by the virulent Delta and the attenuated BA.1, but also to show an increased virulence of other Omicron sub-lineages.
The pipeline used in this study, enabled us to show that BA.5, BA.2.75 and EG.5.1 have acquired an increased virulence phenotype compared to BA.1, as also suggested in other recent studies [6,23,27].We saw no differences instead in virulence between XBB and BA.1, also in keeping with another recent study [8].Interestingly we found a tendency for BA.2 to induce a higher number of macrophage infiltrates and type 2 pneumocyte hyperplasia compared to BA.1, although differences did not reach statistical significance.Other published studies found no major differences in virulence between BA.1 and BA.2 [8,48], in contrast to a previous study which used recombinant viruses with either BA.1 or BA.2 spike [64].
We propose that the approach described in this study can be used to quantify moderate differences in virulence between variants, although animal group sizes may need to be adjusted to address specific experimental questions.Indeed, most studies focusing on SARS-CoV-2 virulence in hamsters use experimental groups between 4 and 6 animals (and often males only), but these may not be sufficient to reveal relatively modest differences in virulence between variants and the inherent individual variability (considering also that scores are set with values between 0 and 3 or 4).It may be argued that modest differences in variant virulence observed in hamsters could have limited biological significance in human patients.Indeed, a limitation of these type of studies is that the intrinsic virulence of variants in naïve hamsters is difficult to compare to their virulence in the "real world" in the human population (with pre-existing immunity derived from vaccination or previous infections).However, the intrinsic virulence of newly emerging variants remains a key phenotype to monitor to determine whether adjustments to public health measures are needed.For example, the emergence of a variant with increased virulence may require different vaccination policies from the existing ones.
We use equal number of both male and female hamsters in all our experiments.As expected, we observed some gender-independent variability in the extent of the lesions caused by the same variant between individual hamsters both within and between different experiments.A standard reference virus could be used in multiple experiments assessing different variants to normalise the data between experiments (S7 Fig) .This approach would also enable meaningful comparison of data between different laboratories.
In conclusion, the digital pathology pipelines developed in this study provide a framework for the quantitative assessment of virulence of SARS-CoV-2 variants in experimentally infected hamsters.The identification of virus standards with defined pathogenicity and sharing of protocols for the quantification of pulmonary lesions, can allow comprehensive comparison of in vivo data between different laboratories.We also developed an online repository to share more effectively with the research community histopathological images derived from scanned whole organ tissue sections.We argue that the reader will be able to appreciate better histopathological micrographs contained in the main body of manuscripts if these were supplemented by images of whole scanned slides similarly to those contained in our virtual microscope.Due to space constraints, micrographs contained in figures of standard manuscripts show only individual areas of interest at single or two magnifications.Virtual microscopes allow the reader to view all areas of a given section and at multiple magnifications, providing therefore a comprehensive spatial context of the data.

Ethics statement
Animals were maintained at the University of Glasgow under specific pathogen free conditions in accordance with UK Home Office regulations (Project License PP0271643), in accordance with the Animals (Scientific Procedures) Act 1986, and approved by the University of Glasgow ethics committee.All animal research adhered to ARRIVE guidelines.

Cells and viruses
Calu-3 cells (ATCC HTB-55) are human lung adenocarcinoma epithelial cells.African green monkey kidney cells (Vero E6) expressing the human Ace2 receptor [65] were used to propagate the viruses.All cell lines were maintained at 37˚C and 5% CO 2 in DMEM (ThermoFisher) supplemented with 10% foetal bovine serum (FBS) (ThermoFisher), except for Calu-3 cells where RPMI-1640 medium (ThermoFisher) was supplemented with 20% FBS.BA.5 virus isolate was obtained as passage P2 from Greg Towers (University College London, UK) after initial access to passage P1 obtained from Alex Sigal (AHRI, South Africa).Other variants used in this study included the Delta variant (B.1.617.2,GISAID accession number EPI_ISL_1731019), BA.1, BA.2, BA.2.75, BQ.1.18,and XBB all isolated from clinical samples obtained either at the CVR or at the Imperial College London.Viruses were isolated from a clinical nasopharyngeal swab sample collected in virus transport medium.Samples were then resuspended in serumfree DMEM supplemented with 100 units mL-1 penicillin-streptomycin, 10 ug mL-1 gentamicin and 2.5 ug mL-1 amphotericin B to a final volume of 1.5 mL.Calu-3 cells were then inoculated and incubated overnight at 37C, 5% CO2.The next day, cell culture medium was replaced with fresh complete medium before further incubation.Infected cells were monitored for signs of CPE and the presence of viral progeny in supernatant by RT-qPCR.All working stocks were propagated in VeroE6 cells.

Animals
Male and female golden Syrian hamsters (HsdHan:AURA) were bred and maintained by Envigo; Wyton, United Kingdom).Animals were shipped as required and acclimatised at the Veterinary Research Facility (VRF) at the University of Glasgow prior to transfer to Containment Level 3 (CL3) suite.Animals were maintained in individually ventilated cages (IVCs) on a 12-hour light/dark cycle and provided with food and water ad libitum.Both male and female hamsters between 8-12 weeks old were used for experiments.

Experimental infections
Animals were randomised to treatment groups, however, due to the nature of this work blinding was not possible.Animals were handled within a Class I microbiological safety cabinet (MSC) in a CL3 suite throughout the experiment.Hamsters were anaesthetised with oxygen (1.5 L/minute) containing 5% isoflurane and intranasally dosed with 50 μl of Dulbecco's Modified Eagle Medium (DMEM) containing 3.75x10 6 genome copies (the equivalent of approximately 1x10 4 PFU) of SARS-CoV-2.We chose to use genome copies for our infection studies to control for the differences in infectivity and plaque formation that have been observed with the BA.1 variant [16].Control animals received DMEM only.Animal weights and temperatures were recorded daily, and animals were monitored twice daily for signs of clinical disease including piloerection, hunching, abnormal breathing and reduced peer interactions.Animals received a disease score based upon weight loss and the presence and severity of these signs.Animals were culled at the end of the experiment via a rising concentration of CO 2 .

Throat swabs
Animals were restrained and swabs (MWE) were inserted into the mouth and rotated five times on each tonsil.The swabs were then placed in 2 ml DMEM (ThermoFisher) containing 2% FBS (ThermoFisher) for 1 minute.For RNA RT-qPCR, 250 μl of media was added to 750 μl TRIzol LS (ThermoFisher).For live virus assays, the media was frozen at -80˚C.

RNA extractions
Viral RNA was extracted from culture supernatants using the RNAdvance blood kit (Beckman Coulter Life Sciences) following the manufacturer's instructions.Tissue samples, approximately 20 mg in size, were collected in homogenisation tubes containing 2.8 mm metal beads (Stretton Scientific) and 1 ml TRIzol reagent (ThermoFisher Scientific).Blood was collected in TRIzol LS reagent (ThermoFisher Scientific) at a 1:4 ratio.Tissue samples were homogenised (6500 rpm, 4x 30s cycle, 30s break, room temperature) using a Precellys Evolution Homogeniser (Bertin Instruments) and the homogenate was mixed with 200 μl chloroform (Merck-1L).Each sample was mixed thoroughly and centrifuged at 12,000 x g for 15 minutes at 4˚C.The aqueous phase was removed and mixed with 230 μl 100% ethanol (Merck).The samples were then added to columns from RNeasy Mini Kits (Qiagen) and processed as per manufacturer's instructions.On column DNase digestion (Qiagen) was performed on all samples for 15 minutes as per manufacturer's instructions.

Digital pathology: Standard signal quantification
For image analysis, slides were scanned with an Aperio VERSA 8 Brightfield, Fluorescence & FISH Digital Pathology Scanner (Leica Biosystems) at 200 x brightfield magnification [66].
To quantify the positive signal of Fast-red (in situ-hybridisation, ISH) and diamino-benzidine (DAB)-brown (immunohistochemistry, IHC) in the tissues (nose, turbinates, larynx, trachea, and lung), we optimised image analysis pipelines based on the signal detected using the QuPath software (version 0.3.2 or later).The initial step was to upload and set the appropriate image type [Brightfield (other) for ISH; Brightfield (H-DAB) for IHC] and manually outline the area of interest using the 'wand and 'brush' tools.Staining obtained for SARS CoV-2 spike RNA (ISH), IBA1 (IHC), MX1 (IHC) as well as various ISGs including IFIT1, MX1, RSAD2 and OAS1 (ISH), we utilised the pixel classifier feature within QuPAth.This specific (pixel classifier) was selected for analysis as, in these stainings, the virus (spike RNA) signal was very abundant, with multiple positive cells being clustered together or in patches.Conversely, for the ISGs, the positive signal was dispersedly distributed with a punctuate pattern.In detail, for each Pixel detection (%), the image type was set as described above.To create the thresholder, the 'Resolution' was set to 1.09 μm/pixel or higher.For the 'Channel', we used 'Red' (ISH) or 'DAB' (IHC); the 'Prefilter' was always 'Gaussian' while 'Smoothing sigma' as well as the 'Threshold' have been tuned for each set of experiments to optimise the correct detection (based on each individual batch of slides with their individual staining intensities and background).The 'Above threshold' option was always set to 'Positive' whereas the 'Below threshold' was always set to 'Negative'.The 'Region' to be analysed was set to 'Any annotation ROI'.The readout is a percentage of positive pixels per total pixels in the annotated area.
The 'Positive cell' detection feature has been used instead for the detection of cell membrane-associated CD3-positive cells.The settings for the 'Requested pixel size' were adjusted to 0.5 μm or smaller and the 'Nucleus parameters' included a 'Background radius' of 8 μm.The 'Median filter radius' was set to 0.8 μm or smaller, 'Sigma' to 1.5 μm or smaller while the 'Minimum' and 'Maximum' areas were adjusted to 10-30 μm 2 and 400 μm 2 respectively.The 'Intensity parameters' included a 'Threshold' of 0.1 or 0.2 and a 'Max.background intensity' of 2.
The 'Cell parameters' included a 'Cell expansion' of 5 μm.The readout is a percentage of positive cells detected per all cells in the annotated area.

Digital pathology: Quantification of alveolar epithelial hyperplasia
We used the algorithm 'AI' within the software HALO (Indica Labs) to detect clusters of TTF1 + nuclei within the lung corresponding to alveolar epithelial hyperplasia, while ignoring individual TTF1 + normal type 2 pneumocytes, and the positively stained normal or hyperplastic bronchial epithelium.The readout is therefore the percentage of positive TTF1-stained area in the lungs with rosette-like shapes.First, the software was trained to detect positive areas of interest.The files with the trained algorithm to detect alveolar epithelial hyperplasia are shared in the Enlighten repository (https://doi.org/10.5525/gla.researchdata.1513) and can be imported to any HALO software (HALO version 3.4.2986.230)by the users.The file ".classifier" can be read by a text editor and includes metadata (resolution, class names, etc.).The ". params" and ".trainer" files are specific to the deep learning framework mxnet.The.params files MXNet's specific file format for storing model weights.The.trainer files are MXNet's specific file format for the optimizer.This keeps track of training information such as iteration, learning rate and optimizer states.A zip file with these three files can be loaded directly into HALO.
The HALO algorithm (Densenet AI V2 plugin) was trained on 3 sets of independent experiments.Each group of animals were handled as separate batches and infected, processed, stained and scanned in 3 different times as independent experiments.The samples were processed in the same way across multiple experiments using the same machine learning algorithm.We selected at least 30-40 slides across 3 different experiments containing delta and omicron-infected lung samples to train the software and optimize its performance.First, regions of interest (ROIs) were selected on each slide.The DenseNet AI plugin classifier was used, and classes were defined for positive (TTF-1 positive and proliferating type 2 pneumocytes) and negative tissue (isolated type 2 pneumocytes, in addition to normal or hyperplastic bronchial epithelium).Example regions were drawn within each class at high magnifications, and any significant TTF1-negative cells within the rosette structures were excluded.After example regions of epithelial hyperplasia were drawn, the classifier was trained over multiple days.Raw data is expressed as the total area covered by the rosette structures (TTF-1 positive nuclei, alveolar epithelial hyperplasia), in comparison to the overall tissue area.Slides with artefacts were excluded from the analysis as well as artefacts on the slides, which have been manually excluded.Additional examples of the performance of the software algorithm developed in this study to detect alveolar epithelial hyperplasia are shown in S6 Fig. Additionally, we also compared data related to quantification of alveolar epithelial hyperplasia from two different lung sections (approximately 100μm apart) from each animal experimentally infected with either BA.1, BA.2 or BA.2.75 (S9 Fig) The data obtained with either section (termed for simplicity "section 1" and "section 2") were very similar.Therefore, we used a single section per animal in subsequent experiments.
All photomicrographs have been captured with Aperio ImageScope (Leica Biosystems).Accuracy of the training outputs was confirmed by a Board certified veterinary pathologist.For some experiments, the data obtained with the digital pathology pipeline was further compared to semiquantitative scoring of HE-stained lung sections by three board-certified pathologists in a blinded fashion according to previously published studies [39].Briefly, lesions (amount of lung affected, respiratory epithelial cell or type II pneumocyte hyperplasia, alveolar damage, perivascular inflammation, damage in bronchi/bronchioles) were scored from 0 to 4 for each lung (n = 6 hamsters) by each pathologist (n = 3).Data are also shown as a cumulative PLOS PATHOGENS score containing all 5 parameters from all pathologists.In this case, results are shown as the sum of the average value scored by each pathologist for each animal.

Online digital pathology tool
Representative whole scanned images are available online at https://covid-atlas.cvr.gla.ac.uk.The "CVR Virtual Microscope," is an online tool where users can zoom in and out of digital images of scanned tissues, accessing therefore the same context and information experienced on a microscope.

RNA sequencing (Bulk RNAseq)
Total RNA extracted from lung and blood samples from uninfected and infected animals was quantified using Qubit Fluorometer 4 (Life Technologies; Q33238), Qubit RNA HS Assay (Life Technologies; Q32855) and dsDNA HS Assay Kits (Life Technologies; Q32854).RNA integrity number was determined on a 4200 TapeStation System (Agilent Technologies; G2991A) using a High Sensitivity RNA Screen Tape assay (Agilent Technologies; 5067-5579).Before library preparation, haemoglobin RNA was removed from 1 μg of RNA extracted from blood samples using the GLOBINclear-Mouse/Rat Globin mRNA Removal Kit (Thermo Fisher Scientific; AM1981) following the manufacturer's instructions.Total RNA (500 ng) was used to prepare libraries for sequencing using the Illumina TruSeq Stranded mRNA Library Prep kit (Illumina; 20020594) and SuperScript II Reverse Transcriptase (Invitrogen; 18064014) according to the manufacturer's instructions.The PCR amplified dual indexed libraries were cleaned up with Agencourt AMPure XP magnetic beads (Beckman Coulter; A63881), quantified using Qubit Fluorometer 4 (Life Technologies; Q33238) and Qubit dsDNA HS Assay Kit (Life Technologies; Q32854).Their size distribution was assessed using a 4200 TapeStation System (Agilent Technologies; G2991A) with a High Sensitivity D1000 Screen Tape assay (Agilent Technologies; 5067-5584).Libraries were pooled in equimolar concentrations and sequenced using high output cartridges with 75 cycles (Illumina; 20024911) on an Illumina NextSeq 550 sequencer (Illumina; SY-415-1002).A Q score of �30 was presented in at least 95% of the sequencing reads generated.

Differential expression genes analysis
After the alignment to the Mesocricetus auratus genome, FeatureCount [68] was used to calculate the mapped reads counts.In this paper, we observed the differential expression genes (DEGs) of mock vs Delta and mock vs BA.1 (2/6 days) on both lung and blood cells.The DESeq2 [69] in Generalized linear models (GLMs) with multi-factor designs (here the factors are gender and condition of the samples) was used for differential expression genes analysis.FDR P-value < 0.05 was used as the cut-off of significant differential expression genes.We analysed the differential expressed gene sets corresponding to molecular pathways of the Reactome database [70].

Statistical analysis
All graphs and statistical analyses were produced using GraphPad Prism 7 (GraphPad Software Inc., San Diego, CA, USA) as indicated in each figure legend.P values < 0.05 were deemed to be significant.In the graphs, values of male animals are displayed as triangles, while circles are used for females.

Fig 1 .
Fig 1. Transcriptomic response of lungs and blood of SARS-CoV-2 infected hamsters.(A) Volcano plots indicate the number of significantly upregulated (red arrows) or downregulated (blue arrows) genes in lungs or blood of hamsters infected with the indicated variant at each timepoint (compared to mock-infected controls).(B-C) Pathway analysis highlighting the enriched ontologies of differentially expressed genes in the lungs (B) and blood (C) of infected hamsters at 2-and 6-days postinfection (dpi).The number of genes involved are indicated beside each pathway.(D-E) Scatterplots representing the relative expression of interferon stimulated genes (ISG) in lungs (D) and blood (E) of hamsters experimentally infected with either the Delta or BA.1 variant at the indicated timepoints.ISG with relative higher expression in BA.1 compared to Delta (FDR<0.05)are shown in red.ISG with relative higher expression in Delta compared to BA.1 (FDR<0.05)are shown in blue.Data derived from n = 8 hamsters (4 females and 4 males per group), infected in two independent experiments.https://doi.org/10.1371/journal.ppat.1011589.g001

Fig 4 .
Fig 4. Lung histopathology of SARS-CoV-2 experimentally infected hamsters.Low (left) and high (right) magnification of lung sections stained with haematoxylin and eosin of hamsters experimentally infected with either the Delta or BA.1 (3.75x10 6 genome equivalent copies per animal), or mockinfected controls.At 2 days post-infection (dpi), Delta causes a higher level of infiltration of macrophages in and around bronchi, while BA.1 causes only mild infiltrations (white arrow) compared to mock-treated hamsters.On the same day, Delta-infected animals showed a moderate vasculitis (inset, right panel) and a moderate sloughing of bronchial epithelial cells (asterisk).Vascular and bronchi pathology was minimal in BA.1-infected animals at 2 dpi.At 6 dpi the lungs of Delta-infected animals show a severe infiltration of macrophages, lymphocytes, plasma cells and neutrophils/heterophils as well as a severe alveolar epithelial hyperplasia covering large areas of the lung lobe (black arrows; see also S4 Fig).BA.1-infected animals show in some cases the same lesions, including type 2 pneumocytes hyperplasia (asterisk) but covering only limited amounts of the lung lobes (black arrows).(Empty scale bars: 2 mm; filled scaled bars: 100 μm).https://doi.org/10.1371/journal.ppat.1011589.g004

Fig 5 .
Fig 5. Quantification of tissue histopathology in Delta or BA.1-experimentally infected hamsters.Images of whole lung sections of hamsters experimentally infected with either Delta or BA.1 (3.75x10 6 genome equivalent copies per animal) or mock-infected controls, and culled at either 2-or 6-days post-infection (dpi).Expression of CD3 (A), IBA1 (B) and TTF1 (C-D) was assessed by immunohistochemistry. Animals were culled 2-or 6-days post infection (dpi) and the lungs were collected for histological analysis.(A-C) Representative photomicrographs for each marker and associated quantification are shown.(D) Mock-infected animals show a positive signal for TTF1 only in bronchial epithelium and in single cells scattered throughout the lung parenchyma.(E) Lungs of Delta-infected hamsters shows severe proliferation of TTF1 + positive cells (arrows).The inset shows a higher magnification of type 2 pneumocytes hyperplasia surrounding centrally located bronchi (asterisk).(F) Same image as in E after analysis with the HALO software algorithm trained to detect alveolar epithelial hyperplasia.The software shows negative areas in green and positive ones in yellow; only proliferating alveolar epithelial cells are identified as positive areas, while TTF1 + type-2 pneumocytes and bronchial epithelial cells are excluded by the software.CD3-and IBA1-positive cells in experimentally infected animals were quantified using QuPath (Version 0.3.2 or newer), while HALO was used to quantify alveolar epithelial hyperplasia.Statistical analysis was performed using a Two-Way ANOVA, *<0.05, **<0.01,***<0.001,****<0.0001.Data were obtained from two independent experiments using 8 hamsters (4 females and 4 males infected in two independent experiments).Black: uninfected; red: Delta-infected; blue: BA.1-infected (males, triangles; females: circles).Blue scale bar: 1 mm; green scale bar: 400 μm.https://doi.org/10.1371/journal.ppat.1011589.g005

Fig 6 .
Fig 6.Quantification of Omicron sub-lineages virulence in experimentally infected hamsters.(A) Photomicrographs of lung sections collected at 6 days post-infection (dpi) from hamsters experimentally infected with either Delta, BA.1, BA.5 (virus load equivalent to 3.75x10 6 genome copies per hamster) or mock-infected.Expression of TTF1 + hyperplastic epithelial cells and infiltrating macrophages (IBA1 + ) was assessed by in situ-hybridisation.Expression of MX1 was instead assessed by immunohistochemistry. Data were quantified as in Fig 5. (B) Differences in the weights between these animals were recorded daily.(C) Daily recorded weight changes in animals mock- Fig 6G).All variants induced limited levels of alveolar epithelial hyperplasia (Fig 6H), in line with the limited ability of these viruses to invade the lung parenchyma.Finally, to visualise and compare all the different variants used in this study, we normalised the data using values obtained in BA.1-infected hamsters as unit of reference between different experiments (S7 Fig).Graphs displayed in S7 Fig show, as expected, a gradient of virulence between variants.XBB induced little or no lung pathology similarly to BA.1.BA.5, BA.2 and BQ.1.18were more virulent than BA.1, while BA.2.75 was clearly more virulent than BA.1 but not as virulent as Delta.The use of a standard virus of reference could therefore enable comparisons between variants used in different experiments (and different laboratories).