NEW INSIGHTS INTO THE PHYLOGENY AND SPECIATION OF KUMQUAT (FORTUNELLA SPP.) BASED ON CHLOROPLAST SNP, NUCLEAR SSR AND WHOLE-GENOME SEQUENCING

Highlights

● Fortunella genus consists of two populations: cultivated kumquat and wild Hong Kong kumquat.

● Artificial selection might involve in the origin of cultivated Fortunella species.

● A hypothesis for the differentiation and speciation of Fortunella species is proposed.

Abstract

Kumquat (Fortunella spp.) is a fruit and ornamental crop worldwide due to the palatable taste and high ornamental value of its fruit. Although Fortunella is classified into the economically important true citrus fruit tree group together with Citrus and Poncirus, few studies have been focused on its evolutionary scenario. In this study, analysis of five chloroplast loci and 47 nuclear microsatellites (nSSR) loci from 38 kumquat and 10 citrus accessions revealed the independent phylogeny of Fortunella among citrus taxa, and that Fortunella mainly comprises two populations: CUL, cultivated Fortunella spp. (F. margarita, F. crassifolia and F. japonica); and HK, wild Hong Kong kumquat (Fortunella hindsii). Genomic analysis based on whole-genome SNPs indicated that the allele frequency of both pupations deviated from the neutral selection model, suggesting directional selection was a force driving their evolutions. CUL exhibited lower genomic diversity and higher linkage strength than HK, suggesting artificial selection involved in its origin. A high level of genetic differentiation (Fst = 0.364) was detected and obviously asynchronous demographic changes were observed between CUL and HK. Based on these results, a new hypothesis for the speciation of Fortunella is proposed.

Graphical abstract

Keywords

Citrus / Fortunella / kumquat / phylogenetics

1 INTRODUCTION

Broiler bodyweight is a key indicator in the commercial rearing process. Monitoring the weight of a broiler flock and integrating the data in farm management benefit the sustainability and poultry welfare improvement^[1]. The obtained real-time weight (i.e., bodyweight) can be compared with expected growth curve and for determining the best finishing day to avoid unnecessary feed consumption and reduce cost^[2,3]. Farmers can also respond to abnormal weight gain and uniformity change events by promptly investigating the causes and accordingly improve management by changing feeding practices, enhancing the environmental control and conducting disease prevention^[4,5]. The weight monitoring of broiler flocks will therefore benefit feeding management and the maintenance of broiler health as a guide for best management^[6].

The most common procedure for weight measurement is a manual process, where chickens are manually caught one by one and weighed on electronic scales with a precision of ± 20 g^[1]. The protocol requires a suitable timing and location for sampling, and a sample size of about 2% of chicken population. This method also requires the chicken to be left on a scale for a period before a stable weight can be read. In the cage-free environment, chickens are mobile, random and difficult to catch^[7–9], which makes the procedure time-consuming and is a source of stress to chickens^[10,11]. In comparison, automatic weighing systems can accurately monitor weight in real-time while avoiding stress responses. Therefore, there is a clear need to develop an automatic scale system for monitoring weight of broiler chickens in commercial. Automatic weighing systems are now increasingly used to assess average weight in studies such as flock growth control^[12], weight gain prediction^[13] and finishing weight prediction^[14].

Automatic chicken weighing research started in the 1980s and achieved an accuracy approaching 97.5%, but the data pretreatment and analysis was insufficient, and wide testing trials for different rearing stages and conditions were not sufficiently conducted^[15–17]. Mollah et al.^[18] in 2010 used digital images acquired with a standard camera to determine top-view surface area, and established a linear regression model between the area and predicted weight. Mortensen et al.^[19] in 2016 used a depth camera for 3D modeling to calculate weight by automatically capturing images of chickens in commercial flocks, but the results had an error of 7.8%. The machine vision method has a large error in measurement results due to the small sampling size and variable body posture of chickens as well as the complex rearing environment. Fontana et al.^[20] in 2015 used sound analysis to automatically determine the bodyweight of broilers and to establish a regression model between broiler calls and bodyweight but it is not technically mature at present.

In response to the shortcomings of the existing technology, this study developed an intelligent system for the weight monitoring of cage-free chickens based on chicken perch habit, and systematically evaluated the performance of the system in various commercial broiler houses. In a floor rearing environment, chickens tended to stand or sit on a weighing platform. The system was based on a weighing platform equipment integrated with weight data collection, transmission, storage, analysis and display. The equipment automatically collected the weight data at a sampling frequency of 5 s and uploaded the data to a cloud server in real-time through a multilevel wireless network. The data pretreatment and analysis module regularly processed the raw data to removed outliers and leave effective data, and analyzed the data to obtain the daily mean weight and uniformity.

2 MATERIALS AND METHODS

2.1 Automatic weighing system

