Genomic signals of migration and continuity in Britain before the Anglo-Saxons

The purported migrations that have formed the peoples of Britain have been the focus of generations of scholarly controversy. However, this has not benefited from direct analyses of ancient genomes. Here we report nine ancient genomes (∼1 ×) of individuals from northern Britain: seven from a Roman era York cemetery, bookended by earlier Iron-Age and later Anglo-Saxon burials. Six of the Roman genomes show affinity with modern British Celtic populations, particularly Welsh, but significantly diverge from populations from Yorkshire and other eastern English samples. They also show similarity with the earlier Iron-Age genome, suggesting population continuity, but differ from the later Anglo-Saxon genome. This pattern concords with profound impact of migrations in the Anglo-Saxon period. Strikingly, one Roman skeleton shows a clear signal of exogenous origin, with affinities pointing towards the Middle East, confirming the cosmopolitan character of the Empire, even at its northernmost fringes.

A ncient genomics has the power to anchor the emergence of modern genetic patterns to archaeological events but, to date, no such genome-scale data have emerged for the Romano-British world, or indeed from any era in the British past. Extensive surveying of modern genomic variation in the British Isles has produced divergent interpretations of the migratory history of the islands. An east-west gradient of Y chromosome, autosomal and mtDNA allele frequencies has been interpreted as reflecting the genetic legacy of substantial Anglo-Saxon invasions following the Roman period [1][2][3][4] . However, it is difficult to distinguish the effects of this much-debated event from other migratory influences from northwest continental Europe, whether these are, for example, Germanic elements in the Late Roman army predating the Anglo-Saxon migrations or Scandinavian settlers arriving some centuries afterwards 5 .
At its maximum, the Roman Empire stretched from Atlantic Europe to the Near East and from Northern Britain to the Sahara, incorporating an advanced transport infrastructure, which would have enabled previously unprecedented levels of mobility 6 . Evidence for the presence of foreigners in Britain has been based on epigraphic sources, material culture and, more recently, bioarchaeological (isotopic) data [7][8][9] . However, there is no way of knowing how representative the people mentioned in inscriptions might be; artefactual imports indicate contact beyond the province but may not point to movement of people; and isotopic baseline values in British data overlap considerably with numerous other regions, including much of Western Europe and the Mediterranean littoral 10 . Thus, measuring who moved about the Roman Empire, and on what scale, remains challenging.
To investigate the genetic identity of Britain in late BC and the early centuries AD we report shotgun genomic sequencing of nine human genomes to coverage depth of B1 Â and analysis of these in the context of extensive genome-wide data from modern populations. Seven ancient genomes are sampled from a cemetery in Roman York dated between the second and the fourth century AD, one from an earlier Yorkshire Iron-Age burial (210 BC-40 AD) and one from a later neighbouring Anglo-Saxon burial (650-910 AD).

