A Population Genomics Analysis of the Native Irish Galway Sheep Breed

The Galway sheep population is the only native Irish sheep breed and this livestock genetic resource is currently categorised as ‘at-risk’. In the present study, comparative population genomics analyses of Galway sheep and other sheep populations of European origin were used to investigate the microevolution and recent genetic history of the breed. These analyses support the hypothesis that British Leicester sheep were used in the formation of the Galway. When compared to conventional and endangered breeds, the Galway breed was intermediate in effective population size, genomic inbreeding and runs of homozygosity. This indicates that, although the Galway breed is declining, it is still relatively genetically diverse and that conservation and management plans informed by genomic information may aid its recovery. The Galway breed also exhibited distinct genomic signatures of artificial or natural selection when compared to other breeds, which highlighted candidate genes that may be involved in production and health traits.


INTRODUCTION
Sheep were domesticated more than 10,000 years ago and have since been bred for a variety of uses including meat, milk and wool production (Taberlet et al., 2011;Larson and Fuller, 2014;MacHugh et al., 2017). During the last 50 years, the focus of the global sheep industry on only a subset of the 1,400 recorded sheep breeds with enhanced productivity and high-quality outputs has resulted in many locally adapted (local) breeds becoming endangered or extinct (Taberlet et al., 2008;Kijas et al., 2009;Kijas et al., 2012). These breeds are generally considered independent genetic units because crosses are usually not used for further reproduction (Taberlet et al., 2008). Local or heritage livestock breeds are important because they constitute reservoirs of biological diversity different to the major production breeds and that may be important genetic resources for domestic animal species in the face of climate change and increased food requirements in the future (Taberlet et al., 2008;Bowles, 2015). To address these future challenges, it will be possible to use targeted genome editing technologies in livestock. Consequently, functionally important natural sequence variants (NSVs) identified in the genomes of locally adapted native and heritage breeds may become increasingly important for genetic improvement programmes (Wells, 2013;Petersen, 2017;Van Eenennaam, 2017).
The local sheep breeds on the periphery of Northern Europe are recognised as heritage livestock populations that should be conserved and represent important sources of novel genetic diversity accumulated over centuries of microevolution and adaptation to marginal agroecological environments (Tapio et al., 2005). In this regard, the Galway sheep breed is the only surviving sheep breed native to Ireland (Curran, 2010); it was once the principal lowland sheep breed in Ireland but is now considered at-risk by the Food and Agriculture Organization (Food and Agriculture Organization, 2019). The Galway breed therefore represents a useful reservoir of genetic variation for domestic sheep and should be conserved.
The Galway breed is thought to have originated as a composite of indigenous and imported sheep populations, present in Ireland in the mid-19th century, through the breeding endeavours at that time, which were concerned mainly with improved wool production (Hanrahan, 1999). Sheep breeds in Ireland during this period include the important Dishley or New Leicester foundational breed developed by Robert Bakewell (Wykes, 2004). However, it was not until 1923 that a formal Galway herd book was established (Curran, 2010;Food and Agriculture Organization, 2019). Therefore, the range of sheep populations ancestral to the Galway breed in the 18th and 19th centuries, coupled with the possibility of more recent gene flow, poses questions concerning the genetic distinctiveness and admixture history of the breed. In addition, the Galway breed has declined from a peak population size in the 1960s when it was the focus of lowland sheep farming in Ireland (Martin, 1975a;Raftice, 2001;Curran, 2010). By 1994, as defined by the UK Rare Breeds Survival Trust, the Galway breed had reached 'critical' status for sheep breeds with only 300 pedigree breeding ewes registered (Curran, 2010). Since being classed as endangered by the Irish Government in 1998, the number of pedigree Galway sheep has increased due to conservation efforts; however, the breed population size is currently decreasing, raising concerns regarding remaining genetic diversity and the overall viability of the population (Curran, 2010;Food and Agriculture Organization, 2019).
As a local breed with a low census population size, the main threat to the long-term survival of the Galway breed is replacement by more productive commercial breeds, which would further reduce the population size, reduce genetic diversity and increase inbreeding. Other challenges faced by threatened local livestock breeds include poor animal husbandry and management, deliberate or inadvertent crossbreeding and geographical isolation, which increases the risk of extinction (Taberlet et al., 2008;Allendorf et al., 2013). In recent years, with the availability of increasingly powerful genomics technologies, a conservation programme for Galway sheep has been proposed that would leverage molecular genetic information (McHugh et al., 2014). McHugh and colleagues also propose that genome-enabled breeding (genomic selection) could be used in threatened livestock populations to improve production, health and reproduction traits, thereby decelerating replacement by modern breeds (Biscarini et al., 2015). Another strategy could leverage multi-breed or across-breed genomic prediction (Iheshiulor et al., 2016). This approach can increase the accuracy of genomic estimated breeding values for small populations such as the Galway breed, since accurate genomic selection requires large numbers of phenotyped and genotyped animals (Iheshiulor et al., 2016).
To provide information that may be relevant to genetic conservation of the Galway sheep breed, we performed highresolution population genomics analyses in conjunction with 21 comparator breeds of European origin. These analyses included multivariate analyses of genomic diversity, phylogenetic network graph reconstruction, evaluation of genetic structure and inbreeding, modelling of historical effective population sizes and functional analyses of artificial and natural selection across the Galway sheep genome.