The automatic weighing system was originally designed by the joint research group of the authors. It had been tested for the performance of proof-of-concept prototypes, and was modified and moved to the pilot test of engineering prototypes in this study. The hardware of the weighing platform (CZ-PYJ20PM-03 by Hefei Shenmu Information Technology Co., Ltd., Hefei, Anhui, China) included a weighing platform and communication terminal (Fig.1). Two strain-gauged sensors were installed under each weighing platform for weight sensing. The equipment had a weighing range of about 200 g to 20 kg, powered by 220 V AC with a rated power of 5 W and IP54 degree of protection. During the weighing procedure, weight sensors automatically read at 10 Hz and recorded a reading when a relative error of less was than 2% of the average value of the most recent 10 readings. The difference between two consecutive stable readings was considered a weight data point generally yielding 10 data points every minute. After a predetermined period of weight recording (24 h in this study), the raw weight data were subject to a protocol termed “Processing of Original Records of Weight Information” (PORWI) to detect and remove repeating weight data and outliers as data preprocessing (beyond a threshold of the pre-estimated average weight of a chicken breed on particular days; the thresholds were kept confidential). The rationales were that the weight of chickens was normally distributed, and most chickens were close to the average weight, the seriously lighter or heavier readings were considered to be noise or not representative, and were removed to reduce error, and a weight of two chickens or more deviated from the distribution with a weight gap. The detailed algorithm design was available in an earlier work of the joint research group^[21].

Fig.1 Weighing platform (a) and communication terminal (b) of the broiler weighing system.

Full size|PPT slide

The software (Fig.2) supported a real-time data transmission network and a user interfaces (smart phone WeChat applet and other software). Monitored weight data was transmitted to the software through the real-time built-in transmission network. This real-time transmission network consisted of self-assembled sensor network, local area network and remote transmission network, which realized data collection of multiple scales within monitoring areas. The self-assembled sensor network was able to communicate with a WeChat applet to provide practical outcomes such as configuration and on-site inspection of monitoring nodes by smart phones. After acquiring data, the data processing module of the software calculated and verified the data, and provided it to users through various functional modules. The WeChat applet is able to interact with the backend of the software through the network so that users can complete various queries and modifications of weight data and system management data.

Fig.2 Weighing system workflow (a) and screenshots of data display on the user interface (b).

Full size|PPT slide

2.2 Manual weighing using the live-bird sales system

Average weight of the birds in a flock obtained by automatic weighing system needed to be compared to standard average weight of the birds (as established by the established live-bird sales system). For the latter, the standard weight was manually collected as follows. Chickens were captured from the floor of production houses, transferred to transportation cages, and transported to the live-bird sales system, where chickens were transferred to standard weighing cages with about eight birds in each. Depending on customer demand, varying numbers of standard weighing cages were manually weighed before transaction. For comparative purposes, the weight data and numbers of chickens were recorded and converted to average weight, and designated as the standard chicken weight. The detailed process for the mean and variance estimation was presented in Section 2.4.1 below.

2.3 Weight data collection in commercial production facilities

From July, 2019 to April, 2020 across the four seasons of a year, 124 trials of chickens of different growth rate types and sexes were tested, and 101 of these trials collected sufficient valid data (over 100 data points) after data preprocessing. Of these, there were 29 trials of medium-growing (finished at about 8−10 weeks of age) and 72 trials of slow-growing (finished at about 12 weeks of age) chickens according to their growth rate. There were 50 trials with cocks and 51 trials with hens. One weighing platform was usually installed per 5000 chickens. The weighing platforms were installed away from feed and water lines to avoid interference from the activity of the chickens. The communication terminal was fixed at about 1.2 m high. To ensure the horizontal stability of the main body of the weighing platform, the equipment was placed on a flat and solid ground (Fig.3).

Fig.3 The installation of weighing platform (a) and communication terminal (b) in real applications.

Full size|PPT slide

2.4 Data analysis

2.4.1 Weight indicators

For the trials with automatic weighing platform, the mean weight (w_a) and standard deviation of the population were estimated based on each weight data point representing the weight of one chicken after PORWI protocol, w₁, w₂, ..., w_n, where n was about 2% to 5% of the total number of farm chickens in a trial. The weight uniformity (u_10%) was defined and calculated as the ratio of the number of chickens falling within ± 10% of the average weight of the sampled chickens to the total number of sampled chickens. Those indicators were calculated as:

Mean weight

(1)$w_a=\frac{\sum _{i=1}^n{w}_i}{n}$

Standard deviation

(2)${\sigma }_a=\sqrt{\frac{\sum _{i=1}^n{\left(w_i-w_a\right)}^2}{n}}$

Coefficient of variation

(3)$c_v=\frac{{\sigma }_a}{w_a}$

Weight uniformity (among w₁, w_2, ..., w_n, assuming there were in total m points of weight data within ± 10% of the mean weight $w_a$)

(4)$w_a\times \left(1-\text{10\%}\right)⩽w_{m_1}\text{,}{w}_{m_2}\text{,}\cdots \text{,}{w}_{m_m}⩽w_a\times \left(\text{1}\text{+}\text{10\%}\right)$

(5)$u_{\text{10\%}}=\frac{m}{n}\text{×}\text{100\%}$

The mean weight calculated from the information collected by the live-bird sales system was defined as standard average weight (w_s), calculated by dividing the total weight by the number of chickens. For each trial, the absolute and relative errors were calculated as:

(6)$E_a=w_a-w_s$

(7)$E_r=\frac{w_a-w_s}{w_s}\text{×}\text{100\%}$

