Introduction
Malaria, caused by parasitic protozoans of the Plasmodium genus, remains a significant global health concern today. In 2023, an estimated 263 million malaria cases and 597,000 deaths occurred worldwide (World Health Organisation, 2024). Most of these are caused by P. falciparum (malaria tropica), the dominant species in Africa, whereas outside of Africa, P. vivax predominates. Although the latter is generally associated with reduced mortality, it can cause severe illness and death, particularly in children under 5 years and pregnant individuals (Singh et al., Reference Singh, Kumar, Rana, Thakur and Singh2013; Naing et al., Reference Naing, Whittaker, Wai and Mak2014; Howes et al., Reference Howes, Battle, Mendis, Smith, Cibulskis, Baird and Hay2016; Phyo et al., Reference Phyo, Dahal, Mayxay and Ashley2022; World Health Organisation, 2024). The remaining human-infecting species, P. ovale, P. malariae, and P. knowlesi, are generally rare and infrequently associated with severe disease (Yaro, Reference Yaro2014; Amir et al., Reference Amir, Cheong, de Silva, Liew and Lau2018), although the number of P. knowlesi cases is rising, particularly in South-East Asia (World Health Organisation, 2024).
Malaria parasites are transmitted to humans via the bites of female Anopheles mosquitoes and exhibit a complex life cycle involving both sexual and asexual stages within the mosquito vector and human host. Following transmission, sporozoites migrate to the liver, where they undergo replication before invading erythrocytes. The subsequent rupture of infected red blood cells leads to the release of new parasites in the bloodstream and it is this process that is responsible for the majority of clinical manifestations (Meibalan and Marti, Reference Meibalan and Marti2017). Infected individuals typically experience non-specific, flu-like symptoms such as fever, diarrhoea, myalgia and headache (Yaro, Reference Yaro2014; Loomis, Reference Loomis2018). Despite their general nature, these symptoms occur in a characteristic periodic pattern, every 2–3 days, reflecting the synchronised erythrocytic cycle and release of merozoites into the bloodstream (Loomis, Reference Loomis2018). Chronic infections often result in haemolytic anaemia and may progress to severe complications, the most notable being cerebral malaria (Yaro, Reference Yaro2014).
Although eradicated from most of the North-western hemisphere today, malaria historically had a far broader geographic distribution. In England, numerous references to malaria-like illnesses occur from the 16th century onward (Dobson, Reference Dobson1989; Gowland and Western, Reference Gowland and Western2012; Kendall, Reference Kendall2014; Perry and Gowland, Reference Perry and Gowland2022), and similar accounts are found throughout northwestern Europe (Knottnerus, Reference Knottnerus, Wefer, Berger, Behre and Jansen2002; Hulden and Hulden, Reference Hulden and Hulden2009; Deschutter and Seys, Reference Deschutter and Seys2016). In these more temperate regions, P. vivax was the dominant malaria species, particularly prevalent in low-lying wet areas, which are the ideal habitats for its vector, Anopheles atroparvus. The Netherlands are no exception to this. From the 17th century onward, historical records describe symptoms consistent with the disease. Sylvius de le Boë, a famous Dutch anatomist, reported an epidemic of cyclical fevers in Leiden in 1669. In 1826, a severe outbreak, apparently linked to North Sea flooding, claimed thousands of lives in the coastal provinces. In Amsterdam alone, approximately 2,400 of 200,000 inhabitants died from ‘malarial disease’ (Bruce-Chwatt and de Zuleta, Reference Bruce-Chwatt and de Zuleta1980). While malaria is thought to have spread northward and eastward with Roman expansion, introducing Plasmodium parasites to previously unexposed mosquito populations via infected human carriers (Knottnerus, Reference Knottnerus, Wefer, Berger, Behre and Jansen2002; Packard, Reference Packard2007; Newfield, Reference Newfield2017; Loufouma-Mbouaka et al., Reference Loufouma-Mbouaka, Gamble, Wurst, Jäger, Maixner, Zink, Noedl and Binder2021), recent DNA studies show Plasmodia endemicity in these areas that dates back thousands of years (Michel et al., Reference Michel, Skourtanioti, Pierini, Guevara, Mötsch, Kocher, Barquera, Bianco, Carlhoff, Coppola Bove, Freilich, Giffin, Hermes, Hiß, Knolle, Nelson, Neumann, Papac, Penske, Rohrlach, Salem, Semerau, Villalba-Mouco, Abadie, Aldenderfer, Beckett, Brown, Campus, Chenghwa, Cruz Berrocal, Damašek, Duffett Carlson, Durand, Ernée, Fântăneanu, Frenzel, García Atiénzar, Guillén, Hsieh, Karwowski, Kelvin, Kelvin, Khokhlov, Kinaston, Korolev, Krettek, Küßner, Lai, Look, Majander, Mandl, Mazzarello, McCormick, de Miguel Ibáñez, Murphy, Németh, Nordqvist, Novotny, Obenaus, Olmo-Enciso, Onkamo, Orschiedt, Patrushev, Peltola, Romero, Rubino, Sajantila, Salazar-García, Serrano, Shaydullaev, Sias, Šlaus, Stančo, Swanston, Teschler-Nicola, Valentin, Van de Vijver, Varney, Vigil-Escalera Guirado, Waters, Weiss-Krejci, Winter, Lamnidis, Prüfer, Nägele, Spyrou, Schiffels, Stockhammer, Haak, Posth, Warinner, Bos, Herbig and Krause2024).
From historical sources, it becomes clear that malaria was an important health concern throughout most of the (early) modern period in the Netherlands. Although genomic evidence confirms its presence (Michel et al., Reference Michel, Skourtanioti, Pierini, Guevara, Mötsch, Kocher, Barquera, Bianco, Carlhoff, Coppola Bove, Freilich, Giffin, Hermes, Hiß, Knolle, Nelson, Neumann, Papac, Penske, Rohrlach, Salem, Semerau, Villalba-Mouco, Abadie, Aldenderfer, Beckett, Brown, Campus, Chenghwa, Cruz Berrocal, Damašek, Duffett Carlson, Durand, Ernée, Fântăneanu, Frenzel, García Atiénzar, Guillén, Hsieh, Karwowski, Kelvin, Kelvin, Khokhlov, Kinaston, Korolev, Krettek, Küßner, Lai, Look, Majander, Mandl, Mazzarello, McCormick, de Miguel Ibáñez, Murphy, Németh, Nordqvist, Novotny, Obenaus, Olmo-Enciso, Onkamo, Orschiedt, Patrushev, Peltola, Romero, Rubino, Sajantila, Salazar-García, Serrano, Shaydullaev, Sias, Šlaus, Stančo, Swanston, Teschler-Nicola, Valentin, Van de Vijver, Varney, Vigil-Escalera Guirado, Waters, Weiss-Krejci, Winter, Lamnidis, Prüfer, Nägele, Spyrou, Schiffels, Stockhammer, Haak, Posth, Warinner, Bos, Herbig and Krause2024), the limited availability of written records means that very little is known about the disease’s distribution in medieval and earlier times. As a result, malaria is rarely considered in discussions of pre-modern health here, leaving us with an incomplete understanding of past societal well-being.
In recent decades, bioarchaeological studies have begun to address this gap, applying macroscopic and biomolecular analyses to skeletal remains from various periods and regions, providing new insights into the prevalence and impact of malaria in past populations (see Schats (Reference Schats2023) for a detailed review). While malaria does not cause any pathognomonic skeletal lesions, the extensive destruction of red blood cells during the erythrocyte stage of an active malaria infection often causes severe haemolytic anaemia (Yaro, Reference Yaro2014). Because of this, several studies have used a skeletal marker often associated with anaemia, termed cribra orbitalia (Figure 1), as a proxy for the disease in the archaeological record (e.g., Gowland, Reference Gowland and Eckardt2010; Gowland and Western, Reference Gowland and Western2012; Schats, Reference Schats2015, Reference Schats2022; Smith-Guzmán, Reference Smith-Guzmán2015a; Loufouma-Mbouaka et al., Reference Loufouma-Mbouaka, Gamble, Wurst, Jäger, Maixner, Zink, Noedl and Binder2021). Although this orbital porosity is non-specific (see Schats (Reference Schats2023) for extensive discussion), the consistent and statistically significant high prevalence of cribra orbitalia in malaria endemic areas suggests that the disease is an important contributing factor to its prevalence (Gowland, Reference Gowland and Eckardt2010; Gowland and Western, Reference Gowland and Western2012; Schats, Reference Schats2015, Reference Schats2022, Reference Schats2023; Smith-Guzmán, Reference Smith-Guzmán2015a). Thus, although cribra orbitalia cannot be used to diagnose malaria at the individual level, its spatial and population-level distribution can still provide meaningful insights.
Figure 1. Cribra orbitalia (different expressions) in the eye orbits of three different individuals from Vlaardingen, the Netherlands: (a and b) Non-adult (age range indeterminate), (c) Child 3–12 years, (d) Child 3–12 years (figure reproduced from Schats, Reference Schats2023).
Spatial epidemiological approaches have shown considerable promise in this regard (e.g., Gowland and Western, Reference Gowland and Western2012; Smith-Guzmán, Reference Smith-Guzmán2015a, Reference Smith-Guzmán2015b; Schats, Reference Schats2021a, Reference Schats2022). By correlating the geographic distribution of skeletal lesions with environmental, archaeological, and historical indicators of malaria risk, these studies have identified regional patterns of disease prevalence. Particularly important for the current study is the work by Gowland and Western (Reference Gowland and Western2012), who conducted a spatial analysis of cribra orbitalia in relation to distribution patterns of indigenous malaria, geographic variables, and the vector’s habitat. The authors found clear spatial clustering of cribra orbitalia in areas that would be high risk for a malaria infection in Anglo-Saxon England. The present study applies a similar approach to the medieval Netherlands, by assessing the occurrence of cribra orbitalia across skeletal collections nationwide and analysing these data spatially in relation to environmental predictors as well as historic mosquito density. In doing so, it investigates whether the spatial clustering of this orbital pathology observed in other regions is also present in the medieval Dutch context, and if such patterns can be linked to specific environmental variables or historic vector patterns. Ultimately, the findings may contribute to a better understanding of the distribution and historical burden of malaria in this region.
Materials and methods
Materials
A total of 2239 individuals from 28 archaeological sites in the Netherlands were analysed for the presence of cribra orbitalia (Table 1 and Figure 2). Sites were selected based on availability, dates between 500 and 1600 CE, and location in the Netherlands. Each site’s geographic coordinates (latitude and longitude, WGS84) were recorded based on modern geolocations corresponding to their historical location. Individuals for whom at least one orbital roof could be observed were included in this study.
Figure 2. Map of the Netherlands showing the archaeological sites included in this study.
Table 1. Demographic composition, cribra orbitalia prevalence, and mosquito density in studied sample
a After mosquito density as estimated by the number of mosquitoes found in stables in 1938 (Swellengrebel and De Buck, Reference Swellengrebel and De Buck1938, 71). High: more than 400 mosquitoes, Medium: 100–400 mosquitoes, Low: 1–100 mosquitoes (see Figure 4).
Methods
Skeletal analysis
As cribra orbitalia is generally more common in non-adults (Schats, Reference Schats2021b), due to the distribution of red marrow in the skeleton (Brickley, Reference Brickley2018, Reference Brickley2024), age-at-death was deemed an important variable in this study. The individuals in the dataset were therefore placed in one of two age groups: non-adults (≤ 20 years at death) and adults (≥ 21 years at death). Age was estimated based on the dental eruption (Ubelaker, Reference Ubelaker1979), permanent tooth formation (Moorrees et al., Reference Moorrees, Fanning and Hunt1963), epiphyseal fusion (Schaefer et al., Reference Schaefer, Black and Scheuer2009), and long bone length (Maresh, Reference Maresh and McCammon1970), as well as the development and morphology of the pubic symphysis (Brooks and Suchey, Reference Brooks and Suchey1990), auricular surface of the os coxa (Lovejoy et al., Reference Lovejoy, Meindl, Pryzbeck and Mensforth1985; Buckberry and Chamberlain, Reference Buckberry and Chamberlain2002), the sternal rib ends (Işcan et al., Reference Işcan, Loth and Wright1984, Reference Işcan, Loth and Wright1985), and on the degree of ectocranial suture closure (Meindl and Lovejoy, Reference Meindl and Lovejoy1985).
Cribra orbitalia was assessed macroscopically, occasionally with a 10 × magnifying glass. The lesion was considered present when porosity appeared in one or both orbital roofs, following criteria defined by Brothwell (Reference Brothwell1981), Rivera and Lahr (Reference Rivera and Mirazón Lahr2017), and Stuart-Macadam (Reference Stuart‐Macadam1985), with each individual classified as either showing or not showing cribra orbitalia (i.e., true prevalence rate). Only the presence or absence of lesions was recorded. Given the unclear relationship between lesion expression and clinical disease status (Grauer, Reference Grauer and Buikstra2019), severity or healing stages were not included in this study. Porosity on the orbital roof that was clearly the result of new proliferative bone deposition was excluded, as these lesions are unlikely to be related to the marrow hypertrophy associated with anaemia (Brickley, Reference Brickley2018, Reference Brickley2024; Brickley and Morgan, Reference Brickley, Morgan and Grauer2022). Cribra orbitalia prevalence was calculated per site as the percentage of individuals showing lesions relative to the total number observed. Prevalence was calculated separately for all individuals, as well as for adult and non-adult subsets, to assess potential age-related effects. Descriptive statistics and spatial visualisations were produced in R (Version 4.5.1) using the dplyr, sf and tmap packages. Site-level prevalence data were converted into a spatial object using the sf package with the recorded latitude and longitude for each site.
Spatial analysis: clustering
Spatial analyses were conducted to examine the geographic distribution of cribra orbitalia and identify potential clustering. Understanding spatial clustering is important because it can reveal non-random patterns in lesion prevalence, suggesting that environmental or disease-related factors may influence the distribution of cribra orbitalia (Pfeiffer et al., Reference Pfeiffer, Robinson, Stevenson, Stevens, Rogers and Clements2008). Analyses were conducted in R (Version 4.5.1) using the packages sf, spdep, tmap, ggplot2, ggrepel, and spatstat for spatial data handling, mapping, nearest neighbour analysis, and spatial autocorrelation. Points were projected to a planar coordinate system (Amersfoort/RD New, EPSG:28992) to allow distance-based analyses. All spatial analyses in this study are based on the digital GIS reconstruction of the Dutch palaeogeography for 1250 CE (Figure 3, Vos et al., Reference Vos, Meulen, Weerts and Bazelmans2025). This reconstruction was selected as it reflects a landscape after the onset of large-scale peat reclamation and human modification in the Netherlands (from c. 900 CE onwards), while also being broadly applicable to the few earlier medieval sites included in this study.
Figure 3. Palaeogeographic map (1250 CE) (adapted from Vos et al., Reference Vos, Meulen, Weerts and Bazelmans2025).
First, the spatial distribution of burial sites was analysed to investigate if sites were clustered, dispersed, or randomly distributed, regardless of cribra orbitalia prevalence. This provides a baseline measure of site distribution, helping to distinguish clustering of sites themselves from clustering of lesions (cf., Gowland and Western, Reference Gowland and Western2012). This was done using an Average Nearest Neighbour (ANN) approach, which calculates the ratio of the observed average distance between sites to the expected distance under complete spatial randomness using the shapefile derived from the GIS files associated with the 1250 CE palaeogeographic map (Figure 3, Vos et al., Reference Vos, Meulen, Weerts and Bazelmans2025). Ratios less than 1 indicate clustering, while ratios greater than 1 indicate dispersion (Siljander et al., Reference Siljander, Uusitalo, Pellikka, Isosomppi and Vapalahti2022). Statistical significance is set at P = < 0.05.
Following this, global spatial autocorrelation was assessed using Moran’s I to determine whether cribra orbitalia prevalence was clustered across the Netherlands (i.e., high or low prevalence values tend to be near similar values at other sites), with significance evaluated via Monte Carlo permutation tests using 999 iterations. To define spatial relationships among sites, a k-nearest neighbours approach was used. The optimal number of neighbours (k) was determined by examining the relationship between k and the average inter-site distance, identifying the point where the curve began to level off (‘elbow point’), while ensuring all sites were connected (Bivand et al., Reference Bivand, Pebesma and Gómez-Rubio2013; Irawan et al., Reference Irawan, Fahmi and Zamzami2024). A final value of k = 5 was selected for subsequent spatial autocorrelation analyses. Moran’s I ranges from −1 to +1, with a value closer to 1 indicating stronger autocorrelation, which means that sites with similar prevalence are spatially grouped rather than randomly distributed (Pfeiffer et al., Reference Pfeiffer, Robinson, Stevenson, Stevens, Rogers and Clements2008). Statistical significance is assessed by calculating Z-scores and associated p-values (Robinson, Reference Robinson2000). To visualise the spatial autocorrelation, a Moran scatterplot was generated, plotting each site’s prevalence against the average prevalence of its neighbours, allowing identification of High–High and Low–Low clusters as well as spatial outliers. Local Moran’s I (LISA) was calculated to detect local spatial clusters, including hotspots of high cribra orbitalia prevalence and cold spots of low prevalence (Anselin, Reference Anselin1995; Pfeiffer et al., Reference Pfeiffer, Robinson, Stevenson, Stevens, Rogers and Clements2008; Siljander et al., Reference Siljander, Uusitalo, Pellikka, Isosomppi and Vapalahti2022).
Spatial analysis: predictive modelling
To assess if environmental features associated with malaria mosquito habitat could explain any spatial variation in cribra orbitalia prevalence, predictive modelling was performed. Although temperature, rainfall, and elevation are factors commonly included in malaria prediction models (e.g., Dabaro et al., Reference Dabaro, Birhanu, Negash, Hawaria and Yewhalaw2021; Tariq et al., Reference Tariq, Bisanzio, Mutuku, Ndenga, Jembe, Maina, Chebii, Ronga, Okuta and LaBeaud2025), reconstructions of these variables at a resolution appropriate for the medieval Netherlands are not currently available. Thus, for each site, the minimum distance (in metres) to three key landscape features was calculated from the digital GIS file associated with the 1250 CE palaeogeographic map (Figure 3, Vos et al., Reference Vos, Meulen, Weerts and Bazelmans2025): (1) wetlands (which include peatlands, polders, salt marshes, stream valleys and river plains, tidal flats, and low-lying Holocene sediments), (2) coastlines and major open-water bodies, and (3) major rivers and waterways. Linear regression models (OLS) were used to investigate if proximity to these environmental features predicted differences in cribra orbitalia prevalence. Predictors were tested for independence to avoid using strongly related variables in the same model (Bivand et al., Reference Bivand, Pebesma and Gómez-Rubio2013). Analyses were conducted separately for the full dataset (all individuals) and for adults only. Because OLS models assume that residuals are spatially independent, residual spatial autocorrelation was assessed using Moran’s I with 999 permutations to check model assumptions. Spatial autocorrelation in the residuals suggests that sites are not independent. For example, neighbouring sites may share a similar environment, leading to similar lesion rates that are not statistically independent. Ignoring this clustering can result in inflated significance levels and ‘false positive’ results (Type 1 errors) (Lichstein et al., Reference Lichstein, Simons, Shriner and Franzreb2002). If significant spatial autocorrelation was detected in the residuals, the assumptions of OSL were violated, and spatial regression models were applied. Specifically, spatial lag models (which account for direct influence between neighbouring sites) and spatial error models (which account for shared unobserved environmental variables) were fitted using the spatialreg package, and model fit was compared using Akaike’s Information Criterion (AIC), where the lowest value indicates the better-fitting model (Bivand et al., Reference Bivand, Pebesma and Gómez-Rubio2013, 298; Xu and Kennedy, Reference Xu and Kennedy2015).
Historical mosquito density correlation
In addition to modelling the cribra orbitalia data in relation to palaeogeographic environmental features, lesion prevalence has been compared to A. atroparvus density, as recorded in 1938 by counting the number of mosquitoes found in stables in 1938 (Figure 4) (Swellengrebel and De Buck, Reference Swellengrebel and De Buck1938). While there is a considerable temporal gap between the medieval period and 1938, and the landscape likely underwent changes, this dataset provides an independent historical approximation of the ecological distribution of A. atroparvus prior to large-scale vector control and eradication efforts. The 1938 density map is used to assess if the spatial patterning observed in cribra orbitalia prevalence corresponds with areas that were historically suitable for the primary malaria mosquito in the Netherlands. The 1938 map was georeferenced using QGIS (version 3.42.0) with the outline of the modern Netherlands as a spatial reference. Considering the coarse and schematic nature of the source, the georeferenced map was used for broad regional comparison rather than detailed spatial analysis. Archaeological sites were plotted on the map (Figure 4) and based on their location, assigned to 1 of 3 ordinal mosquito density categories (Low (1-100), Medium (100-400), High (>400), or No data) (Table 1). To assess the correlation between historical mosquito density and medieval cribra orbitalia prevalence, Spearman’s rank correlation tests were conducted on both the complete and the adult-only datasets. The strength of association is expressed as rs, ranging between −1 and +1, reflecting a perfect negative or positive monotonic relationship. Following Akoglu (Reference Akoglu2018), correlations were interpreted as: weak (<0.4), moderate (0.4 –0.7), or strong (>0.7). Additionally, the spatial location of LISA-identified clusters was compared descriptively to mosquito density to assess if lesion clustering corresponded to areas of higher or lower historical vector presence.
Figure 4. Map of the Netherlands indicating mosquito density as estimated by the number of mosquitoes found in stables in 1938 with the sites in the current study. Black: more than 400 mosquitoes, narrow lines: 100–400 mosquitoes, wide lines: 1–100 mosquitoes, white: no data (adapted from Swellengrebel and De Buck, Reference Swellengrebel and De Buck1938, 71).
Results
Cribra orbitalia prevalence
Site-level prevalence values of cribra orbitalia for all individuals, adults, and non-adults are summarised in Table 1. Among all individuals, cribra orbitalia prevalence ranged from 6.5% to 33.9% with an average prevalence of 17.9%, while in the adult-only dataset it ranged from 0% to 32.7% (average 14.4%). Prevalence among non-adults demonstrated greater variability (0–70%) and a substantially higher average (37.3%). However, as 13 of the 28 sites have fewer than 10 non-adults, the results are skewed both to the high and low ends of the range. Yet, when removing those 13 sites, the average prevalence remains close to the overall average at 37.8%, confirming the higher prevalence compared to adults. Overall, the distribution of cribra orbitalia prevalence shows clear visual geographic variability (Figure 5), supporting the need for further analysis of spatial patterning. To minimise the influence of uneven age distributions and limited non-adult samples in some of the sites, these analyses were conducted on the combined sample (all individuals) and the adult subset, and not independently for non-adults.
Figure 5. Cribra orbitalia prevalence (in percentage) in the studied sites plotted on map of the Netherlands: (a) all individuals, (b) non-adults, (c) adults.
Spatial analysis: clustering
The Average Nearest Neighbour analysis results show a tendency towards clustering (R = 0.83, Z = –1.76, P = 0.078) indicating that the sites are more grouped than expected under complete spatial randomness. While not statistically significant, this result suggests that site locations were not strongly dispersed and may reflect mild regional clustering.
To investigate if cribra orbitalia prevalence demonstrated spatial patterning, global Moran’s I was calculated. When all individuals were included, results showed a strong and significant positive spatial autocorrelation (I = 0.437, P = 0.001). The adult-only dataset showed weaker spatial autocorrelation, yet the results remained positive and statistically significant (I = 0.160, P = 0.04). The Moran scatterplots (Figure 6) illustrate this: sites in the upper-right quadrant (High–High) are locations with high prevalence surrounded by similarly high-prevalence neighbours (hotspots), whereas sites in the lower-left quadrant (Low-Low) indicate low-prevalence clusters (cold spots). Locations in the remaining quadrants (High–Low and Low–High) reflect spatial outliers where site-level prevalence is different from that of the surrounding sites. In the complete dataset, which includes all individuals, most sites fall within the High–High and Low–Low quadrants, suggesting consistent clustering of similar values, with only a few outliers. In contrast, the adult-only dataset has several Low–High outliers, fitting with the observed weaker global spatial autocorrelation found for this group.
Figure 6. Moran scatterplots: (a) all individuals, (b) adults.
In the all-individual dataset, local Moran’s I (LISA) identified several significant clusters (Figure 7a). Five sites formed significant High–High clusters, and five Low–Low clusters were detected. Two High–Low outliers were present, representing high-prevalence sites bordered by lower-prevalence neighbours. In contrast, the adult-only dataset showed a weaker local structure, with only two Low–Low clusters and no High–High clusters or significant spatial outliers (Figure 7b).
Figure 7. LISA cluster maps: (a) all individuals, (b) adults.
Spatial analysis: predictive modelling
Linear regression models were fitted using distance to wetlands, rivers and the coastline as predictors to study if these environmental factors played a role in the spatial distribution of cribra orbitalia. For both the full dataset and the adult-only subset, none of the environmental variables were found to be significant predictors of cribra orbitalia prevalence (P > 0.19), and the models explained little of the observed variation (all individuals: R 2 = 0.14; adults: R 2 = 0.05). Residuals from the adult model showed no spatial autocorrelation (Moran’s I = 0.08, P = 0.12). Yet, residuals from the model created for the complete dataset were significantly autocorrelated (Moran’s I = 0.29, P = 0.004). Therefore, spatial regression was conducted on this group. Both a spatial error model (SEM) and a spatial lag model (SLM) improved model fit (AICOLS = 203.3; AICSEM = 198.9; AICSLM = 197.7), with strongly significant spatial parameters (SEM λ = 0.60, P < 0.001; SLM ρ = 0.58, P < 0.001). However, even after accounting for spatial dependence, none of the environmental distance predictors were significant (P > 0.33).
Historical mosquito density correlation
Spearman’s rank correlation tests were conducted to assess the relationship between historical mosquito density and medieval cribra orbitalia prevalence. The analysis of the complete dataset showed moderate, statistically significant positive association between cribra orbitalia prevalence and mosquito density (rs = 0.517, P = 0.014). When considering the adults only, the correlation is weaker and not statistically significant (rs = 0.362, P = 0.098).
The LISA analysis identified five High–High clusters in the complete dataset, four of which are located in areas with medium to high historical mosquito density. Only Delft falls in an area where fewer than 100 malaria vectors were recorded in 1938 (see Figure 4). The Low–Low cluster sites are all situated in regions for which mosquito density data are unavailable. Of the two High–Low outliers, Kampen is located in a low-density area, while Amersfoort lies on the border of a medium-density region. The analysis of the adult-only dataset only identified two Low–Low clusters, located in an area without density data.
Discussion
This study investigated (1) the spatial distribution of cribra orbitalia prevalence across 28 archaeological medieval sites in the Netherlands, and (2) if this distribution can be explained by underlying environmental features that would be consistent with malaria. Overall, the results show spatial clustering of cribra orbitalia for both analysed datasets, yet the patterning is not explained by broad-scale environmental factors such as distance to the coast, rivers, and wetlands. The sections below will focus on interpreting the observed spatial patterns and investigating the role malaria may have played in determining cribra orbitalia prevalence, taking into account the limitations of the sample.
Cribra orbitalia prevalence showed high variability across the studied populations, particularly for the non-adults. While the variation in non-adult prevalence rates can in part be explained by sample limitations, the adults-only dataset, for which sample sizes are more consistent, also demonstrated a broad prevalence range that, visually, appears to cluster in certain areas.
Before assessing the spatial patterning of cribra orbitalia prevalence, clustering irrespective of the orbital lesions was investigated. While the ANN results showed that the sites are not completely randomly distributed across the Netherlands, as can be expected due to settlement history and excavation and sampling bias, the clustering is mild and not statistically significant and is therefore unlikely to be a driver in the subsequent spatial analyses.
The results of the global Moran’s I showed positive and significant spatial autocorrelation for both datasets (all individuals and adults-only). This means that archaeological sites where cribra orbitalia prevalence is high are located near other high-prevalence sites, and likewise for low-prevalence sites. The analysis of the complete dataset showed the strongest spatial patterning with only a few visual outliers in the Moran scatterplot. The adult-only subset demonstrated weaker, but still significant, spatial autocorrelation, but with quite a few more Low–High outliers (sites with low prevalence surrounded by high-prevalence sites). Thus, when the non-adults are removed from the analysis, several of the sites become outliers. This observation could be explained by the fact that several of the adult outliers (e.g., Blokhuizen, Rijsterbos, Scheemda, Vlaardingen) have high non-adult cribra orbitalia prevalence and a relatively low adult prevalence compared to surrounding sites. This result may point to the (paradoxical) possibility that these are sites with high health stress, resulting in high childhood mortality and thus fewer adults with lesions (Wood et al., Reference Wood, Milner, Harpending, Weiss, Cohen, Eisenberg, Hutchinson, Jankauskas, Cesnys, Česnys, Katzenberg, Lukacs, Mcgrath, Roth, Ubelaker and Wilkinson1993; Gowland and Western, Reference Gowland and Western2012; DeWitte and Stojanowski, Reference DeWitte and Stojanowski2015). While small non-adult sample size is an issue for many of the collections, the patterns may point to a non-random distribution of (childhood) disease risk factors, potentially malaria. This will be further discussed after the predictive modelling.
Local Moran’s I (LISA) results allow for better insights into where clusters and outliers occur. The results for the complete dataset support the Moran’s I scatterplot trends, demonstrating a significant High–High cluster in the Southwest of the Netherlands, particularly in the current province of Zeeland, confirming that this area represents a regional hotspot of cribra orbitalia. While the Moran’s I scatterplot also identified High–High sites in the North of the Netherlands (Alkmaar, Hoogwatum and Groningen), these are not part of a statistically significant local cluster, indicating less spatial coherence (i.e., unclustered high values) with regard to prevalence. Low–Low clusters are observed to the East and South of the Netherlands, in line with the Moran’s I result, indicating a regional cold spot. Even though not all sites are part of a cluster, these LISA results of the complete dataset demonstrate regional structuring in cribra orbitalia prevalence and clear separation of the sites in the West and East of the Netherlands. The LISA results for the adult-only dataset only identified two Low-Low clusters, in the East of the Netherlands, supporting the pattern observed in the complete dataset. The lack of High-High clusters in the adults is in line with the Moran’s I results: the observed spatial structure, especially in the Southwest, is largely driven by non-adult individuals.
While the Moran’s I and LISA results demonstrate spatial clustering of cribra orbitalia prevalence, with a separation between the East and the West, none of the environmental variables (i.e. proximity to the coast, rivers, and wetlands) were significant predictors. This indicates that these broad-scale geographic factors alone are unlikely to be the main determinants of the observed differences in prevalence. If malaria was a contributor to the cribra orbitalia prevalence, a relationship between these variables was expected, since low-lying marshy areas with slow-flowing or stagnant water provide ideal breeding conditions for the Dutch malaria vector, A. atroparvus (Jetten and Takken, Reference Jetten and Takken1994). However, what the predictive model fails to take into account is the salinity regime of the water in the different regions. Wetlands and rivers occurred throughout the Netherlands, and while these would all be excellent mosquito habitats, the A. atroparvus was mainly found in areas with a NaCl content of 800–2500 mg/l. Although this malaria vector can also breed in fresh water, it was outcompeted there by other species, particularly by the better-adapted A. messeae (a non-malaria mosquito) (Jetten and Takken, Reference Jetten and Takken1994; Van der Kaaden, Reference Van der Kaaden2003), meaning that fresh-water wetlands are unlikely to be associated with malaria transmission risk. A study by Jacob Zwarteveen in the Noordoostpolder from 1943 supports this. He looked at the relative abundance of A. messeae and A. atroparvus and found that the percentage of A. atroparvus larvae increased with higher NaCl concentrations (Zwarteveen, Reference Zwarteveen1948; Van der Kaaden, Reference Van der Kaaden2003).
When zooming in on the clusters identified in the LISA analysis in relation to the palaeogeographic map, Figure 8 shows that both clusters occur in areas identified as wetland within the predictive model. However, the High–High sites (red) are located in an area where water salinity is expected to have been high, with frequent flooding of the land (Henderikx, Reference Henderikx, Brusse, Henderikx and Beekman2012), whereas the Low–Low clusters (blue) occur near what is most likely fresh water. Thus, general distance to wetlands or rivers is likely not the best environmental predictor for medieval malaria in the Netherlands, particularly when the breeding ecology of its local vector is taken into account.
Figure 8. Palaeogeographic map with significant LISA clusters.
This, however, does not explain why the ‘distance to the coast’ variable was not a significant predictor for cribra orbitalia prevalence in this study. If malaria caused by the A. atroparvus is a contributing factor, proximity to the coast, as it is likely to impact water salinity to some extent, was expected to correlate with prevalence. Yet, it is important to realise that distance to the coast likely only works as a salinity proxy if this salinity decreases gradually with distance inland. For the Netherlands, this is very unlikely to have been the case, particularly in the medieval period. The coastal areas of the Netherlands, especially the province of Holland, became increasingly characterised by artificial ecosystems. Large-scale peat reclamations, starting in the 10th century AD, created a new, human-made landscape (Borger, Reference Borger and Verhoeven1992; Hoppenbrouwers, Reference Hoppenbrouwers, De Nijs and Beukers2002; Van Bavel, Reference Van Bavel2010; Hoffmann, Reference Hoffmann2014), which may have resulted in saline intrusion further inland or freshwater influx in areas which would have been more brackish before. Thus, distance to the coast is unlikely to reliably reflect salinity, and thus A. atroparvus vector ecology, in the medieval Netherlands and this may be an explanation for why it was not a significant predictor for cribra orbitalia prevalence in this study, even if malaria was present. This suggests that the palaeogeographic information used to model lesion prevalence is likely too coarse to adequately reflect the micro-environments important for A. atroparvus vector ecology and is therefore not accurately predicting potential malaria transmission risk.
This is further supported by the comparison with historical mosquito density as recorded in 1938. While there is a considerable time jump from the medieval period to 1938, it is clear that mosquito presence in the 20th century showed marked heterogeneity, with high density more inland and low-density areas at the coast, illustrating the existence of micro-environments, some of which were more or less conducive for A. atroparvus breeding. Moreover, although the blank area on the map does not have any mosquito data, the strong East–West divide observed in cribra orbitalia prevalence is also notable here. The moderate positive association between cribra orbitalia prevalence and historical mosquito density in the complete dataset supports the notion that malaria may have contributed to skeletal anaemia in medieval populations. The weaker association observed in adults is likely a consequence of lower adult prevalence at several sites within historically high-density areas, possibly reflecting high childhood mortality or selective survival, as suggested before. Most High-High LISA clusters coincide with regions of medium to high mosquito density, reinforcing the link between areas conducive to breeding and elevated cribra orbitalia prevalence. Delft is the exception here, a High–High site located in a historically low-density area, which could point to landscape changes that may have impacted the malaria transmission risk over time, potentially highlighting one of the limitations of using the 1938 map for comparison. Interestingly, one of the High-Low outliers (Amersfoort) is located in an area which had high mosquito density in 1938, which may potentially explain the higher-than-expected cribra orbitalia prevalence if medieval ecological conditions were similar to the 20th century.
While the limitations of the comparison with a 1938 map cannot be ignored, the results suggest that the archaeological clustering in cribra orbitalia prevalence aligns closely with historical vector distribution, providing supporting evidence that malaria contributed to this orbital lesion in certain regions of the Netherlands.
Conclusion
This research has used spatial epidemiology to investigate patterning in cribra orbitalia in the medieval Netherlands and to study if this could be linked to environmental and historical indicators of malaria risk. Although differences between age groups were observed, statistically significant clustering was present, with hotspots in the Southwest and cold spots in the East. None of the environmental predictors explained the observed clustering, likely because the variables did not fully capture the complex A. atroparvus vector ecology nor the micro-environments created by medieval human engineering, such as peat reclamation and drainage. The identified high-prevalence clusters are located in low-lying, likely brackish areas, whereas the low-prevalence clusters occur in freshwater zones where the malaria vector was absent. Moreover, cribra orbitalia prevalence correlates with historical mosquito density. Taken together, the spatial clustering in relation to local environmental conditions and the correlation with historical mosquito presence all suggest that malaria was a contributing factor to cribra orbitalia prevalence in the medieval Netherlands.
While these results provide evidence for the role of malaria, interpretation is constrained by small non-adult sample sizes at several sites, as well as uneven preservation and excavation biases across the Netherlands. The comparison with historical mosquito density data is valuable, but the coarseness of the map, the considerable time gap, and potential landscape changes through time should be taken into account. Future research combining more extensive bioarchaeological sampling, refined environmental and palaeogeographic reconstructions, and potentially biomolecular approaches could help further clarify the impact malaria may have had on past populations.
Data availability statement
Data underlying this research can be found at 10.5281/zenodo.17803944.
Acknowledgements
The author wants to express their gratitude to all the archaeological depots that have granted them access to the skeletal material. Thanks are due to Milco Wansleeben for help with determining appropriate spatial analysis tools in the initial phases of the project. Marijke Langevoort assisted the author with QGIS work for this paper. The author acknowledges the use of ChatGPT (OpenAI, accessed September 2025 through January 2026) as a tool to assist in coding and executing the statistical analyses in R.
Author contributions
R.S. conceived and designed the study, conducted data gathering, performed statistical analyses and wrote the text for the article.
Financial support
This work was financially supported by the Dutch Research Council (NWO) (Grant reference: 016.Veni.195.195) and Elise Mathilde Fund/Leiden University Fund (Grant reference: W18357-6-EM).
Competing interests
The author declares there are no conflicts of interest.
Ethical standards
For working with archaeological human remains, the author strictly followed the Ethics Guidelines prepared by the Laboratory for Human Osteoarchaeology, Faculty of Archaeology, Leiden University as well as those compiled by the Dutch Association for Physical Anthropology (NFVA). The individuals were treated with utmost respect, and the analyses were performed with explicit permission of the relevant stakeholders.
Open access
