A Deeper Insight into Evolutionary Patterns and Phylogenetic History of ORF Virus through the Whole Genome Sequencing of the First Italian Strains

Orf virus (ORFV) is distributed worldwide and is the causative agent of contagious ecthyma that mainly occurs in sheep and goats. This disease was reported for the first time at the end of 18th century in Europe but very little is currently known about the temporal and geographic origins of this virus. In the present study, the use of new Italian whole genomes allowed for better inference on the evolutionary history of ORFV. In accordance with previous studies, two genome types (S and G) were described for infection of sheep and goats, respectively. These two well-differentiated groups of genomes originated for evolutive convergence in the late 1800s in two different areas of the world (Europe for S type and Asia for G type), but it was only in the early 1900s that the effective size of ORFV increased among hosts and the virus spread across the whole European continent. The Italian strains which were sequenced in the present study were isolated on the Mediterranean island of Sardinian and showed to be exclusive to this geographic area. One of them is likely representative of the early European forms of ORFV which infected sheep and became extinct about one century ago. Such an ancient Sardinian strain may have reached the island simple by chance, where it quickly adapted to the new habitat.


Introduction
Orf virus (ORFV; family: poxviridae) is the etiological agent of the zoonotic disease contagious ecthyma (CE or ORF), also known as contagious pustular dermatitis (CPD). This virus belongs to the genus Parapoxvirus (PPV) [1] that, to date, includes four further recognized PPV species based on the classification of the International Committee on Taxonomy of Viruses (ICTV): bovine papular stomatitis virus (BPSV), grey sealpox virus, pseudocowpox virus (PCPV), and red deerpox virus (RDPV).
Contagious ecthyma was reported for the first time in sheep by Steeb in 1787 [2], and as early as the end of the 19th century, was described by Hansen (1879) [3] in goats and humans as contagious dermatitis, ORF, crusta labialis, carbuncle of the coronary band, or coupar angus [4]. It was in 1929 that Howarth [5] showed that the disease in lambs and kids in California (USA) was due to a virus which was better described for specimens from Texas (USA) by Boughton and Hardy, in 1932 [6]. The patterns of lesions connected to ORFV infection can occur in a clinical spectrum ranging from mild papular and pustolar to severe proliferative. These types of ORFV lesions are generally reported in sheep and goats, but similar kinds of lesions have also been found in other animal species [7] and humans [1,8], leading to CE being considered as a worldwide-diffused zoonotic disease [9]. In ruminants, the lesions that are provoked by ORFV infection usually involve mouth, muzzle, nostrils, gums and tongue, and occasionally feet, udders, and sporadically, gastrointestinal tract and respiratory apparatus [10]. This disease usually affects lambs and kids, thus causing their death when severe lesions avoid milk sucking from mothers. The level of morbidity from this disease is very high and, although mortality is generally rare, it can rate up to 10% in lambs [11] and to 93% in kids [12,13]. In humans, lesions are reported mainly on hands and undergo a spontaneous benign resolution in a few weeks. However, malignant cases have been also described with atypical proliferative lesions, particularly in immunocompromised individuals [7,14,15].
ORFV's genome is a linear double stranded DNA which is 135 kb-long and encodes for 132 genes [7,16]. The central region contains highly conserved genes involved in viral replication and in the assembly of the viral structure, while terminal regions are more variable and include genes for virulence and immunomodulation [7,17]. The occurrence of high levels of nucleotide variation for the terminal regions of ORFV may be connected to the phenomenon of reinfection typical of this virus that also appears as a product of the regulation acted by viral genes on the host innate immune response [14]. In such a context, many studies focused their attention on the immunomodulatory action of specific ORFV genes with the aim to shed light on the complex processes involved in viral pathogenesis [18,19].
Only two studies, which were performed to describe novel ORFV genomes isolated in China [20] and France [21], provided phylogenetic inferences based on ORFV whole genomes. The first study [20] was performed on nine genomes and evidenced the occurrence of a well-supported genetic structuring between the strains isolated from sheep and those from goats. Consistently, results obtained in the second study [21], where twelve genomes from almost all the continents were used, evidenced a genetic differentiation depending on whether the virus host was a goat or a sheep. This latter study also described, for the first time, an ORFV genome isolated from a human after infection from sheep. A genetic affinity between genomes from human infections was also found in this study. These authors suggested the importance to analyze a larger number of ORFV whole genomes to confirm the possible occurrence of two types of ORFV (from sheep vs. from goats) and shed new light of the level of genetic variation associated with the host species.
In such a context, in the present study, seven new ORFV whole genomes from the Mediterranean island of Sardinia (Italy) were isolated from sheep and goats and merged with all of the ORFV genomes currently available on public databases to perform a highresolution phylogenetic analysis of the virus.
A similar phylogenetic deep inference was previously performed for samples from all over the world based on the ORFV gene, encoding the dsRNA-binding protein (VIR) [22]. Results evidenced a high worldwide mutational viral evolutionary rate, along with a well-supported genetic divergence between the viral strains isolated from sheep and those from goats.
The aim of this study is to provide new insight into the evolutionary history of ORFV based on the analysis of whole genomes, to better understand the temporal origin of its strains with corresponding patterns of distribution, and to make inferences on the possible occurrence of two different genome types of ORFV associated with sheep and goat infections.

