SNP Identification
After the Illumina HiseqTM sequencing data (Raw Data) is decommissioned, the decommission data is quality controlled, and the low-quality data is filtered to obtain high-quality data (Clean Data). A total of 823.64G of clean data were obtained based on RAD sequencing of 198 samples. A total of 2,852,333 SNP remains for the calculation of PIC, SHI, HE, HO, Tajimas' D, Pi, FST, and Nm. To avoid calculation errors caused by linkage of SNP, the markers used when calculating the inbreeding coefficient (FIS) using PLINK (PLINK-indeep-pairwise 100 10 0.1) should be uncorrelated (i.e., all in linkage equilibrium). Therefore, the obtained 2,852,333 SNPs were subjected to LD filtering, and the resulting 177,447 SNP were used for FIS calculation.
Genetic diversity analysis of 9 populations
Due to the small sample size of C. duclouxii and C. ovata, there may be deviation in their genetic diversity index, so they are not included in the comparison of genetic diversity.The 169 samples of Huangxin were divided into groups by sampling route, and a total of 9 groups were divided (Fig. 2).
The observed heterozygosity (HO) and expected heterozygosity (HE) are important indicators for population genetic diversity analyzing. In the Huangxin populations, the highest HO value was pop6 (0.3545), followed by pop5 (0.3524), with a mean value of 0.3408. The HE values ranged from pop6 (0.2530) to pop7 (0.2823), with an average value of 0.2736. Under the standard neutral evolutionary model, the theoretical value of Tajimas' D should be 0, and when greater than 0, the group observation heterozygosity is greater than the expected heterozygosity, and the group is in equilibrium (Shameek B& Joshua MA, 2006, Hideki I& Wolfgang S,2000).
The Huangxin as a whole was in equilibrium with a mean Tajima’D value of 0.722 and a minimum value of 0.534 (pop9). The polymorphic information content (PIC) value displays the polymorphism of each locus (Botstein et al. 1980). The highest PIC value continues to be pop8 (0.2911), followed by pop7 (0.2899). The highest inbreeding coefficient (FIS) was pop3 (-0.09511) and the lowest was pop5 (-0.15566). In summary, the highest level of genetic diversity in the Huangxin population was pop7. The initial inference may be that the genetic diversity is a little higher because the topography of the pop7 range is more complex than that of other populations.
Table 1
Genetic diversity indicators for 9 populations of Huangxin
Pop
|
N
|
HO
|
HE
|
PIC
|
Tajima’D
|
FIS
|
1
|
48
|
0.3360
|
0.2731
|
0.2761
|
1.042
|
-0.14734
|
2
|
20
|
0.3508
|
0.2790
|
0.2865
|
0.802
|
-0.15367
|
3
|
18
|
0.3180
|
0.2715
|
0.2796
|
0.658
|
-0.09511
|
4
|
19
|
0.3352
|
0.2775
|
0.2853
|
0.737
|
-0.12229
|
5
|
16
|
0.3524
|
0.2772
|
0.2865
|
0.746
|
-0.15566
|
6
|
5
|
0.3545
|
0.2530
|
0.2820
|
0.534
|
-0.14343
|
7
|
20
|
0.3428
|
0.2823
|
0.2899
|
0.804
|
-0.13006
|
8
|
14
|
0.3482
|
0.2804
|
0.2911
|
0.656
|
-0.13703
|
9
|
9
|
0.3292
|
0.2683
|
0.2849
|
0.515
|
-0.1004
|
Mean
|
|
0.3408
|
0.2736
|
0.2847
|
0.722
|
-0.1317
|
(Note: N, samples number, HO, Observed heterozygosity, HE, Expected heterozygosity, PIC, Polymorphic information content, Tajimas’D, Tajimas’D Median test, FIS, Inbreed coefficient.) |
The coefficient of genetic differentiation (FST) > 0.25, represents very significant genetic differentiation, and 0.25 > FST>0.15 represents significant genetic differentiation, below 0.15 is very little differentiation (Wright S, 1965). As we can see from Table 2, the FST for the entire Huangxin populations is very low(The largest FST values, pop1 and pop6, were only 0.05973) and even appears to have a negative value (in the calculation of FST values, a negative value indicates almost no differentiation and is usually regarded as zero).
This suggests that our approach of dividing populations by sampling route to calculate their genetic diversity is not reasonable, as it is possible that two populations in close proximity are the same asexual line that led to this result (tentatively due to human influence). Therefore, we performed genetic diversity analysis of Huangxin in populations delineated exclusively by SNP data.
Table 2
Pairwise comparison of genetic differentiation (FST) among 9 populations.
|
pop2
|
pop3
|
pop4
|
pop5
|
pop6
|
pop7
|
pop8
|
pop9
|
pop1
|
0.00043
|
0.01636
|
0.02258
|
0.03910
|
0.05973
|
0.02948
|
0.00981
|
0.03606
|
pop2
|
|
0.00960
|
0.01237
|
0.02367
|
0.04008
|
0.01834
|
0.00277
|
0.02377
|
pop3
|
|
|
0.00998
|
0.01281
|
0.03333
|
0.00884
|
0.00149
|
0.01137
|
pop4
|
|
|
|
0.00301
|
0.01101
|
0.00404
|
0.00057
|
0.00950
|
pop5
|
|
|
|
|
0.00418
|
0.00590
|
0.00739
|
0.01062
|
pop6
|
|
|
|
|
|
0.01909
|
0.02277
|
0.02127
|
pop7
|
|
|
|
|
|
|
0.00553
|
-0.00365
|
pop8
|
|
|
|
|
|
|
|
0.01057
|
Phylogenetic analysis
Based on SNP, a neighbor-joining method of RAxML software (Stamatakis, 2006) constructed a population evolutionary tree of 198 samples (Fig. 3).
The five individuals of C. ovata, CO001-CO005 (outgroup), are genetically the furthest from other individuals and are located at the bottom of the evolutionary tree, indicating that these five individuals are less closely related to the Huangxin population. In contrast, the C. duclouxii population (the green part) is closely related to the Huangxin population and is clustered together on a branch (Fig. 3A). However, it is noteworthy that the branches of Huangxin and C. duclouxii are completely separate, which indicates that they are two separate species in agreement with what is referred to as the "type" in the application.
To clarify the genetic relationship among the Huangxin population, a non-branched long version of the phylogenetic tree was constructed based on Fig. 3A. One branch of 24 individuals in C. duclouxii were separated from other samples, DQ264 and DQ262 were first separated from other C. duclouxii individuals, and the remaining C. duclouxii individuals were divided into two branches as a whole. The individual HX246 in the Huangxin group was first separated from other individuals, then HX125 and HX113 were separated from other individuals, and the remaining individuals were roughly divided into 5 subgroups (Fig. 3B).
In the group of Huangxin, individual HX246 separated from the other individuals first, followed by HX125 and HX113 separating from the others. The remaining individuals were roughly divided into five subgroups: Subgroup 1 (63 individuals), Subgroup 2 (10 individuals), Subgroup 3 (10 individuals), Subgroup 4 (38 individuals), and Subgroup 5 (42 individuals).
Population genetic structure analysis
In this study, 198 samples (including Huangxin, C. duclouxii, C. ovata) were genetically diverse using RAD-seq, and 7 individuals with insignificant genetic structure were excluded. Filtered SNP data were analyzed using ADMIXTURE software for population genetic structure analysis. Bayesian cluster analysis of samples (Fig. 4A), when K = 7, the ratio of Cross-Validation (CV) errors reaches the lowest value, and the second lowest value is K = 10. Therefore, it can be inferred that the optimal group size is either 7 or 10. When K = 7, 191 individuals were divided into seven clusters, and their genetic information was independent of each other. The five subgroups of Huangxin correspond to the five clusters in K = 7, and in combination with the sampling route distribution map, Subgroup 1 distributed within pop1, pop2, pop8, Subgroup 2 location within pop7 and pop9, Subgroup 3 distributed in the pop4, 5, 7, and pop8 range, the distribution of Subgroup 4 is scattered, except for the pop6 range where it occurs in other populations, and Subgroup 5 scattered throughout the range of other populations except pop1,2.
The remaining five individuals were between subpopulations 1 and 2, and from Fig. 4B, their genetic information was the same as six ancestors in the seven clusters (without outgroup). The possible reasons that may lead to the occurrence of this phenomenon: populations may have intermixed through historical migrations or other events, resulting in individuals who carry the genetic makeup of multiple ancestral populations. These factors can complicate population structure analysis because they can make it difficult to accurately assign individuals to specific ancestral populations.
When K = 10, the population structure remains unchanged, but there are internal changes within the population. The genetic information of Subgroup 1, 4, and 5 have become mixed: according to Fig. 4C, all three subgroups contain three ancestors, with two of them being shared but in different proportions. The two shared ancestors are relatively more prominent in subgroup 1, while Subgroup 4 and 5 are more similar in their proportions. The ancestral origins of the remaining subgroups are clear and independent. It is tentatively speculated that the existence of this phenomenon may be due to the proximity of the three subgroups' geographic locations or other factors that lead to gene flow between them.
Principal component analysis
PCA plots can visualize genetic structures in a multivariate space, the results showed that the contribution rate of the first three principal components to genetic variation is 39.66%, 19.64%, and 11.07%, respectively (Fig. 5). PCA shows that 191 samples are divided into 7 groups.
According to Fig. 5A, it is evident that there is a considerable distance between the C. duclouxii and C. ovata outgroup populations and the Huangxin population. Within the Huangxin population, the Subgroup 2 and 3 are very close in distance, and the individuals that are located between these subgroups but are not included in either are genetically mixed as shown in Fig. 4. An interesting phenomenon is that individual HX246 is located between the Huangxin and C. duclouxii populations (Fig. 5B, C), and in Fig. 4, it can be seen that most of its genetic information comes from C. duclouxii. Therefore, it can be inferred that HX246 may be an early variant of C. duclouxii. Based on Fig.s 5A, and D, it is clear that Subgroup 1, 4, and 5 within the Huangxin population are distinct and independent. This result is consistent with phylogenetic trees (Fig. 3) and population structure (Fig. 4).
Genetic diversity analysis
The results of the genetic diversity analysis of the 7 groups are shown in Table 3. Among the Huangxin population, Tajimas' D values range from 1.860 (Subgroup 2) to 3.343 (Subgroup 1) which means they're both in equilibrium. The average expected heterozygosity (HE, 0.2935) was less than the observed heterozygosity (Ho, 0.44014). Subgroup 5 had the highest HO value (0.4495), followed by Subgroup 1 (0.4476), and Subgroup 2 had the lowest HO value, at 0.4271. The HE value of Subgroup 5 is still the highest (0.3011), followed by Subgroup 4 (0.3002).
Shannon’s information index (SHI) reflects the heterozygous nature of the population, the mean of the SHI is 0.429, with a range of 0.412 (Subgroup 2) ~ 0.440 (Subgroup 5). The overall polymorphism of Huangxin is low, with a PIC value < 0.25, which ranged from 0.221 (Subgroup 2) to 0.236 (Subgroup 5), with an average of 0.2298. Pi is a concept used in molecular genetics to measure the degree of polymorphism within a population (Liu et al.,2020). The Pi is 0.001153 (Subgroup 2) to 0.001178 (Subgroup 5). The inbreeding coefficient (FIS) reflects the degree of differentiation between individuals, the mean of FIS is -0.255 and ranges from − 0.209 (Subgroup 2) to -0.285 (Subgroup 5), the FIS values of each subpopulation are negative, indicating that there is no inbreeding in the Huangxin groups. Based on all the indicators mentioned above, Subgroup 5 has the highest genetic diversity, while Subgroup 2 has the lowest.
Table 3
Genetic diversity indicators for seven groups.
Pop
|
N
|
HO
|
HE
|
SHI
|
PIC
|
Pi
|
Tajima’D
|
FIS
|
Subgroup1
|
63
|
0.4476
|
0.2999
|
0.439
|
0.235
|
0.001169
|
3.343
|
-0.273
|
Subgroup2
|
10
|
0.4271
|
0.2828
|
0.412
|
0.221
|
0.001153
|
1.860
|
-0.209
|
Subgroup3
|
10
|
0.4342
|
0.2837
|
0.413
|
0.222
|
0.001156
|
1.886
|
-0.241
|
Subgroup4
|
38
|
0.4423
|
0.3002
|
0.439
|
0.235
|
0.001176
|
2.950
|
-0.265
|
Subgroup5
|
41
|
0.4495
|
0.3011
|
0.440
|
0.236
|
0.001178
|
3.084
|
-0.285
|
C. duclouxii
|
24
|
0.1056
|
0.1064
|
0.161
|
0.086
|
0.000417
|
1.345
|
0.695
|
Outgroup
|
5
|
0.0479
|
0.0409
|
0.061
|
0.033
|
0.000103
|
0.180
|
0.860
|
(Note: Pi, Nucleotide polymorphisms, Tajimas’D, Tajimas’D Median test.) |
Genetic diversity analysis
The mean FST between groups was 0.2966, the gene flow (Nm) was 0.70945, and the degree of genetic differentiation between groups was high (FST>0.25). The five subgroups of Huangxin had the highest degree of genetic differentiation with the C. duclouxii. group (FST=0.399 ~ 0.420). The degree of genetic differentiation between the five subgroups was high, with FST values ranging from 0.181 to 0.222(mean = 0.1983). Among them, Subgroup 4 and Subgroup 1 have the highest degree of genetic differentiation, followed by Subgroup 5 and Subgroup 1 (FST =0.212). The degree of genetic differentiation between Subgroup 2 and Subgroup 3 is the lowest.
Wright (1943) determined that Nm values greater than 1 were sufficient to overcome the effects of genetic drift. There was less gene exchange between different groups, ranging from Nm = 0.125 to 1.134 (mean 0.7094), with Subgroup 2 and Subgroup 3 having the closest genetic exchange of 1.134, followed by Subgrou p4 and Subgroup 5 (Nm = 1.124).
Table 4
Genetic differentiation of seven groups (lower diagonal: pairwise Nm, upper diagonal: pairwise FST).
|
C.duclouxii
|
Outgroup
|
Subgroup1
|
Subgroup2
|
Subgroup3
|
Subgroup4
|
Subgroup5
|
C. duclouxii
|
0
|
0.666
|
0.420
|
0.399
|
0.400
|
0.405
|
0.412
|
Outgroup
|
0.125
|
0
|
0.313
|
0.303
|
0.305
|
0.310
|
0.311
|
Subgroup1
|
0.346
|
0.549
|
0
|
0.200
|
0.205
|
0.222
|
0.212
|
Subgroup2
|
0.377
|
0.574
|
0.998
|
0
|
0.181
|
0.201
|
0.198
|
Subgroup3
|
0.375
|
0.569
|
0.969
|
1.134
|
0
|
0.192
|
0.190
|
Subgroup4
|
0.367
|
0.555
|
0.877
|
0.991
|
1.053
|
0
|
0.182
|
Subgroup5
|
0.356
|
0.555
|
0.929
|
1.015
|
1.060
|
1.124
|
0
|
(Note: FST, Genetic differentiation coefficients, Nm, gene flow) |