Bioarchaeological insights into the last plague of Imola (1630–1632)

The plague of 1630–1632 was one of the deadliest plague epidemics to ever hit Northern Italy, and for many of the affected regions, it was also the last. While accounts on plague during the early 1630s in Florence and Milan are frequent, much less is known about the city of Imola. We analyzed the full skeletal assemblage of four mass graves (n = 133 individuals) at the Lazaretto dell’Osservanza, which date back to the outbreak of 1630–1632 in Imola and evaluated our results by integrating new archival sources. The skeletons showed little evidence of physical trauma and were covered by multiple layers of lime, which is characteristic for epidemic mass mortality sites. We screened 15 teeth for Yersinia pestis aDNA and were able to confirm the presence of plague in Imola via metagenomic analysis. Additionally, we studied a contemporaneous register, in which a friar recorded patient outcomes at the lazaretto during the last year of the epidemic. Our multidisciplinary approach combining historical, osteological and genomic data provided a unique opportunity to reconstruct an in-depth picture of the last plague of Imola through the city's main lazaretto.

Imola is a city located in the Italian province of Bologna and the region Emilia-Romagna. The city, which was formerly also known as Forum Cornelii during Antiquity, was founded around 82 BCE along the Via Emilia by the river Santerno and was an important urban centre during the Italian Renaissance. Imola was hit by a pestilence during the Italian Plague of 1629-1631, which had reached Northern Italy via soldiers engaged in the Thirty-Year War 1,2 . The plague epidemic started in Milan, where first cases of plague were reported in mid-October 1629. Shortly after, French troops also introduced plague to the Piedmont region 1 . The disease then proceeded to ravage through Northern and Central Italy 3 . In 1630 alone, 34 Italian cities were affected by the disease 3 . In fact, the seventeenth century marks one of two periods of demographic regression for Italy, with the second being after the Black Death during the 1300s 4 . Throughout the Northern Italian plague epidemics of 1629-1631 and the Southern Italian plague epidemics of 1656-1657, a total estimated 20% of the Italian population was lost 4 . Northern Italy was hit particularly heavily, with mortality rates ranging from 20 to 60% (see Fig. 1). Preceding typhus epidemics and famines in 1628 and 1629 are speculated to have weakened the population's immune defense prior to the arrival of plague 1 .
First reports of an epidemic in Imola were recorded on the 29th of June 1630 2 , when the disease reached the city from Bologna despite the implementation of dedicated quarantine and safety measures. While contemporary historical texts describe buboes in Imola 2 , no molecular evidence for the presence of the plague pathogen Yersinia pestis during these epidemics has been available thus far. As plague does not leave identifiable skeletal lesions on the bones of infected individuals, ancient DNA (aDNA) analysis is the only available method to confirm the presence of Y. pestis in archaeological samples unequivocally. We apply a multidisciplinary approach to study the plague outbreaks in Imola between 1630 and 1632 by analysing the remains discovered in four mass graves excavated at the Lazaretto dell'Osservanza and by integrating data from an unpublished register detailing the outcomes of patients, which were treated at the lazaretto at the time of the outbreak.

