Fine mapping and identification of a candidate gene for a major locus controlling maturity date in peach

Background Maturity date (MD) is a crucial factor for marketing of fresh fruit, especially those with limited shelf-life such as peach (Prunus persica L. Batsch): selection of several cultivars with differing MD would be advantageous to cover and extend the marketing season. Aims of this work were the fine mapping and identification of candidate genes for the major maturity date locus previously identified on peach linkage group 4. To improve genetic resolution of the target locus two F2 populations derived from the crosses Contender x Ambra (CxA, 306 individuals) and PI91459 (NJ Weeping) x Bounty (WxBy, 103 individuals) were genotyped with the Sequenom and 9K Illumina Peach Chip SNP platforms, respectively. Results Recombinant individuals from the WxBy F2 population allowed the localisation of maturity date locus to a 220 kb region of the peach genome. Among the 25 annotated genes within this interval, functional classification identified ppa007577m and ppa008301m as the most likely candidates, both encoding transcription factors of the NAC (NAM/ATAF1, 2/CUC2) family. Re-sequencing of the four parents and comparison with the reference genome sequence uncovered a deletion of 232 bp in the upstream region of ppa007577m that is homozygous in NJ Weeping and heterozygous in Ambra, Bounty and the WxBy F1 parent. However, this variation did not segregate in the CxA F2 population being the CxA F1 parent homozygous for the reference allele. The second gene was thus examined as a candidate for maturity date. Re-sequencing of ppa008301m, showed an in-frame insertion of 9 bp in the last exon that co-segregated with the maturity date locus in both CxA and WxBy F2 populations. Conclusions Using two different segregating populations, the map position of the maturity date locus was refined from 3.56 Mb to 220 kb. A sequence variant in the NAC gene ppa008301m was shown to co-segregate with the maturity date locus, suggesting this gene as a candidate controlling ripening time in peach. If confirmed on other genetic materials, this variant may be used for marker-assisted breeding of new cultivars with differing maturity date.


Background
Fruit ripening is a complex process that involves the coordinated regulation of many metabolic pathways which influence numerous traits such as colour, aroma and flavour. During ripening, softening of the fruit tissues occurs in parallel with the accumulation of sugars, acids, and volatile compounds. Together, these traits contribute to increased palatability. Other characters, such as colour, size, and maturity date, have been selected to offer products with improved consumer appreciation. For practical reasons, fruits are harvested before they are physiologically ripe and are thus often perceived by the consumer as of poor quality. Consequently, the ideotype pursued by modern peach breeding is characterized by an increase of fruit quality associated to an easy harvest (fruit with a firmer texture or very slow softening), in order to reach optimal fruit quality at consumption. Fruit maturity is an important trait and early fruit ripening could allow market growth by extending the length of the production season. Studies in a wide range of plant species have provided insights into the genetic mechanisms that mediate fruit ripeningrelated processes, such as pigment synthesis, cell wall and sugar metabolism [1][2][3][4].
Climacteric fruits such as peach are characterized by an increase in respiration and an autocatalytic burst of ethylene production late in fruit development which is essential for normal fruit ripening. Blocking ethylene synthesis or perception prevents ripening [5][6][7]. Some tomato mutants (e.g. rin, nor, cnr) are unable to ripen even when treated with exogenous ethylene, although not impaired in the hormone signal transduction pathway [8; and references therein]. Molecular analyses of these mutants and identification of the corresponding genes indicate that they are necessary for the expression of both ethylene-dependent and -independent genes during ripening, acting upstream of and possibly in parallel to ethylene [8]. In banana, another climacteric fruit, NAC transcription factors have been shown to physically interact with a downstream component of ethylene signaling, ethylene insensitive 3 (EIN3)-like protein, which was down-regulated during ripening [9]. Thus, the control of NAC TFs on ripening is important as that regulated by the hormone.
In peach, QTLs controlling MD have been mapped on different chromosomes, with major QTLs located on linkage group (LG) 4 and 6 [10][11][12][13]. By analyzing the Contender x Ambra (CxA) F 2 population we have shown that the QTL detected on LG4 behaves as a Mendelian trait and has pleiotropic effects masking the effects of QTLs for different fruit traits [10]. In the present work we refer to this locus as qMD4.1. Recently, a major QTL was also identified in the collinear region of apricot (P. armeniaca L.) and sweet cherry (P. avium L.) suggesting that a common mechanism may control fruit maturation in related Prunus species [13].
In the present study, qMD4.1 was genetically dissected to gain insight into the mechanisms controlling MD in peach and to better understand the genetic interactions of this locus with other fruit traits. To this end, a fine map of the qMD4.1 locus was constructed, increasing the number of markers and individuals in the CxA population and analyzing another F 2 population where MD is segregating as a Mendelian trait. Analysis of the region located candidate genes for the MD QTL and a scan of publicly available re-sequencing data for the CxA F 1 parent identified polymorphisms in the region harbouring qMD4.1 [14]. A sequence variant in a NAC candidate gene was shown to co-segregate with the MD trait in these two populations and can be used in selection for early/late maturity genotypes.