Galway and Irish Suffolk Sheep DNA Sampling
The Galway and Irish Suffolk sheep DNA samples used for the current survey were generated from peripheral blood samples collected in standard heparinised Vacutainer blood collection tubes (Becton-Dickinson Ltd., Dublin, Ireland). High-quality genomic DNA was then purified from 200 µl of blood from each animal using standard laboratory methods (Howard, 2008). The 49 Galway sheep were sampled from 14 different flocks and pedigree information was consulted to minimise relationship among the animals sampled. The sample size breakdown across the 14 flocks in order of decreasing size is: 6, 6, 5, 5, 5, 5, 3, 3, 3, 3, 2, 1, 1, 1. The flocks were geographically dispersed across County Galway in western Ireland (Howard, 2008). The 55 Irish Suffolk sheep were sampled in approximately equal numbers from two experimental flocks maintained by University College Dublin and Teagasc, the Agriculture and Food Development Authority of Ireland (Howard, 2008).

Additional SNP Data Sources and Data Filtering
Medium-density SNP data were obtained from the International Sheep Genomics Consortium Sheep HapMap Project and consisted of 2,819 sheep from 74 breeds genotyped for 49,034 evenly spaced SNPs using the Illumina ® OvineSNP50 BeadChip (Kijas et al., 2012). To focus on the Galway breed, a core sample set of 11 breeds, including the Galway breed, was selected for the primary population genomic analyses (n = 615 animals). This included populations previously examined and known to be more closely related due to their shared European origins (Howard, 2008;Kijas et al., 2012). These comparator populations also included widely used breeds, such as the Merino (MER), and at-risk heritage breeds, such as the Dorset Horn (DSH), Soay (SOA) and Wiltshire (WIL) (Food and Agriculture Organization, 2019). Figure 1 and Supplementary Table 1 provide further information on the geographical origins of the 11 breeds used for the core sample set analyses. In addition, Supplementary Table 1 provides information on an expanded sample set of 22 European and Asian breeds, including the core sample set, used for the phylogenetic tree and network graph reconstructions (n = 1,003).
The initial data set had already been filtered to remove SNPs with <0.99 call rate, assay abnormality, MAF <0.01, discordant genotypes and inheritance problems (Kijas et al., 2012). The core and extended sample genome-wide SNPs data sets for this study were filtered using PLINK v1.07 (Purcell et al., 2007) to remove SNPs lacking positional information, SNPs unassigned to any chromosome, or SNPs assigned to the X and Y chromosomes (Patterson et al., 2006;Purfield et al., 2012). The final filtered data set was composed of 47,412 SNPs with a total genotyping rate of 99.7%.

Principal Component Analysis
Principal component analysis (PCA) was performed using 47,412 genome-wide SNPs and SMARTPCA from the EIGENSOFT software package (version 4.2) (Patterson et al., 2006). The number of autosomes was set to 26 and breed names were included. The number of outlier removal iterations was set to 0 since outliers could flag individual animals that were the result of crossbreeding. PCA plot visualisations were generated using ggplot2 (Wickham, 2016).

F ST Analysis
Pairwise F ST values (Weir and Cockerham, 1984) were calculated for each pair of breeds using 47,412 genome-wide SNPs and PLINK v1.9 (Chang et al., 2015). Weighted values were chosen to account for different sample sizes for each breed (Weir and Cockerham, 1984).