B.
A. Osteology. We analyzed all 133 skeletons stemming from the mass graves excavated at the site "Osservanza".
We identified 48 female and 40 male individuals (Table 1), but were unable to confidently determine the biological sex for the remaining 45 individuals, of which 43 were subadults. Age of death could not be assigned within the defined categories for eight adult individuals, but it was possible to determine the age of death for all but one subadult individual. All age classes were represented (Table 1), as can be expected for epidemic death assemblages 14 . The least represented age classes were infants (8.1%) and old adults (5.6%).
Our paleopathological analysis recorded signs of multiple infectious diseases. Lesions associated with tuberculosis were identified on five individuals: one subadult (8-10 y), three middle adult females and one male classified as young adult. Skeletal changes were frequently recorded on vertebrae, but lesions were also observed on ribs and sternum. Affected vertebrae and ribs showed signs of lytic lesions and cavitation on their vertebral body and on their proximal ends, respectively. However, these lesions are not specific to tuberculosis and could have been caused by other conditions 15 . One young adult male could be a potential case of syphilis, as it displayed signs of periosteal reaction with prominent periosteal striations and cortical thickening, as well as a saber malformation of the tibia, which is a common observation for the treponemal disease 15,16 . One child showed signs of poliomyelitis; the femurs were asymmetrical, and the bone showed signs of malformation due to the premature fusion of the femoral head. We also found evidence for metabolic diseases. We identified potential cases of rickets or osteomalacia on the skeletons of two young adults, which showcased characteristic malformations on their lower limbs. We also recorded a potential case of malignant neoplasia in a middle adult female. We observed a thickening of the individual's diploe, as well as increased porosity of the internal cranial vault and lytic lesions on the external cranial vault. Similar lytic lesions were also found on the postcranium, particularly on the left clavicle, on some ribs, on the sternum, the femoral neck, and on the os coxae. We could not identify any perimortem traumas, but observed skeletal evidence for antemortem trauma, i.e. trauma exhibiting signs of remodelling, on three adult males and one old adult female. Ancient DNA. We sampled teeth from 15 skeletons (Table S1). Of these, IMO2 and IMO7 amplified aDNA using primers for the Y. pestis genes caf1 and pla respectively and IMO3a was screened via shotgun sequencing. Based on the PCR results, IMO2a and IMO7a were selected for further metagenomic exploration via shotgun sequencing. Using Kraken2 17 , we found that most reads assigned to Yersinia (Table 2) were classified at the genus level or to the so-called "Yersinia pseudotuberculosis complex", which encompasses all Y. pestis and Yersinia pseudotuberculosis species. Contrary to Kraken2, Metaphlan2 18 yielded no Y. pestis hits for all shotgun datasets (Table 2 and Fig. S2). This is not surprising considering the low number of observed hits using Kraken2. Kraken2, which implements a k-mer based approach, is less sensitive, but incorporates a much bigger database www.nature.com/scientificreports/ than Metaphlan2, which is limited to unique clade-specific marker genes. This distinction is especially important when considering the misincorporations found in aDNA. Shotgun reads for all three samples were also mapped against the Y. pestis CO92 reference sequence (Table 3 and Fig. S5). While all samples had comparable numbers of reads mapping to the reference sequence, sample IMO2, which also yielded the most hits using Kraken2, had the highest number of reads mapping to CO92 above a mapping quality of 30 and ca. 14% of reads displayed the characteristic C > T aDNA deamination at their 5' end ( Fig. S3). Across all three samples, most reads mapped to unspecific regions of the genome. Additionally, IMO2 showed coverage for all plasmids (Table 3). We could not identify clear signs for the presence of Y. pestis or any other pathogen within the shotgun datasets from IMO3 and IMO7.
In order to estimate the overall aDNA preservation across all three samples, we mapped the shotgun data to the rCRS build of the human mitochondrion as well as the Hg38p.12 build of the human allosomes. We recovered 40-100% of the human mitochondrion at mean depths of coverage between 1.4 and 114.7X (Table 4) across all three sequenced samples. The reads mapping to the human mitochondrion also showcased clear signs of deamination (Fig. S4), which combined with the small insert sizes (mean 48.59-69 bp) observed across all samples, is consistent with aDNA. Finally, a thermal age analysis performed for Imola also indicates that the recovery of good quality DNA was possible at the site (λ = 0.0047).
Following the screening of shotgun metagenomic datasets, we enriched libraries of sample IMO2 for the plague pathogen Y. pestis, using a custom RNA baits design. Unique reads mapping to the Y. pestis CO92 reference genome (MQ > 30) were extracted and taxonomically classified using both Kraken2 and Metaphlan2 (Table 2, Fig. S2). We saw a clear increase in hits for the pathogen with up to 81.7% of identified reads being classified as Yersinia or to a lower taxonomic rank using Kraken2 and all identifiable reads being assigned to the species Y. pestis by Metaphlan2. While we saw an increase of Y. pestis hits when compared to the shotgun dataset, the amount of Y. pestis sequences recovered via target enrichment still proved inadequate to reconstruct sufficient parts of the Y. pestis genomes at a viable depth of coverage. We observed C > T deamination rates of ca. 9.2% at position 1 for the capture mapping (-n 0.01, MQ > 0) ( Table 5). Our mapping had a sequence coverage of 16.8% Table 2. Results of the metagenomic analysis of shotgun and target enrichment datasets using Kraken2 (Wood et al. 17 ) and Metaphlan2 (Segata et al. 18 ). Kraken2 records the number of hits for each taxon and Metaphlan2 the abundance of identified reads assigned to each taxon. Mapped data denotes reads extracted from a mapping to the CO92 reference genome.  . While low when considering that multiple enriched libraries were used for this assembly, the coverage across plasmids seems to increase with copy number and we could confirm the presence of the Y. pestis specific plasmids pPCP1 and pMT 19 .

