Large-scale recent expansion of European patrilineages shown by population resequencing

The proportion of Europeans descending from Neolithic farmers ∼10 thousand years ago (KYA) or Palaeolithic hunter-gatherers has been much debated. The male-specific region of the Y chromosome (MSY) has been widely applied to this question, but unbiased estimates of diversity and time depth have been lacking. Here we show that European patrilineages underwent a recent continent-wide expansion. Resequencing of 3.7 Mb of MSY DNA in 334 males, comprising 17 European and Middle Eastern populations, defines a phylogeny containing 5,996 single-nucleotide polymorphisms. Dating indicates that three major lineages (I1, R1a and R1b), accounting for 64% of our sample, have very recent coalescent times, ranging between 3.5 and 7.3 KYA. A continuous swathe of 13/17 populations share similar histories featuring a demographic expansion starting ∼2.1–4.2 KYA. Our results are compatible with ancient MSY DNA data, and contrast with data on mitochondrial DNA, indicating a widespread male-specific phenomenon that focuses interest on the social structure of Bronze Age Europe.

C ontroversy has surrounded the origins and antiquity of the people of Europe, focused on the proportions descending from Neolithic farmers originating B10 thousand years ago (KYA), or from earlier Palaeolithic hunter-gatherers. Early studies observed a European SE-NW cline in classical gene frequency data which was ascribed to demic diffusion of farmers 1 , or, in an alternative view, to the first Palaeolithic colonization 2 . More recent autosomal genome-wide SNP data sets reflect current population structure 3,4 and admixture during the last 3,000 years 5 , but have provided little insight into older population processes.
Most debate on European prehistory has been stimulated by analyses of uniparentally-inherited markers. Spatial patterns in maternally-inherited mitochondrial DNA (mtDNA) are nonclinal, with age estimates of haplogroups (hg) taken to suggest a major Palaeolithic contribution 6 . Analyses of diversity in the male-specific region of the Y chromosome (MSY) show significant frequency clines in major lineages 7 , and geographical distributions and dates based on short-tandem repeats (STRs) have led to interpretations of both Palaeolithic 8 and Neolithic 9 major components. The most frequent western European lineage, hg R1b-M269, was originally believed to have originated in the Palaeolithic 10 , but in more recent analysis was assigned a Neolithic origin 11 , a claim challenged in turn 12 on the basis of STR choice and sample ascertainment. In general, dates based on STRs are problematic because of uncertainty about appropriate mutation rates, and possible long-term mutation saturation due to their stepwise mutation processes 13 . Palaeolithic dates for the major lineages are challenged by scanty ancient MSY DNA data, which suggest a marked discontinuity between 5-7 KYA and the present 14 .
A major cause of the controversy about MSY evidence is that unbiased estimates of diversity and time depth have until recently been impossible to obtain in large samples. Next-generation sequencing (NGS) generally offers unbiased ascertainment of MSY SNPs, providing phylogenies in which topologies inform about past demography, and branch lengths are in principle proportional to time, avoiding dating problems associated with STRs. Some insights have emerged from recent work 15,16 , but no systematic population-based NGS study across Europe has yet been undertaken.
Here, we use targeted NGS of European and Middle Eastern populations to show that Europe was affected by a major continent-wide expansion in patrilineages that post-dates the Neolithic transition. Resequencing at high coverage of 3.7 Mb of MSY DNA, in each of 334 males comprising 17 population samples, defines an unbiased phylogeny containing 5,996 highconfidence single-nucleotide polymorphisms (SNPs). Dating indicates that three major lineages (I1, R1a and R1b), accounting for 64% of the sampled chromosomes, have very recent coalescent times, ranging between 3.5 and 7.3 KYA. In demographic reconstructions 17 a continuous swathe of 13/17 populations from the Balkans to the British and Irish Isles share similar histories featuring a minimum effective population size B2.1-4.2 KYA, followed by expansion to the present. Together with other data on maternally inherited mtDNA 16,18 and autosomal DNA 19 , our results indicate a recent widespread male-specific phenomenon that may point to social selection, and refocuses interest on the social and population structure of Bronze Age Europe.