Construction of Phylogenetic Trees and Ancestry Graphs
Maximum likelihood (ML) phylogenetic trees with ancestry graphs were generated for the core and extended sample data sets using 47,412 genome-wide SNPs and the TreeMix (version 1.12) software package. For the core sample set, the Italian Comisana breed (COM) (Ciani et al., 2014) was used as an outgroup and five migration edges were used for TreeMix visualisation (Pickrell and Pritchard, 2012). The analysis was repeated using the extended sample set of 21 European breeds (Supplementary Table 1) and the Indian Garole breed (GAR) was used as an outgroup, again with five migration edges for TreeMix visualisation.

Genetic Structure and Admixture history
Genetic structure and admixture history was investigated for the core sample set of the Galway and 10 other breeds using 47,412 genome-wide SNPs and fastSTRUCTURE (version 1.0) (Raj et al., 2014) as described previously by us (Browett et al., 2018). The analysis was performed with the model complexity, or number of assumed populations, K = 2 to 11. The simple prior approach described by Raj et al. (2014) was used, which is sufficient for modelling population/breed divergence. The 'true' K-value for the number of ancestral populations was estimated using a series of fastSTRUCTURE runs with pre-defined K-values that were examined using the chooseK.py script (Raj et al., 2014). Outputs from the fastSTRUCTURE analyses were visualised using the DISTRUCT software program (version 1.1) with standard parameters (Rosenberg, 2004).

Modelling of Current and historical effective Population Size
Current and historical effective population size (N e ) trends were modelled with genome-wide SNP linkage disequilibrium data from 47,412 genome-wide SNPs for the core sample set using the SNeP software tool (version 1.1) (Barbato et al., 2015) implementing the method for unphased SNP data as described previously by us (Browett et al., 2018). Graphs used to visualise trends in N e were generated using ggplot2 (Wickham, 2016).

Analysis of Genomic Inbreeding and Runs of homozygosity
Analysis of genomic inbreeding based on the inbreeding coefficient (F) estimated from SNP heterozygosity data was performed using 47,412 genome-wide SNPs and the PLINK v1.07 -het command (Purcell et al., 2007) since comparable inbreeding results have been observed using pruned or unpruned data for a SNP data set of similar size (Binns et al., 2012). Runs of homozygosity (ROH) are continuous tracts of homozygosity that most likely arise due to inbreeding and can be identified through surveys of genome-wide SNP data in populations (Curik et al., 2014;Peripolli et al., 2017). Individual animal genomic inbreeding was evaluated as genome-wide autozygosity estimated from the SNP data using runs of homozygosity (ROH) values generated with PLINK v1.07 (Purcell et al., 2007) and the F ROH statistic introduced by McQuillan et al.
(2008) with methodologies previously described in detail by Purfield et al. (2012) and Browett et al. (2018). The F ROH statistic represents the proportion of each individual animal's genome covered by ROH, which is generally a consequence of historical inbreeding. Statistical analysis was carried out in R and graphs used to visualise F, F ROH and ROH distributions were generated using ggplot2 (Wickham, 2016;R Core Team, 2018).

Genome-Wide Detection of Signatures of Selection and Functional enrichment Analysis
The composite selection signal (CSS) method (Randhawa et al., 2014) was used to detect genomic signatures of selection as previously described (Browett et al., 2018). The CSS approach combines the fixation index (F ST ), the directional change in selected allele frequency (ΔSAF) and cross-population extended haplotype homozygosity (XP-EHH) tests into one composite statistic for each SNP in a population genomics data set (Randhawa et al., 2014). For the present study, we used 47,412 genome-wide SNPs genotyped in 49 Galway sheep (GAL) and a sample of 50 randomly selected sheep (5 selected at random from each of the other 10 breeds in the core data set). To mitigate against false positives, genomic selection signatures were only considered significant if at least one SNP from the set of the top 0.1% genome-wide CSS scores was flanked by at least five SNPs from the set of the top 1% CSS scores.
As described previously (Browett et al., 2018), the Ensembl BioMart data mining resource (Smedley et al., 2015) was used to identify genes within ±1.0 Mb of each selection peak (Ensembl release 85, July 2016). Ingenuity ® Pathway Analysis (IPA ® : Qiagen, Redwood City, CA, USA; release date July 2016) was then used to perform an overrepresentation enrichment analysis with this gene set to identify canonical pathways and functional processes of biological importance. The total gene content of Ensembl release 85 version of the OAR3.1 ovine genome assembly (Jiang et al., 2014) was used as the most appropriate reference gene set for these analyses (Timmons et al., 2015).