Discussion
The results of this study identify the causative pathogen of the epidemic, which affected Imola between 1630 and 1632, as Yersinia pestis. In the light of these results and using a multidisciplinary analysis, we evaluated unpublished historical records 11 in combination with our osteological analysis to provide insights into the epidemic event represented by the four excavated mass graves at the Lazaretto "Osservanza". The comparison of our osteological analysis with the register of Francesco Da Gazzo showed some demographic discrepancies, which can probably be attributed to fluctuation in mortality throughout the epidemic. While the number of identified males and females was comparable in the excavated mass-graves (see Tables 1 and  6), the friar's records show a clear prevalence of female patients and female deaths at the lazaretto at the end of the epidemic. Similarly, we detected a higher-than-expected prevalence of subadult individuals in our assemblage (27%) (see Table 1), relative to the low numbers of child deaths reported in the register (4.9%) (see Table 6). It should also be pointed out that one of the difficulties when comparing the data of the register with our analysis was the definition of the word "putto", which designated children in the register. It is unclear how broad the age category falling under "putto" is, based on the social norms of the time. The age distribution observed in our data, with a high percentage of children, adolescents and young adults, is typical for a catastrophic assemblage. However, the mortality rates we observed for infants in the mass graves of Osservanza are lower than expected. While this can often be observed in plague assemblages 14,[20][21][22] , the incredibly small number of individuals below the age of three identified in Imola "Osservanza" could also be explained by a separate burial place or grave pit for infants, for example outside of the city, or the existence of another location for the care of sick children. However, the records from the lazaretto of 1632 also show that only 4.9% of admitted patients were children, which seems to coincide with the low incidence of infants identified during our analysis and demonstrates that children were treated at the lazaretto. In light of the bad bone preservation observed in the mass-graves, it is also possible that additional subadult skeletons, which are more likely to be affected by diagenetic processes, were deposited, but not recovered during the excavation. We also cannot exclude the possibility that the mass-graves contained individuals, which had not been treated at the lazaretto.
Overall, our palaeodemographic results are consistent with those recorded for other plague sites 20,23,24 . Apart from a few individuals for which we recorded pathogenic markers, the analysed individuals were generally in good health (Table 7). We found little to no evidence of physical trauma, which is typical for epidemic www.nature.com/scientificreports/ assemblages. In addition to the effects of comorbidities, further analyses on the "frailty" of the individuals recovered from the mass graves of Imola "Osservanza'' should be considered, as they could provide additional evidence for the overall health of the plague victims in Imola by examining the presence of metabolic markers, such as porotic hyperostosis, cribra orbitalia and enamel hypoplasia 25 .
The mass-graves excavated at the site of Imola "Osservanza" were dug during the plague epidemic of 1630-1632 to accommodate the burial of large amounts of plague victims, which was confirmed by our metagenomic analysis, which provides evidence for the presence of the plague pathogen Y. pestis in the teeth of deposited individuals. The graves excavated during the survey campaign in 2007 seem far from the described 100 organized graves, planned to house no more than ten individuals each. In fact, there is little order within the excavated burials that each contained between 19 and 65 individuals, which were seemingly deposited without particular care in between layers of lime (Fig. 2C). Considering the archaeological records, which certainly only represent a fraction of the plague pits to be expected at the site of "Osservanza", we could speculate that graves of trench II and III 7 were dug in 1632, the final and most lethal year of the epidemic. However, the recorded discrepancy in child deaths between the Da Gazzo register and our skeletal assemblage (Tables 1 and 6) probably place the time of death for the individuals discovered in the mass graves prior to June-August 1632 and thus before Da Gazzo started recording the activity at the Lazaretto dell'Osservanza.
The use of lime to cover the cadavers is noteworthy (Fig. 2). Lime was believed to have disinfectant properties, to limit the scavenging of graves by predators and to reduce the odour of decay 8 . It has been documented in a number of known plague pits (eg. Barcelona 26 and Marseille 20 ) and has been shown to considerably slow down the decomposition of cadavers 27,28 . The use of lime has also been shown to have a desiccating effect on low-depth tissue 27 , such as the sampled anatomical region, the cranium. Interestingly, lime was used in multiple layers rather than in one final layer at the site of Imola "Osservanza", leaving some cadavers effectively protected from the surrounding soil and its biome (Fig. 2). This could have influenced the diagenesis of the buried bodies and might explain some of the variation in aDNA preservation observed within grave 3. Both skeletons IMO3/ n2 and IMO7/l2 were recovered from the very bottom of the grave, and were, therefore, completely exposed to the surrounding soil from beneath. These skeletons yielded little to no evidence of the plague pathogen. IMO2/ u2, however, was recovered in between layers of lime and on top of other skeletons and yielded low, but detectable, amounts of Y. pestis aDNA.
The overall preservation of the skeletal human remains recovered from the Lazaretto "Osservanza'' was poor, with some individuals missing all, but the densest bones of their skeleton upon discovery. However, our mapping to the human mitochondrion shows that while aDNA was not conserved equally amongst the sequenced samples, the overall aDNA preservation was adequate for the identification of bacterial sequences. The sample IMO2, which carried traces of the plague pathogen, also carried the highest amounts of human DNA with its mitochondrial mapping showing a depth of coverage of more than 15X across 99.9% of the sequence. Yet, IMO7a, a sample from which no clear sign for the presence of Y. pestis could be recovered, also yielded a high coverage assembly of the human mitochondrion with a good deamination signature, indicating that the overall aDNA preservation of the sampled tooth was good. Beyond environmental and taphonomic factors, the course of disease and the timing/occurence of bacteraemia could have played a role as well.
Much like today, plague probably manifested in various forms throughout an epidemic during the Second Plague Pandemic. However, the low mortality rates recorded in Imola 2,11 , especially compared to other cities affected by plague at the time (Fig. 1), tentatively allows us to exclude the much more lethal pneumonic plague 29 , as the predominant form of the disease spreading through the city and its rural surroundings during the early 1630s. Hence, the population of seventeenth century Imola was very likely affected by bubonic plague. This is further corroborated by descriptions of buboes found in historical text documenting the epidemic 2 . Notwithstanding, other historical sources 2 hint at possible cases of primary septicaemic plague, as they feature reports of plague victims without buboes.
The third wave of plague in 1632 is characterised by a shift in seasonality. Instead of running from June to November/December, as it had in previous years, the plague hit Imola in March and ended in August. The city of Imola also recorded its highest yearly mortality in 1632, which was 6.7 and 1.7 times higher than in 1630 and 1631, respectively. As Imola was one of only a handful of cities still affected by plague past 1631 30 , there is little available data for comparison. Florence, which also experienced plague in 1632, recorded cases of the disease for the entire year, starting with sporadic cases in January and a peak in mortality between the end of September and the end of December 31 . We investigated if the increase in mortality in Imola for the year 1632 could be related to an increase in birth rates in the previous years, but found that there had been no significant change to the birth rates recorded for the city over the last 5 years 2 . Trends for the region of Emilia-Romagna from 1530 to 1632 also Table 7. Frequency of pathologies identified on individual from the excavaed mass-graves at the Lazaretto dell'Osservanza (ND: not determined). www.nature.com/scientificreports/ did not show any significant increase except for a peak in 1624 before the region was hit by famine and typhus in 1628-1629 32 . We investigated recorded changes in temperature and humidity and found that, based on a weather diary recorded in Lisbon (December 1630-January 1633) 33 , the spring of 1632 was the wettest period on record with 57 days of rain, resulting in particularly favourable conditions for arthropod vectors 34 .