Results
Samples and approach. To address the lack of a systematic population-based NGS study of European MSY diversity, we assembled a collection of 20 randomly chosen male DNA samples from each of 17 populations from Europe and the Middle East (Supplementary Table 1). We used a sequence-capture approach to successfully generate 3.7 Mb of unique MSY sequence from 334 of the total set of 340 males. Mean read-depth of 51 Â allowed us to call a set of 5,996 high-confidence SNPs (Supplementary Data 1), with a minimum 6 Â read-depth. SNP calls were validated using publicly available Y-chromosome sequence information from genome-wide data (Supplementary Table 2).
SNPs have been deposited in NCBI dbSNP database and ss numbers can be found in Supplementary Data 1.
Phylogeography of European MSY lineages. We constructed a maximum-parsimony tree displaying the phylogenetic relationships between SNP haplotypes ( Fig. 1a; Supplementary Fig. 1), rooted by reference to two MSY sequences 13 from the basal haplogroups A and B. Our sequenced regions cover many previously known SNPs, which allowed us to apply established haplogroup names 20 to clades. Figure 1b shows the geographical distribution of these haplogroups in our samples, which is consistent with previous studies of specific SNPs using larger perpopulation sample sizes 10 . As expected, the commonest haplogroup is R1b-M269 (43.1%), with highest frequency in the north-west, followed by I1-M253 (13.8%), I2-P215 (9.0%), R1a-M198 (7.5%) and J2-M172 (7.5%). Some clades show geographically-restricted distributions, with hg N1c-M178 being most frequent in the Saami, and sub-lineages of haplogroups E, G and J prevalent in the Mediterranean area.
The shapes of different clades within the tree (Fig. 1a) vary greatly. Haplogroups E1b-M35, G2a-L31, I2-P215, J2-M172, L-M11 and T-M70 contain long branches with deep-rooting nodes, whereas I1-M253, N1c-M178, R1a-M198 and R1b-M269 show much shallower genealogies. Haplogroup R1b-M269 is particularly striking, containing a remarkable star phylogeny within which 44 terminal branches (13.2% of the total), found in 13 of the 17 sampled populations, descend as a multifurcation from a single node without any sub-structure whatsoever, despite the extensive nature of the sequencing carried out. These qualitative features of the phylogeny are supported by values of the average number of mutations from the ancestral node to branch tips, and also by estimates of time-to-most-recent-common-ancestor (TMRCA) ( Table 1) derived by two different methods. Considering haplogroups R1b-M269, R1a-M198 and I1-M253, and the 95% highest posterior density intervals of their TMRCAs, 64% of the MSY sequences sampled in our study descend from three ancestors who each lived more recently than B7.3 KYA.
Inferences on demographic history. Although it offers a tool to formulate hypotheses, the phylogeographic analysis above is limited in its power to illuminate demographic history. To further understand past demography, we applied a population approach: Table 2 shows diversity parameters for the 17 populations. When we consider the diversity from a molecular perspective (as the number of polymorphic sites) the highest diversity is in Turkey and Greece, closely followed by other southern populations (including Palestinians), and the lowest in Saami and Orkney. Consistent with this, there is a significant correlation of decreasing diversity both from south to north, and east to west (Supplementary Table 3), a clinal pattern that might be compatible with a model of demic diffusion from the Middle East. Considering instead the distribution of haplotypes, assessed as median number of singletons (here defined as variants that appear only once within a given population), by far the highest diversity is seen in Turkey and Greece, with the lowest diversity in the Saami and Palestinians-probably reflecting the effect of recent isolation and drift. Neither this measure, nor its standard deviation, shows any significant correlation with either latitude or longitude (Supplementary Table 3).
Bayesian skyline plots (BSPs) (Fig. 2) reveal the variation of effective population size with time 21 . The plots are consistent with patterns seen in the relative numbers of singletons, described above, in that the Saami and Palestinians show markedly different demographic histories compared with the rest, featuring very recent reductions, while the Turks and Greeks show evidence of general expansion, with increased growth rate around 14 KYA. A different pattern is seen in the remaining majority (13/17) of populations, which share remarkably similar histories featuring a minimum effective population size B2.1-4.2 KYA (considering the 95% confidence intervals (CIs) reported in Supplementary  Table 4), followed by expansion to the present. Considering only these 13 populations, the only significant geographical correlation is of decreasing diversity in the number of polymorphic sites from east to west (Supplementary Table 3); notably, there is no significant correlation between the age at which effective population size was at a minimum before expansion (Supplementary Table 4), and either latitude or longitude (Supplementary Table 3). Taken together, the very recent age of the demographic shift and its lack of geographical pattern suggest that its origin is distinct from that of the diffusion of agriculture.
In contrast to BSPs, approximate Bayesian computation (ABC) offers an alternative model-based approach to understanding complex population histories, based on coalescent simulations 22,23 . We analysed several demographic models including single population expansion, reduction or   Supplementary Fig. 2). Despite the high-molecular resolution of our data, the fact that they originate from a single locus did not allow precise estimation of demographic parameters, but the ABC generally confirmed the population size dynamics reconstructed by the BSP analysis (with some exceptions: Supplementary Tables 5-7). The proportion of the parameter variance explained by the summary statistics (R 2 ) is in most cases higher than 10% (and hence generally considered as a good estimation), but the 95% credible intervals are wide. This is particularly true for T1, the time of the start of the demographic change (reduction or expansion), thus preventing us from drawing any conclusion about the timing of these events (whether post Neolithic, or more ancient) from the ABC analysis.
In general, the non-parametric genealogical approach represented by BSPs better explores the variation found in our data, compared with a more conservative ABC analysis based on a single locus. We note that both analyses assume panmixia, and that population structure might influence effective population size estimates. However, it seems improbable that such an effect would extend to so many populations.