Results and discussion
Fine mapping of the qMD4.1 locus The objective of this study was to fine map qMD4.1, the MD locus described by Eduardo et al. [10] on LG4 of the CxA F 2 population. This locus had been mapped using 169 individuals and a genetic map comprising 31 SSRs and two phenotypic markers, with flanking markers M12a and EPPISF032 defining a 6.6 cM interval, corresponding to 3.56 Mb on the peach reference genome. In order to fine map qMD4.1, the number of individuals of this population was increased to 306 and 14 SNPs were added around and into the original map interval. The SNPs were manually selected from the publicly available genome re-sequencing data of the CxA F 1 self-pollinated parent (NCBI Sequence Read Archive, biosample SRS335631). In parallel, we also analysed another F 2 population of 103 individuals derived from the cross between PI91459 (N.J. Weeping; W) and Bounty (By), also segregating for MD and genotyped with the 9K SNP peach array v1 [15] (Da Silva et al., unpublished).
In both populations MD showed a trimodal distribution ( Figure 1) and very high correlation between years (0.92) ( Table 1). This distribution suggested a Mendelian behaviour. Thus, a genome-wide QTL analysis for MD was carried out for confirmation using both populations independently, with the purpose of excluding the possible segregation of other MD QTLs (data not shown). Although segregating as a Mendelian trait in both the CxA and the WxBy F 2 populations, MD behaves as a quantitative trait in other populations such as the F 1 population derived from the cross between Bolero and OroA [10]. Figure 2 shows the region where qMD4.1 is located in both maps, CxA and WxBy, as compared to the data from Eduardo et al. [10]. In both maps high segregation distortion was observed in the target region ( Table 2). The analysis of the alleles showed that under-represented alleles came from the C genome in CxA and from W genome in WxBy, in both cases corresponding to the late ripening allele.
To define the target region for candidate gene analysis and to anchor the two maps we used the Sequenom platform to genotype the CxA F 2 progeny for 2 SNPs that were also included in the 9K SNP peach array v1 (i.e., SNP_IGA_385272 and SNP_IGA_415301; Figure 2).
Considering the region between the anchor markers SNP_IGA_385272 and SNP_IGA_415301 present in both populations, in the CxA map, the region spans 33.99 cM and includes 2 SSRs and 6 SNPs, whereas in the WxBy fine map this region includes 16 SNPs and spans 31.12 cM. The full map from WxBy will be described elsewhere and was used here to run a QTL analysis for the MD trait, showing no other QTL for MD segregating in this population (data not shown). Fine mapping analysis showed that in CxA, qMD4.1 is localized in a window of 2.13 cM between the markers PTP4 and PTP8 spanning the chromosome region between 10.87-12.09 Mb, whereas in the WxBy map qMD4.1 is located in a window of 0.98 cM between markers SNP_IGA_411637 and SNP_ IGA_412338 spanning the chromosome region between 10.97-11.19 Mb (Table 2). Moreover, in the CxA progeny LOD score in qMD4.1 region reaches a value of more than 220 (97% of explained variance) while in WxBy LOD score is about 50 with more than 90% of explained variance ( Table 2). Therefore this is the region where we focused the search for candidate genes.
Although the CxA F 2 population has a higher number of individuals, the WxBy F 2 population has a higher number of segregating genotyped markers, which provided a better resolution in the mapping of qMD4.1 ( Figure 3). These results showed the utility of the new genomic tools, such as the 9K SNP peach array v1 [15], to improve the efficiency and resolution in peach mapping studies. In addition, the use of the WxBy F 2 population allowed us to identify a relatively small region harbouring a limited number of candidate genes compared to the CxA F 2 population ( Figure 3).
QTLs for MD in a region overlapping with qMD4.1 were already identified by other groups [11,16]. For example, Dirlewanger et al. [13] mapped a major QTL for MD on LG4 using an F 2 population derived from the cross Ferjalou Jalousia x Fantasia (JxF): based on available data, qMD4.1 seems to overlap with their confidence intervals in different years, as supported by the SSRs UDP97-402 and AMPA 103, that in JxF defined this region and are located on the peach genome sequence in scaffold 4 at position 10,486,180 and 13,509,355, respectively (approximately 3.02 Mb apart). In contrast to the Mendelian behaviour of qMD4.1 in our populations, in these previous studies MD segregated as a quantitative trait with QTLs mapped on different chromosomes. Interestingly, these studies showed that the region spanning qMD4.1 hosted a cluster of QTLs associated with other traits such as fructose and sucrose content, titrable acidity, pH, sorbitol, maltose, citric acid, and quinic acid content [11]. Recently, an association mapping study showed that another locus mapped at the beginning of chromosome 4, associated with SSR marker CPPCT028 (at 2.1 Mb) was strongly associated with variation for the MD trait [17]. This indicates that other loci controlling MD can be identified using different segregating populations or association mapping. In addition, strong associations between MD, SSC, flavonoids, sorbitol, and total sugars levels were also found in this study [17].