It is important to note that on the live-bird sales system the study could not control the sales process such as the number of birds each customer would purchase, the standard deviation (σ_s) was calculated according to the protocol of multiple samplings and pooling standard deviation. Among the chickens weighed in cages on the live-bird sales system, and the weight data were divided into p categories according to the number of chickens weighed at once, as l₁, l₂, ..., l_i, ..., l_p (l_i ≠ l_j for i ≠ j). The number of weighing events for each of p category was a₁, a₂, ..., a_i, ..., a_p, therefore, with the weight of a single weighing event for weighing category of l_i defined as $r_{l_i\text{,}\text{1}}$, $r_{l_i\text{,}\text{2}}$, …, $r_{l_i,p}$.

Average weight of a chicken of a weighing event of a category l_i

(8)$w_{l_i,n}=\frac{r_{l_i,n}}{l_i}$

The mean value of the average weight of a chicken of a category l_i

(9)$\overline{w_{l_i}}=\frac{\sum _{n\text{=1}}^{a_i}{\text{w}}_{l_i\text{,}n}}{a_i}$

Variance of the sampling mean of the weight of a chicken of a category l_i

(10)$s_{l_i}^2=\frac{{\left(w_{l_i\text{,}1}-\overline{w_{l_i}}\right)}^2+{\left(w_{l_i\text{,}\text{2}}-\overline{w_{l_i}}\right)}^{\text{2}}+\text{…}\text{+}{\left(w_{l_i,a_i}-\overline{w_{l_i}}\right)}^{\text{2}}}{a_i-\text{1}}\text{=}\frac{\sum _{n=1}^{a_i}{\left(w_{l_i,n}-\overline{w_{l_i}}\right)}^2}{a_{\text{i}}-\text{1}}$

The estimated overall variance of the sampling weight of chicken flocks based on $s_{l_i}^{\text{2}}$

(11)${\hat{\sigma }}_{population,l_i}^{\text{2}}=l_i\text{×}{s}_{l_i}^{\text{2}}=\frac{l_i\text{×}\sum _{n=\text{1}}^{a_i}{\left(w_{l_i,n}-\overline{w_{l_i}}\right)}^{\text{2}}}{a_i-1}$

The estimated overall variance of the weight of chicken flocks with the method of pooled variance based on ${\hat{\sigma }}_{population,l_i}^{\text{2}}$

(12)$\begin{array}{rl}& {\hat{\sigma }}_{population}^2={\hat{\sigma }}_{pooled}^2\\ & =\frac{\left(a_1-1\right)\times {\hat{\sigma }}_{population,l_1}^2+\left(a_2-1\right)\times {\hat{\sigma }}_{population,l_2}^2+\dots +\left(a_p-1\right)\times {\hat{\sigma }}_{population,l_p}^2}{(a_1-1)+(a_2-1)+\cdots +(a_p-1)}\\ & =\frac{\sum _{i=1}^p\left[l_i\times \sum _{n=1}^{a_i}{\left(w_{l_i,n}-\overline{w_{l_i}}\right)}^2\right]}{\sum _{j=1}^p{a}_i-p}\end{array}$

2.4.2 Statistical analysis

Linear regression was used to model the relationship between absolute error, standard average weight and estimated weight by fitting a linear equation according to least-squares method. Analysis of variance (ANOVA) was used to analyze the differences among means of errors in terms of rearing season, chicken growth rate type and sex. A two-tailed t-test was used to test the equality of mean weight of methods. A two-tailed F-test was used to determine the equality of the variances of chicken weights from the automatic weighing platforms and the live-bird sales system.

3 RESULTS AND DISCUSSION

3.1 Performance of the automatic weighing platforms

3.1.1 Weight comparisons

Fig.4 shows the automatic weighing platform data as histograms before and after PORWI preprocessing protocol for three trials. The scales sensed and collected raw data at high frequency and may include noises induced by environmental disturbance such as chickens touching the weighing platform and pecking activities^[22,23]. Therefore, a high proportion of the raw data clustered in the first peak of the histograms. The PORWI protocol effectively removed those noises and repeating readings from the raw data for further calculation.

Fig.4 Histograms of automatic weighing platform data before (a–c) and after (d–f) POWRI preprocessing protocol for three broiler houses (a and d, b and e, and c and f for the same house, respectively).

Full size|PPT slide

For the three trials, weight estimation according to the raw data, data after PORWI, and data from the live-bird sales system were compared (Fig.5). It also shows the PORWI protocol was necessary and effective as the estimates of the average weight and standard deviation were not in accordance with the other two groups. For each trial, the estimates of average weight based on the automatic weighing platform and the live-bird sales system were close, i.e., 3085 vs 3077 g for the first trial, 2812 vs 2788 g for the second trial, and 2250 vs 2205 g for the third trial. The relative errors for the three trials were only 0.24%, 0.83%, and 2.0%, and all the absolute errors were less than 50 g. Results of t-test further confirmed the equality as it failed to reject the null hypothesis that the means were the same. Therefore, the automatic weighing platform can be used as an appropriate tool to monitor and estimate the average flock weight at the live-bird sales system, i.e., the standard average weight.

Fig.5 Average weight estimation according to the raw data, data after PORWI, and data from the live-bird sales system. Note: the error bars represent the estimated standard deviations.