Analyses of Breed Divergence, Genetic Differentiation and Admixture
The results of multiple population genomics analyses support the genetic distinctiveness of the Galway sheep population as a discrete breed. The PCA results plotted in Figure 2 demonstrate separation of the majority of breeds into distinct population clusters, with the notable exceptions of the Australian Merino (MER) and Scottish Blackface (SBF). However, it is important to note that the PCA plot visualisation shown in Figure 2 did not include the 110 samples from the Soay breed (SOA). A long history as a relatively small isolated island population (Berenos et al., 2016) has led to a marked pattern of genetic differentiation from other breeds, which is evident in the first principal component (PC1) of Supplementary Figure 1. Consequently, when the Soay breed is included in a PCA, PC3 is required to separate the Galway breed from the other populations (Supplementary Figure 2). Otherwise, the Galway breed clusters with the Scottish Texel breed (STX) and is located close to the Border Leicester breed (BLR). This result supports the documented role for the foundational New Leicester breed in the formation of the Galway and Texel breeds (Porter et al., 2016) and is compatible with the results of a previous study using autosomal microsatellites (Howard, 2008).
The PCA plot shown in Figure 2 also demonstrated that a number of individual sheep do not cluster closely with other animals from their breeds. This is likely due to recent unacknowledged or inadvertent crossbreeding between animals from different populations (Patterson et al., 2006) or, alternatively, potential mislabelling of particular samples. For example, the 2D and 3D PCA plots shown in Supplementary Figures 1 and 2 indicate that one of the Irish Suffolk animals (ISF25) was most likely a mislabelled Scottish Texel sample as it emerged within the main Texel cluster for PC1, PC2 and PC3. Consequently, this sample ISF25 was removed from all subsequent analyses.
The PCA results are supported by the interpopulation weighted F ST values for each pair of breeds shown in Supplementary  (Curran, 2010;Porter et al., 2016;Food and Agriculture Organization, 2019).
The ML phylogeny and ancestry graph in Figure 3 shows that the Galway breed groups closely with sheep populations of English and Dutch origin, particularly the Border Leicester (BRL) and the Scottish Texel (STX) breeds. This observation is concordant with previous population genomics studies (Kijas et al., 2012;Fariello et al., 2013) and known breed histories due to the shared historical input of the foundational New Leicester breed (Curran, 2010). The ML phylogeny and ancestry graph generated with additional European breeds and shown in Supplementary Figure 5 also supports the close relationship among the Galway, BRL and STX breeds. The arrows (graph edges) on Figure 3 indicate gene flow modelled between populations with the colour scale representing the weight of each migration event.
Results of the genetic structure analysis for individual animals grouped by population are shown in Figure 4. Model complexity or numbers of assumed populations (K) ranging from 2 to 11 are visualised to explain the structure in the data and to maximise the marginal likelihood. These results demonstrate that the 11 breeds can be considered discrete populations, thereby supporting interpretation of sheep breeds as separate genetic units (Taberlet et al., 2008) Table 2) (Porter et al., 2016).  Table 4 shows that the modelled historical trends in N e for the 11 breeds analysed decline towards the present. However, the GAL breed are intermediate between the breeds with large census populations (FIN, ISF, MER, ROM, SBF and STX) and at-risk breeds with relatively small census populations (BRL, DSH, SOA, WIL) breeds. In addition, the most recent modelled N e value for the GAL breed is 184 animals 13 generations ago, which is comparable to some of the breeds (e.g. ISF and STX with 178 and 150 animals, respectively). These modelled N e values, which are based on linkage disequilibrium, may be underestimates due to the physical linkage between many SNPs (Hall, 2016).