Identification and analysis of candidate genes
Based on overlapping map positions of the qMD4.1 locus in the CxA and WxBy F 2 populations, our hypothesis was that allelic variation at the same gene could explain the segregation of the MD trait in both populations. Thus, we searched for genomic variations within the region 10.97-11.19 Mb identified by two WxBy recombinants ( Figure 3). Twenty-five genes are annotated in this region of the peach genome (Table 3) and we focused our attention on ppa007577m and ppa008301m, both predicted to encode NAC transcription factors (TFs). NACs constitute one of the largest plant TF families and are key regulators of developmental programs and stress response [18-20, and references therein]. Initially, we Sanger re-sequenced ppa007577m in the four parents C, A, W and By trying to identify polymorphisms within the gene, as well as in the 1.5 kb upstream and downstream flanking sequences. Sequence analysis did not show polymorphisms within the coding region, but revealed a 232 bp deletion, compared to the reference genome (Lovell), at 579 bp upstream of the start codon. PCR amplification of this region showed that C is homozygous for the reference Lovell allele, W homozygous for the deletion, whereas A and By are both heterozygous (data not shown). In addition, PCR analysis of the CxA F 1 parent for the 232 bp deletion showed that this individual was homozygous for the reference Lovell allele (C haplotype). Hence, no segregation was expected in the F 2 CxA population. To further verify this hypothesis, inspection of this deletion was conducted screening some individuals of the three phenotypic classes (early, middle and late ripening) of the CxA F 2 population. As expected from the parents' analysis, the deletion of 232 bp was absent in the analysed individuals. In contrast, the 232 bp deletion co-segregates with the phenotype in the selected individuals of the WxBy F 2 population being the F 1 parent heterozygous at this locus. Altogether, these results led us to search for others candidate genes for the MD trait because ppa007577m does not segregate with the MD trait in both populations. We thus focused our attention on the second NAC candidate gene ppa008301m. The Sanger re-sequencing of the four parents and of the two F 1 individuals uncovered an in-frame 9 bp insertion (compared to the reference genome) resulting in a tandem duplication of three amino acids in the C-terminal domain ( Figure 4). A, C, and By are heterozygous for this variant, whereas W is homozygous without the 9 bp insertion. Both CxA and WxBy F 1 individuals were heterozygous for this polymorphism and further analyses on the F 2 individuals were conducted to verify its correct segregation with the MD trait. To this end, the 9 bp insertion/deletion (INDEL) was scored using the genotyping protocol used for microsatellites analyses (data not shown). The allelic variants co-segregated with the MD of all the individuals in both F 2 populations, with the early ripening individuals showing the 9 bp insertion and the late ripening individuals carrying the Lovell reference allele (Additional file 1). To further verify the map position, we added these data to the previous genotyping dataset used for map construction and QTL analysis. We found that this INDEL co-mapped with the MD phenotypic marker (Additional file 2). Finally, repeating the QTL analysis in both populations showed co-location of this INDEL with qMD4.1 (Additional file 2).
As an additional confirmation that this 9 bp INDEL found in the last exon of the NAC gene ppa008301m is the most probable polymorphism within the qMD4.1 region, we searched the publicly available paired-end whole-genome next-generation sequencing data of the CxA F 1 parent for sequence variants within the original target region. Table 4 reports the results of this analysis. Most of the identified polymorphisms were located in intergenic regions. Only two polymorphisms were found within or near coding regions: the first in the 5′-UTR of a β-glucosidase 18 gene (ppa025660m) involved in carbohydrate and phenylpropanoid metabolism, and the second being the 9 bp insertion in the NAC candidate gene ppa008301m. Even if polymorphisms in intergenic regions could cis-regulate the expression of target genes, the presence of the 9 bp insertion in the coding sequence  of this NAC candidate gene supports the hypothesis that the NAC TF ppa008301m is a strong candidate gene for the control of MD in the CxA and WxBy populations. A phylogenetic analysis of NACs identified in Arabidopsis, rice, wheat, potato, banana and tomato indicates that NAC genes ppa007577m and ppa008301m belong to two distinct clades ( Figure 5). NAC ppa007577m clustered with tomato NOR (non-ripening) known to be involved in fruit ripening (GeneBank accession n°AY573802). First mapped on chromosome 10 by Giovannoni et al. [21], the nor gene seems to have a global effect on fruit ripening, probably acting upstream of ethylene synthesis [22]. Related Arabidopsis genes AtNAC2 and AtNAM (namely NARS1 and NARS2, respectively) ( Figure 5) control embryogenesis by regulating the development and degeneration of ovule integuments [23]. The clade of ppa007577m includes also rice ONAC10 and wheat NAC TFs NAM_B1 and NAM_B2 ( Figure 5). NAM-B1 and NAM-B2 have been shown to be involved in delaying seed ripening with pleiotropic effects on grain protein, zinc, and iron content [24]. In contrast, functional analysis of ONAC10, initially selected by the authors as the most similar rice NAC TF to the wheat NAM_B1, excluded a role in seed ripening while showing that it is essential for correct anther development [20]. Altogether, genes from the ppa007577m clade seem to have undergone a functional diversification, playing roles in distinct biological processes. The clade of NAC TF ppa008301m included genes mainly involved in stress responses as Arabidopsis ANAC072, ANAC055 and ANAC019 [25]. Recently, it has been shown that the closely related ANAC019 and ANAC055 play different roles during developmental senescence [26]. In fact, ANAC019 may be involved in triggering senescence by controlling and activating flavonoid and anthocyanin biosynthesis, whereas ANAC055 seems to be involved in the response to chitin [26]. In addition, these genes showed opposing roles in the regulation the jasmonic and salicylic acid pathways   Figure 4 Sequence alignment of the two variants of the NAC TF ppa008301m. ppa008301m_WT indicates the protein sequence without the deletion (identical to the reference genome) whereas ppa008301m_MUT indicates the three amino acids insertion found in this study (red box). The black box indicates the NAC domain. [26]. Unfortunately, literature about the involvement of NAC TFs in ripening is at present very scanty. In banana gene expression profiles in fruit with four different ripening characteristics revealed that MaNAC genes are expressed differentially in peel and pulp during ripening [9]. Additionally, MaNAC1 and MaNAC2 physically interact with MaEIL5 (involved in ethylene signalling), suggesting that these TFs could be involved in banana fruit ripening via interaction with ethylene signaling components.
Our results also support the involvement of ppa008301m in the control of peach fruit ripening, indicating that genes in this clade may have evolved to fulfill different biological functions as proposed for the ppa007577m clade. While orthologs are often assumed to retain equivalent functions in different organisms, deviations from this situation are not uncommon [9,27]. Being NACs one of the largest TF families, indication of the precise function of each member should be experimentally validated. The 9 bp INDEL in candidate gene ppa008301m results in the duplication of a threonine-aspartic acid-proline stretch with possible impact on protein function. Simple amino acid repeats and regions rich in serine and threonine, proline and glutamine, or acidic residues are frequent in the C-terminal region of NAC proteins [28][29][30][31][32][33]. While few studies were focused on NAC post-translational modifications, lysine 268 and 269 in tomato SlNAC1 C-terminal region contribute additively to its degradation, whereas the substitution of both lysine residues with two aspartic acids significantly stabilized the SlNAC1 C-terminal protein sequence [34]. In addition, the C-terminal regions of several NAC proteins function as transcriptional activation domains [19; and references therein]. Domain deletion analysis on the Arabidopsis NAC gene ATAF1, revealed that the transactivation activity was conferred by the C-terminal domain and did not require the N-terminal domain [35]. Hence, the three amino acid INDEL in ppa008301m, possibly modifying the C-terminal domain structure, might influence the protein stability or its transcriptional activation capability and, consequently, the transcriptional activation of target genes.
Interestingly, we found in both populations that the early ripening individuals possessed the three amino acid insertion in the protein sequence, whereas delayed ripening was associated with the reference allele.

