Three NILs obtained in the spring cv. Bowman (“Bowman From Fargo,” NGB22812)¹ genetic background was used in the study. The first one is the i:BwBlp1 line (NGB20470; hereafter: BLP – black lemma and pericarp) carrying the Blp1 locus mapped to chromosome 1HL (Druka et al., 2011) and characterized by melanin accumulation in hulls and the grain pericarp. The second NIL is i:BwAnt1Ant2 (NGB22213; hereafter: PLP – purple lemma and pericarp), which is characterized by anthocyanin accumulation in the grain pericarp and leaf sheath and carries complementary genes Ant1 and Ant2 located on chromosomes 7HS and 2HL, respectively (Druka et al., 2011). The third NIL is the i:BwAnt1Ant2Blp1 line (hereafter: BP – black and purple), which was developed previously via marker-assisted selection on the basis of lines BLP and PLP and is characterized by the accumulation of both anthocyanins and melanins in the grain (Glagoleva et al., 2022). The plants were grown in a greenhouse of the Institute of Cytology and Genetics (ICG) SB RAS (Novosibirsk, Russia) under a 12 h photoperiod in a temperature range of 20–25°C.

To analyze the phenylpropanoid content of the grains, three NILs (BLP, PLP, and BP) and cv. Bowman (parental line) were collected at the early dough stage (BBCH-83) of spike maturity. The whole grains were frozen in liquid nitrogen and dried in vacuum. Then, the dried grains were ground up and used for phenolic-compound extraction. Phenolic compounds (phenolic acids, flavonoids, and anthocyanins) were extracted from 300 mg of mature seeds with 2.5 mL of a methanol:2% formic acid mixture (7:3, v/v) (Shoeva et al., 2016) with ultrasonication at 40°C. Three independent extraction procedures were performed for each genotype under study. After centrifugation, the extracts were analyzed with Agilent 1100 liquid chromatography system (Agilent Technologies, Santa Clara, CA, United States) equipped with a quaternary pump, online degasser, autosampler, and diode array detector. Separation was implemented on an Eclipse SB-C18 (4.6 mm × 150 mm, 5 μm) column in a binary solvent system consisting of 0.1% CF₃COOH in water and methanol. The gradient of methanol was linear from 2 to 98% for 40 min, the injection volume was 4 μL, the column temperature was 23 ± 2°C, the flow rate was 0.8 mL/min. The high-performance liquid chromatography (HPLC) runs were monitored at 254 nm for benzoic acids, at 320 nm for cinnamic acids, at 360 nm for flavonoids, and at 520 nm for anthocyanins. Identification and quantification of benzoic and cinnamic acids, flavonoids, and anthocyanins were based on comparisons of the UV data and HPLC retention times of standard compounds. The limit of detection (LOD) was 0.2 μg/g in sample and the limit of quantification (LOQ) was 0.5 μg/g in sample.

The following reagents were used as a standard compounds: 3,4-dihydroxybenzoic acid, catechin, dihydromyricetin, vanillic acid, trans-caffeic acid, 4-hydroxycinnamic acid, dihydroquercetin, trans-ferulic acid, dimethoxybenzoic acid, 3-hydroxycinnamic acid, dihydrokaempferol, hesperidin, rutin, trihydrate, 2-hydroxycinnamic acid, cinnamic acid, luteolin, apigenin, quercetin, naringenin, 3-O-glucoside cyaniding, 3-glucoside malvidin (Merck, Darmstadt, Germany).

For a comparative high-throughput RNA sequencing (RNA-seq) analysis, total RNA was isolated from the lines BLP, PLP, BP, and parent cv. Bowman at three stages of spike development. The stage 1 is the booting stage, when first awns are visible (BBCH-49); stage 2 is the late milk stage (BBCH-77), when the start of anthocyanin pigmentation is observed; and the stage 3 is the early dough stage (BBCH-83), when melanin pigmentation appears. At the stage 1, RNA was extracted from fresh whole flowers, which was separated from the spike and homogenized in liquid nitrogen. At the stages 2 and 3, pericarp and hulls were peeled from fresh grains of the NILs and homogenized in liquid nitrogen. Spikes of near-isogenic lines under study are presented in Figure 2. RNA was extracted in three biological replicates using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany); each replicate combined material from three plants (i.e., nine plants of each line were analyzed in total as three biological replicates). After the removal of DNA traces with the Turbo DNA free kit, RNA quality was assessed using Bioanalyzer 2100 with the RNA Nano kit. An mRNA fraction was isolated, and barcoded RNA-seq libraries for the Illumina system were generated by means of TruSeq Stranded mRNA Library Preparation Kit A. The resulting 36 barcoded libraries (3 biological replicates × 4 barley lines × 3 stages of development) were sequenced on an Illumina NextSeq 550 instrument with the NextSeq 500 High Output v2 Kit (75 cycles).

