{'fraction_non_alphanumeric': 0.032360339587045386, 'fraction_numerical': 0.04501624567655382, 'mean_word_length': 4.930846930846931, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Recovering high-quality host genomes from gut metagenomic data through genotype imputation

Metagenomic data sets of host-associated microbial communities often contain host DNA that is usually discarded because the amount of data is too low for accurate host genetic analyses. However, if a reference panel is available, genotype imputation can be employed to reconstruct host genotypes and maximise the use of such a priori useless data. We tested the performance of a two-step strategy to input genotypes from four types of reference panels, comprised of deeply sequenced chickens to low-depth host genome (~2x coverage) data recovered from metagenomic samples of chicken intestines. The target chicken population was formed by two broiler breeds and the four reference panels employed were (i) an internal panel formed by population-specific individuals, (ii) an external panel created from a public database, (iii) a combined panel of the previous two, and (iv) a diverse panel including more distant populations. Imputation accuracy was high for all tested panels (concordance >0.90), although samples with coverage under 0.28x consistently showed the lowest accuracies. The best imputation performance was achieved by the combined panel due to the high number of imputed variants, including low-frequency ones. However, common population genetics parameters measured to characterise the chicken populations, including observed heterozygosity, nucleotide diversity, pairwise distances and kinship, were only minimally affected by panel choice, with all four panels yielding suitable results for host population characterization and comparison. Likewise, genome scans between the two studied broiler breeds using imputed data with each panel consistently identified the same sweep regions. In conclusion, we show that the applied imputation strategy enables leveraging insofar discarded host DNA to get insights into the genetic structure of host populations, and in doing so, facilitate the implementation of hologenomic approaches that jointly analyse host genomic and microbial metagenomic data. Author summary We introduce and assess a methodological approach that enables recovering animal genomes from complex mixtures of metagenomic data, and thus expand the portfolio of analyses that can be conducted from samples such as faeces and gut contents. Metagenomic data sets of host-associated microbial communities often contain DNA of the host organism. The principal drawback to use this data for host genomic characterisation is the low percentage and quality of the host DNA. In order to leverage this data, we propose a two-step imputation method, to recover high-density of variants. We tested the pipeline in a chicken metagenomic dataset, validated imputation accuracy statistics, and studied common population genetics parameters to assess how these are affected by genotype imputation and choice of reference panel. Being able to analyse both domains from the same data set could considerably reduce sampling and laboratory efforts and resources, thereby yielding more sustainable practices for future studies that embrace a hologenomic approach that jointly analyses animal genomic and microbial metagenomic features.


49
The large molecular data sets generated through shotgun DNA sequencing usually contain useful 50 information to characterise taxa, functions and structures beyond the primary aim of the study. 51 This is especially true in metagenomic data sets that often present mixtures of DNA from 52 eukaryotic, prokaryotic and viral origin (1,2). While primarily used for characterising the genomic 53 architecture of microbial communities, metagenomic data generated from gut contents or faeces 54 can also be used for extracting useful genomic information of the animal host (3). In fact, 55 hologenomic approaches that entail joint analysis of animal genomes along with metagenomes of 56 associated microorganisms to study animal-microbiota interactions, can benefit from such 57 optimisation strategies (4,5). 58 However, mining host genomic data from metagenomic data sets entails a number of challenges. 59 The fraction of host sequences in the metagenomic mixture is often unpredictable, and can range 60 from a negligible proportion (<5%) to an almost complete representation (>95%) of the sample (6), 61 even within a single taxon and sample type (7). Hence, a given amount of metagenomic 62 sequencing effort does not ensure that the desired depth of host DNA sequencing will be reached. 63 When the host DNA fraction in the metagenomic mixture is low, achieving the desired sequencing 64 depth requires increasing sequencing effort, with its respective economic burden. In consequence, 65 the amount of host DNA sequences generated is often insufficient for accurate variant calling. 66 Reference samples 153 Internal and external high-quality genome sequence data was used to create the reference panels. 154 The internal reference data were generated from ileum content samples of 12 randomly selected 155 individuals included in the target population (7 Ross and 5 Cobb), following the same procedures 156 as explained above. In contrast to caecum samples, ileum contents contain a very large fraction 157 (>90-95%) of host DNA, and a small representation of microbial DNA. Hence, in order to generate 158 a comparable amount of microbial data to that of the caecum, ileum samples were sequenced 159 aiming 100 GB/sample. This sequencing effort yielded about 90 GB of host DNA (ca. 80-90x depth 160 of chicken genome), which enabled generating a high-quality internal reference panel from a 161 subset of the studied population. In addition, chicken DNA sequence data of 40 broilers, 20 layers 162 and 20 red junglefowls (RJF) generated by Qanbari et al. (2019) from blood samples were used 163 as external reference data ( Fig 1A). 164 Composition of reference panels 165 We used different combinations of the internal and external reference samples to create the four 166 reference panels used to evaluate imputation accuracy and impute the target population: (i) The 167 internal panel comprised 12 animals from our target population (7 Ross and 5 Cobb), (ii) the 168 external panel comprised 40 animals from two broiler breeds (different to our target population), 169 (iii) the combined panel combined the previous two panels, and (iv) the diverse panel contained 170 more distant populations (Fig 1B). The four panels varied in size and genomic diversity in order to 171 see whether the composition of the reference panels affected imputation accuracy. With the 172 internal panel, we tested if a small subset of the target population was enough for a proper 173 imputation in low-quality host sequence data derived from metagenomic samples. Assembly accession GCF_000002315.6) was conducted with BWA-MEM (v0.7.16a) (29). We 187 employed default parameters except for the minimum seed length (-k), which was increased to 25 188 in order to reduce the number of incorrectly aligned read pairs. We added the flag -M, which was 189 used to mark shorter split hits as secondary mappings. Aligned reads were sorted and converted 190 into sample-specific BAM files before filtering out the metagenomic fraction (unmapped) using 191 SAMtools view (v1.11) (30) with "-b" and "-F12" flags. Mapping statistics including depth and 192 breadth of coverage as well as percentage of mapped reads were calculated using SAMtools' 193 depth and flagstat functions. 194 Pure genomic data (with no microbial fraction) generated by others (31) was downloaded from the 195 EMBL-EBI ENA database, and mapped to the same chicken reference genome using BWA-MEM 196 with -k default value and -M flag. 197