Conclusions
In this study we report the fine mapping of a major locus controlling maturity date trait in peach, using two different populations. The strategy allowed us to restrict The NAC gene ppa008301m was identified as a strong candidate within this region and a consistent 9 bp insertion in its last exon was proposed as the variant possibly causing early ripening. Further experiments are needed to explore the functional significance of this INDEL. In any case, the marker developed on this sequence polymorphism provides a convenient molecular tool to discriminate early vs. late ripening individuals within the CxA, WxBy and likely other breeding populations. Considering the proposed pleiotropic effect of qMD4.1 on other important fruit traits, our results provide a basis to better define breeding programs for fruit quality improvement.

Plant material and phenotyping
Two F 2 populations were used in this study: a population of 306 individuals derived from the cross between Contender (C) and Ambra (A), and a second population of 103 individuals derived from the cross between N.J. Weeping (W) and Bounty (By). Trees were located in the ASTRA orchards in Castel S. Pietro and Tebano (Emilia Romagna, Italy). Trees of CxA and WxBy progenies were planted on their own roots with a spacing of 1 m within and 4 m between rows and trained as slender spindle (one stem with short lateral scaffolds). Pruning was performed yearly and standard cultural practices were applied. The fruits were thinned before pit hardening to a load of only 30-40 fruits per tree according to vigour, in order to allow a full expression of fruit size not limited by competition. Twenty fruits per tree were harvested at commercial maturity based on visual colour change and manual evaluation of firmness and the date was recorded. Maturity date (MD) phenotyping was obtained in years 2007-2008 and 2010-2011 for CxA and WxBy, respectively. The MD is defined in Julian days at harvest.

