Images of grains from the various accessions of the recombinant inbred lines (RILs) of the mapping population of bread wheat (Triticum aestivum L.) from the International Triticeae Mapping Initiative were used. The ITMI mapping population was obtained by crossing the T. aestivum spring wheat cultivar Opata 85 and the synthetic hexaploid spring wheat W7984 (Arif et al., 2022b). Plants of each genotype were grown in the 2003, 2004, 2009 and 2014 seasons. After harvest, the seeds were stored at IPK genbank with an -18 ± 2 °C and 8 ± 2% seed moisture content. Images were taken from previous works (Afonnikov et al., 2022; Arif et al., 2022b). The images represent grains on a white background, next to the ColorChecker calibration palette, which was used to determine the scale (x-rite ColorChecker^® Classic Mini, https://xritephoto.com/camera; accessed on 20 January 2022). Examples of seed images are shown in Supplementary Figure S1 ( Supplementary File 1 ).

To analyze the diversity of traits in the population and identify QTL, grain images of 114 RILs harvested in 2014 were used (Afonnikov et al., 2022; Arif et al., 2022b). Two images per RIL were obtained for this sample, each counting 15 and 5 grains, respectively. Our preliminary analysis demonstrated that this splitting does not affect the estimation of the seed traits.

Seeds from 44 lines harvested in 2003, 2004, 2009, and 2014 were used to analyze the relationship between texture traits and harvest year (Afonnikov et al., 2022). To analyze the relationship between grain germination and texture characteristics, images of seeds from 19 lines harvested in 2003, 2004, 2009, and 2014 were used (Afonnikov et al., 2022). Images for these samples included from 17 to 20 grains.

These samples were used to perform various types of analysis of grain coat texture traits as summarized in Supplementary Figure S2 ( Supplementary File 1 ): statistical relationship between traits, population diversity analysis, relationship with storage duration in genbank and germination rate, QTL identification, and gene prioritization.

Digital image processing was performed using the SeedCounter application (Komyshev et al., 2017), a desktop PC version supplemented with a color characteristics calculation module (Afonnikov et al., 2022). Areas corresponding to grains were identified in the image, and their size, shape, and color characteristics were evaluated as described earlier. In this study, we took 12 color characteristics, the average values of the color components of the pixels in the grain area for the RGB, HSV, Lab, and YCrCb color spaces (Afonnikov et al., 2022). Their list is given in Supplementary Table S1 ( Supplementary File 1 ).

To evaluate texture properties, we used 16 second-order characteristics determined based on gray level co-occurrence matrices (GLCM) and gray level run-length matrices (GLRM) (Haralick et al., 1973; Galloway, 1975; Majumdar and Jayas, 1999). These characteristics allow us to quantitatively describe the features of micro-relief or surface color inhomogeneities of grains, which are not always distinguishable by visual inspection. GLCM describes the spatial distribution of the brightness of neighboring pixels by evaluating the frequency of co-occurrence of certain combinations of gray-level values. GLRM evaluates sequences of pixels with the same brightness, providing information about the spatial organization of texture elements. This allows characterizing such surface properties as uniformity, contrast, and texture complexity. The list of texture features is given in Supplementary Table S1 ( Supplementary File 1 ), and the definition is given in Supplementary Tables S2 – S4 ( Supplementary File 1 ).

When calculating texture characteristics as both GLCM (for neighboring pixels) and GLRM (for series) matrices, eight main directions are distinguished: up, down, left, right, and four diagonally. In this work, texture characteristics were determined based on grayscale images summed across all 8 directions (omnidirectional).

A preliminary analysis of the images was conducted to exclude outliers from further consideration, i.e., grains with texture characteristics whose values deviated from the mean by more than 3 standard deviations.

To evaluate the Pearson correlation r of texture, size, shape, and color characteristics, grains from the 2014 harvest were analyzed (114 RILs). Based on this, the distance d = 1−∣r∣ between pairs of characteristics was calculated, and then a tree of similarity of characteristics was reconstructed using the UPGMA method. To assess the diversity of texture traits in the ITMI population for grains from the 2014 harvest, the principal component analysis (PCA) method was used based on a correlation matrix.

Grains from 44 RILs harvested in 2003, 2004, 2009 and 2014 were used to assess the contribution of genetic factors (RIL) and harvest year to the variability of texture characteristics implementing one-way analysis of variance (ANOVA). The contribution of a factor was considered significant at p < 0.05.

The linear correlation between the trait value and the harvest year for each grain in this sample was assessed based on the approach proposed earlier (Afonnikov et al., 2022). Data for 3460 grains were used for correlation analysis. The harvest year was coded for each grain in three ways: binary (Year01, values 0 were assigned to the years 2003 and 2004; values 1 were assigned to the years 2009 and 2014); numerical (Year, numerical values of the year were used); rank (YearRank, values 1, 2, 3, and 4 were assigned to the years 2003, 2004, 2009, and 2024, respectively). The significance of correlation between the trait and the harvest year in the three encodings was independently verified using 2000 replicates of permutation and bootstrap tests (randomization of texture trait values was used) in the sample of 3460 grains. The relationship between the trait and the harvest year was considered significant if the correlation coefficient was less than the minimum or greater than the maximum values in both randomization tests.

The statistical relationship between the trait and germination was assessed for grains obtained from 19 RILs harvested in 2003, 2004, 2009 and 2014 (1279 grains in total) using Pearson’s correlation coefficient, as was done previously (Afonnikov et al., 2022). Preliminary evaluation of the contribution of harvesting year and genotype to the germination variance demonstrated that the year (but not genotype) has significant effect. Therefore, to eliminate this effect mean germination values for the corresponding year were subtracted from the each genotype and year values. Some outlier values were removed after that in the germination data (Afonnikov et al., 2022). Before the analysis, each trait values were standardized so that the means were equal to 0 and the standard deviations were equal to 1 for all genotypes. The significance of correlation between the grain coat texture trait and the germination rate was assessed by randomization tests as described above.

Statistical data processing was performed using Python 3.10 software (SciPy, sklearn, pandas libraries).

Grains from 114 RILs harvested in 2014 were used for QTL analysis. Mean values of the grain texture characteristics for each RIL served as input data.

Experimental procedures to obtain SNPs in ITMI plants for QTL analysis were described in (Arif et al., 2021).

To capture the variance explained by the molecular markers such as SNPs mapped to any genome, a refined method known as “inclusive composite interval mapping” was used as implemented in the QTLIciMapping 4.2.53 (http://www.isbreeding.net/(latest released in September 2019). This method currently considered as the most modern method of QTL detection (Arif et al., 2021). It was used successfully to detect several QTLs for Fusarium head blight (Sgarbi et al., 2021) and seed longevity (Arif et al., 2022a) in wheat and germination related traits in tobacco (Agacka-Mołdoch et al., 2021) applying the QTLIciMapping tool. Therefore, we convened the IciMapping 4.2.53 to detect the putative additive QTLs of the traits under consideration applying the inclusive composite interval mapping (ICIM) command where 1.0 cM was the walking speed. An LOD score of > 2.0 ≤ 3 was applied to detect QTLs as significant and > 3.0 as highly significant (Meng et al., 2015).

In order to discover digenic epistasis QTLs to find clues for latent variation, the ICIM-EPI command was used where LOD was kept 5.0 cM. Here, the epistasis QTLs with LOD ≥5 and explaining ≥ 5% phenotypic variance were reported. All QTLs were assigned names according the rules set out in the Catalog of Gene Symbols (McIntosh et al., 2008). Epistasis QTLs were visualized using “circlize” package in R (Gu et al., 2014).

In addition to 16 texture traits for QTL identification, we used their two linear combinations represented by PC1 and PC2 from the PCA analysis for seeds harvested in 2014. (see above). The QTL locations for texture traits were compared with QTLs for traits such as size, shape, and color of grains in images from our previous work (Arif et al., 2022b).

Gene identification in QTL regions, their functional annotation, and prioritization were performed according to previously described procedures (Arif et al., 2022b). The sequences of markers delimiting QTLs were aligned on IWGS 2.1 wheat genome sequence (Zhu et al., 2021). Genome sequence and annotation data were obtained from URGI (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v2.1; accessed on 10 January, 2022). Only ‘high confidence’ gene annotations were considered. Marker sequences were obtained from reference (Wang et al., 2014) and Gramene marker Database (https://archive.gramene.org/markers/; accessed on 10 January, 2022) (Tello-Ruiz et al., 2021). Marker sequences were aligned using blastn of the BLAST+ package (Camacho et al., 2009) using e-value=1e-17 (other parameters were set by default). Marker locations were selected by choosing appropriate chromosome and highest sequence identity with the reference. Search for genes and their functional analysis were performed only for QTLs which have both left and right markers mapped on the reference genome.

Genes located within marker borders were selected by expression level in the grain (TPM>=1). For this purpose, wheat gene expression data from the expVIP database (Borrill et al., 2016) were used. Data in text format were downloaded from URGI (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.1/iwgsc_refseqv1.1_rnaseq_mapping_2017July20.zip; accessed on 10 January, 2022). We used data from RNA-seq experiments in which the column ‘High level tissue’ contains ‘grain’ term. Additional conversion was performed between annotation ver. 2.1 (genome) and 1.2 (transcriptome) gene IDs.

Since there was no prior knowledge about possible molecular mechanisms related to seed texture characteristics in wheat, two approaches for functional annotation and gene prioritization were used. First, full list of QTL related genes expressed in seeds was analyzed by DAVID web service, https://davidbioinformatics.nih.gov/, accessed on 12 April, 2025 (Sherman et al., 2022). Functional clusters and functional charts of genes were obtained. Clusters and functional categories were selected using p-values corrected for multiple hypothesis testing (Benjamini correction and FDR < 0.05). Second, sequences of selected genes were used to search for KEGG Orthology (KO) annotation by BlastKOALA and GhostKOALA web-services (Kanehisa et al., 2016). List of KO IDs was compared with orthologous groups of genes related to seed development in Arabidopsis and rice according to literature data, see details in (Arif et al., 2022b).

The values of Pearson’s correlation coefficients for pairs of traits are given in Supplementary Table S5 ( Supplementary File 1 ). The table shows that significant correlation coefficients are more often observed for pairs within groups of color, size/shape, and texture traits, with several texture traits having significant correlation coefficients with color traits. The tree of grain traits similarity based on these coefficients is shown in Figure 1 . Eight clusters are distinguished in the tree (at a clustering threshold of d = 1.5). The first cluster (top to bottom, light green) includes four grain texture features, three of which are based on GLRM (GLRMsr, GLRMlr, GLRMrr), and one is calculated based on the GLCM matrix (GLCMh). The second cluster (gray) includes two texture features (GLCMu, GLCMmp). The third cluster (pink) includes five grain color features that characterize color (Lab_mb, YCrCb_mCb, Lab_ma, HSV_S, YCrCb_Cr). The next cluster (brown) includes texture features based on GLRM (GLRMglnu, GLRMrlnu, GLRMe). The next two clusters, lilac and red, include shape and size features, respectively. The green cluster includes five texture features based on GLCM calculation. Finally, the orange cluster includes nine features. Of these, two outlying features characterize texture (GLCMc, GLCMi). The remaining features form a tight cluster and include three components of the RGB space; the remaining ones, except for HSV_mH, characterize the lightness/brightness of pixels.

Thus, the clusters on the dendrogram correspond to several interrelated groups of characteristics that describe the shape, size, color, lightness, and texture of the grains. The assignment of characteristics to clusters reflects their common biological nature. Note that texture characteristics are grouped into several clusters. The exceptions are GLCMc and GLCMi, which fall into the brightness trait cluster but are simultaneously quite distant from it. This means that grain texture traits reflect specific surface characteristics that are not related to shape or size and are to some extent related to the lightness of the grain shell ( Figure 1 ).

We analyzed the diversity of grains in the ITMI population based on texture characteristics using PCA. The results are shown in a scatter plot for the two principal components ( Figure 2 ). These two components account for 78% of the total variance. The diagram demonstrates the wide variability of the ITMI population in terms of texture, with no clusters standing out.

The first principal component (PC1) explains half of the total variance and shows a high positive correlation with features such as GLRMglnu, GLCMc, and GLRMe. This component also shows a high negative correlation with the GLCMi feature ( Figure 2 ).

It should be noted that high values of this component are observed for the ITMI_62 line ( Figure 2A ), whose grains appear to be the smoothest ( Figure 2B ). Opata grains also have positive values for this component, and their surface also appears smooth with small wrinkles. Conversely, for the ITMI_11 line, the value of this component is the lowest and negative ( Figure 2A ). In the image, these grains appear to be the roughest ( Figure 2B ). It can be assumed that the first component reflects the smoothness of the grain: the higher its value, the smoother the grain surface; the lower the value, the rougher the grain surface. Note that smooth grains appear lighter in color than rough ones (ITMI_62, ITMI_11, Figure 2B ), which may partly explain the high correlation coefficients between some texture and lightness features in Supplementary Table S5 ( Supplementary File 2 ) and Figure 1 .

The second principal component (PC2) explains a quarter of the total variance and shows a high positive correlation with features such as GLCMcp, GLCMcs, GLCMv, and GLRMrlnu. A negative correlation with this component is observed for the GLCMe feature. Interestingly, a high positive value for this component is observed for the Synthetic_W7984 sample, whose grains appear wrinkled, with wrinkles extending along the grain ( Figures 2A, B ). Low values for this component are characteristic of the ITMI_74 line ( Figures 2A ). Small transverse folds are observed for its grains ( Figure 3B ). It can be assumed that the second component reflects the folding of the grain surface. At the same time, its values probably characterize the direction of the grain surface wrinkles: low values correspond to wrinkles directed across the grain, and high values correspond to wrinkles directed along the grain.

Figure 2 demonstrates the high diversity of grain texture in the ITMI population and allows to distinguish two main types of variability: smoothness/roughness and wrinkling along and across.

Data for various harvest years were presented for 44 lines. ANOVA tests were used to estimate the contributions from both genetic and harvest year factors to the variability of traits related to grain texture in the image. The results are presented in Supplementary Table S6 ( Supplementary File 1 ). Results demonstrated that the variability of all 16 traits is due to a significant contribution of both genotype and harvest year. The p-values were significantly less than 5% for all 16 traits and both factors. The highest p-values were observed for the “genotype” factor and the GLRMglnu and GLRMe traits (0.00047 and 0.00043, respectively). This ensures the reliability of the contribution of the two factors, even when taking into account the correction for multiple comparisons.

Thus, the results of the analysis confirm that the diversity of grain texture traits in the studied wheat samples is influenced by both genetic and environmental (harvest year) factors.

For 44 ITMI lines, trends in grain texture variability depending on the harvest year were assessed based on Pearson’s correlation coefficient with three variants of numerical representation of the harvest year: binary Year01, rank YearRank, and numerical Year. For each representation of the harvest year, we conducted two randomization tests with 2,000 replicates, permutation and bootstrap. With their help, we determined the minimum and maximum confidence limits of the correlation coefficient.

Results are shown in Table 1 . Only two of the 16 texture characteristics do not have a significant correlation between their values and the harvest year. The correlation coefficients for the GLCMv and GLCMe features are less than 0.07 in magnitude. In addition, for the GLCMcs feature, the correlation coefficient with the YearRank parameter (0.07) only slightly exceeds the threshold value (0.06) in magnitude.

In other cases, the correlation coefficient of the trait with the YearRank variable is equal to or greater than 0.1. The highest positive correlation coefficients (greater than 0.28) are observed for traits such as GLCMc, GLCMh, and GLRMlr. These results indicates that these traits are greater for grains from a later harvest year (or shorter storage period in the gene bank). Interestingly, these three traits on the principal component diagram ( Figure 2 ) are collinear with the first principal component, which can be interpreted as an increase in grain smoothness ( Table 1 ).

The lowest correlation coefficients (< -0.26) are demonstrated by GLRMsr, GLRMrr, and GLCMi. The latter has the highest correlation coefficient among all of them with the YearRank feature (-0.44). Negative correlation coefficients mean that the lower the harvest year (and therefore the longer the storage time of the grains), the greater the value of the feature. Note that the three features indicated in the principal component diagram ( Figure 2 ) are opposite in relation to the PC1 component (associated with grain smoothness). Thus, an increase in the value of these features can be interpreted as an increase in grain roughness.

Table 1 shows that the correlation coefficients of the features with the harvest year, determined by different year encodings, are generally consistent with each other. For example, for GLCMh, the correlation coefficient with Year01 was r = 0.269, with Year r = 0.292, and with YearRank r = 0.333. At the same time, the YearRank coding generally shows the lowest absolute values of correlations with features. Therefore, selecting the threshold for randomization tests based on it gives more conservative estimates of significance.

Examples of the relationship between the magnitude of certain traits and the harvest year are shown as box plots in Figure 3 (for other traits, they are shown in Supplementary Figure S3 , Supplementary File 1 ). These graphs clearly show trends in trait variability depending on the year if the correlation coefficient estimate differs significantly from 0. For example, for the GLCMh trait ( Figure 3A ), there is a steady increase depending on the harvest year. This is consistent with the high correlation coefficient values ( Table 1 ).

For the GLCMi trait, on the contrary, the opposite trend is observed: as the harvest year increases, the value of the trait decreases ( Figure 3B ). This is consistent with the negative value of the GLCMi correlation coefficients with the harvest year ( Table 1 ). It should be noted that higher GLCMi values correspond to greater grain roughness ( Figures 2A, B ) and are higher for earlier harvest years, i.e., for longer grain storage periods ( Figure 3 ).

For the GLCMe trait ( Figure 3C ), the values for different harvest years differ, but no trend with increasing harvest year is observed. This is consistent with the results in Table 1 : there is no significant statistical relationship between this trait and the harvest year. In Figure 3D , the values of the GLCMcp trait for the harvest years 2003 and 2004 are slightly higher than the values for 2009 and 2014. There is a noticeable downward trend with increasing harvest year, but it is less pronounced than for the GLCMh and GLCMi traits. This is also consistent with the data in Table 1 : the absolute value of the correlation coefficient with the harvest year for GLCMcp is less than for GLCMh and GLCMi.

Summarizing the results presented, it can be assumed that the correlations between texture characteristics and harvest year that we have identified reflect, in general, an increase in the roughness of wheat grain coat as the storage period in the genbank increases. Thus, the results show that storing grains in a genbank leads to changes in their texture. These changes may reflect structural or metabolic changes in the grain shell.

We assessed the relationship between grain texture characteristics and germination rates. The results are presented in Supplementary Table S7 ( Supplementary File 1 ). They show that only one characteristic, GLCMc (GLCM correlation), meets the criteria for a significant deviation from 0. At the same time, its correlation coefficient with germination (0.098) only slightly exceeds the threshold obtained for the permutation (0.087) and bootstrap (0.097) tests. We note another feature, GLCMi (GLCM inertia). For it, the correlation coefficient with similarity was -0.102, which is slightly less than the lower threshold for the bootstrap test (-0.095), but exceeds the lower threshold for the permutation test (-0.104). Thus, this feature satisfies the criterion of a significant deviation from 0 based on the results of only one randomization test. Interestingly, both GLCMc and GLCMi correlate with grain lightness features (see Figure 1 ). At the same time, the GLCMi was interpreted as characteristics of grain shell roughness ( Figure 2 ).

For the remaining texture parameters, the correlation coefficients between the normalized values of the trait and germination are within the ranges obtained from randomization tests. Thus, it can be concluded that the statistical relationship between grain texture traits and their germination is weak.

For seed texture traits, a total of 36 texture related QTLs were discovered on chromosomes 2D, 3B, 3D, 4D, 5B, 5D, 7A and 7B where majority were overlapped with other QTLs ( Supplementary Table S8 , Supplementary File 2 ; Figure 4 ). For example, there were four QTLs (Q.GLCMe-2D, Q.GLCMcs-2D¹ , Q.GLCMcp-2D¹ and Q.GLRMsr-2D) on chromosome 2D. These QTLs were related to GLCMe, GLCMcs, GLCMcp and GLRMSr where the log of odds (LOD) values ranged from 2.5 to 6.1 and the phenotypic variation explained (PVE) varied from 6.19 to 17.52%. In addition, the QTLs for GCLMCs and GCLMcp overlapped. Likewise, on chromosome 3B, there were six QTLs (Q.GLCMu-3B² , Q.GLCMmp-3B² , Q.GLRMsr-3B³ , Q.GLRMlr-3B³ , Q.GLRMrr-3B³ and Q.GLRMrlnu-3B) linked with GCLMu, GCLMmp, GLRMsr, GLRMlr, GLRMrr and GLRMrlnu. Here the LOD ranged from 2.5 to 3.30 and the PVE ranged from 8.58 to 11.87% whereas all the QTLs except Q.GLRMrlnu-3B were nearly overlapping with each other ( Figure 4 ).

On chromosome 3D there were nine QTLs (related to Q.GLCMu-3D⁴ , Q.GLCMmp-3D⁴ , Q.GLCMc-3D⁴ , Q.GLCMh-3D⁴ , Q.GLCMi-3D⁴ , Q.GLRMsr-3D⁴ , Q.GLRMlr-3D⁴ , Q.GLRMrr-3D⁴ and Q.GLCMe-3D) GLCMu, GLCMmp, GLCMc, GLCMh, GLCMi, GLRMsr, GLRMlr, GLRMrr, and GLCMe and all of them except GLCMe were at the same location. The PVE by these QTLs ranged from 8.5 to 33.23% and the maximum LOD was 9.56. On chromosome 4D, there were three QTLs (Q.GLCMc-4D⁵ , Q.GLCMi-4D⁵ and Q.GLCMh-4D) related to GLCMc, GLCMi and GLCMh that explained 4.48 to 10.24% variation and the LOD remained between 2.67 and 3.41. On chromosome 5B there was one single QTL Q.GLCMh-5B) related to GLCMh responsible for 9.31% variation with an LOD value of 3.40. On chromosome 5D, there were three QTLs (Q.GLCMu-5D⁶ , Q.GLCMmp-5D⁶ and Q.GLRMglnu-5D⁶ ) at the exact location related to GLCMu, GLCMmp and GLRMglnu responsible for > 10% variation and their LOD ranged from 2.53 to 3.52. There were two separate QTLs (Q.GLCMc-7A.1 and Q.GLCMc-7A.2) related to GLCMc on chromosome 7A with LOD values of 9.89 and 14.36 explaining 15.14 and 24.38% variation. Finally, on chromosome 7B, we detected three overlapping QTLs (Q.GLRMsr-7B⁷ , Q.GLRMlr-7B⁷ and Q.GLRMrr-7B⁷ ) for GLRMsr, GLRMlr and GLRMrr where the LOD value was >3 and the PVE ranged from 8.76 to 10.84.

To capture additional variation, we also used the first two principal components as trait values and performed the QTL mapping. Interestingly, we detected five QTLs (three with PCI and two with PC2) ( Supplementary Table S8 , Supplementary File 2 ). The QTLs of PC1 on chromosomes 3B, 3D and 5D (Q.PC1-3B² , Q.PC1-3D⁴ and Q.PC1-5D⁶ ) overlapped exactly with the texture related QTLs. On the other hand, the two QTLs with PC2 (Q.PC2-2D and Q.PC2-5B) on chromosomes 2D and 5B did not overlap with other QTLs. These QTLs explained > 17% phenotypic variance and their LOD values were also > 5.

Additionally, epistatic analyses further detected a total of eight pairs that further explained up to 52.66% variation individually ( Supplementary Table S9 , Supplementary File 2 ; Figure 5 ). From trait perspective, there were four epistatic pairs of QTLs detected for GLCMm on chromosomes 1D-3A, 3B-7A, 4A-7A and 4D-7B. The variation explained by these pairs ranged from 9.79 to 21.96%. One pair for GLCMe was detected on chromosomes 5B-5D responsible for 17.95% variation. Another pair was detected for GLRMglnu on chromosomes 4A-7D causing 21.92% variation in trait expression. Further, another pair was detected for GLrMrlnu on chromosomes 4A-6B. This pair was responsible for 20.58% variation in the trait. We also detected an epistatic QTL pair for PC1 on chromosomes 3B-3D that explained > 50% variation in PC1 ( Figure 5 ).

The results of comparing the positions of QTL for texture traits on chromosomes with QTL for grain size/shape and color traits in images are shown in Supplementary Table S10 ( Supplementary File 2 ). This table shows that QTL for texture traits do not coincide with any of the QTL for grain size and shape traits in the wheat genome. However, there is an overlap with a large number of grain color traits in two locations. The first region is located on chromosome 3B at positions 299.179–300.179 cM (delimited by markers Excalibur_c5309_286 and BS00085434_51).

In this location, three QTLs for texture characteristics (Q.GLRMsr-3B³ , Q.GLRMlr-3B³ , Q.GLRMrr-3B³ ) and 14 QTLs for color traits were found. Second region is located on 3D chromosome at 100–102 cM (delimited by SNPs CAP12_c2615_128 and BS00067163_51). In this location, 9 QTLs for texture characteristics (Q.GLCMu-3D⁴ , Q.GLCMmp-3D⁴ , Q.GLCMc-3D⁴ , Q.GLCMh-3D⁴ , Q.GLCMi-3D⁴ , Q.GLRMsr-3D⁴ , Q.GLRMlr-3D⁴ , Q.GLRMrr-3D⁴ and Q.PC1-3D⁴ ) and 34 QTLs for color traits were observed.

Other notable QTLs for texture characteristics (Q.GLCMu-3B² , Q.GLCMmp-3B² and Q.PC1-3B² ) were located on 3B chromosome at 297.179 cM (delimited by SNPs Excalibur_c36725_96 and Ku_c24974_674). This region is close to seven QTLs for the color traits (position 298.179 cM, SNPs Ku_c24974_674 and Excalibur_c5309_286) and has common SNPs Ku_c24974_674 for these loci.

The remaining QTLs for texture traits do not share common locations on chromosomes with QTLs for color traits.

To prioritize genes based on highly significant QTLs (LOD>3), we selected 27 QTLs out of 36. Of these, 23 QTLs were related to texture traits, two to their linear combination PC1 (on chromosomes 3B and 3D), and two to PC2 (on chromosomes 2D and 5B). In the wheat genome, these QTLs correspond to 13 regions. For nine of them, we determined the coordinates of the left and right markers in the IWGS 2.1 genome sequence. The sizes of the regions ranged from 2 to 80 Mb. Using genome annotation, 1,688 wheat genes were identified that are localized within the boundaries of these regions. The number of genes per region ranged from 15 (QTLs on chromosome 3B at position 297.179) to 526 (15 QTLs on chromosome 4D at position 164.488). Lists of selected QTLs, their position in the genome sequence, and lists of corresponding genes are provided in Supplementary Table S11 ( Supplementary File 2 ). Of the 1,688 wheat genes, 1,213 genes were selected based on their expression level in grain. Their list is provided in Supplementary Table S12 ( Supplementary File 2 ).

These genes were analyzed for functional enrichment using the DAVID service. Table 2 presents information on two functional clusters of genes identified by us that demonstrate a statistically significant association with functional annotation terms. The first cluster includes genes whose function is associated with the INTERPRO domains IPR000490 (Glycoside hydrolase family 17; 13 genes) and IPR044965 ‘Glycoside hydrolase family 17, plant’ (detected for 12 genes), as well as the term GO:0004553 ‘hydrolase activity, hydrolyzing O-glycosyl compounds’ (detected for 21 genes).

The second cluster contains three significant terms. All of them are related to ceramide biosynthesis: ‘ceramide biosynthetic process’ (GO:0046513, 7 genes), ‘Lag1/Lac1-like’ domains (IPR016439, 5 genes), sphingosine N-acyltransferase activity (GO:0050291, 5 genes).

Supplementary Table S13 ( Supplementary File 2 ) presents annotation terms that are significantly represented in the sample of genes localized in regions of highly significant QTLs according to DAVID data. In addition to the terms presented in Table 2 , the table contains Uniprot annotation terms ‘Disordered region’ (detected for 635 genes). Another group of genes is characterized by the annotation of UP_SEQ_FEATURE sequences (COMPBIAS, Basic residues).

Thus, the most significantly represented genes we identified for QTL loci are primarily associated with the glycoside hydrolase family and ceramide biosynthetic process functions.

The list of genes associated with glycoside hydrolase function and ceramide biosynthetic process, together with the locus to which they belong, is given in Supplementary Table S14 ( Supplementary File 2 ). This table shows that 40% of genes related to glycoside hydrolase activity (9 out of 23) are located in the QTL region on chromosome 3D, at position 102 cM. This position was previously noted for its overlap with a large number of QTLs for the color trait ( Supplementary Table S10 , Supplementary File 2 ). This region is associated with the largest number of texture traits, nine, which indicates its high significance. According to Figure 2 , some traits associated with roughness are among these QTLs: Q.GLCMi-3D⁴ (GLCM inertia), Q.GLCMmp-3D⁴ (GLCM max probability), Q.GLCMc-3D⁴ (GLCM correlation), Q.PC1-3D⁴ (first principal component). At the same time, two traits, GLCMi and GLCMc, are closely related to grain shell lightness traits ( Figure 1 ; Supplementary Table S5 , Supplementary File 1 ). Eight of the 23 glycoside hydrolase genes were localized in the QTL region on chromosome 4D, position 164.488. Interestingly, this region corresponds to two QTLs also associated with grain roughness/smoothness traits: Q.GLCMс-4D⁵ (GLCM correlation) and Q.GLCMi-4D⁵ (GLCM inertia), which are associated with grain shell lightness traits.

Most of the genes associated with the ceramide biosynthetic process are located on chromosome 5B, position 75 cM (4 out of 7). Two genes are located on chromosome 4D and one on chromosome 3D.

KEGG orthologous groups were identified for 531 of the 1,213 genes. Analysis of the KEGG annotation showed that a significant proportion of the genes (almost half) belong to the functional category “Genetic information processing.” The next categories are ‘Carbohydrate metabolism’, ‘Signaling and cellular processes’, ‘Environmental information processing’, and ‘Lipid metabolism’. Based on KO identifiers, 35 genes associated with seed development processes were identified. Their list is given in Supplementary Table S15 ( Supplementary File 2 ). They are represented in all loci we detected, except for chromosome 3B. Most (14 genes) are located in the 4D region of the chromosome, at position 164.488 cM. The genes associated with seed development include transcription factors (MADS-box, EREBP-like, MYB, HD-ZIP, AP2-like, NFYC, 19 genes in total), while the remaining genes encode various enzymes, kinases, translation regulators, and a number of others. Based on these data, it is difficult to identify any functional group specific to these genes.