Full size|PPT slide

The similarity of the standard deviations estimated from the automatic weigh platform and the live-bird sales system were compared according to the strict F-test as a tool of homogeneity of variance test. Of the three trials, the first and third passed the test, and the second failed. When looking into the calculated uniformity values of the second trial from the two approaches, u_10% was 52.5% based on the automatic weighing platform, and 29.0% (very low uniformity) based on the live-bird sales system. As the farmer did not observe or complain about the chicken grades of this particular flock, the value 52.5% seemed closer and realistic. The above analysis indicates that the automatic weighing platform provided more accurate estimates for uniformity of chicken weight compared to the data from the live-bird sales system. It may be attributed to the fact that the live-bird sales system did not control consumer need on chicken numbers which had an impact on standard deviation estimation; therefore, it yielded less accurate data. The automatic demonstrated its suitability and superiority in the estimation of standard deviation. The accurate estimation of standard deviation indicates accurate estimation of uniformity and coefficient of variation of chick flocks, both of which are important management indicators in chicken farming^[24].

3.1.2 Regression analysis

Among the 101 trials, the error between the estimated weight by automatic weighing platform and the standard average weight ranged between –109 and 94 g and averaged –10.3 g (σ = 39.5 g). The relative error ranged between –6.21% and 4.88% and averaged –0.54%, i.e., an average accuracy of 99.5% (σ = 2.3%). Pearson’s correlation analysis (Tab.1) further shows a strong (r = 0.990) and highly significant linear relationship between standard average weight and estimated weight. A small but significant correlation between the absolute value of error and the standard average weight, and a slight correlation between the absolute value of error and the estimated weight were found. This may be attributed to the sensitive range of weight sensors, or different activity patterns of small and large chickens^[25]. The scatter plot of standard average weight compared to the estimated weight by automatic weighing platform showed clear agreement (Fig.6).

Tab.1 Pearson’s correlation of standard average weight, estimated weight, and absolute value of error

	Standard average weight	Estimated weight	Absolute value of error
Standard average weight	1	0.990**	0.201*
Estimated weight	0.990**	1	0.155
Absolute value of error	0.201*	0.155	1

Note: **, correlation is significant at the 0.01 level; *, correlation is significant at the 0.05 level.

Fig.6 Scatter plots of standard average weight vs estimated weight (a), absolute value of error vs standard average weight (b), and absolute value of error vs estimated weight (c).

Full size|PPT slide

3.2 Stability of the automatic weighing platform

3.2.1 Monitoring of a growing period

The growth curve of a rearing cycle (Ross 308 broiler hens) was monitored starting from the age day seven when chicks were imported to cage-free house (Fig.7). The number of effective readings after PORWI was between 61 to 1613 per day, indicating an adaption period for the first few days. After that period, the number of daily readings shows a clear linear decrease (R² = 0.7393) that is in accordance with an early study observing a decrease in the number of roosting chickens^[16], which can be attributed to a decreased mobility with growth or the fact that heavier chickens less visited the scales^[26]. Nevertheless, even for the final day (166 data points), the number of readings from a single automatic weighing platform was enough for weight estimation for 8300 chickens guaranteeing the commonly required sampling proportion of 2% of the population, within the limit of the recommendation of one weighing platform per 5000 birds. A normal growth curve and daily weight gain can be derived by weight monitoring, indicating its potential use for forecasting the finishing weight and monitoring abnormal growth events.

Fig.7 Number of effective weight readings (a), growth curve (b), and daily weight gain (c) of a broiler flock monitored over 29 days.

Full size|PPT slide

3.2.2 Influence of season, growth rate type and sex

The growing seasons, chicken growth rate types and chicken sexes may have an influence on the accuracy of an automatic weighing platform. The samples were divided into four groups according to seasons (spring, summer, autumn and winter), two groups according to the growth types (medium-growing and slow-growing chickens), and two groups according to the sexes of chickens (hens and cocks), to determine if there were any significant difference in the absolute errors of the weight data in each group. To test the stability of scale performance, the weight data of 101 rearing trials growing at different seasons, growth rate types and sexes were statistically analyzed for effects (Tab.2). The absolute error of the equipment was less than 10 g and the absolute value of relative error was less than 0.62% in all four seasons, with the largest absolute and relative errors in summer. The absolute and relative errors of the equipment for medium-growing chickens were –27 g and –1.45%, respectively, larger than those for slow-growing chickens. The absolute and relative errors of the equipment for cocks were –10 g and –0.62%, and –4 g and –0.23% for hens. Among all the categories of variables, the absolute error was less than 27 g and the relative error less than 1.45%, showing clear agreement between the estimated average weight using the automatic weighing platforms and the standard average weight. As shown by the ANOVA (Tab.2), there was no significant difference for season and sex for the absolute error of the weights. There was a significant difference in the absolute error of weights between medium- and slow-growing chickens at α = 0.001, but the effect size was small and does not indicated any substantive implications. The automatic weighing platform is therefore suitable for varying categories of these variables of growing seasons, growth rate types and sexes.

Tab.2 The standard average weight, estimated weight, absolute error and relative error in different production conditions of rearing season, breed and sex