FASTQC v.0.11.9 was employed to evaluate the libraries’ quality (Andrews, 2010). Filtering of the libraries was performed in Trimmomatic software v.0.39 with the following parameters: “LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:50” (Bolger et al., 2014). This approach results in the removal of unidentified bases (N) or bases with a Phred quality score below 20 from both 3′ and 5′ ends of a read, in the removal of reads with length less than 50, and in the trimming of 5′ ends of the read when mean quality in the sliding window of length 4 drops below 30. DART tool v.1.4.2 (Lin and Hsu, 2018) was used to align the filtered reads to barley genome assembly IBSC v.2 (Mascher et al., 2017) release 47 from the Ensembl Plants database². The reads aligned to each gene were counted with the help of the featureCounts function in the Subread software (Liao et al., 2014). Based on the obtained counts, principal component analysis was conducted using the DESeq2 function of plotPCA (Love et al., 2014). Raw read counts data before normalization and genes expression after CPM (counts per million mapped reads) normalization are presented in Supplementary Data Sheets 1, 2.

Differential gene expression analysis was performed by means of the edgeR package for R (Robinson et al., 2010). Weakly expressed genes were eliminated using the “filterByExpr” function. Differential expression between samples was detected via the generalized linear model approach. Differential expression was calculated between time stages within each barley line and between barley lines within each time stage. Thus, 26 comparisons were performed. Genes with a false discovery rate (FDR) of <0.05 and an absolute value of a log₂-transformed fold change greater than 2 (| log₂(FC)| > 2) were considered differentially expressed.

For functional annotation of differentially expressed genes (DEGs), Gene Ontology (GO) enrichment analysis and metabolic pathway involvement analysis were undertaken. Singular enrichment analysis of DEG lists was performed using the AgriGO v.2 toolkit (Tian et al., 2017). For each comparison performed for the differential expression analysis, lists of upregulated and downregulated genes were examined separately. GO terms with adjusted p-value < 0.05 were regarded as significantly enriched. KEGG pathways enriched in the DEG lists were identified in BlastKoala v.2.2 (Kanehisa et al., 2016). In this case, amino acid sequences of DEGs were subjected to a BLAST search against a non-redundant dataset of pangenome sequences at the genus level and then annotated in KEGG Ontology (KO) terms by the KOALA algorithm. Barley genes involved in the phenylpropanoid biosynthesis pathway were detected via a homology search in the BioMart software using Arabidopsis thaliana genes participating in the phenylpropanoid biosynthesis pathway (KEGG id: KO00940) as a query. Heatmap figures were generated using the seaborn v.0.11.2 library in Python.

Nine genes with different expression patterns were chosen randomly for the qRT-PCR verification of the RNA-seq data: HORVU1Hr1G011930, HORVU2Hr1G103040, HORVU2Hr1G103000, HORVU3Hr1G077790, HORVU4Hr 1G090870, HORVU6Hr1G075900, HORVU3Hr1G077960, HORVU7Hr1G098280, and HORVU4Hr1G014010. The qRT-PCR was performed on a QuantStudio 5 machine (Applied Biosystems, Waltham, MA, United States)³ with the HS-qPCR Lo-ROX SYBR kit (Biolabmix, Novosibirsk, Russia) in a 15 μL reaction mixture. The number of amplification cycles and annealing temperature were optimized for each primer pair (Supplementary Table 1). Three technical replicates of each reaction were run. Gene expression levels were calculated by the relative-standard-curve method and normalized to the geometrical mean of actin and ubiquitin expression (Himi et al., 2005; von Zitzewitz et al., 2005).

A total of 36 mRNA libraries were sequenced that were derived from grain envelopes of the three NILs (BLP, PLP, and BP) and parental cv. Bowman at three stages of spike development: stage 1, booting; stage 2, late milk; and stage 3, early dough (three biological replicates for each type of sample). In total, 801,184,069 single-end 75 bp raw reads were obtained. After filtering, 741,311,996 (92.51%) reads were retained. On average, 92.43% (71.11% unique) reads were successfully mapped to the reference barley genome in the DART software (Supplementary Table 2).

