Footprints of domestication revealed by RAD-tag resequencing in loquat: SNP data reveals a non-significant domestication bottleneck and a single domestication event

Sequencing data

SNP calling and its distribution pattern

Population genetic structure and genetic differentiation

Phylogenetic analysis

Unrooted phylogenetic tree of 19 samples showed that nine samples of wild group clustered together at one end of the tree, and 10 cultivated samples were clustered at the other end of the tree (Fig. 1). In the wild group, three samples from Bijie county, Anlong county, and Meitan county, Guizhou province clustered into close proximity to each other, and the sample from Shimian county, Sichuan province was also close to the aforementioned accessions, indicating that the genetic distance of wild loquat population had some regional relationships. However, the wild samples from Hubei and Shaanxi province were more close to cultivated loquat samples and the closest one was TS_W. Interestingly, the cultivar, Sijipipa, exhibited the closest relationships with the wild loquat population, and the cultivars Dawuxing, Yunan and Biano were clustered far from the wild loquat population.

PCA analysis

The results of principal components analysis of the 19 genotypes showed that the first principal component could explain 5.09% genetic variation of whole loquat population, and the second principal components can explain 3.38% genetic variation of whole loquat population (Fig. 2). The first and the second principal components could clearly separate wild and cultivated loquats, and the results showed that the wild and cultivated loquat have entirely different evolutionary trends.

Population structure

Population structure helps us to understand the evolutionary process of a population or species through the association between genotype and phenotype, and to determine the groups or subgroups of different individuals/populations. When k was set as 2, Al_W, MT_W, BJ_W, SM_W, LC_W and BD_W were clustered into one wild group, and nine cultivated loquat cultivars were clustered into cultivated group, and HZ_W, WF_W, TS_W and SJPP_C were assigned to aforementioned two groups. However, HZ_W, WF_W and TS_W have a higher probability to be included in wild group, and the SJPP_C has a higher probability to join cultivated group. When k was set as 3, Al_W, MT_W and BJ_W from Guizhou province constituted a group, MBC_C, Biano_C, YN_C and DWX_C were clustered into one group, and BD_W, HZ_W, WF_W TS_W and SJPP_C were assembled into another group, while SM_W, LC_W, BY_C, ZZ6H_C, MM_C and Algeria_C were clustered into two separate groups. When k was set as 4, three samples collected from Guizhou province were further divided into two groups, and the clustering pattern of other samples was almost consistent with K = 3; When k was set as 5, the wild samples were divided into three groups with three samples in each group, and the cultivated loquat cultivars was clustered into two groups (Fig. 3).

The representativeness of samples

According to the extensive survey, the wild loquat (Eriobotrya japonica Lindl.) populations have been found in Hubei, Sichuan, Yunnan, Guizhou, Guangxi and Guangdong provinces of China. Among them, Hubei and Guizhou provinces are the major distributing areas of wild loquat [48]. In the present study, a total of 9 wild accessions and 10 loquat cultivars were selected for sequencing. Among 9 wild loquat samples, four samples were collected from different regions of Hubei province, three samples from different regions of Guizhou province, one from Sichuan province, and one from Shanxi province. All these nine samples were regional representatives of species distribution range and separated from each other by far geographical distance. Among 10 loquat cultivars, six are the main cultivars of different loquat producing areas of China, such as Zhejiang province, Fujian province, Jiangsu province, Sichuan province, Yunnan province and Guangdong province, and other four are the main cultivars of Japan, Italy, Spain and the United States of America. Molecular pedigree analysis showed that those cultivars almost represent different genetic clusters of current loquat cultivars around the world [25, 26].

Genetic bottleneck in cultivated loquat during domestication

In the process of crop domestication, the genetic bottleneck is common, which indicate that the effective size of crops population were significantly fewer than its corresponding wild progenitor population [49]. Compared to the nearest wild progenitors, crops population appeared to have lower genetic diversity, higher level of linkage disequilibrium and lower SNP frequencies, such as rice [38, 39], soybean [11, 50], cucumber [51], peach [52], and tomato [53]. However, there was no clear genetic bottleneck in several crops, such as grape [54]. In this study, we found that the SNP frequencies in wild loquat population were only 2.4% higher than that in cultivated loquat population. Meanwhile, AMOVA analysis showed that 85.55% of genetic variants existed among- and inter-individuals, while only 14.45% happened between wild and cultivated populations. The independent samples t-test exhibited non-significant differences between cultivated and wild loquat groups (P > 0.05) for the frequency of total number of SNPs, homozygous SNPs and heterozygous SNPs. Above all results suggested that a non-significant genetic bottleneck appeared during the domestication of loquat. We deduced following major reasons for these results: 1) the domestication history of loquat is not so long; 2) loquat domestication is not yet complete, because there are non-significant differences for physiological and morphological traits between cultivated and wild loquat except that the fruit size of cultivated loquat is bigger than that of wild loquat.