Modelling historical effective Population Size
To examine these historical trends in N e more systematically, the data for each breed were shown to be not normally distributed using the Kolmogorov-Smirnov test (Supplementary Table 3).
Therefore, the non-parametric general Kruskal-Wallis test followed by pairwise Wilcoxon rank sum tests for all population/breed comparisons with adjustment for multiple statistical tests performed with the Bonferroni correction. This analysis demonstrated that the GAL historical N e trend is significantly different only from the MER breed (P adj. = 0.006; Supplementary Table 5). Livestock populations tend to exhibit lower N e values than comparable wild mammal populations (Waples et al., 2016). Notwithstanding this, from a conservation perspective, it is reassuring that the most recent estimated N e value of 184 for the GAL is above the critical threshold of 100 animals considered essential for the long-term survival of livestock populations (Meuwissen, 2009). This 'demographic fingerprint' (Barbato et al., 2015) is most likely a consequence of the widespread use of the Galway breed for lowland sheep production in Ireland up until the 1980s (Raftice, 2001;Curran, 2010).

Genomic Inbreeding and Runs of homozygosity
The recent N e of each of the sheep breeds modelled in Figure 5 will have been substantially influenced by their inbreeding histories. In this regard, the genomic inbreeding coefficient (F) values estimated for individual animals across all breeds range up to 0.389 for a single Dorset Horn (DSH) animal (Figure 6). The majority of F values for individual animals in each breed were not normally distributed based on Shapiro-Wilk test results (Supplementary Table 3); therefore, the median F values were generated and evaluated for each breed (Supplementary Table 6). The breeds with the highest median F values were the SOA (0.308) and the WIL (0.299) and the two breeds with the lowest median (SOA) have existed as a relatively small and isolated population on the island of Soay for hundreds of years while the Wiltshire breed (WIL) has recently experienced a dramatic decline in census population and is considered at risk by the FAO (Food and Agriculture Organization, 2019). From a genetic conservation perspective, except for a single outlier (GAL26), it is encouraging that the Galway breed (GAL) exhibits an intermediate median F value calculated using genome-wide SNP data. A systematic analysis of F value distributions using the nonparametric Kruskal-Wallis test indicated there were significant differences among breeds (H = 477.33, df = 10, P < 0.001). An analysis of all pairwise breed comparisons using the non-parametric Wilcoxon rank sum test (with Bonferroni correction) was then performed (Supplementary Table 8). These results showed that the majority of pairwise comparisons were highly significant, again reflecting the distinct demographic histories of each breed.
Overall, comparable results to those obtained using the genomic inbreeding coefficient (F) were observed for inbreeding coefficients estimated using ROH (F ROH ) (Figure 7, Supplementary Tables  3, 6, 7 and 9). However, there were some notable differences; in particular, the lower median F ROH value of 0.101 for the Soay breed (SOA) is likely due to their longer geographical isolation and a Frontiers in Genetics | www.frontiersin.org October 2019 | Volume 10 | Article 927 consequence of early historical inbreeding that produced ROH tracts, which have broken down due to recombination (Barrett, 2012;Purfield et al., 2012). It is also notable that the Galway breed contains several individual animals with higher F ROH values (GAL15, GAL16, GAL18, GAL26 and GAL36) indicating that this statistic is useful for identifying animals that should not be prioritised for conservation programmes. With regards to historical inbreeding in the Galway breed (GAL), inbreeding coefficients have previously been calculated using pedigree information for the population in 1969 (F = 0.019; Martin, 1975b), 1999 (F = 0.020; Raftice, 2001) and 2012 (F = 0.023; McHugh et al., 2014). These results indicate that the general trend in inbreeding has been relatively moderate, which may also be reflected in the results obtained using genomic information reported in the present study. It is important to note that monitoring of inbreeding for genetic conservation and management of potentially deleterious recessive genomic variants can be greatly informed through evaluation of ROH parameters using SNP data (Peripolli et al., 2017). The mean sum of ROH for different length categories varies among the breeds (Figure 8); however, none of the breeds exhibit large mean values for the total length of ROH in the 1 to 5 Mb category. This is because the SNP density on the OvineSNP50 BeadChip is too low to accurately detect ROH in this size range and may not accurately estimate F ROH when short segments are included (Supplementary Table 7) (Purfield et al., 2012;Ferenčaković et al., 2013). Notwithstanding this limitation, patterns of ROH, which reflect both recent and older inbreeding histories, are evident. For example, the Wiltshire breed (WIL) has large mean total ROH lengths for the other categories, presumably reflecting both historical and recent inbreeding. Other breeds, such as the Australian Merino (MER), have smaller mean total lengths of ROH in all categories, an observation that is concordant with the results of the genomic inbreeding and the analysis of N e estimates. This is because individual animals from breeds with larger effective population sizes-such as the Australian Merinoare less likely to be the result of inbreeding and are therefore less likely to contain large ROH segments in their genomes (Curik et al., 2014;Peripolli et al., 2017). The converse of this is true for breeds with lower N e values and large ROH tracts in their genomes, such as the endangered Wiltshire breed. In terms of mean total Frontiers in Genetics | www.frontiersin.org October 2019 | Volume 10 | Article 927 FIGURe 5 | Trends in effective population size (N e ) estimated using 47,412 genome-wide SNPs. length of ROH, the Galway breed emerges between these extremes, reflecting an intermediate effective population size and history of moderate inbreeding (Figure 8). In conjunction with the other analyses of genomic diversity, these results are also encouraging for genetic conservation and the long-term viability of the breed.