Variables		Standard average weight mean (SD) (g)	Estimated weight mean (SD) (g)	Absolute error (g)	Relative error (%)
Season (p = 0.658)	Spring	1883 (259)	1879 (260)	–4	–0.21
	Summer	1623 (322)	1613 (305)	–10	–0.62
	Autumn	1561 (134)	1554 (134)	–7	–0.45
	Winter	1751 (191)	1750 (188)	–1	–0.06
Growth rate type (p < 0.001)	Medium-growing	1864 (376)	1837 (370)	–27	–1.45
	Slow-growing	1644 (202)	1642 (203)	–2	–0.12
Sex (p = 0.442)	Hen	1616 (153)	1606 (153)	–10	–0.62
	Cock	1738 (305)	1734 (305)	–4	–0.23

Note: p-values indicate the result of one-way ANOVA for absolute errors of two weighing methods.

3.2.3 Comparison to established systems

Broiler weighing is a useful process in chicken rearing to evaluate the growth performance of a flock. The obtained average weight and uniformity can reflect the daily growth rate, feed-to-meat ratio, health conditions, and finishing day prediction. The standard protocol is to manually sample and weigh a certain ratio of a flock one by one (e.g., 2% of the flock population or 50 birds, whichever is larger). Technology development is steadily replacing time-consuming and inaccurate manual weighing protocols, focusing on machine vision-based approach^[19,27,28], electronic scale algorithm optimization^[29–31], or hybrid systems^[32] for the purpose of better accuracy and stability under varying conditions. Established mature technologies in commercial contexts are mostly electronic scales (Tab.3) manufactured by Fancom, Big Dutchman, Veit and the like (Fig.8). Most products require installation suspended from ceilings, which tends to create an installation difficulty for workers. Those products usually have high precision ranging between 97% and 99.9%, but do not consistently compare to live-bird sales system data for the equality of mean and variance of methods or tested in practical production conditions.

Tab.3 Comparison of existing products

Products	Lumina 47	Fancom 747	Swing 20	Swing 100	DWS-3-ZW	DWS-4-ZW	PS1 and BAT2
Manufacturer	Fancom	Fancom	Big Dutchman	Big Dutchman	Hotraco Agri	Hotraco Agri	DACS and Veit
Poultry types	Broilers, turkeys	Broilers, turkeys	Broilers, ducks, turkeys	Turkeys	Turkeys	Broiler hens and roosters	Broilers, turkeys
Weight range	Up to 150 kg	Up to 150 kg	Up to 20 kg	Up to 100 kg	Up to 100 kg	Up to 12 kg	Up to 50 or 100 kg
Installation	Consisting of a control computer with 2 scales	Consisting of a control computer with a maximum of 8 scales	Consisting of a load cell and a platform made of plastic material; adjustably suspended from ceiling	Consisting of a 1 m × 1 m plastic plate that is attached directly to the load cell by suspension ropes	Consisting of a 1 m × 1 m × 1.5 m cage and suspended in the poultry house	Consisting of hooks, chains, a load cell and a platform	Memory capacity up to 1 year for the number of weighed heads, the average weight, the standard deviation, CV, uniformity, daily increment and sex
Accuracy	97%	97%	Not available	Not available	Not available	Not available	99.9%

Fig.8 Common chicken weighing scales: (a) Lumina 47; (b) turkey weighing scales (a part of Lumina 47 and 747 weighing system); (c) turkey weighing scales in a house for turkey production; (d) Swing 100 in a house for turkey production; (e) Swing 20 in a broiler house; (f) Swing 20 in a house for duck production; (g) DWS-3-ZW; (h) DWS-4-ZW; and (i) BAT 2 GSM.

Full size|PPT slide

4 CONCLUSIONS

This study developed a stress-free, real-time, remote and precise automatic weighing platform for obtaining the average weight and uniformity of broiler flocks based on the biological perching habit of chickens. The performance and stability of this weighing platform were systematically evaluated in commercial chicken production systems. With the adoption of the PORWI protocol, the value of average error between the estimated weight using the automatic weighing platform and the standard average weight assessed using the live-bird sales system was 5.4 g among 101 trials, and the value of relative error averaged 0.54% (an average accuracy 99.5% with a standard deviation 2.3%). The t-test of three monitored trials confirmed the equality of the average weight of birds in broiler flocks estimated using the automatic weighing platform and the live-bird sales system, and the F-test confirmed the similarity of the standard deviations, which indicated accurate estimation of average weight, uniformity and coefficient of variation of chicken flocks using the automatic weighing platforms. This novel weighing systems also had superior stability in terms of growth stage, rearing season, growth rate type (medium- and slow-growing chickens) and sex. The new automatic weighing platform can be used for daily cage-free broiler weight monitoring, potentially benefiting breeding, feeding management, growth monitoring and finishing day prediction, therefore enhancing the sustainability of the poultry industry with this smart poultry farming tool.

References

[1] HuangC. Flora of China. Beijing: Science Press, 1997 (in Chinese)

[2] SwingleW T, Reece P C. The botany of Citrus and its wild relatives. In: Reuther W, Webber H J, Batchelor L D, eds. The Citrus Industry. Berkeley: University of California, 1967, 190–430