After the removal of genes with low expression and CPM (counts per million mapped reads) normalization, there were 30,969 genes with an expression level higher than the significance threshold. Then, differential expression was calculated between stages within each line and between lines within each stage (Supplementary Table 3). An average, the number of upregulated and downregulated DEGs between the NILs at the stage 1 were 23 and 54, respectively, at the stage 2 – 283 upregulated and 130 downregulated DEGs, and at the stage 3 – 677 upregulated and 304 downregulated DEGs between the NILs under study were revealed (Supplementary Figures 1–3). In a comparison of numbers of DEGs between stages 1 and 2 within the lines, on average 1,788 upregulated and 2,415 downregulated DEGs that are common to all lines were found. The average number of line-specific DEGs was approximately an order of magnitude less than the number of common genes (Figure 3). In contrast, 76 upregulated and 15 downregulated DEGs that are common for all lines were revealed between stages 2 and 3, whereas the number of line-specific DEGs was substantially greater in all lines characterized by pigment accumulation but not in the uncolored parental line (Bowman). It is possible that at the early stages of spike development, genes responsible for plant primary metabolism that are common to all the tested lines were activated regardless of grain pigmentation; on the contrary, during spike maturation, genes of special or secondary metabolism may be activated, which is associated with pigment synthesis in the barley grain.

The lists of DEGs between different lines and different stages were annotated by singular enrichment analysis of GO terms to gain a better understanding of the DEGs functions. In total, 49.4% of upregulated genes and 52.2% of downregulated genes were annotated successfully. Full lists of DEGs along with GO term enrichment are presented in Supplementary Table 4.

At the booting stage (stage 1), no significantly enriched GO terms were found between the non-pigmented Bowman cultivar and the lines characterized by accumulation of melanins and/or anthocyanins. In contrast, there were differences between the lines in metabolic processes in conjunction with pigmentation at the later stages. In the lines accumulating melanins (BLP and BP), the lists of upregulated DEGs at stages 2 and 3 were found to be enriched with GO terms related to fatty acid biosynthesis, lipid biosynthetic process, and organic acid biosynthesis when compared with Bowman. On the other hand, the upregulation of genes involved in the L-phenylalanine metabolic process and ammonia-lyase activity was detectable only in the BLP line at the later stage of spike maturation (stage 3) in comparison with Bowman (Figure 4). In the list of downregulated DEGs, GO terms “photosynthesis,” “thylakoid,” and “photosystem” were found for the BLP line at stage 2, but they were absent in the other pigmented lines in comparison with Bowman. Moreover, the genes involved in cell wall biogenesis and assembly, lignin metabolic process, and cell growth turned out to be downregulated in lines BLP and BP when compared with Bowman at stage 3 but not in the PLP line. In the comparison of the PLP line with Bowman, there were no line-specific GO terms that were not detectable in the other lines.

Barley genes encoding the main enzymes of the phenylpropanoid biosynthesis pathway were identified in the KEGG database via homology with A. thaliana genes (Supplementary Table 5). Most of these enzymes are encoded by gene families in the barley genome; therefore, we used mean levels of gene expression for heatmap construction (Figure 5A). According the expression patterns, the genes can be subdivided into two groups: early and late genes. The early genes are activated at the booting stage (stage 1) in all four lines regardless of pigmentation. These include the following genes: PAL, C4H, 4CL, C3′H, CSE, CCoAMT, COMT, F5H, CCR, CHS, CHI, and FLS. The late genes are transcribed at later stages of spike development (stages 2 and 3), and their expression proved to be line-specific and depended on the pigment type accumulated. In lines characterized by the presence of anthocyanins (PLP and BP), specific activation of the flavonoid biosynthesis pathway at the later stages was observed (Figure 5A). The pathway includes enzymes F3H, F3′H, and ANS. Furthermore, secondary activation of genes encoding CHS and CHI, which are implicated in the key steps of flavonoid biosynthesis, was registered in lines PLP and BP. It should be noted that in the BP line, these genes were found to be activated earlier, at stage 2, whereas in BP’s parental line PLP, they are highly transcribed only at stage 3. Moreover, activation of some specific genes from the monolignol branch of the biosynthesis pathway was observed in lines PLP and BP, including gene HORVU3Hr1G056560 encoding CCR, whereas the expression of gene HORVU3Hr1G099990 encoding CCoAMT was noticed specifically in the PLP line (Figure 5B).

The expression of regulatory genes Ant1, Ant2, and WD40 was investigated too; they are necessary for MBW complex assembly and anthocyanin biosynthesis in the grain pericarp. Specific activation of Ant1 and Ant2 was detectable in anthocyanin-accumulating lines; moreover, the expression pattern of the Ant2 gene was the same as that of the late genes of flavonoid biosynthesis pathway, and activation timing was different between lines PLP and BP. By contrast, the expression of WD40 was constant in all the lines at all examined stages (Figure 5A).