Sampling
Samples were collected, between May 2017 and March 2021, in five different Lacaune and Sarda breed sheep flocks and in one Sarda goat herd (see Table 1 for details). ORFV infection evidenced by clinical signs in individuals was confirmed by virological and molecular analyses of biological samples. Clinical CE-outbreaks have been detected and reported by private practitioners and veterinary public health components involved in the present research, during their diagnostic activity. Sampling collection protocols were as reported by Coradduzza et al., 2021 [22]. Viral DNA was extracted from lesions isolated in five infected sheep and two infected goats (see Table 1 for details on the samples).

Virus Isolation
ORFV was isolated on Vero cells from homogenates tissues of ovine and caprine mild and severe-type lesions according to the methodologies reported by Coradduzza et al., 2021 [22]. In particular, a total of 0.5 g of each tissue sample was homogenized in 5 mL (10% w/v) of DMEM medium with the use of the following antibiotics: 400 UI/mL penicillin, 400 µg/mL streptomycin, 300 µg/mL gentamicin, and 2.5 µg/mL amphotericin B. The suspension obtained was centrifuged at 1000× g for 15 min and used for the infection of Vero cells using cells from the 100th to 120th passage (BSCL86 ATCC, Istituto Zooprofilattico Sperimentale della Lombardia e dell'Emilia Romagna, Italy). Twelve-well plates of cells in DMEM were prepared with the addition of 10% (v/v) FBS and then incubated for 1 h in a thermostat at 37 • C in an atmosphere of 5% CO 2 . After 18-24 h (80-90% confluence), the medium was removed and the cells were incubated with 0.5 mL of tissue homogenate under the same incubation conditions. Furthermore, 3 wells without the virus, used as a negative control, along with 3 wells containing the virus, were also prepared; after incubation, the cells were washed three times with 1× PBS and the new DMEM with antibiotic, and fetal bovine serum was added (final concentration 100 UI/mL penicillin, 100 µg/mL streptomycin, and 0.5 µg/mL amphotericin B with 2% (v/v) FBS). The cytopathic effect (CPE) was checked daily, when it was detected, or on the 5th day of culture, the material was freeze-thawed three times, collected, and centrifuged at 200× g for 10 min. The collected supernatant was stored at 80 • for future use.

Viral DNA Extraction, Sequencing and Genome Assembly
After centrifugation at 10,000× g for 1 min at 4 • C, viral DNA was extracted from cell culture supernatant using a QIAmp UltraSens Virus Kit (Qiagen, Hilden, Germany), as described in Fiori et al. [23] to perform the whole genome sequencing using Illumina platform.

Phylogeny, Molecular Dating and Evolutionary Rate
Seven whole genome sequences were obtained for ORFV in the present study. The whole dataset (n = 26) was aligned using the Clustal Omega package [31] (available at https://www.ebi.ac.uk/Tools/msa/clustalo/; accessed on 15 April 2022) after editing by means of Unipro UGENE v.35 [29], and subsequently deposited in GenBank (see Table 2 for accession numbers). These Sardinian sequences were included in a large dataset containing all the ORFV whole genomes available on GenBank to date (see Table 2 for details and references). This dataset included a total of 26 sequences; 7 from the present study (Italy), 2 from Germany, 1 from France, 10 from China, 1 from India, 1 from New Zealand, and 4 from United States of America (USA) (see Table 2 for details).
To identify potential genetic clusters within the dataset and to determine the dissimilarity represented by the genetic variability among genomes, a PCoA (principal co-ordinate analysis) was performed using GenAlEX 6.5 [38]. PCoA reconstruction was based on a genetic pairwise p-distance matrix.
The jModeltest 2.1.1 software [39] was used to find the best probabilistic model of genome evolution with a maximum likelihood optimized search. MrBayes 3.2.7 [40] was used to carry out a Bayesian phylogenetic analysis, setting nst = 6, rates = invgamma, and ngammacat = 4 as model parameters. Two independent runs, with 4 Metropoliscoupled Markov chain Monte Carlo (MCMC) chains (1 cold and 3 heated chains), were run synchronously for 5 million generations, sampling trees every 1000 generations. The first 25% of sampled trees were discarded for burn-in. The average standard deviation of split frequencies (which should approach 0) was verified checking the convergence of chains [40] and the potential scale reduction factor (which should be approximately 1) [41] following Scarpa et al. [42].
When available, the collection date of the samples (at least month and year were necessary) were used to set molecular dating which was performed by means of a Bayesian approach under the MCMC algorithm as implemented in Beast 1.10.4 software [43]. Only the genomes whose sample collection date was known were included in this analysis. Both strict and uncorrelated log-normal relaxed clock models were tested with fast runs of 100 million generations to identify the best clock model to perform the dating analyses. Selection was made comparing Bayes factor values using Tracer 1.7 software [44]. All available demographic models (both parametric and nonparametric) were also tested. The phylogenetic time-scaled (ultrametric) trees and the evolutionary rates were co-estimatedafter the selection of the Bayesian skyline demographic model under the uncorrelated log-normal relaxed clock model-running 500 million generations, with samplings repeated every 50,000 generations. The resulting log files were inspected using Tracer 1.7 software [44]. Only values of ESS (effective sample size) ≥200 were accepted. The maximum clade credibility tree was drawn, visualized, and annotated by means of the software TreeAnnotator (Beast package) and FigTree 1.4.1, respectively. Beast software was also used to perform further runs under the coalescent Bayesian skyline demographic model [45] to estimate the evolutionary rates for ORFV strains isolated in sheep and goats. All phylogenetic/phylodinamic runs were carried out using the CIPRES Science Gateway (Cyberinfrastructure for Phylogenetic Research available at https://www.phylo.org/index.php/; accessed on 8 May 2022) [46] with a total effort of about 5000 CPU hours of computation.

Results
Assembly performed for the Italian ORFV strains with the reference genome (NC_005336) provided seven new whole genomes. See Tables A1 and A2, provided in Appendix A, for details regarding each new genome.

Phylogeny, Molecular Dating, and Evolutionary Rate
In the present study, a dataset including 26 ORFV whole genomes that were 137,151 bp long, was analyzed. The sequences were isolated in all continents with the only exception being Africa (see Table 2 for details). The nucleotide frequency analysis carried out on the whole dataset reveals that conserved sites amount to 118,671 (86.5%), variable sites amount to 18,478 (13.5%), parsimony informative sites amount to 11,801 (8.6%), and singletons amount to 6677 (4.9%).
The phylogenetic mid-point tree analysis evidenced the occurrence of two fullysupported genetic clusters (cluster S and cluster G in Figure 1). The Bayesian phylogenetic time-scaled maximum clade credibility tree (ultrametric tree), used for the molecular dating, is provided in Appendix A (as Figure A1) with the confidence interval (C.I.) at 95% of the highest posterior density (HPD) indicated for each coalescence time estimate. Strains whose phylogenetic position was puzzling in the phylogenetic tree (GenBank # HM133903, MN648218) were excluded from this latter analysis. The phylogenetic mid-point tree analysis evidenced the occurrence of two fullysupported genetic clusters (cluster S and cluster G in Figure 1). The Bayesian phylogenetic time-scaled maximum clade credibility tree (ultrametric tree), used for the molecular dating, is provided in Appendix A (as Figure A1) with the confidence interval (C.I.) at 95% of the highest posterior density (HPD) indicated for each coalescence time estimate. Strains whose phylogenetic position was puzzling in the phylogenetic tree (GenBank # HM133903, MN648218) were excluded from this latter analysis.  Furthermore, only genomes whose sample collection date was available were used to perform the molecular dating (two genomes were removed from the dataset used for molecular dating: GenBank # HM133903, DQ184476).
Results obtained from the ultrametric tree ( Figure A1) are consistent with those obtained from the phylogenetic tree (Figure 1), with a few discrepancies due to the lack of three genomes in the dataset analyzed, as reported above. The coalescence time at each node of the phylogenetic tree was inferred according to the molecular dating estimates obtained with the ultrametric tree.
The common ancestor of all the ORFV strains present in the analyzed dataset, which corresponds to the root of the tree, dates back to 257.59 years before 2021 (i.e., the end of the year 1763). The cluster G of the tree includes only ORFV genomes isolated from goats with a single exception. Indeed, within this group of sequences, a genome isolated in a sheep from Germany (GenBank # HM133903) is also present and sets in a basal position, as an ancestral and quite divergent strain, outside the ingroup of genomes from Italy, India, China, and USA. For this strain, the date of collection is unknown (probably before year 2010), and for this reason it was not possible to estimate its divergence time (indeed, this genome was not included within the dataset used to construct the ultrametric tree, see Figure A1). Within the ingroup, which dates back to 255.93 years before 2021 (i.e., the middle of the year 1765), the genomes from Sardinia (Italy) (two kids) that were isolated in 2019 and 2020, set as an external sub-cluster that originates 10.39 years before 2021 (i.e., the final part of the year 2010). The other genomes of cluster G were grouped in a sub-cluster The second sister clade of cluster S dates back to 70.62 years before 2021 (i.e., the end of the year 1950) and includes only sequences from Sardinia (Italy) collected between 2017 and 2020. In particular, within this Sardinian clade, a sequence isolated in an adult from the northern part of the island in 2020 sets in a basal position, while all the other strains (also isolated in the northern of the island), collected from a lamb and two adults, between the 2017 and 2020, grouped together in a sub-cluster that originates 57.62 years before 2021 (i.e., the end of the year 1963). Within this sub-cluster, the two sequences collected in 2019 and 2020 belong to a group that dates back to 19.19 years before 2021 (i.e., the end of the year 2001).
The PCoA (Figure 2) performed to evidence the occurrence of genetic clusters among the genomes analyzed, was able to explain 61.14% of the variability (PCoA1/X axis: 41.21%, PCoA2/Y axis: 19.93%). Results were consistent with those obtained from the phylogenetic tree analysis, thus evidencing the occurrence of divergence between ORFV strains isolated from sheep and humans and ORFV strains isolated from goats (cluster S and cluster G, respectively, in Figure 2). Furthermore, it is interesting to note that in the PCoA graphic, the genome isolated in a sheep from Germany, which sets in a basal position of cluster G (goats) in the phylogenetic tree, is included as an outlier within the variability of strains isolated from sheep. In accordance with the phylogenetic tree analysis, the genomes from Sardinia (Italy) isolated in goats, set as outliers within the variability of strains isolated from goats. The two genomes from China isolated from a sheep and a goat, which are included as a highly divergent group within cluster S of phylogenetic tree, could be considered as divergent strains outside the variability of the two groups of genomes evidenced by PCoA.
analysis, the genomes from Sardinia (Italy) isolated in goats, set as outliers within the variability of strains isolated from goats. The two genomes from China isolated from a sheep and a goat, which are included as a highly divergent group within cluster S of phylogenetic tree, could be considered as divergent strains outside the variability of the two groups of genomes evidenced by PCoA. Under the Bayesian skyline lognormal uncorrelated relaxed clock model, the viral evolutionary mean rate calculated for the ORFV whole genome dataset was estimated to be 1.324 × 10 −6 substitution/site/year, with a C.I. 95% HPD of 7.488 × 10 −10 -5.439 × 10 −6 (1.324 × 10 −6 [7.488 × 10 −10 -5.4396 × 10 −6 ]). The viral evolutionary mean rate was also calculated for clusters S and G of phylogenetic tree and principal co-ordinates analyses; the strains whose sampling collection dates were unknown, along with those whose phylogenetic position was puzzling (GenBank # HM133903, MN648218), were excluded from the datasets. The two values were estimated to be 4.

Discussion
The present study allowed for a 37% increase in the number of whole genomes available for ORF virus from all over the world. In particular, we provided the first Italian sequences for this virus, thus expanding the number of European countries for which ORFV genomes are known. According to Chi et al. [20], who were the first to report a great genetic divergence between goat ORFVs and sheep ORFVs, and Andreani et al. [21], who suggested that more indepth research was needed to understand worldwide ORFV distribution, the present research expands the knowledge on the evolutionary history of ORFV, thus providing the first hints of the temporal and geographic origins of the virus and evidencing the occurrence of two highly-differentiated types of genomes. They are Under the Bayesian skyline lognormal uncorrelated relaxed clock model, the viral evolutionary mean rate calculated for the ORFV whole genome dataset was estimated to be 1.324 × 10 −6 substitution/site/year, with a C.I. 95% HPD of 7.488 × 10 −10 -5.439 × 10 −6 (1.324 × 10 −6 [7.488 × 10 −10 -5.4396 × 10 −6 ]). The viral evolutionary mean rate was also calculated for clusters S and G of phylogenetic tree and principal co-ordinates analyses; the strains whose sampling collection dates were unknown, along with those whose phyloge-

Discussion
The present study allowed for a 37% increase in the number of whole genomes available for ORF virus from all over the world. In particular, we provided the first Italian sequences for this virus, thus expanding the number of European countries for which ORFV genomes are known. According to Chi et al. [20], who were the first to report a great genetic divergence between goat ORFVs and sheep ORFVs, and Andreani et al. [21], who suggested that more indepth research was needed to understand worldwide ORFV distribution, the present research expands the knowledge on the evolutionary history of ORFV, thus providing the first hints of the temporal and geographic origins of the virus and evidencing the occurrence of two highly-differentiated types of genomes. They are the lineages of ORFV associated with the infection of goats (type G genomes) and sheep (type S genomes). The strong genetic divergence found between these two genome types could be related to the general lack of transmission of the ORFV infection between sheep and goats that would have prevented the occurrence of viral recombination and genetic similarity, among types S and G strains. Indeed, although this virus can be transmitted among hosts by direct contact with tissue lesions or fomites, thus infecting animals through skin cuts and abrasions, in general it is reported to be transmitted from wild to domestic goats, but not to domestic sheep [47]. Furthermore, in accordance with Chi et al. [20], who found that ORFV genomes isolated from goats are more similar to each other than those isolated from sheep, the evolutionary rate calculated for type G genomes in the present study was lower than that calculated for type S genomes, suggesting the occurrence of a stronger selective pressure acting on the viral strains that infect sheep along with a higher capability of type S genomes to quickly adapt to environmental changes. A possible explanation for this finding may involve the immune response of sheep that could be highly reactive towards ORFV as a consequence of extensive vaccination campaigns, thus prompting the evolutive rate of the virus with the occurrence of frequent viral recombination among type S genomes.
Furthermore, the evolutionary mean rate calculated for the ORFV type G genomes was similar to the value obtained in the present study when the strains from all hosts (goats, sheep and humans) were considered together. However, the evolutionary mean rate calculated for the ORFV type S genomes was higher and more similar to the value obtained by Coradduzza et al. [22] for the ORFV VIR gene which was isolated in hosts from all continents. This latter finding suggests a possible connection (that must be supported by further specific genetic studies) between the occurrence of contagious ecthyma in sheep and the expression of the gene encoding the dsRNA-binding protein in the viral strains that infect these animals.
In accordance with Andreani et al. [21], in the present study, different ORFV strains isolated after a human infection from sheep were closely related to one another and belonged to an external heterogeneous genetic sub-group within the type S genomes clade. Interestingly, considering that human-to-human transmission has never been reported [48], with only few exceptions which involve direct contacts with other human lesions or fomites [49], and according to Andreani et al. [20], who suggested the use of many genomes to evidence the occurrence of ORFV variants related to the strains isolated from humans that may be able to again cross species barriers; in the present study an ORFV strain (GenBank # MN454854) from a cell culture infected with contagious ecthyma vaccine in Texas (USA) in 2019 might represent a viral strain that can easily produce new spillover between sheep and humans. These events are frequent during the whole evolutive history of zoonoses and represent the main reason for which infections originating in wildlife require constant monitoring (see  and references therein) [50]. In such a context, although contagious ecthyma is considered as an occupational disease in "human risk populations", such as shepherds, sheep shearers, butchers, and veterinary surgeons, studies on a large number of ORFV genomes isolated from animals and humans are needed to better infer the evolutionary patterns occurring among types S lineages and to better describe the characteristics of the viral strains that infect sheep and can be easily transmitted to humans in specific conditions (i.e., direct or indirect contact with infected animals). In fact, considering that the ORFV strains isolated in farmers infected by sheep in the present study resulted in genetically divergence from those that are common among sheep, there is a strong possibility that only specific viral strains, with high levels of virulence, are able to infect humans even when environmental conditions make the possibility of ORFV transmission extremely low.
In the phylogenetic tree (Figure 1) of the present study, an ORFV strain isolated from a sheep in Germany (GenBank # HM133903) before 2010, which was considered as a puzzling outgroup genome by Andreani et al. [21], was included within the G (goats) cluster, and thus confirmed the peculiar nature of this viral allelic variant (as evidenced in the previous studies) that could be representative of an ancestral form of ORFV. Notably, this same viral strain in the PCoA (Figure 2) of our research clustered within the type S genomes isolated from sheep evidencing that this latter analytical approach could provide a more representative description of the genetic structuring among viral genomes. Indeed, PCoA, which was reconstructed on a p-distances matrix (without information on evolutionary models) is able to depict divergences or similarities based on differences in the genomic nucleotide composition of strains.
The common ancestor to all viral strains included in the present study (see also Figure 1 for molecular dating of the tree clusters), which could be considered as a proto-ORFV form not yet well-differentiated in the present types S and G ORFV, dates back to the end of 1763. This value represents a molecular dating which is consistent with the first description of contagious ecthyma in European sheep some years later by Steeb in 1787 [2], and in European goats and humans in 1879 [4]. In the present study, types S and G ORFV genome groups are almost coevals (1796 and 1765, respectively), and the rise of these two well-differentiated and contemporary genetic clades may represent a typical case of multiregional geographic origin by evolutive convergence of viral strains in different host species (indeed, identical pathological phenotypes/patterns of clinical lesions are reported for sheep and goats) (see e.g., Sackman et al., 2017 and references therein) [51]. In particular, regarding the type G strains, they may represent the first group of lineages that differentiated, thus suggesting that the early form of ORFV associated with the development of contagious ecthyma originated from viruses infecting goats. Their geographic origin could be tentatively placed in Asia, which currently represents the continent where the highest and the oldest (year 1831) level of genetic variability was evidenced in our analyses.
Furthermore, considering that contagious ecthyma was described for the first time in goats from European farms in 1879 [4], ORFV's type G strains may have taken about one century (from 1765 to 1879) to spread from Asia toward Europe, thus increasing its population size among hosts. This virus may have reached Europe, not only transported by live hosts involved in animal trade, but also with goods of animal origin along the commercial routes between the two continents. Indeed, ORFV remains viable on the wool of animals (infected and recovered) and contaminated materials for long periods, enabling indirect transmission to new susceptible hosts [52,53]. In general, this virus remains viable on the host's wool and survives for about one month after the lesions have healed, but it can also be carried by clinically normal sheep. Furthermore, ORFV was reported to survive in laboratory lesion tissue samples for up to 12 years [48], and for up to 17 years in natural environments with a dry climate [14] as the trade commercial routes between Asia and Europe likely were in the 1800s.
Within the type G clade, two sub-lineages are exclusive to the USA and Italy (Sardinia). In particular, the Italian allelic variants isolated in Sardinia may be considered as representative of a group of lineages exclusive to the island that originated in loco in 2011. Indeed, contagious ecthyma was first reported on this island in the early 1990s, according to a booklet edited by the Zooprophylactic Institute of Sardinia (Italy) as part of a health education project, and these recently-differentiated strains might have diverged from the first founders under the selective pressure of the new island habitat. The molecular dating (2011) provided for the Sardinian type G ORFV genomes cluster is consistent with the value (years 2009 and 2012) obtained by Coradduzza et al. [22] for strains isolated from Sardinian goats, on the basis of the ORFV VIR gene.
Our results suggest that the ORFV type S genomes clade is 31 years younger than type G and differentiated in the year 1796. This finding is consistent with the first report of contagious ecthyma in European (maybe Danish) sheep by Steeb in 1787 [2], suggesting that the early form of ORFV, associated with the development of contagious ecthyma in sheep, originated in northern Europe. However, with the only exception being a Sardinian, likely ancestral strain (GenBank # ON691522), all the type S genomes included in the present study share a common ancestor which dates back to 1828. This latter molecular dating may correspond to the beginning of the ORFV expansion across Europe and Asia and it matches with the first detailed description of contagious ecthyma in European sheep provided as early as 1890 by Wallay [9]. In this case, considering the young age of the Asian (Chinese) cluster (year 1953) of type S genomes, ORFV may have taken about a century (or even less) to spread from Europe toward Asia.
This evolutionary scenario in consistent with the extinction of most of the early European ORFV strains, as only those strains which performed better adaptative variation to the habitat are the ancestors of the genomes which expanded at the end of 19th century and, later in the 1920s, spread across the whole continent as suggested by several reports of contagious ecthyma in many countries [52][53][54][55].
The Sardinian ancestral type S strain (GenBank # ON691522) was isolated in 2021 from a lamb belonging to the Sarda sheep breed that lived on a farm in the southern part of the island. In particular, in accordance with personal communications obtained by these authors from the farm's owners, new heads were never introduced into their livestock and parent animals always either lived on the same farm or came from neighboring areas. In such a context, the Sardinian Sarda sheep breed was selected on the island in 1927/28 (in particular, in 1928 the first aries and the first sheep were registered in the Herd Book, they were owned by the itinerant Chair of Agriculture of Cagliari directed by Professor Francesco Passino) from heads arriving from the European mainland. Therefore, the ancestral ORFV Sardinian strain isolated in the present study could be a direct descendant of the early ancient ORFV strains that probably originated in northern Europe at the end of 18th century and then became extinct. This likely uncommon form of ORFV may be present in Sardinia simply by chance according to a typical effect of genetic drift and may have survived on the island as, after its arrival, it became well-adapted to the new environment. In such a context, this ORFV strain might be considered as a relic of the first ORFV variants that originated early on European mainland and is worthy of properly investigation to better understand the origins of this virus.
Conversely, the other Sardinian ORFV strains included within the type S genomes clade of the present study were isolated in the northern part of the island and they belong to a divergent genetic group that dates back to 1950 and includes viral allelic variants that originated in 1963 and 2001. The strains that have differentiated in the 21st century may be exclusive to the island, specifically the northern area. Another point of interest is that the molecular dating provided by the present research for the common ancestor of the ORFV genome strains isolated from sheep in Sardinia (year 1950) is quite consistent with the estimation reported by Coradduzza et al. [22] based on the ORFV VIR gene (year 1925). The small discrepancy could be related to the different molecular datasets used to perform molecular dating (whole genome 137,151 bp vs. VIR gene 382 bp). In general, the results obtained in the present study for ORFV strains isolated from heads living in northern and southern Sardinia, suggest that different strains of this virus arrived on this Mediterranean island during distinct periods of the 20th century. These highly divergent allelic variants were introduced in different areas of the island where various kinds of selective pressure might have originated new private and endemic ORFV variants. As a consequence, and as previously reported by Coradduzza et al. [22], the enclosed nature of Sardinia may have prompted the growth of a strong genetic differentiation among viral strains from different geographic areas. In the future, a higher number of ORFV whole genomes isolated from Sardinian hosts could help to shed further light of the evolutionary history of ORFV on the island and to provide evolutive predictive models of the distribution of this virus for areas of the world with similar environmental and selective conditions.

Conclusions
In conclusion, further analyses on many ORFV whole genomes are needed to corroborate the multiregional geographic origin in the late 18th century for ORFV strains infecting sheep (type S) and goats (type G), in Europe and Asia, respectively, as it is suggested in the present study.
Institutional Review Board Statement: Ethical review and approval were waived for this study, as this study did not involve any animal experiments. Samples were collected from sheep and goats using standard procedures by the Sardinian Veterinary Services A.T.S. and submitted to the Experimental Zooprofylactic Institute of Sardinia for ORFV testing. Special authorization for sampling activities was not necessary; this action is regulated by the Italian Ministry of Health and performed in the case of infectious diseases.

Informed Consent Statement:
Written informed consent has been obtained from the patient(s) to perform molecular analyses on biological tissue samples.

Data Availability Statement:
The sequences of the ORFV whole genomes obtained during the present study are openly available in the GenBank nucleotide sequence database. Accession numbers: ON691519-ON691525. See Table A2 in Appendix A for more information on the data accessibility.

Conflicts of Interest:
The authors declare no conflict of interest. * The contigs information row refers to de novo assembly. ** The percentage of identity has been inferred by using the blast tool implemented in the NCBI Virus portal. ON691519 SAMN28405493 SRR19732916 S10 ON691520 SAMN28405494 SRR19732915 S15 ON691521 SAMN28405495 SRR19732914 S19 ON691522 SAMN28405496 SRR19732913