[3] ZhouK, Ye M. Fruit Trees in China. Citrus Volume. Beijing: China Forestry Publishing House, 2009 (in Chinese)

[4] Barreca D, Bellocco E, Laganà G, Ginestra G, Bisignano C. Biochemical and antimicrobial activity of phloretin and its glycosilated derivatives present in apple and kumquat. Food Chemistry, 2014, 160 : 292–297

[5] Sadek E S, Makris D P, Kefalas P. Polyphenolic composition and antioxidant characteristics of kumquat (Fortunella margarita) peel fractions. Plant Foods for Human Nutrition, 2009, 64( 4): 297–302

[6] Wang Y W, Zeng W C, Xu P Y, Lan Y J, Zhu R X, Zhong K, Huang Y N, Gao H. Chemical composition and antimicrobial activity of the essential oil of kumquat (Fortunella crassifolia Swingle) peel. International Journal of Molecular Sciences, 2012, 13( 3): 3382–3393

[7] Nagahama K, Eto N, Shimojo T, Kondoh T, Nakahara K, Sakakibara Y, Fukui K, Suiko M. Effect of kumquat (Fortunella crassifolia) pericarp on natural killer cell activity in vitro and in vivo. Bioscience, Biotechnology, and Biochemistry, 2015, 79( 8): 1327–1336

[8] Lou S N, Ho C T. Phenolic compounds and biological activities of small-size citrus: kumquat and calamondin. Journal of Food and Drug Analysis, 2017, 25( 1): 162–175

[9] DengX, Peng S. Citrus Science. Beijing: China Agriculture Press, 2013 (in Chinese)

[10] Fu X Z, Gong X Q, Zhang Y X, Wang Y, Liu J H. Different transcriptional response to Xanthomonas citri subsp. citri between kumquat and sweet orange with contrasting canker tolerance. PLoS One, 2012, 7( 7): e41790

[11] YeY M. Kumquat germplasms in China. China Seed Industry, 1983, 4: 2−5 (in Chinese)

[12] Zhu C, Zheng X, Huang Y, Ye J, Chen P, Zhang C, Zhao F, Xie Z, Zhang S, Wang N, Li H, Wang L, Tang X, Chai L, Xu Q, Deng X. Genome sequencing and CRISPR/Cas9 gene editing of an early flowering mini-citrus (Fortunella hindsii). Plant Biotechnology Journal, 2019, 17( 11): 2199–2210

[13] Tanaka T. Species problem in Citrus: a critical study of wild and cultivated unites of citrus, based upon field studies in their native homes (Revisio Aurantiacearum IX). Japanese Society for the Promotion of Science, 1954, 1–141

[14] DengX, Peng C, ChenZ, DengZ, XuJ, LiJ, LiuY, TangX, Zhong G. Citrus Varieties in China. Beijing: China Agriculture Press, 2008 (in Chinese)

[15] Barkley N A, Roose M L, Krueger R R, Federici C T. Assessing genetic diversity and population structure in a citrus germplasm collection utilizing simple sequence repeat markers (SSRs). Theoretische und Angewandte Genetik, 2006, 112( 8): 1519–1531

[16] Cheng Y, de Vicente M C, Meng H, Guo W, Tao N, Deng X. A set of primers for analyzing chloroplast DNA diversity in Citrus and related genera. Tree Physiology, 2005, 25( 6): 661–672

[17] Pang X M, Hu C G, Deng X X. Phylogenetic relationships within Citrus and its related genera as inferred from AFLP markers. Genetic Resources and Crop Evolution, 2007, 54( 2): 429–436

[18] Yasuda K, Yahata M, Komatsu H, Kunitake H. Phylogeny and classification of Fortunella (Aurantioideae) inferred from DNA polymorphisms. Bulletin of the Faculty of Agriculture Miyazaki University, 2010, 56 : 103–110

[19] Garcia-Lor A, Curk F, Snoussi-Trifa H, Morillon R, Ancillo G, Luro F, Navarro L, Ollitrault P. A nuclear phylogenetic analysis: SNPs, indels and SSRs deliver new insights into the relationships in the ‘true citrus fruit trees’ group (Citrinae, Rutaceae) and the origin of cultivated species. Annals of Botany, 2013, 111( 1): 1–19

[20] Carbonell-Caballero J, Alonso R, Ibañez V, Terol J, Talon M, Dopazo J. A phylogenetic analysis of 34 chloroplast genomes elucidates the relationships between wild and domestic species within the genus Citrus. Molecular Biology and Evolution, 2015, 32( 8): 2015–2035

[21] Curk F, Ollitrault F, Garcia-Lor A, Luro F, Navarro L, Ollitrault P. Phylogenetic origin of limes and lemons revealed by cytoplasmic and nuclear markers. Annals of Botany, 2016, 117( 4): 565–583

[22] Xu Q, Chen L L, Ruan X, Chen D, Zhu A, Chen C, Bertrand D, Jiao W B, Hao B H, Lyon M P, Chen J, Gao S, Xing F, Lan H, Chang J W, Ge X, Lei Y, Hu Q, Miao Y, Wang L, Xiao S, Biswas M K, Zeng W, Guo F, Cao H, Yang X, Xu X W, Cheng Y J, Xu J, Liu J H, Luo O J, Tang Z, Guo W W, Kuang H, Zhang H Y, Roose M L, Nagarajan N, Deng X X, Ruan Y. The draft genome of sweet orange (Citrus sinensis). Nature Genetics, 2013, 45( 1): 59–66