The origin center of cultivated loquat

Discovering the origin center of a crop is not only an interesting scientific question, but also an important social science issue, because the ancient civilization of mankind is actually a farming civilization, and the action of crops domestication represents the height of the agricultural civilization to a certain extent. In fact, most crops were thought to have originated in the more developed areas of ancient civilizations [2, 9]. The valley of Qingshuijiang river of South-West Hubei province or valley of Daduihe river of South-West Sichuan province have been inferred as the original place of cultivated loquat according to the field investigations of morphological traits and geographical distribution [55, 56]. In this study, the sample, TS_W, collected from Hubei province was the most closest to the cultivated loquat according to the phylogenetic tree. Structure analysis also showed that this genotype was clustered with some cultivated loquat cultivars, such as SJPP_C (K = 2, 3, and 4). However, the wild sample, SM_W from Sichuan province, was far from the cultivated loquat based on above mentioned analysis. According to the archaeological work, in the 1970s, when people dug an ancient tomb that was constructed in Han dynasty, which is located in Jiangling county of Hubei province and near to Tongshan county, the carbonization relic of loquat seed was found among many funeral items, and this was the earliest archaeological evidence of loquat cultivation [13]. So, based on molecular and archaeological evidence, we deduced that the Hubei province, not Sichuan province, might be the origin center of cultivated loquat.

A single domestication origin of cultivated loquat

A crop that comes from a single domestication event or multiple domestication events is often present dispute based on the ancient book record and archaeology, which was used in corresponding research in the past [57–59]. However, modern molecular phylogeny, molecular population genetics and phylogeographical studies showed that most crops were evolved from single domesticated event, such as einkorn wheat [60], potato [61], soybean [50], maize [62], japonica rice [10], and cucumber [63]. Multi-domesticated events were not so common and it was observed in common bean [64] and barley [65].

China has about 2000 years long history for loquat cultivation until the twelfth century. Loquat was introduced into Japan from China and then from China and Japan to other parts of the world [19]. In the process of artificial selection and ecological adaptability, cultivated loquats differentiated into different cultivars or strains. Liu et al. [66] revealed three groups of cultivated loquat by analyzing 100 morphological and biochemical traits. Cultivated loquat of China was divided into two categories, namely as subtropical and tropical groups according to the regional ecological adaptability [67]. Two ecological groups of cultivated loquat were found according to the plant morphology such as leaf color, fruit size, and sugar content [48]. These analyses suggested that cultivated loquat may have been domesticated at multiple places, or may have been domesticated at just one place, and then spread outside along different routes, and differentiated into different groups. In the present study, structure, phylogenetic and PCA analyses have shown that all the cultivated loquat samples were clustered into the same group, which showed that cultivated loquat samples exhibited mono-phylogenic origin. So, we concluded that cultivated loquat was originated from a single domestication event. However, cultivated loquat cultivars exhibited differentiation by structure analysis (k = 5). One group, SJPP_C, RTBS_C, BY_C, ZZ6H_C and MM_C, might be comprised of early domesticated cultivars, and another group involving Biano_C, YN_C and DWX_C might be domesticated recently, while Algeria_C and MBC_C clustered between these two groups.

Plant sampling and DNA extraction

The leaves of ten cultivated loquat cultivars were sampled from Horticultural germplasm conversation center of South China Agricultural University (SCAU), where most of the major cultivated loquat cultivars of the world have been planted, and nine wild loquats were sampled from native wild loquat populations representing their natural distribution range (Table 6; Fig. 4). The DNA was extracted according to CTAB method [68].