Discussion
Our approach has led to the confident identification of many MSY sequence variants in European population samples and a highly resolved phylogeny, but our conclusions are also influenced by a more contentious factor, the choice of mutation rate. We chose a rate (1.0 (0.92-1.09) Â 10 À 9 per bp per year, considering a 30-year generation time) based on the observation of 609 MSY mutations (excluding palindromic regions) in Icelandic deep-rooting pedigrees 24 . The point estimate of this rate is the same as an earlier pedigree-based estimate in which only four mutations were observed 25 , and which we applied in our broader MSY-phylogeny study 13 .  We note that the rate we have used is higher than the estimate of 0.76 (0.67-0.86) Â 10 À 9 per bp per year based on counting the 'missing' mutations in the genome of the Ust'-Ishim male 26 , radiocarbon dated to B45,000 YBP. Other studies 27,28 have inferred slower mutation rates (0.62 or 0.64 Â 10 À 9 ) based on scaling the genome-wide de novo rate to account for male-specific transmission, though this has been criticized 29 , and is not consistent with phylogenetic mutation rates estimated from human-chimpanzee MSY comparisons 30,31 . Some have chosen to calibrate the pedigree mutation rate against external (for example, archaeological) data 15,32 , but we have rejected this idea, firstly because of uncertainty over how archaeological date estimates correlate with demographic changes, and secondly because we have used a coalescent-based dating method that itself models genealogies. Despite recent advances, mutation rate remains a difficult issue, and more data are needed.
The recent and rapid continent-wide demographic changes we observe suggest a remarkably widespread transition affecting paternal lineages. This picture is confirmed in an independent analysis of MSY diversity in the pooled HGDP CEPH panel European samples 16 , and is compatible with current (n ¼ 98) ancient DNA data for MSY ( Fig. 3; Supplementary Table 8), in which hgs R1a, R1b and I1 are absent or rare in sites dating before 5 KYA, whereas hgs G2a and I2 are prevalent.
Analyses of ancient autosomal sequence data 19 demonstrate discontinuity 7-5 KYA between western European huntergatherers, tending to have genetic affinity to northern Europeans, and farmers, resembling southern Europeans. Consideration of the genomic ancestry of modern Europeans 33 reveals ancestry from these two groups, but also from north Eurasian hunter-gatherers 33 . Recent analysis 34 better defines this latter component, supporting a two-migration model into a hunter-gatherer substrate, involving an early-Neolithic (7-8 KYA) arrival of farming populations from the Near East, followed by a late-Neolithic (4.5 KYA) migration of pastoralists from the steppe region north of the Caspian Sea, whose genomic contribution is ubiquitous in modern Europeans. Ancient MSY sequences 34 show that hgs R1a and R1b are present in the steppe much earlier than observed in any European sites (Supplementary  Table 8), making this region a likely source for these MSY expansion lineages.
Ancient mtDNA data 18 also indicate large-scale population discontinuity since the Neolithic transition, with a massive shift in haplogroup composition B7.5 KYA between Central European hunter-gatherers (carrying exclusively hg U lineages) and farmers (a much broader range of hgs), followed by later fluctuations. Demographic inference from whole mtDNA sequences 16 , however, does not show recent and sudden expansion. This suggests that the recent events responsible for shaping modern MSY variation were male specific.
The period 4-5 KYA (the Early Bronze Age) is characterized by rapid and widespread change, involving changes in burial practices that might signify an emphasis on individuals or kin groups, the spread of horse riding, and the emergence of elites and developments in weaponry 35 . In principle male-driven social selection 36 associated with these changes could have led to rapid local increases in the frequencies of introgressing haplogroups 34 , and subsequent spread, as has been suggested for Asia 37 . However, cultures across Europe remain diverse during this period; clarifying the temporal and geographical pattern of the shift will rely heavily on additional ancient DNA data.

