Characterization of genomic variations in CRI12 pedigree

We employed a combination of Nanopore long reads, Illumina short reads, and Hi-C data to assemble the CRI12 genome, capturing a total of 2.24 Gb sequences (contig N50 = 11.2 Mb) (Table 1; Additional file 1: Table S1-S3). The assembled genome exhibited high syntenic ratios (ranging from 93.5 to 98.5%) with other upland cotton genomes (Additional file 1: Table S4) [7, 8, 17]. Moreover, we evaluated the completeness of the CRI12 genome using Benchmarking universal single-copy orthologs (BUSCO) and Core Eukaryotic Genes Mapping Approach (CEGMA), achieving completeness of 96.1 and 95.56%, respectively [18, 19] (Table 1) [7, 8, 17]. We annotated a total of 72,262 protein-coding genes after conducting continuity, consistency, and completeness (3C) evaluations of the assembled CRI12 genome (Table 1). Our results suggest that the CRI12 genome assembly is highly competent and can be used as a reference genome for downstream analyses. To investigate genetic pattern of segments from CRI12 in cultivar-driven strategy, we collected 20 CRI12 pedigree members including CRI12’s parents ( Xingtai6871 and Uganda4; defined as Parent#1 and Parent#2 in this study, respectively) and 18 progenies for whole-genome sequencing using ONT long reads to a median of 20 × coverage (Fig. 1a; Additional file 1: Table S5-S6). We identified a total of 13,138 high-quality PAVs (7172 deletions and 5966 insertions) after aligning the long reads of the pedigree members to the CRI12 genome (Fig. 2a; Additional file 1: Table S7). Upland cotton is an allotetraploid, and its genome (AADD) consists of sub-genome A and sub-genome D. In CRI12 pedigree, there is no preference on insertion or deletions between A and D sub-genomes (Fig. 2a; P = 0.55, χ2 test). The pan-core curve of the 13,138 PAVs eventually reached a plateau, confirming the representativeness of the 20 pedigree members we selected (Fig. 2b).

Table 1 BUSCOs and CEGMAs indicated ratios of benchmarking genes in 2 pipelines were correctly assembled in CRI12 genome. The ratios could represent completeness of assembled genome. Repeat ratio indicated ratio of repeat sequence in CRI12 genomeFig. 2

Characterization of genomic variations in CRI12 pedigree. a Statistic about number of PAVs identified in CRI-12 pedigree (insertions and deletions) within A and D sub-genomes. b Saturation curve of core and pan SVs in the CRI12 pedigree with cultivars added. c Statistic about PAV sequences influenced by transposons. Barplot was the distribution of transposon ratio in transposon-influenced segments. Pie plot represents PAVs of high transposon ratio (≥ 90%) and common transposon ratio (< 90%). d Genomic annotation of transposon-influenced segments. e Statistic about number of transposons of various categories in pedigree

Previous studies have demonstrated that transposon activity drives genomic variation, thereby promoting genetic innovation [10, 20, 21]. In this study, we annotated sequences of 13,138 PAVs, 6076 among them were generated through transposition events (Additional file 1: Table S8). A considerable proportion (54.3%) of the transposon-influenced PAVs displayed a high transposon ratio (> 90%), showing transposons also contribute to genomic variations in pedigree (Fig. 2c). Furthermore, we found that exons and introns exhibited fewer transposon activities (299 in exons and 449 in introns) as compared to interval and promoter, i.e., 3931 in interval and 833 in promoter, indicating genes tend to repel transposon activity (Fig. 2d). Analysis of transposon categories revealed a total of 26 distinct categories found within the 6076 transposon-influenced PAVs. Of these categories, LTR/Gypsy, the most abundant transposon, is highly related with PAV length in exon, promoter, and interval regions (Pearson correlation is 0.81, 0.72, and 0.73, respectively; 5627 of 6076 PAVs) and relatively lowly related with PAV length in intron (Pearson correlation is 0.64; 449 of 6076 PAVs) (Fig. 2e; Additional file 2: Fig. S1a), indicating its important role in SVs proved in G. rotundifolium [20]. Surprisingly, we observed a strong correlation (r2 = 0.85) between PAV length in exons and the number of Simple_repeats, which were present at a notably higher proportion in exons (19%) when compared to interval regions (13%) (Additional file 2: Fig. S1b). These findings suggest that Simple_repeat is highly correlated with exon, and their relationship should be further investigated in future studies. Overall, we constructed a high-quality PAV map which could present pedigree’s variation, providing a platform for functional interpretation on CRI12 pedigree’s large sequence variations.

CRI12 pedigree contains valuable genomic variations

SVs contribute significantly to phenotype during crop breeding. In cotton, the use of an SV-based pan-genome has accurately reflected the breeding history in China. To identify valuable SVs during the formation of the CRI12 pedigree, we constructed a Graphical Pedigree Genome (GPG) by integrating 13,138 PAVs into the CRI12 reference genome.

Using the GPG of CRI12, three previous cohorts (Ma_2018, Wang_2017, and He_2021) were genotyped, and trait-related PAVs were detected through PAV-GWAS for 7 agronomic traits (Fig. 3a; Additional file 1: Table S9-S11; Additional file 2: Fig. S2-S5). After applying strict filters to the detected trait-related PAVs (“Methods”), a functional PAV pool consisting of 623 variations was created (Additional file 1: Table S12). Within this PAV pool, we annotated 131 genes whose exons, introns, and promoters were directly influenced by PAVs (Additional file 1: Table S13). Of these 131 genes, 48 were related to resistance against Verticillium wilt (VW), and only 25 had transcriptional activity during Verticillium dahliae infections (Additional file 1: Table S14). We focused on GhNAC083, a gene previously reported as a repressor for xylem vessel formation in Arabidopsis thaliana, due to its opposite transcription patterns between VW-susceptible and VW-resistant accessions (Additional file 1: Table S15) [22]. As NAC083 is involved in both biotic and abiotic stress responses, we inferred that it may also play an important role in cotton’s VW resistance [23]. In addition, GhNAC083-silenced plants are more susceptible to Verticillium dahliae infection compared to the wild type, indicating that GhNAC083 could enhance VW resistance and be utilized in cotton breeding (Fig. 3b).

Fig. 3

PAV-GWAS for agronomic traits. a PAV-GWAS for Verticillium wilt resistance, and the threshold of P value was set as 0.01. b VIGS assay on GhNAC083. c Genomic annotation for trait-related PAVs. d Concept illustration for consistent and conflict SV-Gene expression-Phenotype triplets. e Display of Parent#1_DEL_2442 and Parent#2_DEL_5283. The alignment of long reads from Parent#1 and Parent#2 were displayed in IGV toolkit. f Fiber length divergence of 2 genotypes compared to CRI12-like genotypes in Wang_2017 cohort (n = 169). Boxplot center, median; bounds of box, lower quartile (Q1) and upper quartile (Q3); minima, Q1 − 1.5 (Q1–Q3); maxima, Q1 + 1.5 (Q1–Q3). The dots represent for original values. P values were calculated by t-test and “**” represent for P ≤ 0.01. g Gene expression of GhKCR1A and GhKCR1D in wang_2017 cohort (n = 169). X-axis was the rank of fiber length (from short to long; divided into quartiles). Y-axes are TPM (transcripts per million mapped reads) of GhKCR1A and GhKCR1D, respectively. The description of boxes is the same as that in f, and “***” represents for P ≤ 0.001 in t-test. h Illustration for prediction of agronomic trait based on functional variations by support vector regression. i Pearson correlation between predicted phenotype and real phenotype

In the functional PAV pool, we discovered that 495 PAVs were located in the interval region, which was typically ignored in previous GWAS analyses, limiting the further utilization of elite genotypes (Fig. 3c) [14,15,16, 24, 25]. However, by utilizing an expression array with a population scale (Wang_2017 and He_2021 cohorts, n = 314 in total), we were able to link these functional non-coding PAVs to genes through eQTL analysis (Additional file 1: Table S10-S11) [15, 26]. Based on the effect of functional PAVs on agronomic traits and gene expression (positively or negatively associated), we classified all 2306 identified SV-Gene-Phenotype (SGP) triplets into consistent and conflict categories (Fig. 3d; Additional file 1: Table S16). Genes within consistent SGP triplets had higher expression with the gradual improvement of fiber quality, while those within conflict triplets showed a decline in gene expression as fiber quality improved (Additional file 2: Fig S6-S7). In the Wang_2017 cohort, we noted that enrichment of genes related to fatty acid metabolism (ko01212) among both consistent and conflict SGP triplets (Additional file 2: Fig. S8). The fatty acid metabolism pathway contained six genes including three genes each from the consistent and conflict SGP triplets (Additional file 1: Table S17). In the two consistent SGP triplets (Parent#1_DEL_2442-GhKCR1A-FL and Parent#2_DEL_5283-GhKCR1D-FL), the transcription abundance of GhKCR1 increased from 1 day post anthesis (DPA) to 8 DPA, indicating its essential role in fiber elongation. Two different elite variations from Parent#1 and Parent#2 linked to homologs of GhKCR1 in the A and D sub-genomes, respectively (Fig. 3e, f). The number of cultivars that contained these two elite variations increased with the improvement of fiber length (Additional file 2: Fig. S9). We noticed that GhKCR1A in long fiber cultivars (rank of fiber length at 50–75% and 75–100%) exhibited higher expression levels than short fiber cultivars (rank of fiber length at 0–25%) with both 6.5% increase (P values were 0.02 and 0.01, respectively) (Fig. 3g). GhKCR1D had similar transcription patterns with GhKCR1A, and long fiber cultivars (rank of fiber length at 50–75% and 75–100%) exhibited 6.5 and 13% increase in expression level (P values were 0.035 and 0.0025, respectively) compared to short fiber cultivars (rank of fiber length at 0–25%), respectively. This result shows that higher expression level of GhKCR1s is related with longer fiber length in a population scale. Interestingly, the expression of GhKCR1 in the A sub-genome was higher than that in the D sub-genome, indicating asymmetric selection on GhKCR1s from different sub-genomes (Fig. 3g). GhKCR1 encodes very-long-chain 3-oxoacyl-CoA reductase, responsible for the first reduction step in synthesis of very-long-chain fatty acids (VLCFAs) which is important for fiber elongation [27,28,29]. Based on gene annotation and the performance of the two SGP triplets in the population, we concluded that GhKCR1s identified from the above two SGP triplets could be used as candidate genes for improving fiber length. To examine the effects of the trait-related PAVs identified in this work, we constructed a set of support vector regression (SVR) models to predict agronomic traits based on genotype data (Fig. 3h). The Pearson correlation coefficient between predicted values and true values of multiple traits ranged from 0.58 to 0.85 (P values range from 1.09e − 52 to 1.38e − 129). These results show that except for fiber uniformity and fiber length, the other agronomic traits could be predicted reliably by models constructed based on functional PAVs (Fig. 3i).

Allelic variations were detected based on an accurate PAV map (482 pairs), and their effects were evaluated using the results of PAV-GWAS (Additional file 1: Table S18). We focused on a pair of non-coding allelic deletions, Xinluzhong7_DEL_3213-Parent#2 _DEL_2661 in the Wang_2017 cohort. On this locus, Parent#2-like genotypes had significantly higher fiber uniformity than Xinluzhng7-like genotypes (P = 0.002) (Additional file 2: Fig. S10a and S10b). However, compared to CRI12-like genotypes, both Parent#2-like and Xinluzhong7-like genotypes did not reach the significant divergence threshold in fiber uniformity (P = 0.07 and P = 0.3 for Parent#2-like and Xinluzhong7-like genotypes, respectively), resulting in their dismissal according to the threshold of PAV-GWAS described in the “Methods” (Additional file 2: Fig. S10a and S10b). Furthermore, we discovered 39 and 7 genes linked to Parent#2_DEL_2661 and Xinluzhong7_DEL_3213, respectively, from eQTL results, indicating that allelic variation might regulate different gene networks and lead to phenotype divergence (Additional file 1: Table S19; Additional file 2: Fig S10c).

Low hereditary stability for functional segments in backbone cultivar

The functional PAV pool of the CRI12 pedigree revealed that valuable variations were generated during its formation. Taking the phenotype effects into consideration, we defined a deleterious/favorable segment in CRI12 genome if the variation in alternative genome happened on corresponding locus could improve/reduce agronomic traits (Additional file 2: Fig. S11). Each functional variation in the pedigree represents either a favorable or deleterious segment in CRI12, and we classified all CRI12 segments into four categories such as Parent#1-inherited, Parent#2-inherited, CRI12-specific, and Biparental, based on the segment source (Fig. 4a). The source of Biparental segments in CRI12 is unclear due to the lack of genomic information about the parents of Parent#1 and Parent#2. Given that the Biparental segments were present in both Parent#1 and Parent#2, we concluded that they were likely inherited by CRI12’s parents from progenitor cultivars. The present/absent information regarding CRI12 segments among pedigree members enabled us to characterize the genetic patterns of the backbone cultivar in cultivar-driven strategy.

Fig. 4

Genetic patterns of functional segments in CRI12 pedigree. a Conceptual illustration for segments of 4 categories. Parent#1-inhrited segments in CRI12 were inherited from Parent#1, presenting a structural variation in Parent#2. Parent#2-inherited segments in CRI12 were inherited from Parent#2, presenting a structural variation in Parent#1. CRI12-specific segments were those created by CRI12, presenting a structural variation in both Parent#1 and Parent#2. Biparent segments were those both present in CRI12’s 2 parents, only presenting a structural variation in progenies. b Relationship between Fusarium wilt resistance and functional segments of 4 categories in cultivars (Wang_2017 cohort, n = 208). X-axis represents for rank of disease index (from low to high; divided into quartiles) and Y-axis represents for segment number in cultivars from each quartile. Orange lines were CRI12-specific segments. Red lines were Parent#1-inherited segments. Blue lines were Parent#2-inherited segments, and turquoise lines were Biparent segments. The dots represent for the median values of segments in cultivars from each quartile. c Top 4 pieplots represent for favorable/deleterious ratio of 4 segment categories. The bottom boxplot was the hereditary stability evaluation of Fusarium wilt-resistant segments. X-axis represents 4 segment categories and y-axis represents the segment number retained in 18 progenies. Boxplot center, median; bounds of box, lower quartile (Q1) and upper quartile (Q3); minima, Q1 − 1.5 (Q1–Q3); maxima, Q1 + 1.5 (Q1–Q3). The dots represent for original values. d P value of t-tests for retained segments of 4 categories. e Distribution of fiber-quality-related segments in pedigree. The left parallel coordinate plots displayed segment number of 4 categories. SV axis represents for variations. Class axis represents for 4 segment categories which have been annotated in figure legends. Progenies axis represents for 18 progenies. Right heatmap were detailed presence of each functional segments in 18 progenies. X-axis represents for segments and y-axis represents for 18 progenies which are consistent with the progenies axis in left parallel coordinate plots. The top parallel coordinate-heatmap illustrates distribution of favorable fiber quality-related segments, while the bottom plot set illustrates deleterious fiber quality-related segments

According to the results of PAV-GWAS on Fusarium wilt resistance (FR), 70 FR-related segments were identified. Among these 70 segments, 4, 7, 14, and 45 of them were assigned as Parent#1-inherited, Parent#2-inherited, CRI12-specific, and Biparental segments, respectively (Additional file 1: Table S20). We evaluated the effects of 70 FR-related segments belong to four categories in the Wang_2017 cohort (Fig. 4b). The cultivars with more CRI12-specific segments exhibited higher resistance to Fusarium wilt (lower disease index), indicating the positive effects of CRI12 in pathogen resistance. The ratio of favorable to deleterious segments in the CRI12-specific category reached 92.9% (13 favorable and 1 deleterious), while the ratios in the other three categories ranged from 42.9 to 53.3% (Fig. 4c). These results indicated that, compared to its parents, CRI12 generated more favorable FR-related segments, contributing to its pathogen resistance. However, we also found that these favorable CRI12-specific segments were less frequent in 18 progenies compared to segments inherited from its two parents and ancestral cultivars (P values ranged from 0.01 to 7e − 10), indicating low hereditary stability of CRI12-specific segments (Fig. 4c, d). The nature of the cultivar-driven strategy is to pass elite genomic segments from the backbone cultivar to its progenies. Our results showed the high ratio of favorable CRI12-specific segments, proving that CRI12 is competent to be a backbone cultivar, while their low hereditary stability suggests that they may be difficult to utilize adequately in a cultivar-driven strategy.

Among the 342 segments that influence fiber quality (FL, FS, FU, and FE), 199 were favorable while the remaining 143 were deleterious (Fig. 4e; Additional file 1: Table S21). The ratio of the four categories did not show any imbalance between the favorable and deleterious segments (P values range from 0.9 to 0.13, chi-square test), implying that the genomic resource for cotton fiber quality improvement in CRI12 pedigree is diverse (Fig. 4e; Additional file 1: Table S22). Similar to genetic pattern of 70 FR-related segments, both favorable and deleterious CRI12-specific segments presented low hereditary stability (P value ranged from 0.0008 to 3.5e − 40) (Fig. 4e; Additional file 2: Fig. S12). Moreover, we found that fiber quality-related biparental segments (both favorable and deleterious) were more stable in heredity than segments of the other three categories (P value were from 0.0012 to 3.5e − 40) (Fig. 4e; Additional file 2: Fig. S12). In the genetic trajectory of biparent segments in pedigree formation, we inferred that at least three rounds of selection for hereditary stability had occurred (from progenitor cultivars to Parent#1|Parent#2, from Parent#1|Parent#2 to CRI12, and from CRI12 to progenies). We concluded that only segments of high hereditary stability could be retained in progenies. This result showed that a cultivar-driven strategy could fix favorable segments of biparental category. However, deleterious segments from progenitor cultivars were fixed simultaneously during pedigree formation, and the accumulation of these highly stable segments could narrow the genetic background of the pedigree, damaging the potential of modern cultivars in cotton breeding.

Geographically specific sub-groups in valuable segments

Cultivars in the CRI12 pedigree have been widely planted in China from the Northwest (Xinjiang Province) to the South (Hubei Province), suggesting the spread of favorable segments in CRI12 across different regions. To cluster functional segments into geographically specific sub-groups (SGs), we used the non-negative matrix factorization (NMF), an algorithm that has previously been used to detect genomic signatures [30].

To perform NMF analysis in the Wang_2017 cohort (195 cultivars retained after trimming foreign-introduced cultivars), we generated a present/absent matrix of 70 FR-related CRI12 segments, and two SGs were classified according to the feature-marker matrix (“Methods”; Fig. 5a; Additional file 2: FigS13a). The 195 retained cultivars were from four regions: the Yellow River region (H, 99 cultivars), the Yangtze River region (Y, 66 cultivars), the Northwest region (NW, 21 cultivars), and the early cotton producing region (N, 13 cultivars) (Additional file 1: Table S23). We noticed that segments in SG1 (45 segments) were preferentially present in cultivars from the N, NW, and H regions, while SG2 (25 segments) was preferentially present in cultivars from the Y region (Fig. 5a; Additional file 1: Table S24). The Fusarium wilt resistance of cultivars from the H, NW, and N regions was higher than that of cultivars from the Y region (P value was from 0.006 to 1.8e − 6), and among the three regions, the Fusarium wilt resistance of cultivars did not show any divergence (P value was from 0.79 to 0.97) (Fig. 5b, c). According to the results of the effect evaluation, we found that 95.6% of the segments in SG1 were favorable, while all segments in SG2 were deleterious (Fig. 5d). For each sample, the ratio of SG1 segments to SG2 segments was negatively related to the disease index, indicating that the ratio of SG1/SG2 segments is related to a cultivar’s Fusarium wilt resistance (Fig. 5d). Among the four categories, only CRI12-specific segments showed a significant divergent ratio between the two SGs (P = 0.01, chi-square test), and all favorable CRI12-specific segments belonged to SG1 (Fig. 5e). However, regardless of SG1 preference (those in the Y region) or SG2 preference (those in the H, N, and NW regions), several CRI12-specific segments were severely absent, indicating inadequate utilization of elite genomic resources (Fig. 5f). In these severely absent segments, a segment (D10: 64,941,033–64,944,924) contained a new gene, CRI12_D10G2690 (Fig. 5g), which encodes a KH domain-containing protein (GhKHCP) which is involved in RNA-binding activities. A KH domain protein was shown to enhance the plant immune response in apple [31]. Cultivars without this segment showed a natural knock-out for GhKHCP (P = 0) and a significant increase in Fusarium wilt susceptibility (P = 0.0004) (Fig. 5h).

Fig. 5

Geographical sub-groups in functional segments. a Cluster result of cultivars from 4 regions. The value in heatmap is the median feature-sample score of cultivars from each region. Four regions were clustered into two clades (red and blue dashed lines represent regions where segments belong to sub-group1 and sub-group2 enriched, respectively) by hierarchical method. b Disease index of cultivars from 4 regions in Wang_2017 cohort (n = 195). c t-test for disease index of cultivars from every 2 regions. d The left pieplots represented for ratio of favorable/deleterious segments in SG1 and SG2. The middle scatter plot was the Pearson correlationship between SG1/SG2 segments ratio and disease index (the disease index was normalized by median disease index). The yellow barplot was the SG1/SG2 segment ratio for each cultivar, while the blue barplot was the disease index for each cultivar normalized by median disease index in Wang_2017 cohort (n = 195). e Statistic about number of segments of 4 categories in SG1 and SG2. f The present/absent information of CRI12-specific segments in cultivars from Yangtze River region and the other 3 regions (Yellow River region, Northwestern region, early cotton producing region). g Illustration for deletions of GhKHCP in CRI12 pedigree. h The top boxplot was the transcription abundance of CRI12-like and Parent#2-like (with GhKHCP deletion) genotypes in Wang_2017 cohort (n = 195). The bottom boxplot was the disease index of 2 genotypes (Wang_2017 cohort; n = 195). P value was calculated by t-test and “***” represents for P ≤ 0.001. Boxplot center, median; bounds of box, lower quartile (Q1) and upper quartile (Q3); minima, Q1 − 1.5 (Q1–Q3); maxima, Q1 + 1.5 (Q1–Q3). i The conceptual illustration for bottleneck effect in cultivar-driven strategy. Orange segments represent for segments with low hereditary stability

This segment contained GhKHCP only appeared in 18 and 43 cultivars from the Y and H regions, respectively (Additional file 2: Fig. S14). Correspondingly, there were only 6 progenies that inherited this CRI12-specific segment, showing its low hereditary stability (Fig. 5g). We calculated the present frequency of the 13,138 segments in the pedigree and large population (733 accessions) (Additional file 1: Table S26). Segments retained by less than 6 pedigree members had lower present frequencies in the large population than those retained by more than 12 pedigree members (P = 6.7e − 14), and we inferred that the “bottleneck effect” had occurred on these segments (Additional file 2: Fig. S15). The cultivar-driven strategy is widely applied in modern cotton breeding, in which cultivars from various pedigrees are utilized as elite gene donors. Thus, favorable segments with low hereditary stability are hard to spread widely among cotton germplasm (Fig. 5i).

We also performed NMF analysis on 59 fiber length-related segments from the Wang_2017 cohorts which were divided into two SGs: 29 segments in SG1 and 30 segments in SG2 and these classified cultivars into two clades (Additional file 1: Table S27; Additional file 2: Fig. S13b; Additional file 2: Fig. S16a). Clade I contained the cultivars from the N and NW regions, while clade II had the cultivars from the H and Y regions. Furthermore, we found that there were only 2 and 3 genes directly influenced by segments in SG1 and SG2, respectively (Additional file 1: Table S28). CRI12_A09G0106 in SG2, which encodes an ISWI chromatin-remodeling complex ATPase, was identified because of an insertion in its last exon by Ekangmian10. Cultivars with this insertion had shorter fiber length (P = 0.007) and a lower expression level of CRI12_A09G0106 (P = 1.33e − 5) (Additional file 2: Fig. S16b). ISWI members have been shown to be essential for the formation of heat stress memory in Arabidopsis thaliana [32]. In cotton, we also found that CRI12_A09G0106 had divergent transcriptional patterns under heat and cold stresses (Additional file 2: Fig. S16c). Interestingly, the transcriptional abundance of CRI12_A09G0106 increased during heat stress, while its expression was repressed in long fiber cultivars, implying its contradictory role in the heat stress response and fiber elongation.

Fingerprint segments in CRI12 pedigree

In the above results regarding “bottleneck effects,” we have identified that segments with low hereditary stability always had low present frequencies in large populations. However, we also found several segments that were retained by most pedigree members but were also rare in the large population. We assigned these pedigree-locked segments as “pedigree FingerPrint Segments” (FPS) (Fig. 6a). To detect the FPS quantitatively, we applied TF-IDF (term frequency–inverse document frequency) algorithm (“Methods”) to calculate the fingerprint score of each segment (Additional file 1: Table S29). Finally, we identified 367 FPS by setting the threshold of the fingerprint score as 2, in which no functional segments were included (Fig. 6b and c; Additional file 1: Table S30). We checked the pedigree-population distribution of both FPS and functional segments, indicating that the TF-IDF algorithm had accurately detected the FPS (Fig. 6d). Most functional segments were broadly introduced to cultivars in the large population, while fingerprint segments were restricted to the pedigree members, presenting a pedigree-lock on them (Fig. 6d). The genomic annotation showed that 265 FPS were in the genomic interval region, hindering the comprehensive interpretation of these pedigree-locked segments (Additional file 1: Table S30). Thus, we focused on those influencing genes directly (overlapped with exons, introns, and promoters), and a total of 126 genes were gained. The KEGG enrichment in 126 genes showed that there were no significant enrichment items, indicating that FPS-influenced genes had no pathway preference and these genes affect multiple traits.

Fig. 6

Fingerprint segments in CRI12 pedigree. a Conceptual illustration for pedigree fingerprint segments. b Fingerprint score of 13,138 segments. c Fingerprint score of functional segments. d Presence/absence of 13,138 segments in CRI12 pedigree and 733 cultivars. e