Infections Trauma Neoplasms Metabolic diseases Joint diseases
In conclusion, although we could not detect high amounts of Y. pestis via our aDNA analysis, the site of Imola and the conducted multi-disciplinary analyses, have given us the opportunity to describe a multi-facetted plague epidemic. The results of our analysis also point to the possibility of differential aDNA preservation in a grave containing multiple layers of lime, which could be valuable information for designing a sampling strategy in similar graves. Although further data will be needed to validate this hypothesis. As we only screened 12% of individuals excavated from Imola "Osservanza" for the presence of plague, further metagenomic exploration of individuals recovered from the mass-graves could yield sufficient data for a phylogenetic analysis in the future.

Material and methods
For this study, a total of 133 individuals were analyzed osteologically. Of these, 15 individuals were chosen for ancient DNA (aDNA) analysis. We extracted two teeth from nine skeletons (IMO1, IMO3-IMO10), one tooth from five subadult skeletons (IMO11-15) and seven teeth from skeleton IMO2 (Table S1). Sampled teeth were either molars or premolars, which were still embedded in alveolar bone and exhibited dental wear equal or below stage 4 as defined in Smith and Kight 35 . Selected teeth also had to be free of dental pathologies and cavities.
We sampled at least two individuals per mass grave, except for grave 3, from which four individuals were selected as it contained the largest amount of skeletons. Individuals were chosen from the top and bottom layers of the graves, ideally separated by at least one layer of lime. The teeth were sampled under clean conditions at the Laboratory of Archaeo-Anthropology and Forensic Anthropology, Department of Neuroscience and Rehabilitation in Ferrara, Italy. They were then sealed and sent off to the aDNA laboratory of the University of Oslo for further analysis. To avoid contamination of the skeletal material, the analysis was performed using protective clothing such as coveralls, masks and gloves at all times.
Osteology. Sex and age-at-death were assessed by morphometric methods. Estimation of the age at death for immature individuals was estimated primarily based on the degree of dental maturation and dental eruption, diaphyseal length of the long bones and stage of ossifications. For adults, we determined the age-at-death by analysing age-related changes to the pubic symphysis and the auricular surfaces of the os coxae. Whenever these anatomical sites were missing, we based our analysis on dental wear and the ossification of the ectocranial sutures of the cranial vault [36][37][38][39][40][41][42][43][44][45][46][47][48][49] . Sex diagnosis was performed only for adults and late adolescents (≥ 15 years of age) using the morphological traits of the skull and the pelvis and metric traits if these bones were not available [50][51][52] . The samples were assigned to the following age classes after Buikstra and Ubelaker 36 : fetus (pre-natal), infant (0-3 yrs old), child (3-12 yrs old), adolescent (12-20 yrs old), young adult (20-34 yrs old), middle adult (35-49 yrs old), old adult (50 + yrs old). Pathologies were also recorded to provide further background information. The morphological examination was based on a macroscopic analysis.
Ancient DNA laboratory work. The aDNA analysis was performed in a dedicated aDNA laboratory at the University of Oslo, Department of Biosciences. The introduction of contamination into the workspace is, amongst other measures, minimized by the use of protective clothing, including Tyvek suits, gloves, masks, helmets/visors and hairnets. Benchtops and other surfaces in the lab are routinely wiped down with soap, bleach and/or DNA ExitusPlus. Additionally, the lab is equipped with filtered ventilation, positive airflow and ceiling UV lights. All reagents and equipment in the aDNA laboratory are dedicated solely to the study of aDNA. No modern samples have ever been processed in the laboratory.
Sample preparation. After the introduction into the laboratory, samples were UV irradiated for 30 min on each side in a dedicated UV cupboard. The outer 1-2 mm surface of the teeth was then mechanically removed using a dental Renfert Basic quattro IS sandblaster. After another cycle of UV irradiation, the samples were milled to powder using a Retsch Oscillating Mill.
Extraction. We used two protocols to extract DNA from the Imola samples. Samples IMO1-10ab and samples IMO11-15a were extracted using protocol A. Samples IMO2d-g were extracted using protocol B.
Protocol A: 0.120 g of tooth powder was extracted using a modified version of the extraction protocols published by Dabney et al. 53 and Bennett et al. 54 . 1 ml of lysis buffer (0.25 mg/ml Proteinase K, 0.45 M EDTA) was added to each sample and incubated on a nutator at 38 ℃ overnight. The lysate was then added to a binding buffer (6X QG Buffer, 4X Isopropanol) and passed through a Qiagen MinElute Spin Column and a Zymo-Spin V 15 ml extension reservoir using a Qiagen Qiavac Vacuum Manifold. The column was washed twice with 1 ml of Qiagen PE Buffer. After a dry spin, the sample was eluted in two steps with a total of 50 µl Qiagen EB Buffer.
Protocol B: 0.20-0.60 g of tooth powder was extracted using a modified version of the protocol published by Brotherton et al. 55 . Samples were digested in 4.3 ml of lysis buffer (0.5 M EDTA, 0.25 mg/mL Proteinase K, 0.5% N-Lauroylsarcosine) on a nutator at 38 ℃ overnight. The lysate was added to 2.5 ml of binding buffer (1X Triton, 0.2 M Acetic Acid, 20 mM NaCl, 13.5 ml of Qiagen QG Buffer) and 120 µl of medium-sized silica suspension and was left to incubate on a nutator for 90 min. The samples were then pelleted, and the silica pellet washed three times with 1 ml of 80% Ethanol. Finally, the dried silica pellets were eluted in 200 µl Qiagen EB Buffer. Genomic and metagenomics analysis of sequencing data. Data preparation and filtering. Raw datasets were demultiplexed at the Norwegian Sequencing Centre, University of Oslo. Using cutadapt, the datasets were trimmed of their adapters and indices. All reads below 30 bp were filtered out, the quality cut-off for read termini was set to 20. Dataset quality was assessed using FastQC 61 .
Metagenomic analysis. To remove redundant and uninformative sequences we filtered our shotgun datasets using prinseq-lite 62 . We removed all exact duplicates, as well as all sequences with a dust score above 5. The paired-end reads were then merged using FLASH (-z -M 125) 63 and analysed using Kraken2 17 using a custom database containing all dusted bacterial, viral, archaeal and protozoan reference sequences, as well as the UniVec_Core database. The datasets were also analysed using Metaphlan2 64 . Reads mapping to the reference sequence CO92 were extracted and converted to fastq format. We then ran Kraken and Metaphlan2 with the same settings as described above. Prior to the extraction of the reads, duplicates had been removed using Picard's MarkDuplicates module 65 . Data used for genomic alignments were not filtered using prinseq.
Thermal age analysis. We performed thermal age analysis after Smith et al. 76 on the portal http:// therm alage. eu/. The parameters were set as follows: age of the samples was based on the year 1632, we chose the available DNA depurination (Bone) model, excavation year was set to 2007, year of analysis was set to 2018 and storage conditions were set to 22 ℃ with an annual variation over a range of 15 ℃.

Data availability
Sequencing data will be made available upon publication at the European Nucleotide Archive (ENA) under the accession number PRJEB27822.