Similarly, the genes with specific activation in the presence of melanin pigments at stages 2 and 3 were identified in the phenylpropanoid pathway (Figure 5B). They are represented by specific genes from gene families encoding enzymes of the general phenylpropanoid pathway such as PAL (gene HORVU2Hr1G089440), and enzymes participating in monolignol synthesis such as COMT (HORVU1Hr1G011930), CSE (HORVU3Hr1G087020), CCR (HORVU3Hr1G056560 and HORVU7Hr1G093260), and CAD (HORVU2Hr1G010830 and HORVU0Hr1G012940). In addition, genes encoding C4H and 4CL proved to be activated at stage 3 in all pigmented lines; what is more, in the case of C4H, there is secondary activation of this gene in response to Ant1, Ant2, and Blp1 presence.

Given that melanin is a product of phenolic-compound polymerization, the expression of PPOs was also analyzed. It was observed that PPO genes are inactive at stage 1, and their expression increases during spike maturation in all the analyzed lines, whereas in BLP, these genes are most transcriptionally active (Figure 5A).

Expression of nine genes was confirmed by qRT-PCR. These genes code for o-methyltransferase (HORVU1Hr1G011930), polyphenol oxidases (HORVU2Hr1G103040, HORVU2Hr1G103000, HORVU3Hr1G077790, and HORVU4Hr1G090870), chloroplastic geranylgeranyl diphosphate reductase (HORVU6Hr1G075900), translation initiation factor SUI1 (HORVU3Hr1G077960), xyloglucan endotransglucosylase/hydrolase (HORVU7Hr1G098280), and chloroplastic photosynthetic NDH subunit of subcomplex B1 (HORVU4Hr1G014010). In vitro expression analysis of these genes was highly consistent with the in silico data (Spearman’s correlation was 0.93), thus indicating the reliability of our RNA-seq data (Figure 6).

Extraction of phenolic compounds with subsequent HPLC analysis was carried out at the 3rd stage of spike development, when the greatest of line-specific differences in gene expression were observed. The quantities of identified phenolic compounds and chromatographic profiles are presented in Supplementary Table 6 and Supplementary Figures 4–7, respectively. The total contents of benzoic acids, hydroxycinnamic acids, flavonoids and anthocyanins are presented in Figure 7A. The total benzoic-acid content was significantly higher in melanin-containing NILs BLP (25.2 ± 0.3 μg/g) and BP (38.6 ± 4.6 μg/g) as compared with PLP (20.6 ± 0.7 μg/g) and cv. Bowman (21.5 ± 0.2 μg/g), whereas in the BP line, it was statistically significantly the highest among all the lines. Among the compounds in this family, protocatechuic, 4-hydroxybenzoic, and vanillic acids were identified. It is noteworthy that an elevated concentration of protocatechuic acid was found to be characteristic of lines PLP and BP, which accumulate anthocyanins; a significantly increased content of vanillic acid was also documented in the BP line, whereas in BLP, the increase was not significant (Figure 7B).

The highest level of hydroxycinnamic acids was registered in PLP (118.2 ± 5.9 μg/g); however, this level was also significantly higher in BP (108.7 ± 3.4 μg/g) than in BLP (78.8 ± 3.9 μg/g) and Bowman (82.2 ± 7.8 μg/g). In this family of compounds, caffeic acid, p-coumaric acid, m-coumaric acid, and ferulic acid were identified. Furthermore, four additional compounds were assigned to the family of cinnamic acids according to their characteristic UV spectra. Two patterns of these acids in the lines under study can be distinguished. The first one concerns caffeic acid accumulation and represents its significantly increased content in NILs with melanin. The second pattern concerns p-coumaric and ferulic acids: greater accumulation in lines PLP and BP (Figure 7B).

The flavonoid content in Bowman and BLP was 41.2 ± 2.3 and 25.3 ± 1.7 μg/g, respectively: dramatically different from the flavonoid content in PLP and BP (404.5 ± 16.5 and 213.8 ± 10.3 μg/g, respectively). In this family of compounds, luteolin, apigenin, and quercetin 3-rhamnoglucoside (rutin) were identified. Besides, a series of glycosides of luteolin and apigenin was identified by means of characteristic UV spectra. In the subfamily of luteolin derivatives, a flavonoid was found, which, according to the UV spectrum and retention time, can be luteolin methyl ester. This compound was not found in Bowman and BLP and is one of the major phenolic compounds in lines PLP and BP. Anthocyanins were not detectable in lines Bowman and BLP, but their presence in lines PLP and BP was confirmed. In the family of anthocyanins, cyanidin-3-glucoside and malvidin-3-glucoside were identified. It was noticed that the total concentration of flavonoids was ∼1.5 times higher in PLP than in BP; however, the total anthocyanin content was significantly greater in BP than in its parental line PLP (26.2 ± 1.7 and 18.0 ± 0.4 μg/g, respectively).