Variant calling and genotyping
198 Variants in the target population were called by chromosome with the mpileup utility of SAMtools 199 using standard parameters (-C 50 -q 30 -Q 20). Variant calling was performed with "-m" and "-v" 200 flags to allow variants to be called on all samples simultaneously. Raw variants were filtered using 201 BCFtools (v. 1.11) (32) commands "-m2", "-M2" and "-v snps" to keep only bi-allelic SNPs. 202 Variants of the internal reference samples were called the same way, but additionally, low quality 203 variants with a lower base quality than 30 (QUAL<30) and variants with a base depth higher than 204 three times the average (DP<(AVG(DP)*3) were removed to ensure only highly reliable variants 205 were retained. 206 207 Since we were solely interested in imputing variants present in our target population, the external 208 reference samples were genotyped by defining variant sites detected in the internal reference 209 samples. Genotyping was performed for all autosomal chromosomes with GATK (v4.1.7.0) (33) 210 HaplotypeCaller using the "--min-base-quality-score 20", "--standard-min-confidence-threshold-211 for-calling 30", "--alleles" and "-L" parameters to obtain calls at all given positions, followed by 212 GATK SelectVariants "--select-type-to-include SNP" to only include SNPs. 213 In preliminary analyses, we also called variants in the external reference panel in order to examine 214 the overlap with the variants present in the internal reference samples. We used the same 215 procedures explained above for GGA1. Genotyping based on the positions of the internal panel 216 and variant calling from scratch were compared by using the 40 broilers from the external reference 217 panel for GGA1 ( Fig 1B). A similar number of variants had been obtained for the genotyped (2.5 218 M) and the variant called VCF files (2.7 M). Moreover, 28% of the variants from the 40 broilers 219 were not present in the internal reference samples ( Fig S2). Thus, we decided to genotype the rest 220 of the samples to reduce possible bias through the high number of variants specific to the external 221 reference for the imputation of our target population. 222 Two-step imputation via genotype likelihood updates 223 We imputed genotypes from the four aforementioned reference panels to the target population 224 using a two-step strategy. Prior to imputation, the reference panels were filtered by excluding 225 variants with missing genotypes to remove any potential noise caused by inference errors, and 226 subsequently phased using SHAPEIT (v4) (34). 227 Imputation was performed in two steps following Homberg et al. (2019) and Hui et al. (2020). First, 228 genotype likelihoods were updated based on one of the reference panels using Beagle 4.1 (35). 229 Beagle 4.1 accepts a probabilistic genotype input with "-gl" mode, and it only updates sites that 230 are present in the input file. Second, missing genotypes in the input file were imputed using Beagle 231 5.1 with "-gt" mode using the same reference panel. Beagle 5.1 only accepts files with a genotype 232 format field, like later versions than Beagle 4.1. Therefore, the latest version cannot be used for 233 both steps. Format field genotype probabilities (GP) were generated in both steps in order to enrich 234 confident genotypes. We required the highest GP to exceed a threshold of 0.99 after both steps 235 using BCFtools +setGT plugin. The rest of the parameters were set to default. Both steps' input 236 and output files were in VCF format. The schematic steps detailed in methods can be found in Fig  237   1C and the scripts in the following link (https://github.com/SofiMarcos/Host-genome-recovery.git). 238 Imputation accuracy using 12 validation samples 239 The accuracy of the imputation using the four reference panels was tested using the 12 individuals 240 for which we generated both low-depth (target population) and high-depth (internal reference 241 samples) sequence data from caecum and ileum contents, respectively, hereafter referred to as 242 validation samples. The low-depth samples of the 12 individuals had a depth of coverage spanning 243 0.05x to 3.73x. For an unbiased evaluation, we employed a leave-one-out cross-validation 244 (LOOCV) approach by excluding each of the 12 validation samples once from the reference panel 245 in each of the different imputation scenarios. Considering the large size-variation of avian 246 chromosomes, a macrochromosome (GGA1, 197.6 MB), a mid-size chromosome (GGA7, 36.7 247 MB) and a microchromosome (GGA20, 13.9 MB) were selected for the test to optimise runtime 248 and computational resources. Concordance between the internal reference samples and imputed 249 genotypes was calculated for each individual chicken using VCFtools, with the"--diff-discordance-250 matrix" option. Precision of heterozygous sites was also calculated, since these alleles are the 251 most difficult to impute correctly. Kruskal-Wallis test was performed to test for differences across 252 chromosomes. A paired sample T-test and F-test were performed for both parameters to verify if 253 the difference in means and variances were significant between reference panels. T-test p-values 254 were adjusted using Bonferrini's correction method. Moreover, imputation accuracy was estimated 255 for variants in different minor allele frequency (MAF) bins to evaluate whether rare and common 256 variants are equally correctly imputed. We thus extracted variant frequencies from the internal 257 panel by analysing precision of heterozygous (het.) sites for the GGA1 in bins of 0-0.05, 0.05-0.1, 258 0.1-0.3 and >0.3. 259 260 We explored the implications of using different reference panels in downstream analyses of 261 population genetic inferences, including population structure, genetic diversity, and genome scans 262 for signatures of selection. 263

Impact of reference panel on population genetics inference
These analyses were run in all but two outlier samples with depths of coverage of 0.07x and 0.05x, 264 which were below the threshold of 0.28x corresponding to the lowest successfully imputed sample 265 in the validation set (genotype concordance of >0.90 and het. sites precision of >0.75, see results 266 below). We thus used 100 samples (53 Ross and 47 Cobb) for which we ran the host DNA recovery 267 pipeline for all the autosomal chromosomes and analysed common population genetics 268 parameters including observed heterozygosity (O.Het), nucleotide diversity (π), pairwise distance 269 as estimated through identity-by-state (1-IBS) and kinship. The same analyses were also 270 conducted for 10 validation samples (for the low-depth and high-depth samples) after excluding 271 two of them, whose respective counterparts in the target populations (with 0.05 and 0.07x depth) 272 were filtered out. The imputed data sets with each of the panels were filtered for missingness 0 273 with PLINK (v1.9) (36). 274 For measuring population genetics parameters, the VCF files were filtered for MAF >0.05. O.Het 275 was calculated for each individual using the command "--het" in PLINK (v1.9). π was calculated in 276 40 kb windows with 20 kb step size across autosomal chromosomes using VCFtools. For the 277 validation samples whole-genome windowed values were averaged to generate a genome-wide π 278 for each individual. For the target population, π was calculated for each breed population. Paired 279 sample T-tests were performed for O.Het and π parameters. Pairwise distance was calculated 280 using "--distance square 1-ibs" in PLINK (v1.9). Kinship was calculated with the command "--make-281 king square" using PLINK (v2). To test the correlation between the resulting matrices from the 282 pairwise distance and kinship analyses using different panels, a Mantel test was performed with 283 the R package ade4 (37). 284 We further tested whether genome scans for selection between the Cobb and Ross population 285 with each of the imputed datasets yielded consistent results. To this end, we calculated population 286 differentiation along the genome using fixation index (FST) between both breeds using each panel. The pipeline required some tests and adjustments to optimise it to our system. The standard 301 alignment (seed length 19) presented an unconventional distribution of reads across the genome, 302 i.e. unspecified read mapping leading to regions being stacked with 80+ reads (Table S1). In order 303 to remove as many remaining microbial reads as possible, we increased the seed length to 25. 304 Standard deviation of the depth of coverage decreased considerably (from 202.79 to 3.66), while 305 the mean depth decreased from 2.78x to 1.73x. The breadth of coverage decreased by 9% (Fig  306   S1). 307 Imputation accuracy of 12 validation samples 308 The internal (n=12), external (n=40), combined (n=52) and diverse (n=92) reference panels were 309 used to study (i) the effect of panel size and diversity and (ii) sample depth of coverage threshold 310 on imputation accuracy in three chromosomes with contrasting dimensions. Variant calling in the 311 internal reference samples detected 2.4 M, 470 K and 182 K putative SNPs in chromosomes 312 GGA1, GGA7 and GGA20, respectively. After genotyping the external reference samples and 313 combining them to create the external, combined and diverse panels, each panel was filtered 314 before being phased. As a consequence, the filtering step decreased the number of SNPs by 315 13.83±1.36% for the external and combined, and by 23.80±0.99% for the diverse panel, which 316 yielded panels with different numbers of SNPs (Fig 2A). More than 96% of the total SNPs in each 317 panel successfully passed the multiple filters of the pipeline, even for samples with less than 1x 318 coverage ( Fig 2B). Furthermore, the proportion of imputed SNPs increased and gained uniformity 319 across samples when the panel was larger but had fewer SNPs. The mean number of imputed 320 SNPs across samples differed between all the panels: internal vs external (t=14.58, p-value < 321 0.001), external vs combined (t=13.56, p-value < 0.001) and combined vs diverse (t=11.63, p-322 value < 0.001). The F-test to compare variances was significant only between the diverse and the 323 rest of the panels: internal vs diverse (F= 30.54, p-value<0.001), external vs diverse (F=24.24, p-324 value< 0.001) and combined vs diverse (F= 11.31, p-value<0.001). Results indicate that the 325 variance across samples for the diverse panel greatly decreased compared to the rest of the 326 panels ( Fig 2B). 327 between panels for each chromosome (Fig 3). Concordance values significantly differed when 344 comparing the internal, external and combined panels (Fig 3A). But no differences were detected 345 between the combined and the diverse panels, indicating that no significant increase in imputation 346 accuracy can be achieved in terms of overall concordance by adding more distant individuals. For 347 precision of het. sites, differences were detected for all panels (Fig 3B), including for the combined 348 and the diverse except for GGA20. This suggests that the heterozygous positions are the most 349 sensitive to the imputation process. In an attempt to further assess imputation accuracy, we classified variants according to their MAF 359 in four bins (0-0.05, 0.05-0.1, 0.1-0.3 and >0.3) and calculated precision of het. sites, and the 360 number of correctly imputed variants for the 12 validation samples for GGA1 (Fig 4). The internal 361 panel, while recovering the largest number of variants, was also the panel with the lowest 362 performance in adequately inferring low-frequency variants, especially for the variants with MAF 363 <0.1 (Fig 4A). Although there was no improvement from the external to the combined panel for the 364 smallest MAF bin, a substantial improvement was seen for the rest of the bins. Some significant 365 differences but not as pronounced were also observed from the combined to the diverse. 366 Therefore, the combined panel showed overall the best results with the highest number of correctly 367 imputed variants in all MAF bins (Fig 4B), while maintaining a very similar number of imputed SNPs 368 as the external panel. The diverse panel inferred fewer low-frequency variants, but did so more 369 effectively (Fig 4). 370 Variants that coincided between imputed (low-depth 12 validation samples) and true (internal 377 reference samples) genotypes were considered correctly imputed variants. Paired T-tests were 378 performed to identify significant differences in means across panels: the following symbols ( "**", 379 "*", "ns") indicate different p-value cut-points (>0.001, 0.001, 0.05). 380 381 Despite the high overall imputation accuracy, the two samples with depths of 0.05x and 0.07x were 382 outliers that did not achieve a sufficiently high concordance (>0.90) and precision (>0.75) with any 383 of the panels and chromosomes (Fig 3). They were thus excluded from the target population, and 384 we refer from now on to 10 validation samples instead of 12. 385 with different panels decreased as more distant individuals were included (Fig 5A). This was due 390 to the missing call rate (MCR) filter during the two-step imputation. Using the internal panel, we 391 recovered 11.7 M filtered SNPs in the target population. These were 30% more recovered variants 392 than when using the diverse panel (8.9 M). Most of the excess variants from the internal panel are 393 low-frequency variants that cannot be confidently recovered (Fig 5B), as seen in the less effective 394 imputation of low-frequency variants with the internal panel ( Fig 4A). Both Ross and Cobb 395 populations showed extreme allele frequencies (peaks at both ends of the distribution, Fig 5B)  imputed and high-depth validation samples (Fig S3). There were some outlier samples (two from 426

Panel choice impact on population genetic inference
Cobb and one from Ross) that presented lower O.Het than the high-depth validation samples (Fig  427   6A). These samples apparently underwent an incorrect imputation process, but it was not 428 necessarily related to a low mapping depth. 429 Nucleotide diversity, on the other side, decreased with increasing panel size and diversity (Fig 6B), 446 which was directly related to the lower number of variants retained in the external, combined and 447 diverse panels compared to the internal. There were significant differences in means except 448 between the external and combined panels for both breeds, most likely because of the similar 449 number of variants both panels share (Fig 5A). When comparing the imputed population with the 450 validation samples, π of imputed samples and of the target population were underestimated for all 451 panels (Figs 6B and S3). 452 Regarding the population genetic interpretation, both populations were very similar, but the 453 imputation tended to accentuate differences between the two populations (Figs 6 and S4). Within 454 population pairwise distance and kinship values did not vary much according to the panel. For 455 pairwise distance, the diverse panel resulted in larger interindividual distances within breeds (Fig  456   S4). Kinship values were lower when computed with the internal panel, since a larger number of 457 SNPs were retained, in particular, low-frequency variants which are typically unique to one or few 458 individuals thus decreasing kinship (Fig S4). Mantel R tests did not show any significant differences 459 for pairwise distance and kinship matrices, giving the same result for all panel comparisons (Mantel 460 statistic, p-value < 0.001). Correlation values for pairwise distance were very similar and close to 461 1 (Fig 6D), even for the validation samples (true genotypes) when compared with any panel (Fig  462   S3). For kinship instead, it seemed that the internal panel differed more from the rest (Figs S3 and  463   6C). In both cases, the 10 validation samples were most correlated with samples imputed with the 464 combined and diverse panels (Fig S3). 465 Whole-genome mean FST values between Ross and Cobb populations were very similar (internal 466 0.071, external 0.071, combined 0.072 and diverse 0.072) indicating overall low differentiation 467 between the breeds. When analysing the putative selective sweep regions using as threshold the 468 99th percentile, 68.2% of the windows coincided across the four panels, but more interestingly, 469 75.9% of the windows were shared across the external, combined and diverse panels. When we 470 raised the threshold to the 99.9th percentile, 77% of windows were identified by the genome-scans 471 regardless of the choice of panel, indicating that the strongest signals are detected with any panel. 472 Yet, there were some regions that only passed the threshold when imputation was performed with 473 a particular panel (Fig 6F). The combined panel did not show specific sweeps when the percentile 474 was set at 99.9, and it was the panel with the lowest panel-specific regions with the 99th percentile 475 as well, potentially indicating the most robust results, i.e. without panel-specific biases. 476 Surprisingly, the diverse panel detected the most panel-specific sweeps after the internal panel 477 (Fig 6F). On the other side, in terms of density of variants in the common windows, the mean 478 number of variants reduced significantly from the internal to the diverse panel ( Fig 6E). This 479 suggests that although in a broad sense the same sweep signals can be detected by all panels, a 480 reduced number of imputed variants might give a smaller chance of detecting causative variants. 481

482
Shotgun metagenomic datasets of host-associated microbial communities often contain host DNA 483 that is usually discarded because the amount of data is too low for accurate host genetic analyses. 484 Here, we introduced an effective and accurate approach to recover high-quality host genomes 485 from gut metagenomic data, which can be used to study host population genetic analyses and 486 ultimately contribute to a better understanding of host-microbiota interactions. 487 Our analyses yielded drastic differences in mapping statistics between caecum samples used to 488 characterise the target population and ileum samples employed to generate the internal reference 489 panel. Although both sample types derived from gut contents, the caecum harbours a very small 490 amount of the host DNA compared to the ileum, because the latter is known to contain fewer 491 bacteria (39), and the higher permeability and a thinner mucus layer of the ileum probably entails 492 higher release of epithelial cells to the lumen (40). Moreover, the low, yet variable, proportion of 493 host DNA retrieved from caecum samples renders sequencing depth adjustment highly 494 unpredictable, as previously reported (7). Notwithstanding, we showed that if a proper reference 495 panel is designed, the low and variable fractions of host DNA recovered from such suboptimal 496 samples, can be used for accurately inferring host genetic features. It must be noted though, that 497 the ratio of host and microbial DNA recovered from chicken caecum and ileum samples can not 498 directly be extrapolated to other host taxa and sample types. 499 The two-step imputation strategy performed efficiently despite the structural (e.g., study design, 500 animal taxa, reference panel size) differences between our study system and the ones the strategy 501 was originally designed for (23,24). First, we used custom reference panels with less than one 502 hundred individuals, while the two-step strategy was originally tested with large reference panels 503 such as the Human 1000 Genomes (12). Nevertheless, our accuracy values were comparable to 504 the previous results, most likely because the individuals in our target population were closely 505 related, as evidenced by the high kinship values. Second, although we had a similar range of target 506 population sample depths (Hui: from 0.05x to 2x and Homburger: 0.54x to 1.76x), our samples 507 consisted of real low-depth sequence data, instead of downsampled sequencing reads from high-508 depth samples. Thus, mapping gaps across the reference genome were unevenly distributed. This 509 is evidenced by the large difference between depth (1.8x) and breadth (50%) of coverage (S1 510 Table), likely hampering accurate computation across the genome. Besides, Hui et al. (2020) 511 documented that the proportion of correctly imputed heterozygous sites started decreasing at 0.5x 512 of depth of coverage, reaching 50% of correctly imputed sites at 0.1x. In our system, >90% of the 513 variants in samples with 0.28-0.5x could be recovered, and accuracy only dropped significantly in 514 samples below 0.1x. Accordingly, we decided to set a mapping depth threshold at 0.28x, but we 515 recommend adjusting it depending on the sample size and quality of the data set, as well as the 516 accuracy needs of each study. 517 The accuracy of low-frequency variants for all panels except for the internal, which showed much 518 lower values, were comparable to previous works (24), most likely owing to the stringent filtering 519 criteria applied in our study (MCR = 0). But the overall accuracy and the accuracy of heterozygous 520 sites depends heavily on variant frequencies, therefore these comparisons should not be decisive. 521 Finally, unlike humans, avian genomes present macro-and micro-chromosomes and the latter 522 frequently undergo interchromosomal translocations (41). However, it seems that the possible 523 interchromosomal translocations of the target population did not affect imputation, since we did not 524 find any significant differences in accuracy between chromosomes, revealing that the strategy 525 worked equally well for large, mid-sized and small chromosomes with potentially different linkage 526 patterns. 527

Effect of reference panel on accuracy statistics 528
Reference panel design depends on data availability as well as computational capacity. It is a 529 common strategy for imputation of inbred populations to resequence a subset of samples with 530 higher resolution in order to optimise imputation performance (35). Based on previous works, we 531 estimated that 12 individuals out of 100 would be sufficient to represent the genetic diversity of the 532 population. For instance, previous chicken studies deep-sequenced 25 individuals to impute 533 approximately 450 chickens genotyped with 600-K SNP arrays (~5% of sample size) (20,42). 534 In terms of panels SNP density, we decided to genotype variants that did appear in our target 535 population rather than calling for specific variants in the rest of the breeds that composed the 536 reference panels. Thereby, we aimed at reducing the noise that the excess of variant density could 537 cause in the imputation process. Nevertheless, as the genetic distance between the selected and 538 our breeds is very small (43), we expected them to share many variants, as we evidenced with 539 preliminary analyses using GGA1 where 72% of variants identified by genotyping or by calling 540 overlapped in the external panel (Fig S2). 541 The internal panel resulted in a larger variance across samples. SNPs with low MAF had the lowest 542 accuracy when imputed with the internal panel. Moreover, incorrectly imputed low-frequency 543 variants can be easily overcome if a strict MAF filter is applied for downstream analysis. Another 544 possible option is to sequence more individuals of the target population to increase the reference 545 panel size. Hence, despite the internal panel only representing a small subset of the target 546 population, and showing lower imputation values than in the external, combined, and diverse 547 panels, for scientists without access to external reference samples, this approach is equally useful 548 as overall imputation accuracy was higher than 90% and biological differences were still visible. In 549 this sense, host resequencing of a small subset of the target population might represent a cost-550 efficient option, especially for researchers working with non-model organisms and inbred 551 populations. Thus, our approach could be useful, for example, to study genome features of 552 endangered populations relying on faecal samples recovered from the environment. 553 Our results showed that the combined panel performed better in terms of overall accuracy, and 554 specifically of minor allele frequency variants, than the internal and the external panels alone. 555 Despite the fact that the external and combined panels had the same number of SNPs, including 556 a subset of individuals from the target population was beneficial. Many studies already mentioned 557 an improvement for the combined option (44,45). Lastly, the diverse panel showed the highest 558 values of concordance and het. sites precision, most probably because of the lower number of 559 SNPs recovered, especially low-frequency variants, which generally yielded lower imputation 560 accuracies. In terms of imputation of low-frequency variants, the combined panel outperformed the 561 diverse one, i.e. it correctly imputed a larger number of variants and tended to improve the 562 precision of het. sites in some MAF bins. A recent large-scale study performed in a Chinese 563 population showed that a population-specific reference panel worked the best compared to 564 European reference panels such as 1000G (21). Imputation was greatly improved when the 565 reference panel contained a fraction of an extra diverse sample, but they obtained a different 566 pattern when the panel size was fixed (21). Thus, taking into consideration our and previous results 567 on selection of imputation panels, it can be concluded that increasing panel size and diversity 568 improves imputation, but a balance has to be found in the composition of the panel. The distance 569 between the panel and the target population has to be taken into account. 570

Effect of reference panel on biological inference 571
Besides crude imputation accuracy statistics, we evaluated the impact of the panels on 572 downstream population genetic parameters and their biological interpretation. As imputation 573 accuracies were generally high with our applied pipeline and the stringent filtering approach, we 574 expected population genetic inferences to follow similarly. 575 Although overall results were in agreement, all the tested parameters showed slight trends 576 according to the used reference panel. O.Het, pairwise distance and kinship values increased 577 while mean FST and π values decreased with panel size and diversity (Figs 6 and S4). Such 578 biases were related to the composition of the panels and the associated number and distribution 579 of recovered SNPs. 580 Imputation worked slightly differently for the two breeds, as Ross population estimations were 581 closer to the true values than for the Cobb population. Thus, accentuating the distance between 582 both breeds. This is most likely due to a smaller representation of Cobb individuals in the reference 583 panels, i.e. 5 Cobb and 7 Ross samples constituted the internal reference panel. Secondly, there 584 were some samples that were incorrectly imputed because of their low O.Het values (Fig 6A). We 585 do not know if there are individuals with lower heterozygosity in our Cobb and Ross populations. 586 For instance, there was a Ross individual from the high-depth validation samples with lower O.Het. 587 Chickens came from two different hatcheries, which might be the reason why some individuals 588 might have slightly different genetic features. We may have under-represented one of the origins 589 in the internal reference samples. Thus, it is necessary to be more cautious for the interpretation 590 of individual genomes. Nevertheless, results appeared to be robust and similar across panels at 591 the population level. The genome scans yielded overall very consistent results with major 592 differentiation signals identified by any of the imputed datasets, likely indicative of a true selection 593 signature between both breeds. However, downstream analyses such genome scans and GWAS 594 must be performed with caution since this method is sensitive to low-frequency variants quality. 595 Both breeds exhibited extreme minor allele frequencies, indicating that the genetic drift due to 596 selection in a closed breeding population has a notable effect. Domestication and breeding history 597 are the two major processes that shape haplotype structure (31,46). Cobb and Ross, together with 598 other commercial breeds, have much smaller effective population size than other chickens (47). 599 Broiler breeding methods are described as a pyramid strategy, in which pure, inbred lines are 600 crossed, then F1 individuals are crossed between each other. In some cases, even a second or a 601 third cross is performed in F2 and F3 generations before raising them for meat (48). Therefore, 602 broilers are highly related populations, but at the same time present high heterozygosity values. 603 Heterozygosity of our studied breeds were much higher O.Het than of local populations (49), but 604 similar to other broiler breeds (46). Similarly, nucleotide diversity and mean fixation index values 605 were comparable to those previously reported (31). 606

607
Our results show that the two-step imputation implemented in this study can be used to 608 successfully reconstruct genotypes and study population genetic properties of hosts from 609 suboptimal metagenomic samples. The comparison among reference panels also demonstrated 610 that this method is versatile and flexible. This approach could be used in many contexts and exploit 611 different data sources to address a variety of research questions. This includes the possibility of 612 mining published metagenomic data sets to recover discarded host DNA sequences. In our 613 particular case, the reconstructed genotypes will be employed in the H2020 project HoloFood to 614 detect interactions with microbial metagenomic features, and thus implement a hologenomic 615 approach to improve animal production (50). Because 'host-contamination' should no longer be 616 considered a problem, but an opportunity. 617