Signatures of Selection in the Galway Sheep Breed
Using defined criteria, five significant peaks of selection were detected with the CSS approach (Figure 9): two on OAR1, one on OAR3 and two on OAR8 (that merge into one peak on the graph). Each selection peak was located in a ROH tract detected in at least three Galway samples, which may indicate reduced genetic diversity in these regions as a consequence of localised selective sweeps (Purfield et al., 2017). Detection of these selection peaks demonstrates that the Galway population has experienced a unique history of both natural and humanmediated selection, presumably because of adaptation to the agroecology of Ireland, a large Northwestern European island with a temperate oceanic climate. The precise locations of the peaks that have clusters of SNPs within the top 0.1% CSS score class are provided with additional information in Supplementary Table 10. The 197 genes within these regions are listed in Supplementary Table 11. Using IPA ® , the top five physiological system development and function pathways enriched for the subset of 119 genes that could be mapped to HGNC symbols were identified and are listed in Table 1 (Krämer et al., 2013).
Of the 119 candidate genes hypothesised to be under selection in the Galway breed, 28 are involved in tissue development and 15 are involved in connective tissue development and function. This is a common observation in studies of selection across the genomes of livestock populations (de Simoni Gouveia et al., 2014;Gutiérrez-Gil et al., 2015;Randhawa et al., 2016). Seven of the 119 genes are involved in hair and skin development and function, which may be explained by the use of Galway sheep in wool production (Curran, 2010). Selection and maintenance of traits that confer resilience to infectious disease is important in domestic animal populations, including many sheep breeds (Bishop and Woolliams, 2014;Bishop, 2015). Thirteen of the 119 genes under the selection peaks are involved in immune cell trafficking, which may be as a result of the climate and unique disease challenges posed by the Irish environment, such as the prevalence of liver fluke (Toolan et al., 2015). A large group of 26 genes enriched for haematological system development and function were also located under the selection peaks; however, a microevolutionary explanation for this is not hypothesised here.

Genetic Conservation of the Galway Sheep Breed
The results of the population genomics analyses presented here are mutually consistent and highlight the utility of Frontiers in Genetics | www.frontiersin.org October 2019 | Volume 10 | Article 927 dense genome-wide marker data for conservation genomics in livestock populations, particularly for at-risk heritage landrace populations such as the Galway breed. Our results show the Galway breed is genetically distinct from other European sheep breeds, emerging in multivariate PCA and phylogenetic tree network graph visualisations as a distinct group but close to the Border Leicester breed (BRL), which has been observed previously (Kijas et al., 2012). In terms of effective population size and genomic inbreeding, the Galway breed emerged as intermediate between non-endangered and endangered sheep breeds. This indicates that there is substantial genetic diversity remaining in the population, which could be managed with a conservation programme that is informed by genomic information.

DATA AVAIlABIlITY STATeMeNT
The Galway sheep (GAL) and additional sheep breed Illumina ® OvineSNP50 BeadChip data are available as part of the International Sheep Genomics Consortium Ovine SNP50 HapMap Dataset (www.sheephapmap.org/download.php).

eThICS STATeMeNT
Animal biological sample collection was conducted under license issued in accordance with Irish and European Union legislation (Cruelty to Animals Act, 1876, and European Community Directive, 86/609/EC) as described previously (Mullen et al., 2013). All animals were managed in accordance with the guidelines for the accommodation and care of animals under Article 5 of the Directive.

AUThOR CONTRIBUTIONS
DEM, DH, MM and JH conceived and designed the project. DH, MM, DAM, ES and JH organised sample collection and