Methods
Samples. DNA donors were recruited with informed consent (University of Leicester Research Ethics Committee reference: maj4-cb66). DNA was extracted from various sources including lymphoblastoid cell lines, peripheral blood and saliva. Samples (340) were included in the design from 17 populations (20 males each) across Europe and the Near East. Samples from Greece, Serbia, Hungary, Germany (Bavaria), Spanish Basque country, central Spain, Netherlands (Frisia), Denmark, Norway, Finland (Saami), England 38 (Herefordshire and Worcestershire), Orkney 38 , Ireland and Turkey were collected by the authors. Twenty random Palestinian male samples were purchased from the National Laboratory for the Genetics of Israeli Populations (www.tau.ac.il/medicine/NLGIP). Finally, samples from two HapMap 39 populations were used, both to supplement the population data set and to provide data on externally analysed samples for validation purposes: the Centre d'Etude du Polymorphisme Humain (CEPH) collection in Utah, USA, with ancestry from Northern and Western Europe (CEU) and the Toscani in Italia (TSI). After the initial analyses, one English and one Spanish individual were identified as females and therefore removed from all downstream analyses in this study, reducing the final number of samples to 338. For further details on samples see Supplementary Table 1.
In this study we focus on the eight X-degenerate regions 31 of the Y chromosome which are likely to yield interpretable sequence data; other captured regions are discussed elsewhere 13 . The total length of targeted regions was B8.6 Mb, and following capture design and the necessary repeat masking, the designed baits covered 2.3 Mb. Coordinates of the eight targeted regions can be found in Supplementary Data 2.
Sequencing and data processing. Genomic DNA (3-5 mg) was used for library preparation and target enrichment using Agilent SureSelect XT Target Enrichment System for Illumina Paired-End Sequencing Library kit (version 1.3). In order to obtain larger insert sizes, DNA samples were fragmented to B250-600 bp without size selection. This resulted in a mean insert size of 330 bp, which increases recovery of sequence data from bait-adjacent regions. Sequencing was done on an Illumina HiSeq 2000 instrument (Illumina, CA, USA) with paired-end 100-bp run to high coverage. Library preparation, target enrichment and sequencing were carried out at the High-Throughput Genomics Centre at the Wellcome Trust Centre for Human Genetics, University of Oxford, UK.
Base calling was done using Illumina Bustard 40 and quality control with FastQC 41 . Sequence data were mapped to the human genome reference (GRCh37) using Stampy v1.0. 20 (ref. 42). Local realignment was done using The Genome Analysis Toolkit (GATK) v2. 6-5 (ref. 43), followed by duplicate read marking with Picard v1.86 (ref. 44) and base quality score recalibration also with GATK. The individual steps and parameters used are listed in Supplementary Table 9.
Variant calling and filtering. Owing to larger insert sizes (see above) and high efficiency of sequence capture, high sequence coverage was obtained not only at baited regions but also at B300 bp flanking the enrichment baits. Therefore, the ARTICLE original bait coordinates were modified by adding 300 bp to either side of each bait followed by merging the overlapping coordinates and increasing the size of the analysed region to 4,433,580 bp. Data on the 338 male samples described above were co-analysed with simultaneously generated data on an additional 117 samples, described elsewhere 13 . Variant calling was done using the samtools mpileup v0.1.19 multi-sample option, calling all samples simultaneously (Supplementary Table 9). In total 19,276 raw SNPs were called from 455 male samples.
Raw variants were filtered using vcftools v0. 1.11 (ref. 45) and in-house Perl scripts. Filters used for the final data are listed in Supplementary Table 9. As well as the two females, seven samples (four from the European population set) were removed from the final data set due to missing 45% of calls. The filtered data set included 13,474 variant sites from 448 samples, from 0 to 643 missing calls per individual, with an average call rate of 99.8%.
To recover as many missing genotypes as possible for subsequent analyses, they were divided into three groups based on read-depth: DP 0-the genotype call was discarded; DP 2-6-the raw call was accepted; all other cases-the sites were re-called using a single-sample approach to obtain the DP4 field in the vcf, the bam file was checked manually, and the most probable allele was inferred by comparing the bam file with the information contained in the DP4 field. After this procedure, 213/13,474 sites still lacked genotype calls, leading to a final number of 13,261 sites for further analyses.
Having applied the above filters to variant sites, it was necessary to apply the same criteria to non-variant sites. We calculated the depth per sample per site using the GATK DepthOfCoverage tool, filtered for base quality 20 and mapping quality 50, and then applied the criterion of Z6 Â coverage in Z95% of samples. This led to a reduction in the figure of base pairs sequenced from 4,433,580 to 3,724,156 bp. Validation. In silico validation of the 13,261 filtered SNP calls was done using two previously published data sets: genomes sequenced to high-coverage with selfassembling DNA nanoarrays by Complete Genomics 46 , and Omni2.5 BeadChip genotype data produced at the Broad Institute as part of the 1,000 Genomes Project 47 . Our samples included 4 and 39 HapMap individuals overlapping with these two data sets, respectively. A Perl script was written to compare the SNP calls in our variant set to overlapping samples and positions in the control sets.
Of the 888 variant sites shared between our data and the Complete Genomics data across four overlapping samples, the false positive and false negative error rates were both 0%. When compared with the Omni data across 241 variant sites and 39 overlapping samples, the error rates were 0.13% for false positives and 1.82% for false negatives. However, all the false calls originated from only 19 variant sites.
To shed light on these comparatively high error rates, we also compared Complete Genomics and Omni data for regions corresponding to our final analysed regions. Across 263 variant sites and 49 overlapping samples, we obtained false positive and false negative rates of 2.85 and 2.36%, respectively. The false calls originated from 30 sites and 15 of those overlap with the sites producing high error rates when comparing our data with Omni. Since the Complete Genomics data set is generally considered to have very high quality then this seems to indicate problems in making correct calls from Omni genotyping data. More detail is provided in Supplementary Table 2.
Phylogenetic inference. PHYLIP v3.69 was used to create a maximum parsimony phylogenetic tree 48 . Three independent trees were constructed with dnapars using randomization of input order (seeds: 4941985, 62529981 and 38185313), each 10 times. Output trees of these runs were used to build a consensus tree with the consense programme included in PHYLIP package.
The tree was rooted using two Y-chromosomes belonging to haplogroups A and B which were sequenced in the complete data set 13 . FigTree v1.4.0 49 was used to visualize the tree (tree.bio.ed.ac.uk/software/figtree/).
Haplogroup prediction. The presence of known markers was checked using AMY-tree v1. 2 (refs 50,51). This software was developed to predict MSY haplogroups from whole-genome sequence data using a list of known markers. Since our data do not cover the whole MSY but only a proportion of it, the software lacks sufficient information for haplogroup prediction. However, it can be used to deduce the presence and allelic states of known MSY markers present in sequence data. The AMY-tree v1.2 conversion file contains a list of 1,453 known Y-SNPs, of which 490 are present in our data. These 490 sites were used to assign a standard haplogroup to all our samples according to the Y Chromosome Consortium phylogenetic tree 20 and its subsequent updates (Supplementary Table 1).
TMRCA and ages of nodes. The TMRCA of the tree and of nodes of interest were estimated via two approaches: BEAST v1. 8 (refs 17,52): Markov chain Monte Carlo (MCMC) samples were based on 25,000,000 generations, logging every 1,000 steps, with the first 2,500,000 generations discarded as burn-in. Three runs were combined for analysis using LogCombiner. We used an exponential growth coalescent tree prior (growth rate prior: uniform(0-0.002)), HKY substitution model, and a strict clock with a substitution rate of 1.0 (95% CI: 0.92-1.09) Â 10 À 9 mutations/nucleotide/year 24 . TMRCAs were estimated in a single run including all 17 European populations and assigning samples to specific clades in agreement with the MP tree shown in Fig. 1. Rho: A Perl script was written to calculate TMRCA and its standard deviation for any given clade within a PHYLIP outfile, using the rho statistic 53,54 . A scaled mutation rate of 268.5 (246.3-291.9) years per mutation was used, based on a published rate of 1.0 (95% CI: 0.92-1.09) Â 10 À 9 mutations/ nucleotide/year 24 and the number of nucleotides in our regions of interest (3,724,156).
Bayesian skyline plots. Bayesian skyline plots were generated using BEAST v1.8 (refs 17,52). MCMC samples were based on 100,000,000 generations, logging every 1,000 steps, with the first 10,000,000 generations discarded as burn-in. We used a piecewise linear skyline model with 10 groups, a HKY substitution model, and a strict clock with a mean substitution rate of 1.0 Â 10 À 9 mutations/nucleotide/ year 24 and a generation time of 30 years, consistent with 55 . For the 13 populations showing a recent expansion in the BSP, the limits of the 95% CI of mutation rate 24 (0.92-1.09 Â 10 À 9 ) were used to define the CI of the time estimate of the minimum effective population size before the expansion (grey shading in Fig. 2).
Intrapopulation diversity and geographical correlation. The number of polymorphic sites per population, Tajima's D 56 , and Fu's FS 57 were calculated using Arlequin 3.5 (ref. 58). The number of singletons was calculated using Vcftools v0.1.11. Correlation tests between measures of genetic diversity and latitude and longitude were run in R 59 with the function cor.test of the package stats.
Approximate Bayesian computation. We generated one million simulated datasets for each tested model ( Supplementary Fig. 2) with the programme FastSimcoal2 (ref. 60), simulating a single haploid locus of 3,724,156 bp. We summarized the data by means of the derived site frequency spectrum (-s -d flags in the command line) considering only categories with at least one observed polymorphic site. Ancestral states in the observed data were defined elsewhere 13 using a custom script.
To compare models we applied the Logistic Regression procedure 61 . Model parameters were estimated by a locally weighted multivariate regression 22 after a logtan transformation 62 of the 10,000 best-fitting simulations from a specific model. To calculate the posterior probabilities for models and parameters we used R scripts from http://code.google.com/p/popabc/source/browse/ #svn%2Ftrunk%2Fscripts, modified by AB and SG.
We also estimated the power of our ABC procedure to correctly recognize the true model calculating for each model the proportion of true positives and false positives. We evaluated 1,000 pseudo-observed data sets generated under each model, counting the number of times a specific model is correctly identified by the ABC procedure (true positives), and the number of times the same model is incorrectly selected as the true model (false positives).
Demographic models and priors. We considered five models depicting different demographic histories, testing each model separately for each population ( Supplementary Fig. 2). M1 is the simplest, in which the effective population size remains constant over time (uniform prior: 20-20,000). In M2 an ancient constant-sized population (uniform prior: 1,001-20,000) starts an exponential reduction T1 generations ago. The reduction spans LEX generations (uniform prior: 5-634), then the population returns to constant size (uniform prior: 20-1,000) T2 generations ago (uniform prior: 0-30). In M3 the reduction is followed by an expansion that starts T2 generations ago (with the effective population size, NER, drawn from an uniform prior: 20-1,000) until the present (uniform prior for the current effective population size, NC: 1,001-20,000). M4 is parameterized in the same way as M2, with an expansion instead of a reduction (NA uniform prior: 20-1,000; NC uniform prior: 1,001-20,000). In M5, the expansion ends at T2 followed by a reduction until present time (NEE uniform prior: 1,001-20,000; NC uniform prior: 20-1,000). We considered a generation time of 30 years. In all the models the Last Glacial Maximum (B20,000 years ago) represents the upper bound of the time for the first demographic change (T1).
In each simulation the per-generation, per-site mutation rate 24 is drawn from a normal distribution with mean 3.01 Â 10 À 8 and 95% confidence intervals 2.77-3.26 Â 10 À 8 . DNA sequences were generated under a finite sites mutational model with no transition/transversion bias.
Perl scripts used in the analysis are available upon request.