[23] Wu G A, Prochnik S, Jenkins J, Salse J, Hellsten U, Murat F, Perrier X, Ruiz M, Scalabrin S, Terol J, Takita M A, Labadie K, Poulain J, Couloux A, Jabbari K, Cattonaro F, Del Fabbro C, Pinosio S, Zuccolo A, Chapman J, Grimwood J, Tadeo F R, Estornell L H, Muñoz-Sanz J V, Ibanez V, Herrero-Ortega A, Aleza P, Pérez-Pérez J, Ramón D, Brunel D, Luro F, Chen C, Farmerie W G, Desany B, Kodira C, Mohiuddin M, Harkins T, Fredrikson K, Burns P, Lomsadze A, Borodovsky M, Reforgiato G, Freitas-Astúa J, Quetier F, Navarro L, Roose M, Wincker P, Schmutz J, Morgante M, Machado M A, Talon M, Jaillon O, Ollitrault P, Gmitter F, Rokhsar D. Sequencing of diverse mandarin, pummelo and orange genomes reveals complex history of admixture during citrus domestication. Nature Biotechnology, 2014, 32( 7): 656–662

[24] Wang L, He F, Huang Y, He J, Yang S, Zeng J, Deng C, Jiang X, Fang Y, Wen S, Xu R, Yu H, Yang X, Zhong G, Chen C, Yan X, Zhou C, Zhang H, Xie Z, Larkin R M, Deng X, Xu Q. Genome of wild mandarin and domestication history of mandarin. Molecular Plant, 2018, 11( 8): 1024–1037

[25] Yang X M, Li H, Yu H W, Chai L J, Xu Q, Deng X X. Molecular phylogeography and population evolution analysis of Citrus ichangensis (Rutaceae). Tree Genetics & Genomes, 2017, 13( 1): 29

[26] Yang X, Zhou T, Su X, Wang G, Zhang X, Guo Q, Cao F. Structural characterization and comparative analysis of the chloroplast genome of Ginkgo biloba and other gymnosperms. Journal of Forestry Research, 2020, 32( 2): 765–778

[27] Wang Y, Zhou T, Li D, Zhang X, Yu W, Cai J, Wang G, Guo Q, Yang X, Cao F. The genetic diversity and population structure of Sophora alopecuroides (Faboideae) as determined by microsatellite markers developed from transcriptome. PLoS One, 2019, 14( 12): e0226100

[28] Cheng Y J, Guo W W, Deng X X. Molecular characterization of cytoplasmic and nuclear genomes in phenotypically abnormal Valencia orange (Citrus sinensis) + Meiwa kumquat (Fortunella crassifolia) intergeneric somatic hybrids. Plant Cell Reports, 2003, 21( 5): 445–451

[29] Bausher M G, Singh N D, Lee S B, Jansen R K, Daniell H. The complete chloroplast genome sequence of Citrus sinensis (L.) Osbeck var ‘Ridge Pineapple’: organization and phylogenetic relationships to other angiosperms. BMC Plant Biology, 2006, 6( 1): 21

[30] Yang X, Li H, Liang M, Xu Q, Chai L, Deng X. Genetic diversity and phylogenetic relationships of citron (Citrus medica L.) and its relatives in southwest China. Tree Genetics & Genomes, 2015, 11( 6): 129

[31] Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Molecular Biology and Evolution, 2016, 33( 7): 1870–1874

[32] Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Research, 2016, 44( W1): W242–W245

[33] Rozas J, Ferrer-Mata A, Sánchez-DelBarrio J C, Guirao-Rico S, Librado P, Ramos-Onsins S E, Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large datasets. Molecular Biology and Evolution, 2017, 34( 12): 3299–3302

[34] Polzin T, Daneshmand S V. On Steiner trees and minimum spanning trees in hypergraphs. Operations Research Letters, 2003, 31( 1): 12–20

[35] Ruiz C, Paz Breto M, Asíns M J. A quick methodology to identify sexual seedlings in citrus breeding programs using SSR markers. Euphytica, 2000, 112( 1): 89–94

[36] Peakall R, Smouse P E. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research —an update. Bioinformatics, 2012, 28( 19): 2537–2539

[37] Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 2005, 14( 8): 2611–2620

[38] Earl D A, Vonholdt B M. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources, 2012, 4( 2): 359–361

[39] Bolger A M, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 2014, 30( 15): 2114–2120

[40] Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009, 25( 14): 1754–1760

[41] Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. The sequence alignment/map format and SAMtools. Bioinformatics, 2009, 25( 16): 2078–2079

[42] Cingolani P, Platts A, Wang L, Coon M, Nguyen T, Wang L, Land S J, Lu X, Ruden D M. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly, 2012, 6( 2): 80–92

[43] Yang J, Lee S H, Goddard M E, Visscher P M. GCTA: a tool for genome-wide complex trait analysis. American Journal of Human Genetics, 2011, 88( 1): 76–82

[44] Alexander D H, Lange K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics, 2011, 12( 1): 246