Code	Genotype	Geographical origin
DWX_C	Dawuxing	Sichuan province, China
BY_C	Baiyu	Jiangsu province, China
YN_C	Younan	Guangdong province, China
SJPP_C	Sijipipa	Yunnan province, China
RTBS_C	Ruantiao baisha	Zhejiang province, China
ZZ6H_C	Zaozhong No.6	Fujian province
Biano_C	Biano	Italy
MBC_C	MBC	America
MM_C	Mogi	Japan
Algeria_C	Algeria	Spain
TS_W	Wild loquat	Tongshan county, Hubei province
WF-W	Wild loquat	Wufeng county, Hubei province
BD_W	Wild loquat	Badong county, Hubei province
LC_W	Wild loquat	Lichuan county, Hubei province
MT_W	Wild loquat	Meitan county, Guizhou province
BJ_W	Wild loquat	Bijie county, Guizhou province
AL_W	Wild loquat	Anlong county, Guizhou province
SM_W	Wild loquat	Shimian county, Sichuan province
HZ_W	Wild loquat	Luding county, Sichuan province

The construction of RAD-tag libraries

Mixed library of RAD-tag (restriction-site associated DNA tags) was constructed in this study by the protocols described by Etter et al. [69], and it is as follows: (A) Restriction enzyme EcoRI digested about 1μg genomic DNA and joint P1 adapter ((P1-FOR-xxxxx: ′ -Phos- AATGATACGGCGACCACCGAGATCTACACTTTCCCTACACGACGCTCTTCCGATCTXxxxxTGC*A-3 ′ P1-REV-xxxxx: 5 ′ -Phos-xxxxxAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGA; TCTCGGTGGTCGCCGTATCAT*T-3 ′) which contains sequences of primers P1-PCR (5 ′-AATGATACGGCCCACCGA-3 ′) for amplification, primer binding sites in Illumina Genome Analyzer and short tags to distinguish samples; (B) Samples with different adapters were mixed to be broken into segments of 300 ~ 700bp using physical method; (C) Joint P2 adapter (P2-FOR:′ -Phos-GATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCAGAACAA-3 ′P2-PE-REV: 5′ -Phos CAAG CAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCG ATC*T-3 ′, including primer sequences: P2-PCR, 5′ -AATGATACGGCGACCACCGA-3 ′), recovery; (D) PCR: 5-10 rounds of amplification target sequence (enrichment and sequencing primers for P1-PCR\P2-PCR).

Sequencing, raw data processing and SNP genotyping

Illumina Solexa sequencing was done by sequencer Hiseq 2000, and sequencing type was SE50. We make sure that clean data were at least more than 0.6 G (0. 8 × C) for each individual. After sequencing, the raw data were processed in three steps: first, we allocated the raw data to its own origin sample according to the labels (4 ~ 8bp) used for sequencing to distinguish samples (4 ~ 8bp); secondly, we filtered label sequences and the joint to exclude the pollution; thirdly, we discard those raw data with low quality base (Q ≤ 5 (E)) number accounted for more than half of the whole reads.

The rest clean and high quality raw data were used for further SNPs calling by software Stacks (http://catchenlab.life.illinois.edu/stacks). Software Stacks uses short-read sequence data to identify and genotype loci in a set of individuals either de novo or by comparing to a reference genome. From reduced representation Illumina sequence data, such as RAD-tags, Stacks can recover thousands of single nucleotide polymorphism (SNP) markers useful for the genetic analysis of crosses or populations [70]. We used SPSS 24.0 (https://spss.en.softonic.com/) to execute the significance test of difference in average frequency of total SNPs, homozygous SNPs and heterozygous SNPs between wild and cultivated loquat by using T test of independent samples. For further analysis, we filtered the genotype data with the loss rate less than 50% (Additional file 1: Table S1) or with no loss (Additional file 1: Table S2).

Population genetics and evolutionary analysis

We used Arliquin 3.11 software [71] to execute the AMOVA and genetic differentiation analysis using SNP dataset with no loss/gap (Additional file 1: Table S2). Neighbor-Joining method [72, 73] was employed to construct the phylogenic relationships between loquat genotypes based on SNP dataset with no loss/gap (Additional file 1: Table S2) by using program MAGE 5.0 [74], and removed all the ambiguous positions for each sequence pair and set bootstrap = 500. The PCA analysis was executed by software EIGENSOFT 3.0 [75] using SNP dataset with the loss rate less than 50% (Additional file 1: Table S1). The software Structure v2.3.4 [76] was used to analyze the genetic structure using SNP dataset with the loss rate less than 50% (Additional file 1: Table S1) with admixture model, and K was set from 2 to 6.