Results
Archaeological samples. York (Eboracum), founded c. AD 71, became the Roman empire's northernmost provincial capital in about AD 200. Its southwest approach road was lined with tombstones and mausolea 11 and excavations conducted there between 2004 and 2005 at Driffield Terrace revealed a cemetery dating from the second to the fourth century AD with a high incidence of decapitated remains 12 . After initial screening of eight individuals chosen for DNA extraction, seven presented superior endogenous DNA content. Near Melton, East Yorkshire, a first century AD late Iron-Age settlement gave human remains predominantly of adult females and non-adults 13 from which five individuals were sampled. From these, the best preserved individual with known date was M1489 (between 210 BC and 40 AD). Norton Bishopsmill, dates to 650-910 AD and was a Christian Anglo-Saxon cemetery excavated in the village of Norton, Teesside, northeast England 14 . We sampled 3 individuals from burials of 100 skeletons and selected the best preserved, NO3423, for the present study (Supplementary Notes 1 and 2 and Supplementary Fig. 1).
Sequencing results and sample contamination. Gamba et al. 15 established that the dense internal centre of the inner ear petrous bone is an excellent source of preserved ancient DNA suitable for high-throughput sequencing. Accordingly, we sampled single petrous bones, extracted DNA and made indexed Illumina sequencing libraries. After preliminary screening for endogenous content, nine samples were chosen for genome-level sequencing. On alignment to the human reference genome, reads showed mismatch patterns typical of archaeological DNA indicating deamination damage ( Supplementary Fig. 6). Contamination estimated from both mtDNA heterozygosity (1.82 ± 0.47%) and X-chromosome contamination in male samples (0.79 ± 0.21%) was low (Supplementary Note 2 and Supplementary Tables 5-7).
Sex determination and uniparental marker analysis. Using the ratio between sequencing reads aligned to the X and Y chromosomes 16 , it was possible to assign biological sex to each individual, confirming skeletal assessments: the Anglo-Saxon and each Roman-period sample were male, whereas the Iron-Age sample was female ( Supplementary Fig. 7). Mitochondrial genomes were retrieved for each sample with between 39 Â and 98 Â coverage and were assigned to known haplogroups (Table 1), which are common in present-day European populations 17 . Y-chromosome haplogroups were determined for each male (Table 1): the majority (6/7) of Driffield Terrace samples belong to sub-lineages of R1b-L52/L11, which reaches its highest frequencies (470%) in Western European countries 18 . Sample 3DRIF-26, on the other hand, despite belonging to the same burial context, presented a lineage consistent with haplogroup J2-L228, which has a modern distribution centred on the Middle East, but which is also present in the Caucasus region, the Balkans and Italy 19 . The Anglo-Saxon (NO3423) sample was assigned to haplogroup I1-S107, which is widespread in Nordic countries 20 .
Affinity with global populations. We called between B210 and B400 thousand single-nucleotide polymorphisms (SNPs) within our ancient samples that had previously been genotyped in a data set of 780 European, West Asian, North African and Middle Eastern individuals 21 . Figure 1a shows a principal component analysis (PCA) where eight of nine ancient genomes cluster close to a collection of northwest European samples. One York Roman 3DRIF-26 gives a clear Middle Eastern signal, with closest neighbours of Palestinian, Jordanian and Syrian origin. This dichotomy is also apparent in maximum likelihood estimation of individual ancestries using NGSadmix (Fig. 1b). In this, when a model of three ancestral populations is imposed across the entire sample, this analysis highlights three major geographical foci: Europe; North Africa; and West Asia/Middle East. The European ancestral component predominates in the majority of ancient samples (which show similar profiles to modern northwestern Europeans), whereas 3DRIF-26 again shows a majority West Asian/Middle Eastern component. Isotopic analyses of the skeletons support this genetic differentiation of 3DRIF-26 from the remainder of the individuals sampled. Strontium isotope ratios ( 87 Sr/ 86 Sr) vary mainly according to geological substrate, while oxygen isotope values (q 18 O), which track locally available drinking water, reflect climatic and geographic variables such as temperature, rainfall levels or distance from the coast 22 . When we compared these ratios in our seven samples with other British Romans, 3DRIF-26 showed both an unusually low 87 Sr/ 86 Sr ratio and an extreme q 18 O p value ( Supplementary Fig. 2).
To maximize resolution of genetic affinity, for each ancient sample we performed pairwise comparisons with each modern sample and calculated the proportion of SNP positions at which these were identical by state (IBS). Only single-SNP alleles were considered at each locus and were randomly sampled from the biallelic genotype. Taking the median IBS score for each modern population sample, we then ranked these for similarity to each ancient genome in turn. Interestingly, the top-ranked modern sample for IBS for each of these ancient British samples was one of the formerly Celtic language-speaking regions of the British Isles, with the single exception of 3DRIF-26, which showed highest IBS with samples from Saudi Arabia. We gauged the sensitivity of this approach by checking whether individual modern samples were assignable to their region of origin. When tested, local individuals were assigned with high frequency (0.97) to the British Isles and also most often to their correct country. The method showed lower sensitivity for Middle Eastern genotypes, with primary assignment to that region in only 39% of instances ( Supplementary Fig. 15). Nevertheless, outside assignments tended to be to Cypriot, Sardinian and Druze, never to Northern Europe. In contrast, specificity of a Middle Eastern assignment was high-only three individuals, from Iran, Tunisia and Morocco, were incorrectly assigned to that region. Thus, assignment of 3DRIF-26 to the Middle East region seems secure, but resolution to an individual population may not be possible. Specificity in assignment to the British Isles was lower, with about half of assignments (0.53) derived from elsewhere, most often from neighbouring populations such as France (0.28) and Norway (0.15). Small sample sizes (B10 per population) render individual scores only weakly informative but, when we compared the six Roman burials after excluding the outlying 3DRIF-26, their rank orders across the geographical sample were highly correlated (Spearman rank correlation coefficients  Tables 11 and 12). This allowed us to consider these together and generate a combined percentile score by scoring each comparator population as the product of the rank percentiles versus each Roman genotype (Fig. 2). The Welsh were most consistently ranked as highly identical by state, followed by Irish and Scottish scores, a result strongly supporting an origin within the British Isles for this Roman sample majority. Interestingly, the modern English sample was ranked only ninth for IBS to the Romans from York, at a level similar to German, Norwegian, Orcadian and Basque samples.
Ancient sample ancestry within Britain. To place our ancient genomes within a detailed British context, we next plotted these in a background PCA using 3,075 published genotypes from British 3 , Irish 23 and southern Netherlands samples 24 . The modern samples were analysed using SNP genotypes at B250,000 loci and projected into a single plot using smartpca (Fig. 3a). As in Burton et al. 3 the first component of the variation was informative for the structure within Britain. Given the close ancestral relationships between these populations and their wellknown history of migrational exchange, a substantial overlap between regional groups was both expected and observed. However, by considering median values, one can see a clear progression from Irish samples at one pole through Scottish, Welsh, English to the Dutch cohort at the other extreme. In this plot the York Romans cluster centrally close to the modern Welsh median value, along with the Iron-Age genome. The local Anglo-Saxon is placed differently, closest to modern East Anglians between the English and Dutch medians. This first component also offers an opportunity to compare within the English sample. Figure 3b shows a boxplot of PC1 values for each subsample and structure is evident, with higher median values in Eastern regions such as East Anglia, East Midlands, intermediate values in the southern and western parts and lower values in the north and northwest. This pattern is more clearly seen in a geographical plot of interpolated values (Fig. 5a). When the York Romans are compared together with each modern cohort, they are most similar to the Welsh distribution of PC1 values and differ significantly from all other regional groups, apart from those from North and Northwest England (Mann-Whitney test; Fig. 3b, Supplementary Note 2 and Supplementary  Table 13). An interesting difference is the marked one between the Driffield Terrace ancient and contemporary Yorkshire samples (P ¼ 0.003), implying regional discontinuity. It is also worth noting that the PC1 coordinate of the Anglo-Saxon individual is closer to the median PC1 value of East Anglians, possibly reflecting a more pronounced contribution of Germanic immigrants to eastern British populations. However, we note the inherent uncertainty in drawing inference from a single sample.
Leslie et al. 4 used a haplotype-based clustering method, implemented in fineSTRUCTURE 25 to deconstruct British populations into geographically and historically meaningful clusters. Accordingly, we used this approach to search for substructure within the modern British genotypes used here and to identify patterns of allele sharing between the ancient samples and the clusters identified. Structure was apparent with separation into subgroups of predominantly Welsh, English and Scottish provenance. The six Romans shared most alleles on average with those clusters consisting primarily of Welsh individuals ( Fig. 4 and Supplementary Fig. 16). Allele sharing patterns also allowed comparison with the other ancients; when median IBS values across the clusters were compared between those for the Iron-Age genome and those for the Roman cohort, these correlated strongly (r ¼ 0.74, P ¼ 0.004), supporting continuity. However, a comparison between the Anglo-Saxon and the Romans showed no correlation (r ¼ 0.06, P ¼ 0.842; Fig. 4).
Genomic change in Yorkshire between the early centuries AD and modern sampling is further illustrated by both the oldest and most extensively typed genetic system (the ABO blood group) and the system known to show maximal differentiation within the British Isles (chromosome Y haplogroup R1b1a2-M269). The plots in Fig. 5 show imputed blood group O and observed chromosome Y haplogroup R1b1a2-M269 frequencies for the Roman genomes (excluding the immigrant outlier) contrasting sharply with interpolated allele frequencies for modern eastern Britain 1,26-28 .
Imputation and phenotype determination. Using a similar approach to that of Gamba et al. 15 , we used phased reference genomes from the 1000 Genomes Project to impute genotypes associated with phenotypic traits. In particular, we inferred genotypes at SNP positions to predict eye and hair    29 . The most common predicted phenotype in the Roman burial samples is brown eyes and black/brown hair. However, one sample, 6DRIF-18, was estimated to have had a distinctive appearance with blue eyes and blonde hair, as did the single Anglo-Saxon individual. We also inferred that blood group O is the most common in the Roman samples (Supplementary Table 17). The Iron-Age sample is also estimated as blood type O and the Anglo-Saxon is likely to have been type B or possibly type A. Five samples returned imputed lactase persistence genotypes: two Roman burials and the Iron-Age individual were likely to have been lactase persistent, while two Romans, 6DRIF-22 and the suspected migrant 3DRIF-26 were homozygous for the ancestral non-persistence variant.