DNA extraction and genotyping
DNA extraction from young leaves was carried out using the DNeasy 96 Plant Kit (QIAGEN). Microsatellite (SSR) amplifications were performed following the multiplexready PCR protocol [10].

SNP selection, design and genotyping of the CxA population
For the CxA F 2 population, based on the confidence interval identified by Eduardo et al. [10], fine mapping was conducted on the 306 individuals with 16 newly developed SNPs spanning the interval 9-12 Mb on LG4. Putative SNPs were manually selected among those available from sequencing of the F 1 CxA and the relevant surrounding sequences (300 bp) were downloaded from the peach Gbrowse available on the IGA website (http:// www.appliedgenomics.org/). The Mass ARRAY Assay Design 3.1 software was used to design multiplex reactions in which the selected SNPs were included. Genotyping was performed using iPLEX Gold technology [36] and Mass ARRAY high-throughput DNA analysis mass spectrometry (Sequenom, Inc). All the 14 selected SNPs gave high quality results and were used for further analyses. Information about the primers used for the genotyping is listed in Additional file 3.

SNP selection, design and genotyping of the WxBy population
The 103 individuals from the WxBy F 2 population were genotyped using the recently developed Illumina 9,000 SNP array v1 for peach [15]. For SNP array genotyping, DNA was extracted with the DNeasy 96 Plant kit (Qiagen), diluted to 50 ng/μl and sent to IASMA Research and Innovation Centre (San Michele all'Adige, Italy) for genotyping. Genotyping was performed following the manufacturer's recommendations as described in Verde et al. [15]. SNP data were scored using GenomeStudio Data Analysis software (Illumina Inc.) using a GenCall threshold of 0.15. SNPs with GenTrain score < 0.6 and those showing severe segregation distortion (χ2 test, p < 10 −6 ) and more than 1% of missing data were excluded.

Genetic linkage maps
Genetic linkage analysis and map construction were performed with JoinMap 4 [37]. CxA and WxBy maps were produced and analyzed as "F2" population in the JoinMap 4 software. The recombination threshold value was set at 0.40 and the Kosambi mapping function was used to convert recombination frequencies into map distances. Markers showing distorted segregation were included in the linkage analysis. All linkage groups of both CxA and WxBy maps were calculated at a minimum LOD score of 5. Linkage maps were drawn using the MapChart 2.1 software [38].

QTLs analysis
QTL analysis was carried out using the software MAPQTL version 6.0 [39]. After a test of 10,000 permutations, LOD thresholds of 3.9 in 2007 and 3.8 in 2008 in CxA map, and 3.6 in 2010 and 3.5 in 2011 in WxBy were assessed to accept QTL significance. A stringent significance level of p = 0.005 was adopted as threshold for the detection of a QTL for the individual test in order to obtain an overall significance level of about p = 0.05, as previously suggested by Van Ooijen [39]. QTLs were drawn using the MapChart 2.1 software [38]. For each plant of both F 2 populations, maturity date phenotypic data and genotypes of the 9 bp insertion/deletion in the gene model ppa008301m are reported in Additional file 1. In order to visualize the genetic map and QTL features, QTL-LOD profiles were rendered in heatmaps produced by Harry Plotter software developed at the Edmund Mach Foundation (http:// genomics.research.iasma.it/download.html).
Discovery and genotyping of the NAC ppa007577m and ppa008301 9 bp INDELs Sanger sequencing for the candidate genes ppa007577m and ppa008301m was conducted spanning the region 11,126,755-11,131,088 and 11,105,608-11,107,907, respectively, with several primer combinations as described in Additional file 4. The 232 bp deletion in gene ppa007577m was identified by amplifying with the forward and reverse primers 5′-CTACTCATACCCGCCAAGGA-3′ and AACGTCGTCATGAGGTACCC, respectively. PCR reactions contained 1-20 ng of genomic DNA, 1x PCR reaction buffer ( Detection of polymorphisms in CxA F 1 parent from whole-genome re-sequencing data Publicly available paired-end whole-genome re-sequencing data of Prunus persica accessions (study SRP013437) were downloaded from the NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra, [41]. For this study we considered biosamples SRS335634 (Lovell Clone PLov2-2 N, run SRR502985), and SRS335631 (F1 Contender x Ambra, run SRR502997). SRA data of each run were dumped in fastq format using the fastq-dump tool of the NCBI sratoolkit v.2.1.16 software (http://www.ncbi. nlm.nih.gov/Traces/sra/sra.cgi?view=software); forward and reverse paired reads were split for each sample into two separate files (option -split-files). Reads were quality trimmed by Trimmomatic v.0.22 [42], trimming leading and trailing bases below a quality threshold of 20, also removing reads having an average quality below 20 (calculated on 8 bp long sliding windows) and trimmed reads shorter than 24 bp. For each sample, only reads passing the quality filtering as matching pairs were retained and aligned to the whole Prunus persica reference genome v1.0 (http://www.rosaceae.org, [14]) using the Burrows-Wheeler Alignment Tool (BWA) v.0.6.2 [43]. The aln (IS linear-time algorithm) and sampe (all default options except -n 25 -N 25) commands were applied, respectively, for constructing suffix array (SA) coordinates of good hits of each individual read and to convert them to chromosomal coordinate and pair the reads. The resulting SAM files were converted to sorted BAM files compliant to the Genome Analysis Toolkit (GATK) format by Picard Tools v.1.77 (http://picard.sourceforge.net/) using, in the order, tools CleanSam, SamFormatConverter and AddOrReplaceReadGroups. GATK-compliant BAM files were submitted to GATK v.2.3-3 [44] for preprocessing procedures, consisting of indel realignment, duplicate removal and base quality score recalibration (BQSR). The data table needed for the recalibration step in BQSR was manually generated upon validated SNP data from the Peach 9 K chip array [15]. Variant discovery procedures were then applied using whole-genome recalibrated alignments across both CxA and Lovell sample simultaneously using GATK HaplotypeCaller tool applying standard hard filtering parameters [45].
Phylogenetic analysis revised the manuscript); LD (genotyping; analysed data; revised the manuscript); GP (analysed data); DB (selected the genetic materials and developed the populations; revised the manuscript); LR (conceived and designed the study; analysed data; helped writing the manuscript). All authors read and approved the final manuscript.