{'fraction_non_alphanumeric': 0.02946810078090467, 'fraction_numerical': 0.02264132409999509, 'mean_word_length': 5.001178897730622, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Targeted genotyping by sequencing: a new way to genome profile the cat

Summary Targeted GBS is a recent approach for obtaining an effective characterization for hundreds to thousands of markers. The high throughput of next‐generation sequencing technologies, moreover, allows sample multiplexing. The aims of this study were to (i) define a panel of single nucleotide polymorphisms (SNPs) in the cat, (ii) use GBS for profiling 16 cats, and (iii) evaluate the performance with respect to the inference using standard approaches at different coverage thresholds, thereby providing useful information for designing similar experiments. Probes for sequencing 230 variants were designed based on the Felis_catus_8.0. 8.0 genome. The regions comprised anonymous and non‐anonymous SNPs. Sixteen cat samples were analysed, some of which had already been genotyped in a large group of loci and one having been whole‐genome sequenced in the 99_Lives Cat Genome Sequencing Project. The accuracy of the method was assessed by comparing the GBS results with the genotypes already available. Overall, GBS achieved good performance, with 92–96% correct assignments, depending on the coverage threshold used to define the set of trustable genotypes. Analyses confirmed that (i) the reliability of the inference of each genotype depends on the coverage at that locus and (ii) the fraction of target loci whose genotype can be inferred correctly is a function of the total coverage. GBS proves to be a valid alternative to other methods. Data suggested a depth of less than 11× is required for greater than 95% accuracy. However, sequencing depth must be adapted to the total size of the targets to ensure proper genotype inference.


Introduction
The global pet care market size (major segments including food, veterinary care and over-the-counter products) was estimated at USD 131.7 billion in 2016 and is expected to reach USD 202.6 billion by 2025, an estimated growth of 4.9% calculated with the Compound Annual Growth Rate (https://www.grandviewresearch.com/press-release/globalpet-care-market, March 2018). Cats are increasingly appreciated as pets because they are known for helping reduce stress and anxiety and for having strong interactions with humans (Hart et al. 2018). In this context, fancy breeds are becoming more and more popular worldwide. In pedigreed cats, studbooks can recommend a DNA-based control of both animal identity and traits of interest for enrolment and selection, together with a permanent electronic identification. The International Society of Animal Genetics (ISAG) fosters the definition and nomenclature standardization of panels of genetic markers for the identification and parentage control of domestic animals, including cats (Lipinski et al. 2007; https://www.isag.us/committees.asp). These panels are used by service laboratories for owners and breeder association requirements. Recently, ISAG introduced single nucleotide polymorphisms (SNPs) as additional markers for the cat genetic profiling and parentage test after more than a decade of exclusively using microsatellites (short tandem repeats (STRs); Lipinski et al. 2007). As novel polymorphisms affecting health and morphology are discovered, the identification of carriers is a growing problem, which requires versatile approaches when regions of interests must be added or changed rapidly. The simultaneous analysis of mutations and neutral markers that define cat identity allows a univocal match between the animal and the biological profile, helping to avoid mistakes and fraud. Following this line of reasoning, and under the expectation of a huge increase in the number of documented trait-associated polymorphisms, a targeted GBS approach was implemented. Such a method can probably be used effectively at reasonable costs.
Genotyping by sequencing has been widely used in plants as a valid alternative to the currently available genotyping technologies for molecular marker discovery, genome-wide associations (Zhang et al. 2015), genomic diversity analyses (Pembleton et al. 2016) and breeding programmes (He et al. 2014;Kim et al. 2016). GBS can also be applied successfully to animal genomes (De Donato et al. 2013;Gurgul et al. 2018). GBS has been used recently in animals, even when no SNP chips were available. Moreover, GBS has been demonstrated to be an effective method for genome-wide SNP discovery and genotyping, opening the door to inbreeding control and genomic selection, for example in Peking duck and Indian cattle (Zhu et al. 2016;Malik et al. 2018).
Genotyping by sequencing is a targeted multiplexed approach exploiting the extremely high throughput of next-generation sequencing (NGS), whereby a high number of regions of interest (for instance containing variants) are simultaneously enriched using specifically designed probes, amplified and then sequenced. After demultiplexing and quality control, bioinformatic analyses are performed to infer the genotypes at the variant loci. In humans, GBS showed high accuracy in the detection of mutations in patients with suspected genetic disorders and has been proposed for the diagnosis of more than 500 Mendelian diseases (Liu et al. 2015). The GBS approach can also be expanded to the monitoring of predisposition to polygenic disease or complex disease (Guan et al. 2015).
In the present work, the design and the data analysis of a customized GBS panel for the genomic profiling of the domestic cat is reported.

Samples
Sixteen genomic DNA Felis catus samples were analysed. Thirteen samples (Cat_1-13) from British Shorthair cats were provided by an Australian service laboratory (Orivet). Two came from the Animal Tissue and DNA Repository at the University of Milan, Italy (Cat_14-15), and one, from an Abyssinian cat (Cat_16), was part of the cat whole genome sequencing (WGS) 99 Lives Cat Genome Sequencing Initiative (NCBIe accession no. SRR6997541).
Genomic DNA had been extracted from tissue samples using a DNeasy Blood and Tissue Kit (Qiagen), according to the company instructions. The DNA was quantified with Qubit (Thermo Fisher Scientific) and adjusted to 25 ng/ll. Approximately 50 ng was used for genotyping each sample.

Genomic target
Cat-specific DNA variants (n = 230) were targeted, partially obtained by screening the references in OMIA (https://omia. org/home/) and relevant reviews (Lyons 2010(Lyons , 2012(Lyons , 2015. SNP1-SNP120 had already been genotyped by Orivet using probes designed based on Felis_catus_6.2 (accession no. GCA_000181335.2; Montague et al. 2014) and the Agena Biosciences iPLEX Gold Genotyping kit, as previously described (Oh et al. 2017) and shown in Table S1. STR loci (no. 19) were included in the 230 variants (Lipinski et al. 2007). The genotypes at all of the target loci were also already available for Cat_16 from its WGS.
DNA sequences of approximately 500 bp and including the target loci were manually downloaded from the genome cat assembly ICGSC Felis_catus_8.0/Felcat8 (https://ge nome.ucsc.edu/cgi-bin/hgGateway; Kent et al. 2002). The genome positions of the target variants were submitted to Agilent Technology for the SureSelect QXT Target Enrichment custom panel design, with a tier of 1-449 kb and a maximum of 57 680 probes (hereinafter the target probes) being requested. The Agilent custom kit included the design and production of five target probes of 120 bp for each target locus. Target probes were then manually controlled for correct design using INTEGRATED GENOME BROWSER (Affymetrix Inc. Research).

Sequencing data output and analysis
Paired reads obtained from the sequencing were combined using FLASH whenever they overlapped (Mago c & Salzberg 2011), setting minimum and maximum overlap to 10 and 150 nucleotides respectively. The resulting reads were aligned to the genome cat assembly ICGSC Felis_catus_8.0/ Felcat8 using BOWTIE2 (Langmead & Salzberg 2012), and the BCFTOOLS suite (Li 2011) was used to obtain the consensus base calling for the covered positions. In detail, BCFTOOLS MPILEUP was used to calculate the likelihood of the alignments present in each position. With BCFTOOLS CALL, using multiallelic caller option, the SNPs were called and the genotypes of the target genome positions were obtained. Informative values associated with the genotypes, such as the likelihood of the second most likely genotype (an indication of the robustness of the reported genotype inference) and the read coverage, were extracted and analysed in R (R Core Team 2013).
Comparisons of the GBS results with those from the Agena Mass Array (hereinafter named Reference) were based on 120 markers and involved 13 common samples (Cat_1-13), whereas the comparison with the WGS involved only the corresponding sample (Cat_16) in the GBS panel and comprised all of the markers. All of the comparisons were performed using custom-made R scripts.

Results
All of the target probes were successfully designed for all 230 regions. The libraries passed the quality controls, and the reads were quality filtered to keep only high-confidence sequences. GBS produced 1 402 870 reads with an average of approximately 87 679 reads per sample (Fig. 1). The length of the Agilent SureSelect amplicons (~120 bp) was not sufficient to generate reads spanning the full repetitive motif of the STRs, and genotypes could not be determined for the majority of the 19 STR markers (16 of 19; data not shown). Therefore, the STRs did not receive further consideration in the analyses. The target loci for the 211 SNPs were identified in each cat.
For the GBS results from the 211 SNPs, the sequencing depth varied among samples, ranging from 56 332 to 172 812 fragments per sample. Additionally, samples differed in the number of read mates combined (from 33 173 to 101 308) or not combined (from 23 159 to 71 504) by FLASH and in the number of fragments mapping to the cat genome (from 39 994 to 129 779) (Fig. 1).
A preliminary analysis with BCFTOOLS of the 211 SNPs allowed inference of the genotypes in the 16 cats. As a measure of significance of the inferred genotype (Table S1), the likelihood of the second most likely genotype was used, as reported by the program (Table S2). This value was chosen because it was informative about the existence of alternative SNPs at the same locus, providing a measure of the reliability of the most likely genotype. BCFTOOLS assigns the significance of the reported genotype a likelihood of 1, without respect to coverage or additional information. Therefore, the negative log likelihood of the second most likely genotype allowed an indirect assessment of the reliability of the reported one: when the value was high, the likelihood was low, which means that a valuable alternative to the called genotype was not detected. The reliability of the inferred genotypes was dependent on the coverage, as a small number of reads was insufficient to identify the correct genotype, especially in heterozygous conditions. Coverage and negative log likelihood of the genotype were positively related, meaning that the genotypes with higher coverage were identified with more confidence. In Fig. 2, sequencing coverage was plotted against the likelihood of the second genotype for all of the SNPs under analysis (data in Table S2). Despite the points being widely dispersed, the trend was clearly positive.
Similarly, by setting a threshold on the negative log likelihood ( Fig. 3; data reported in Table S3) at 15 (likelihood of the second genotype of 10 À15 ), the percentage of genomic positions that pass the threshold increased, making the coverage deeper but reducing the number of loci for which the genotype could be inferred. Briefly, the threshold at 15 was empirically selected based on the number of the genomic positions considered and on the proportion of errors admitted, which minimized the number of incorrectly assigned genotypes and allowed an adequate number of loci to be retained. In Fig. 3, the histogram refers to the left yaxis and shows that the number of SNPs for each interval of coverage decreased at higher coverage for the given threshold. On the right y-axis, the fraction of target loci overcoming the likelihood threshold per coverage bin was also calculated. The chart clearly shows that the fraction of target loci above the coverage threshold increased when the coverage was larger.
The average number of target loci per sample (indicated on the y-axis) with a coverage exceeding the threshold (indicated on the x-axis) is shown in Fig. 4. Combined with the information in Fig. 3, this number highlights the compromise among the fraction of recovered target loci (which was maximal for low coverage thresholds) and the quality of the resulting genotypes (maximal for high coverage).

Genotyping by sequencing
Genotyping results were compared with the available genotype data (Reference), and the accuracy was analysed.
As a first assessment, the work focused on the 120 SNP markers genotyped in 13 samples (1560 genotypes) and included both in the GBS panel and Reference data, with the latter assumed to have 100% accuracy. Of the 120 SNPs, three markers recorded low performances: SNP33, SNP60 and SNP101. These markers showed large failure (missing genotypes) and low concordance in genotyping with the Reference, suggesting a design failure. SNP45, SNP86 and SNP90 also showed low concordance. Furthermore, SNP20, SNP71, SNP75 and SNP115 did not have MassArray data for comparison for most of the samples. All discordant genotypes concerned homozygous loci. This finding was expected, as the probability of observing both bases at heterozygous loci decreases at smaller coverage (Table S4).
The ratio of the number of GBS matching Reference genotypes to the number of genomic positions shared was calculated. For the GBS results, the performances for several coverage thresholds were reported, as the coverage had previously been shown to have a strong impact on the accuracy and completeness of the results (Figs 3 & 4). These thresholds were chosen as representative of shallow to strict filtering, but at least half of the total genomic positions were always included so as not to lose too much information. For each coverage threshold, the genomic positions typed with both methods were considered for comparison. The resulting assessment (Table 1A), indicated good performances for To check if any of the samples had a strong impact on the global accuracies, the accuracy on a per sample basis was calculated, with the results showing that the global performances were reduced by a few samples (Table 1B, &  Table S4 limited to the heterozygous performances). Here, the same approach mentioned above that was used to calculate performances was applied once more. Notably, different pools of genomic positions might be compared depending on the samples being contrasted, which is different from the previous comparison analysis.
The comparison of GBS with WGS (309 average depth, one single cat, n = 16) further indicated the importance of coverage: with a lower threshold, the performance was low compared with the results discussed above, but with a stringent cut-off on coverage, the accuracy increased, even though the genomic position overcoming the coverage threshold was only 25-50% of the original dataset (Fig. 5). This confirmed that the accuracy can be easily increased to the highest values by increasing the depth just above 119.

Discussion
Since 2004, ISAG has organized cat comparisons using anonymous DNA markers and microsatellites and, since 2017, also SNPs to standardize genetic profiles for subject identification and parentage control, so these profiles could be recognized by genetic testing service laboratories worldwide. However, the availability of a unique assay to genotype neutral loci for identity profiling, together with coding loci for traits of interest, could provide all of the possible requested information and protect against mistakes and errors.
DNA chips can be used to type large numbers of SNPs (Alhaddad et al. 2013;Gandolfi et al. 2018), combining identity profiling and coding trait loci, but their customization is expensive in both the production and application phases because they are static tools that must be entirely used or wasted before adding content. The possibility to calibrate the average sequencing coverage, the high barcode sample multiplexing/pooling and the easy addition of new markers all make GBS a flexible, customizable, upgradable and effective alternative to SNP chip or other methods and, overall, make it economically competitive in a short time prospect (Gorjanc et al. 2015). In fact, the use of NGS through intensive multiplexing dramatically increases the captured fraction of the genome. This fraction can be readily adjusted according to the needs of the project, and samples obtained from different species can potentially be pooled for sequencing.
When the present GBS experiment was planned, the cost of the design was reduced as much as possible, with a compromise among the depth of the sequencing and completeness of the results being sought. In the design, in addition to SNP variants, some indels (STRs) were also included for testing. As expected, the short length of the reads and the customized, but standard, design were insufficient to resolve the full repetitive motif of most STR markers, causing the overall genotyping to fail at these loci.
All of the remaining loci (n = 211 SNPs) were successfully genotyped. An assessment of the accuracy of the GBS results was conducted on 13 samples (Cat_1-13) that were previously genotyped with Agena at 120 loci. One additional sample (Cat_16) was compared with the information from its WGS (309 average depth), in this case considering the full set of 211 SNPs. Analyses showed that, by using a relatively shallow sequencing (~119), achieving good performances in genotyping is possible. The range of accuracy was always affected by the coverage levels. That is, the range of accuracy was 89-96% in the WGS sample (Fig. 5), 92-96% in the GBS when total genotypes were considered (Table 1A and Fig. 5) and 82-100% (Table 1B) when the genotypes per sample were considered. More importantly, the present data show that the coverage has a strong effect on the number of genomic positions for which a genotype can be inferred with genomic DNA.
The present work therefore suggests that, when genotype inference for all or most genomic regions of interest is needed, such as in genotyping for service, a deeper sequencing is required. In our analysis, a coverage of 259 provided 100% of SNPs, exceeding the threshold that  minimized the number of incorrectly assigned genotypes, as shown in Table S3. Therefore, we suggest a coverage depth of 259 when maximal accuracy is required. Conversely, when exploratory analysis is performed, such as in variant discovery or population screenings, a shallow sequencing can be enough to identify the genotype at most of the targeted genomic positions.
These results echo other studies in which GBS was used to genotype known variant sites. The coverage level considered varied depending on the studies (minimum two reads per site in Torkamaneh et al. 2016; minimum five reads per site in Ariani et al. 2016) and no consensus was found. However, the present trial showed that such low coverage levels often preclude the determination of genotypes with satisfactory performances. Additionally, Brouard et al. (2017) showed results consistent with those reported in the present work, with filtering by coverage and by genotype significance leading to better performances. This is the first work aiming to identify anonymous (suitable for DNA profiling) and non-anonymous loci at the same time using a targeted NGS approach in cats. The comparatively high throughput, the high sensitivity and the opportunity to combine different types of variants (SNPs and indels) and samples from different species, thanks to the barcoding, together with the exhaustive performance characterization performed in this work, all suggest that GBS can be used proficiently in population screening of mutations and, at deeper coverage, in DNA diagnostic testing.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article.