[45] Danecek P, Auton A, Abecasis G, Albers C A, Banks E, DePristo M A, Handsaker R E, Lunter G, Marth G T, Sherry S T, McVean G, Durbin R. The variant call format and VCFtools. Bioinformatics, 2011, 27( 15): 2156–2158

[46] Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira M A R, Bender D, Maller J, Sklar P, de Bakker P I W, Daly M J, Sham P C. PLINK: a tool set for whole-genome association and population-based linkage analyses. American Journal of Human Genetics, 2007, 81( 3): 559–575

[47] Vilella A J, Blanco-Garcia A, Hutter S, Rozas J. VariScan: analysis of evolutionary patterns from large-scale DNA sequence polymorphism data. Bioinformatics, 2005, 21( 11): 2791–2793

[48] Li H, Durbin R. Inference of human population history from individual whole-genome sequences. Nature, 2011, 475( 7357): 493–496

[49] Meirmans P G, Hedrick P W. Assessing population structure: F(ST) and related measures. Molecular Ecology Resources, 2011, 11( 1): 5–18

[50] Weir B S, Cockerham C C. Estimating F-statistics for the analysis of population structure. Evolution, 1984, 38( 6): 1358–1370

[51] Wang X, Xu Y, Zhang S, Cao L, Huang Y, Cheng J, Wu G, Tian S, Chen C, Liu Y, Yu H, Yang X, Lan H, Wang N, Wang L, Xu J, Jiang X, Xie Z, Tan M, Larkin R M, Chen L L, Ma B G, Ruan Y, Deng X, Xu Q. Genomic analyses of primitive, wild and cultivated citrus provide insights into asexual reproduction. Nature Genetics, 2017, 49( 5): 765–772

[52] LinD, WuF. Distribution and variety of kumquat in China. South China Fruits, 1987, 1(1): 3−5 (in Chinese)

[53] Garcia-Lor A, Luro F, Ollitrault P, Navarro L. Genetic diversity and population structure analysis of mandarin germplasm by nuclear, chloroplastic and mitochondrial markers. Tree Genetics & Genomes, 2015, 11( 6): 123

[54] ChenP. The investigation and genetic diversity evaluation of wild Hong Kong kumquat (Fotunella hindsii swingle) in China. Dissertation for the Master’s Degree. Wuhan: Huazhong Agricultural University, 2011 (in Chinese)

[55] Wu G A, Terol J, Ibanez V, López-García A, Pérez-Román E, Borredá C, Domingo C, Tadeo F R, Carbonell-Caballero J, Alonso R, Curk F, Du D, Ollitrault P, Roose M L, Dopazo J, Gmitter F G, Rokhsar D S, Talon M. Genomics of the origin and evolution of Citrus. Genomics of the origin and evolution of Citrus, 2018, 554( 7692): 311–316

[56] Atahan P, Itzstein-Davey F, Taylor D, Dodson J, Qin J, Zheng H, Brooks A. Holocene-aged sedimentary records of environmental changes and early agriculture in the lower Yangtze, China. Quaternary Science Reviews, 2008, 27( 5−6): 556–570

[57] ZhangJ Z, Chen C F, YangY Z. Origins and early development of agriculture in China. Journal of National Museum of China, 2014, 1: 6−16 (in Chinese)

Supplementary materials

The online version of this article at https://doi.org/10.15302/J-FASE-2021436 contains supplementary materials (Figs. S1–S5; Tables S1–S11; supplementary data set). Resequencing data are available in the CNCB-NGDC database (BioProject PRJCA005010) and NCBI database (BioProject PRJNA736109).

Acknowledgements

This work was funded by the National Key Research and Development Program of China (2018YFD1000106), the National Natural Science Foundation of China (31630065), and Special Project for External Science and Technology Cooperation of Science and Technology Department of Yunnan Province (202003AD150014). We sincerely thank Prof. Jian Li (Department of Agriculture of Fujian Province, China), Prof. Jian’guo Xu (Citrus Research Institute of Zhejiang Province, China), Prof. Shiping Zhu (Citrus Research Institute of CAAS, China), Prof. Wu Wen (Guangdong Academy of Agricultural Sciences, China), Zhenhua Tan (Agricultural Bureau of Yizhang, China), Prof. Guohua Li (Meizhou Academy of Agricultural Sciences, China), and Prof. Min Zhang (Zhejiang A&F University, China) for their assistance in sample collection.

Compliance with ethics guidelines

Chenqiao Zhu, Peng Chen, Junli Ye, Hang Li, Yue Huang, Xiaoming Yang, Chuanwu Chen, Chenglei Zhang, Yuantao Xu, Xiaoli Wang, Xiang Yan, Guangzhou Deng, Xiaolin Jiang, Nan Wang, Hongxing Wang, Quan Sun, Yun Liu, Di Feng, Min Yu, Xietian Song, Zongzhou Xie, Yunliu Zeng, Lijun Chai, Qiang Xu, Chongling Deng, Yunjiang Cheng, and Xiuxin Deng declare that they have no conflicts of interest or financial conflicts to disclose. This article does not contain any study with human or animal subjects performed by any of the authors.

RIGHTS & PERMISSIONS

The Author(s) 2022. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)