Discussion
Combined genomic and isotopic evidence support the inference that the origins and childhood of individual 3DRIF-26 lay far outside Britain. His modern genomic affinities clearly lie with the Middle East. Isotopically, the most plausible suggestion is an arid environment on igneous or limestone geology, which is consistent with the same regions ( Supplementary Fig. 2 and Supplementary Note 2). Hence, although this individual is indistinguishable from the other inhumations in terms of burial practice and osteology, the analyses show that, even in its northernmost provincial capital, the profoundly cosmopolitan nature of the Roman Empire suggested by documentary and epigraphic sources continued to hold sway.   The peoples of Britain show marked genetic structure (Fig. 5), which has been the focus of generations of investigation 30 . A common theme in this research has been a contrast between a southern and eastern lowland zone and a western and northern upland zone. This patterning resembles the geography of Anglo-Saxon settlement in the fifth to seventh centuries AD, inviting the conclusion that the cultural and linguistic change effected by this migration was also reflected by major genetic change 1,2,31 .
Projections from modern data to the past are, however, subject to considerable uncertainties and may be compounded by unknown complexities, which do not feature in their underlying models. Prehistorians point out that the Germanic affinity of eastern Britain could also be a result of earlier communications with the northwest European mainland. For example, there may have been 'Belgic' peoples in Britain at the time of the Claudian conquest in AD43 (ref. 32), and the Roman army that arrived in Britain was composed of recruits from various provinces 33 . Recently Leslie et al. 4 have used haplotype-based statistical methods applied to modern genome-wide SNP genotypes to infer several distinct ancestral influences from migratory events into Britain. This included a major 35% contribution to modern Central and Southern English populations from a German source, which, they surmise, occurred in the century after AD 800, some 200 years or more after archaeological evidence for initial Anglo-Saxon influence. Evidence from direct observations of ancient genomes is required, however, if we are to draw conclusions about genetic exchange that distinguish between closely dated events.
Six of the seven individuals sampled here are clearly indigenous Britons in their genomic signal. When considered together, they are similar to the earlier Iron-Age sample, whilst the modern group with which they show closest affinity are Welsh. These six are also fixed for the Y-chromosome haplotype R1b-L51, which shows a cline in modern Britain, again with maximal frequencies among western populations. Interestingly, these people differ significantly from modern inhabitants of the same region (Yorkshire and Humberside) suggesting major genetic change in Eastern Britain within the last millennium and a half. That this could have been, in some part, due to population influx associated with the Anglo-Saxon migrations is suggested by the different genetic signal of the later Anglo-Saxon genome. Iron-Age, Roman, Anglo-Saxon, Viking and other migrations have all been proposed as contributors to the genetic structure in modern United Kingdom 34 .
The thesis that the mountainous regions of Wales may have held populations that are representative of earlier, more widely dispersed indigenous British genetic strata is not new, yet it finds some support in our analyses. The genomes of modern Scottish and Irish populations diverge from this group of early inhabitants of northern Britain, whereas their Welsh counterparts do not. Modern data for genetic structure among non-Saxon samples from the British Isles are said to deny the existence of a single 'Celtic' population 4 . Our data indicate that differentiation within such groups may have happened before the early centuries AD. By the same token, it lends support for genetic exchange between Scotland and Ireland, as attested in some historical sources and mirrored by linguistic affinity: Irish and Scottish Gaelic are sister P-Celtic languages, whereas Welsh is a divergent Q-Celtic language, akin to that thought to have been spoken throughout pre-Roman Britain 35 .
In the Roman York burials at Driffield Terrace, the majority were adults under 45 years old, male and most had evidence of decapitation 36 . They were slightly taller than average for Roman Britain, displayed a high occurrence of trauma, potentially related to interpersonal violence and evidenced childhood stress and infection (Supplementary Tables 1 and 2 and Supplementary Note 1). This demographic profile resembles the population structure in a recently excavated burial ground of the second and third century AD at Ephesus, which has been interpreted to be a burial ground for gladiators 37 . However, the evidence could also fit with a military context; the Roman army had a minimum height for recruitment 38 and fallen soldiers would match the young adult profile of the cemetery. In this later Roman period increasingly large numbers of soldiers were enlisted locally 33 .
Whichever the identity of the enigmatic headless Romans from York, our sample of the genomes of seven of them, when combined with isotopic evidence, indicate six to be of British origin and one to have origins in the Middle East. This is the first refined genomic evidence for far-reaching ancient mobility and (although from an unusual context) also the first snapshot of British genomes in the early centuries AD, indicating continuity with an Iron-Age sample before the migrations of the Anglo-Saxon period. Read processing and analysis. Next-generation sequencing reads were trimmed with Cutadapt v. 1.3 (ref. 44). Two bases from each side of reads were removed with seqtk (https://github.com/lh3/seqtk). Reads were aligned to the human reference genome (UCSC hg19) and mtDNA (rCRS, NC_012920.1) with Burrows-Wheeler Aligner v.0.7.5a-r405 (ref. 45 Contamination estimates and authenticity. To determine the extent of contamination in the ancient samples sequenced, we calculated the number of mismatches in mtDNA haplotype defining mutations 15 and X-chromosome polymorphisms in samples determined to be male (Supplementary Note 2 and  Supplementary Tables 5-8) 48 . We also used PMDtools 49 to select reads with evidence of deamination and compared sex determination and PCA (Supplementary Note 2 and Supplementary Figs 7 and 8). Finally, we confirmed the presence of aDNA misincorporations by analysing a subset of 1 million reads for each sample with mapDamage 2.0 (ref. 47).

Methods
Sex determination and uniparental lineage determination. We used the method published in ref. 16 to determine the sex of the ancient individuals ( Supplementary  Fig. 7). Y-chromosome lineages of ancient male samples were identified using clean_tree software 50 (http://www.erasmusmc.nl/fmb/resources/cleantree/; Supplementary Table 10). Regarding mtDNA analysis, reads were separately aligned to the revised Cambridge Reference Sequence (rCRS; NC_012920.1) 51 , filtering for base (qZ20) and mapping (qZ30) quality and duplicate reads as above. We then used SAMtools to obtain a consensus sequence in fasta format, which we uploaded to HaploFind 52 , which determines mtDNA haplogroups based  (Figs 1a and 3a and Supplementary Figs 8 and 9) were performed using smartpca from the EIGENSOFT v.5 package 56 and LASER 2.02 (ref. 57; Supplementary Fig. 10). Model-based clustering analysis was done with ADMIXTURE v. 1.23 (ref. 58). We then extracted genotype likelihoods from aDNA data with ANGSD v.0.592 (ref. 48; Supplementary Fig. 11), which we combined with genotype data of present-day populations. We analysed this data with NGSadmix v.32 (ref. 59; K ¼ 3; Fig. 1b and Supplementary Note 2).
Identity-by-state analysis. IBS between ancient and present-day samples was estimated with PLINK v. 1.9 (ref. 60). Median IBS proportions between aDNA samples and European, Middle Eastern and North African populations were obtained and plotted individually on maps ( Supplementary Fig. 12). Then, we selected only the Roman York samples, except for 3DRIF-26 (Middle Eastern affinity), ranked their median IBS score in relation to modern populations and combined these ranks by calculating their product (Fig. 2). Spearman rank correlations were estimated with R (ref. 61; Supplementary Note 2 and Supplementary Table 11). Regarding the WTCCC1 genotypes, we followed the same approach as described above, but present scaled median IBS values to the 0-1 range on Fig. 4. fineSTRUCTURE analysis. We randomly selected 100 individuals from each region of the WTCCC1 data set, excluding SNPs with missing genotypes, which resulted in a total of 431,366 variants and 1,000 samples. We used SHAPEIT2 to phase genotypes 62 and ran the ChromoPainter pipeline 25 with default parameters as implemented by fineSTRUCTURE v.2. For the fineSTRUCTURE analysis, the following settings were used: 3,000,000 burn in iterations, 1,000,000 sample iterations for the Markov chain Monte Carlo (MCMC) and 10,000,000 tree comparisons (Supplementary Note 2). We then called genotypes in our ancient samples for estimation of IBS between these and the inferred fineSTRUCTURE population clusters ( Fig. 4 and Supplementary Figs 16 and 17).
Imputation of phenotype associated loci. To investigate loci associated with selective sweeps, we took a similar approach as in ref. 15, where alleles observed in the 1000 Genomes Project 63 were called with Genome Analysis Toolkit v.2.5, from which we extracted genotype likelihoods and converted to BEAGLE format. BEAGLE v.3.3.2 (ref. 64) was used to phase and subsequently impute genotypes at SNP positions described in the HIrisplex system 29 , loci associated with blood groups 65,66 , lactase persistence 67,68 and pigmentation phenotypes 69,70 . Only posterior genotype probabilities Z0.85 were kept 15 (Supplementary Note 2 and  Supplementary Table 14 -17). We generated interpolated frequency maps of blood group 26 and Y-chromosome frequency data 1,27,28 with ArcMap v.10.1 from the ArcGis suite (Environmental Systems Research Institute) using the default settings of the geospatial analysis plugin (Fig. 4).