Abstract

Background

   Highland barley (Hordeum vulgare L. var. nudum) is a key crop of the
   Qinghai-Tibet Plateau, renowned for its nutritional value and
   exceptional adaptability to high-altitude environments. Induced
   mutagenesis offers a powerful approach to developing new crop varieties
   and elucidating the genetic basis of functional traits.

Results

   In this study, nitrogen ion beam implantation was employed to induce
   mutations in the highland barley cultivar Kunlun 14 (K14), generating
   71 novel mutation materials and enriching the genetic resources for
   barley breeding. Phenotypic trait correlation analysis identified two
   mutation lines exhibiting significant variations: E8-38 with highly
   increased 1000-grain weight, and D7-67 displaying a two-row spike
   phenotype. Comparative transcriptomic analysis was applied to the above
   two mutation materials. It was revealed that the high 1000-grain weight
   of E8-38 was driven by synergistic regulation of phytohormone
   signaling, metabolic pathways, and epigenetic modifications,
   particularly the upregulation of the cytokinins signaling pathway and
   starch metabolism genes. In D7-67, the two-rowed spike phenotype was
   underpinned by the upregulation of VRS1 genes. Adaptive mechanisms to
   high-altitude environments were investigated, revealing upregulation of
   PAL and 4CL genes in phenylpropanoid biosynthesis, which enhances UV
   resistance and antioxidant capacity. Additionally, optimization of the
   photosynthetic pathway may contribute to acclimation under harsh stress
   conditions. Through PPI analysis, BZIP transcription factors were found
   to regulate downstream genes, facilitating adaptation to light changes
   and oxidative stress.

Conclusion

   Nitrogen ion beam implantation was demonstrated to be an efficiency
   method to introduce mutation on the highland barley, generating a set
   of mutation materials with a wide range of genetic variations. Through
   comparative transcriptomic analysis, this study elucidates the
   molecular basis underlying the 1000-grain weight and spike-type
   mutants, as well as the adaptive mechanisms enabling highland barley to
   thrive in high-altitude environments. These findings provide critical
   insights into the genetic and molecular mechanisms driving high-yield
   and stress-resilient traits in highland barley.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12864-025-11856-8.

   Keywords: Highland barley, Ion beam implantation, Mutation, High
   altitude adaptation, Transcriptomic sequencing

Introduction

   Highland barley (Hordeum vulgare L. var. nudum) is the major crop in
   the Qinghai-Tibet Plateau of China, which belongs to the Poaceae
   family. Characterized by its separation of the glumes, which exposes
   the seed, it is a variant of two-rowed barley. With a short growing
   period, strong stress resistance, and wide adaptability, it is
   particularly well-suited for cultivation in the Qinghai-Tibet Plateau
   region. It serves as a typical example of evolutionary adaptation in
   extreme high-altitude ecological environments. Barley is one of the
   earliest domesticated crops [[30]1], with widespread applications in
   brewing, food production, and other industries. Owing to its extensive
   adaptability and excellent stress tolerance, barley is cultivated
   globally and serves as a key model species for molecular research in
   cereal crops. The harsh environmental conditions of the Qinghai-Tibet
   Plateau, coupled with human-induced domestication pressures, have led
   to the evolution of enhanced stress resistance in highland barley,
   making it an ideal model species for studying stress tolerance
   mechanisms in cereal crops.

   Altitude is one of the key environmental factors influencing plant
   growth. As altitude increases, changes in temperature, sunlight, and
   precipitation occur, which subsequently affect plant development. Due
   to the complexity of impact of altitude on plant growth, research on
   plant adaptation to high-altitude conditions is limited. By means of
   genomic sequencing of the highland barley and analysis of the genome
   regions under selective sweeps due to the plateau environment, many
   genes involved in the pathways of plant hormone signal transduction,
   phenylpropanoid biosynthesis, flavonoid biosynthesis, etc., were
   identified to be directly or indirectly associated with environmental
   stress adaptation [[31]2]. In addition, several studies focusing on
   individual factors such as light or drought are more detailed and
   comprehensive. For example, for UV-B adaption of the highland barley,
   several genes of phenylpropanoid pathway were identified to be key
   determinants controlling natural variation of phenylpropanoid content,
   resulting in UV-B protection [[32]3]. In barley, the plastidnucleus
   localized DNAbinding protein WHIRLY1 was identified to be required for
   acclimation to high light [[33]4]. Through comparative transcriptomic
   analysis, several metabolic pathways were identified in response to the
   drought stress in highland barley [[34]5, [35]6]. In another
   high-altitude plant Primula sikkimensis, the genes involved in
   metabolic processes, secondary metabolism, and flavonoid biosynthesis
   as well as photosynthesis and plant hormone signal transduction were
   revealed to be modulated across various altitudes [[36]7]. Furthermore,
   the Arabidopsis UV receptor UVR8 dimer dissociated into monomer once
   perception of UV-B, which was then bound to the E3 ubiquitin ligase to
   form the COP1 complex that accumulated in the cell nucleus to induce
   photomorphogenic responses and UV-B adaptation [[37]8]. Although the
   exact mechanism by which UVR8-COP1 regulates gene expression remains
   unclear, it is known to be involved in the stabilization and activation
   of BZIP transcription factors [[38]9].

   The morphology of the inflorescence is a key determinant of grain yield
   in Poaceae crops. The inflorescence structure of barley differs from
   that of other cereal grasses such as rice, wheat, and maize [[39]10].
   Barley features a specialized, unbranched, spike-like inflorescence
   [[40]11]. The differentiation of two-rowed and six-rowed spike types is
   a unique feature of barley, with VRS1 [[41]12], VRS2 [[42]13], VRS3
   [[43]14], VRS4 [[44]15], and VRS5 [[45]16] identified as genes
   controlling row-type, encoding transcription factors involved in
   lateral organ boundary structures. The gene Int-C (Intermedium-C)
   [[46]17], which encodes a transcription factor, regulates lateral
   fertility in spikelets, controlling spikelet development and
   influencing variations in male fertility and grain development between
   two-rowed and six-rowed spikes. Another class of yield-determining
   factors involves plant hormone biosynthesis and metabolism, including
   cytokinins (CKs), auxin and ethylene. Cytokinin dehydrogenase (CKX)
   catalyzes the degradation of CKs [[47]18]. In rice mutants with the
   OsCKX2 gene knockout, a significant increase in grain number was
   observed, likely due to reduced CKX activity, leading to the
   accumulation of CKs in the inflorescence meristem [[48]19]. In barley,
   silencing two HvCKX genes using RNAi technology also resulted in an
   increase in seed number and spike number [[49]20]. Auxin, the first
   identified plant hormone, plays a crucial role in regulating plant
   growth by promoting cell division and proliferation [[50]21]. It is
   involved in various developmental stages, such as seedling growth,
   gamete formation, and embryo development [[51]22], and also influences
   grain development and maturation through its signaling pathways
   [[52]21]. VILLIN 2 (VLN2), an actin-binding protein (ABP), demonstrates
   conserved actin filament bundling activity in vitro. Mutations in VLN2
   lead to hypersensitive gravitropic responses and accelerated recycling
   of OsPIN2, resulting in altered plant architecture and wrinkled seeds
   [[53]23]. Ethylene, the only gaseous hormone produced in plants,
   regulates diverse physiological processes, including leaf senescence,
   fruit ripening, and seed dormancy [[54]24]. Ethylene receptor 2 (ETR2),
   which encodes a serine/threonine kinase, has been shown to negatively
   influence grain weight when overexpressed [[55]25]. This finding
   strongly suggests a positive correlation between ethylene sensitivity
   and traits such as grain weight and size.

   Barley is one of the most genetically diverse cereal crops, classified
   into spring and winter types, two-rowed and six-rowed varieties, among
   other classifications based on distinct traits [[56]26]. Therefore, the
   high level of genetic diversity provides significant opportunities for
   progress in breeding programs. Early breeding efforts primarily focused
   on high yield and disease resistance. Through hybridization and other
   methods, numerous lines, such as “Kunlun 1” [[57]27] and “Zangqing 320”
   [[58]28], were developed to meet the increasing demand for barley in
   the Qinghai-Tibet Plateau, driven by population growth. These breeding
   efforts helped meet the growing demand for barley in the Qinghai-Tibet
   Plateau region, driven by population expansion. However, as awareness
   of the health benefits of Tibetan barley has increased, breeding
   objectives have diversified. While contemporary barley genomics and
   transcriptomics have elucidated numerous stress-responsive genes and
   pathways, significant challenges remain in generating and
   characterizing novel genetic variations for crop improvement.
   Conventional mutagenesis methods typically produce limited mutational
   spectra, and genome editing predominantly targets known genes. Although
   natural genetic diversity studies provide abundant variation
   information, these variants are often constrained by natural selection
   and environmental factors. In contrast, ion beam irradiation induces a
   broader mutation spectrum, including complex structural variations and
   clustered DNA damage, fundamentally differing from spontaneous
   mutations or targeted editing. This technology also enables rapid
   generation of large mutant populations, facilitating efficient
   screening for agronomically valuable traits.

   Ion beam biotechnology, a novel physical mutagenesis technology,
   induces chromosomal and genetic variations [[59]29]. It offers
   advantages such as low damage, high mutation rate, and a broad mutation
   spectrum, making it widely applied in crop breeding. For example, rice
   seeds subjected to heavy-ion irradiation have produced mutants with
   reduced cadmium accumulation, early maturity, dwarfism, and chlorosis
   [[60]30, [61]31]. Sweet sorghum seeds irradiated with carbon ions,
   followed by multi-generational field selection, led to the development
   of a new high-yield, high-sugar variety, “Jin Tians1” [[62]32].
   Transcriptomic sequencing of the Vlamingh barley variety, induced by
   γ-ray irradiation, revealed the mechanisms of γ-ray-induced mutations
   [[63]33], proving that irradiation mutagens can be used in the
   development of important genetic resources and crop breeding. Previous
   studies have demonstrated the effectiveness of nitrogen ion beam
   irradiation in mutagenesis applications. In wheat, irradiation
   treatment successfully generated a mutant population consisting of 60
   distinct lines, with the M4 generation showing phenotypic variations in
   growth duration, agronomic traits, gliadin content, and SSR
   polymorphisms [[64]34]. Similarly, research in rice has revealed that
   low-energy ion beams can effectively induce mutations while maintaining
   low biological damage, achieving both high mutation frequency and a
   diverse mutation spectrum [[65]35]. Compared to traditional hybrid
   breeding, ion beam biotechnology significantly shortens the breeding
   cycle, making it an ideal method for cultivar selection and
   development. Among various mutagenesis techniques evaluated, nitrogen
   ion beam irradiation demonstrates distinct advantages: it produces
   higher mutation rates with wider mutational diversity while causing
   less physiological damage compared to chemical mutagens [[66]34], and
   achieves high-frequency mutations at low doses. While carbon ion beams
   induce similarly high mutation rates, their higher linear energy
   transfer may cause more extensive genomic damage [[67]36]. Notably,
   nitrogen ion beam’s efficiency varies across species and is influenced
   by seed coat thickness, sometimes resulting in relatively lower
   mutation frequencies than high-energy carbon ions. Nevertheless,
   nitrogen ion beam has successfully generated stable mutants with
   improved agronomic traits in wheat and rice, establishing it as an
   effective and relatively gentle mutagenesis approach particularly
   suitable for crops with complex genomes like highland barley.

   In this study, we employed ion beam implantation technology to induce
   mutations in the elite highland barley variety “Kunlun 14” (K14) by
   treating seeds with nitrogen ion beams of varying energies and doses.
   Through winter propagation and generation acceleration, we rapidly
   selected and constructed a novel mutant population, broadening the
   genetic base and enriching the parental pool for high-quality barley
   breeding. Furthermore, we performed RNA-seq-based differential
   expression genes (DEGs) analysis combined with functional enrichment
   (GO/KEGG) and SNP to elucidate the genetic basis of phenotypic
   variations in two unique mutants. In addition, we investigated the
   adaptation mechanisms of highland barley to high-altitude environments.
   In summary, this study aims to identify mutant lines with desirable
   agronomic traits through mutagenesis breeding and transcriptomic
   analysis, providing valuable genetic resources and theoretical insights
   for the genetic improvement and breeding of highland barley varieties.

Materials and methods

Plant material

   The highland barley cultivar “Kunlun 14” (here referred as K14) was
   used as plant material throughout this study. The K14 variety was
   developed through sexual hybridization with “Bai 91-97-3” as the
   maternal parent and “Kunlun 12” as the paternal parent, which has been
   widely planted in northern Qinghai-Tibet Plateau [[68]37].

Nitrogen ion beam implantation treatment

   Dried seeds of the highland barley cultivar K14 were divided into six
   groups, with each group consisting of about 5,000 seeds. These groups
   were subjected to nitrogen ion beam implantation with various ion beam
   intensities and dosages at the Southwestern Institute of Physics in
   Chengdu, Sichuan, China.

Acquisition of mutation materials and phenotypic data analysis

   The six groups of seeds subjected to ion beam implantation treatment
   were sown in a row-seeding fashion in Luhuo County, western Sichuan,
   China (31°39’ N, 100°67’81” E, altitude 3,229 m, referred as Highland),
   with approximate 4,200 seeds per treatment (25 rows × 170 seeds). These
   seeds generated M1 generation mutation materials. By the end of the
   same year, M1 seeds were sown by a single spike per row in a winter
   breeding field in Chongzhou City, Chengdu (30°40’4” N, 103°43’28” E,
   altitude 540 m, referred as Lowland). Routine field management was
   conducted after germination, and phenotypic observations were made
   periodically throughout the growth period. Based on the phenotypic
   variations and agronomic traits, 71 mutation materials were selected
   and harvested in the following year. For each material, 15–30 spikes
   were collected and designated as M2 generation mutation materials.

   Subsequently, the seeds of each M2 mutation materials were divided into
   two parts and sown at two locations with different altitudes, Chongzhou
   City (Lowland) and Luhuo County (Highland) for field trials. Field
   observations and data collection continued throughout the growth
   period. At harvest, seeds were collected from both locations, which
   were designed as the M3 generation. In addition, 30 individual plants
   were sampled and the key traits, including spike length, number of
   spikelets, and 1,000-grain weight, were measured and analyzed
   statistically by using SPSS 20 software (Version 26.0,
   [69]https://spss.mairuan.com/).

Genomic DNA extraction and simple sequence repeat (SSR) analysis

   During the M3 generation planting in Chongzhou City, about 10 g of
   young leaf tissue from each mutation material was collected at the
   seedling stage and stored at −80 °C for genomic DNA extraction. For
   each material, about 5 g leaf tissue were ground in liquid nitrogen,
   and genomic DNA was subsequently extracted. The quality of the
   extracted DNA was assessed using a Nanodrop 2000 spectrophotometer
   (Thermo Fisher Scientific, Massachusetts, USA).

   By reviewing relevant literatures, 50 pairs of SSR primers with good
   polymorphism in the barley species were selected [[70]38], with
   sequences listed in Table [71]S1. Using the genomic DNA from the
   highland barley mutation materials as the template, PCR reactions were
   carried out with a total volume of 20 µL. The reaction mixture
   contained 50 ng genomic DNA, 2×PCR Mix, and 1 µM each of the forward
   and reverse primers. Electrophoretic band detection was then performed,
   and the electrophoretic bands of the 50 pairs of SSR primers detected
   in each mutation material were statistically analyzed, categorized into
   a presence (1) or absence (0) matrix. Finally, dendrogram construction
   was carried out using the NTSYS 2.1 software [[72]39] to assess
   molecular differences at the genomic level among the different mutation
   materials.

RNA extraction, library construction, sequencing, and quality control

   Spike tissue samples from M3 generation mutant materials (D7-67, E8-38)
   and the control K14 were randomly collected at early heading stage from
   high-altitude (Luhuo) and low-altitude (Chongzhou) sites. Three
   biological replicates were sampled per plant material, flash-frozen in
   liquid nitrogen, and stored at −80 °C. Total RNA was extracted using
   TRIzol reagent according to the manufacturer’s instructions. The RNA
   concentration, integrity and purity were assessed using Nanodrop 2000,
   Qubity 2.0, and Agilent 2100, respectively. High-quality RNA samples
   were sent to Beijing Novogene Bioinformatics Technology Co., Ltd. for
   library construction and sequencing on the Illumina platform. Raw reads
   were processed using Trimmomatic [[73]40] to filter out low-quality
   sequences, and quality was assessed with FastQC [[74]41]. After
   removing adapters and primers, the clean reads were prepared for
   further analysis.

   Clean reads were mapped to the reference genome of Hordeum vulgare
   “Morex” [[75]42] using Bowtie2 [[76]43], with alignment parameters set
   to “-N 0 -L 20” (-N < int > specifies the allowed number of mismatches
   during alignment; -L specifies the length of the alignment). Upon
   completion of the alignment, SAM format files were generated, and the
   mapping statistics for each sample were subsequently analyzed.

Functional annotation and differentially expressed genes (DEGs) analysis

   After mapping the clean reads of the transcriptome samples to the
   reference genome to obtain BAM files, the SAM-to-BAM format conversion
   was performed using SAMtools [[77]44], followed by sorting and index
   file creation for the BAM files. StringTie [[78]45] was then used to
   calculate the TPM (Transcripts per Kilobase Million) value for each
   gene in the transcriptome samples, which was employed to quantify gene
   expression levels. The calculation of TPM values allows for the
   normalization of gene length and differences in total sequencing reads
   across samples. To assess the intra-sample reproducibility of the
   transcriptome data, we analyzed the expression levels of each gene from
   the resulting TPM files using the Pspikeson correlation method from the
   cor function package in R. The pairwise correlation coefficients
   between samples were calculated to evaluate the consistency of the
   transcriptome data within each group. Additionally, the prepDE.py
   script provided by StringTie was used to extract the reads count values
   for each sample from the GTF files generated during the expression
   quantification process. These read counts values were subsequently
   normalized using the edgeR package [[79]46] in R.

   Genes demonstrating an absolute |log2FoldChange| >1 and padj < 0.05
   were categorized as differentially expressed genes (DEGs). The
   identified DEGs were then visualized using the Heatmapper tool
   [[80]47]. To annotate the sequences of DEGs, the blastp function of
   Diamond [[81]48] was employed to align the gene sequences with the NCBI
   Nr (Non-redundant) protein database. The alignment results were
   imported into Blast2GO [[82]49] for functional annotation using the GO
   database [[83]50]. The DEGs were categorized based on their molecular
   functions, cellular components, and biological processes. Furthermore,
   KEGG [[84]51] pathway analysis was conducted using the KEGG Automatic
   Annotation Server (KAAS) [[85]52], where functional annotation of
   metabolic pathways for the DEGs was performed via the bi-directional
   best hit (BBH) method.

Single-nucleotide polymorphism (SNP) analysis

   The transcriptomic data preprocessing began with sorting the BAM files
   using SAMtools [[86]44] (version 1.4.1, default parameters) and merging
   the triplicate BAM files from the same sample. Next, the MarkDuplicates
   module of Picard software was employed to detect and mark duplicate
   reads, thereby minimizing the impact of sequencing-induced duplicates
   on subsequent SNP detection. SNP calling was then performed using the
   HaplotypeCaller tool from GATK (The Genome Analysis Toolkit) [[87]53].
   Following SNP calling, the VariantFiltration tool in GATK was applied
   to filter the SNP sites based on the conditions “QUAL < 20 || DP < 20||
   MQ < 2 0.0”, retaining only those SNPs with a distance greater than 5
   base pairs, ensuring high-quality, reliable variant data. Subsequently,
   SNP sites were annotated to the “Morex” genome using SnpEff [[88]54],
   generating the filtered_variants.vcf file for SNP quantity and type
   statistical analysis. Finally, bcftools isec [[89]55] was used to
   compare the VCF files of two samples, extracting shared and specific
   SNPs along with gene IDs for functional annotation analysis.

Protein-protein interactions (PPI) network analysis of the selected DEGs

   First, the protein sequences of the selected DEGs were extracted using
   TBtools software (version 2.083) [[90]56]. The PPI, including both
   known and predicted interactions, were then identified using the
   Hordeum vulgare reference genome file from the STRING database
   ([91]https://www.string-db.org/). Finally, Cytoscape software (version
   3.8.2) [[92]57] was utilized to further download and visualize the
   interaction network.

qRT-PCR confirmation of RNA sequencing data

   For qPCR validation of VRS1 expression, we utilized the same highland
   barley cDNA samples employed in the RNA-seq analysis with the specific
   primers (Table [93]S1). Quantitative PCR reactions were performed using
   Hieff™ qPCR SYBR Green Master Mix (Yeasen Biotechnology, Shanghai,
   China) with the following thermal cycling protocol: initial
   denaturation at 95 ℃ for 2 min, followed by 40 cycles of 95 ℃ for 10 s
   and 60 ℃ for 30 s. All reactions were conducted in triplicate, and
   relative gene expression levels were calculated using the 2^−ΔΔCT
   method. The housekeeping gene HvActin was amplified using specific
   primers (Table [94]S1) as an internal control for data normalization
   [[95]58].

Results

Statistical comparison of phenotypic variations in the Highland barley
mutation lines

   The highland barley cultivar K14 was selected as parent material for
   mutagenesis by nitrogen ion beam implanting. Six treatments were
   applied to the dried seeds of K14 and the effect on the seed
   germination rate was initially determined by laboratory experiment
   (Table [96]1). It was shown the seed germination rate decreased along
   increasing dosages. The treated seeds were first sowed at highland
   location (Luhuo county) to generate the M1 mutation materials. However,
   phenotypic variation was not observed significantly among the M1
   materials. Therefore, more than 100 single spikes were randomly
   selected from each treatment, and then sowed at Chongzhou (lowland) in
   a row-fashion to generate M2 mutation materials. Phenotypic variation
   such as plant height, spike morphology, etc. was observed carefully
   throughout the entire growth period. Finally, 71 mutation materials
   were selected out from M2, each of which was divided into two portions
   to be sowed at two locations of highland and lowland, respectively.

Table 1.

   Intensity and dose of nitrogen ion beam used to treat seeds of K14
   Treatment Ion beam intensity (keV) Dosage
   (N^+/cm²) Seed germination rate (%) (n = 50)
   T1 15 2 × 10^16 95.3
   T2 15 5 × 10^16 95.6
   T3 15 1 × 10^17 81.3
   T4 30 2 × 10^16 95.3
   T5 30 5 × 10^16 92.2
   T6 30 1 × 10^17 62.0
   [97]Open in a new tab

   Yield is a critical economic trait in crops. In this study, the
   yield-related traits of M3 from highland and lowland were measured for
   the 71 mutation materials (Table [98]S2). Statistical analysis
   indicated that 20, 26, and 37 mutation materials showed significant
   differences compared to the control parent in the number of spikes,
   spike length, and 1000-grain weight, respectively. Among these, 33
   mutation materials, including D1-21, D10-58, and E8-38, exhibited a
   significant increase in 1000-grain weight, while only E3-90 showed a
   notable increase in spike length. No significant increase in the number
   of spikes was observed compared to the control. It was noticeable that
   the mutant E8-38 exhibited the highest 1000-grain weight among all
   mutation materials, significantly surpassing that of K14 by approximate
   42% and 47% at high and low altitudes, respectively. This variant holds
   great potential for further breeding as a high-yield cultivar.
   Interested, D7-67 exhibited only the development and grain filling of
   the two-rowed seeds in the spike, while the remaining four rows were
   sterile.

   Furthermore, the average values, standard deviations, variation ranges,
   and coefficients of variation (CV) for the yield-related traits in the
   M3 generation mutation lines were calculated (Table [99]2). Among the
   different populations in the M3 generation from two regions, the CV for
   1000-grain weight was the most stable, remaining around 18%. In
   contrast, the CV for single-spike yield was the largest, exceeding 20%,
   indicating that this trait may be the most sensitive to nitrogen ion
   beam-induced mutation treatments. For 1000-grain weight, 69.01% of the
   materials had higher mean values than the control parent; for
   single-spike weight, 50.7% of the materials had higher mean values
   compared to the parent; for spike length, 26.76% of the materials
   exceeded the control; while for spike fertility and number of spikes
   per row, only 14.08% and 12.68% of the materials showed values greater
   than the parent, respectively. These results suggest that nitrogen ion
   beam had a significant mutagenic effect on yield-related traits in
   highland barley, with the mutant population exhibiting considerable
   variation in these traits, which is advantageous for selecting positive
   mutations. Notably, 1000-grain weight and single-spike yield traits
   exhibited more positive mutations within the population.

Table 2.

   Data statistics of Highland barley mutation population M3 (n = 72)
   Generation and region Trait Mean Standard deviation Range of variation
   Coefficient of variation (CV)
   M3/Lowland Spike length 11.36 2.02 6.55–23.23 0.18
   Number of rounds 15.79 1.52 12.50–19.00 0.10
   1000-grain weight 40.71 6.64 21.66–66.73 0.16
   Grain number 93.85 11.68 32.58–114 0.12
   Single spike yield 3.84 0.78 0.71–5.72 0.20
   M3/Highland Spike length 11.40 1.63 6.91–15.22 0.14
   Number of rounds 13.27 1.80 9.40-16.45 0.14
   1000-grain weight 37.14 6.55 17.74–62.63 0.18
   Grain number 78.85 12.33 28.57–98.73 0.16
   Single spike yield 2.93 0.71 1.40–4.29 0.24
   [100]Open in a new tab

SSR genetic diversity analysis of Highland barley mutant materials

   In the M3 generation, a variety of phenotypic variations and
   significant differences in yield-related traits were observed,
   suggesting that these variations may be the result of genetic mutation.
   To investigate the genetic differences at the DNA level, SSR analysis
   was conducted on these materials. For the M3 generation materials
   planted in Chongzhou, young leaves were collected and total DNA was
   extracted. A total of 50 pairs of SSR primers, known for their
   polymorphism in highland barley (Table [101]S1), were selected to
   amplify DNA from 72 genomic samples and to analyze the electrophoretic
   banding patterns. The results showed that 38 out of the 50 SSR primers
   could distinguish the variant materials, indicating DNA polymorphism.
   Among these, the primer EBmag0793 exhibited the most significant
   polymorphism across different variant materials, amplifying 8
   polymorphic loci with 7 distinct banding patterns (Fig. [102]S1).
   Statistical analysis and clustering of the electrophoretic bands
   obtained from the 38 SSR primer pairs revealed that the 71 mtataion
   materials and the control (K14) could be grouped into 8 clusters at a
   distance coefficient of 0.89 (at which point the inter-group
   differences became significantly larger) (Fig. [103]1). Four clusters
   included only single materials (D3-72, E1-70, E1-102, or E3-90), and
   two clusters contained two members each (E2-43 and E2-97, or E3-81 and
   E3-79). The remaining 63 mutation materials and the parent were
   assigned to clusters I and VI. Cluster I included 57 materials, such as
   K14, D7-67, and E8-38, with most materials being distinguishable from
   one another. However, nine mutation materials, including D8-107, D9-19,
   and E6-56, despite differences in plant morphology and yield-related
   traits, clustered together. Analysis of the grouping of mutation
   materials in responding to the ion implantation dosages showed that the
   mutation materials selected from the T1, T2, T4, and T5 experimental
   groups were all grouped in cluster I. The seven mutation materials in
   cluster VI were all derived from seeds of the T3 experimental group,
   and except for D3-27, the remaining six materials exhibited late
   heading. Furthermore, materials from the T6 group were found in all
   seven clusters outside of cluster VI, suggesting that higher energy and
   dosage of nitrogen ion implantation (T6) may be more favorable for
   producing mutation in highland barley, resulting in a broader spectrum
   of variation.

Fig. 1.

   [104]Fig. 1
   [105]Open in a new tab

   Dendrogram based on SSR polymorphism in highland barley mutation
   population. I to VIII represent the eight clustered groups

Overview of comparative transcriptome analysis of the heading stage among
three mutation materials

   Highland barley “K14” seeds were treated with nitrogen ion beam
   implantation, and after three generations of planting and selection,
   several materials with distinct phenotypic variations were obtained,
   demonstrating both scientific value and significant economic potential.
   Among them, the 1000-grain weight of E8-38 was the most outstanding
   among all the mutants (Fig. [106]2A, Table [107]S2); while D7-67
   exhibited only the development and grain filling of the two-rowed seeds
   in the spike, whereas the remaining four rows were sterile
   (Fig. [108]2A). In order to explore the genetic basis of the above two
   mutation materials, comparative transcriptomics analysis was performed.

Fig. 2.

   [109]Fig. 2
   [110]Open in a new tab

   Phenotypic characteristics of the highland barley mutation materials
   and statistical analysis of RNA-seq DEGs. Mutation materials of E8-38
   and D7-67 (A). Statistics on the number of DEGs between two-two sample
   groups, the latter is a control group (B). Comparative Venn diagram
   between high and low altitude mutation materials and their parent K14
   (C). A Venn diagram comparing the effect of different altitudes on the
   same barley mutation materials (D)

   A total of 18 RNA samples were extracted from spike tissues of the
   mutation materials D7-67 and E8-38, as well as the parent “K14” at the
   early heading stage growing at different high altitudes. Sequencing was
   performed using the Illumina Hiseq PE150 platform, generating
   780,758,278 reads and 117.31 GB of clean data. The Q20 value for each
   sample exceeded 96%, and the clean data aligned to the “Morex”
   reference genome with a mapping rate of 94% (Table S3), confirming the
   reliability of the transcriptome sequencing results. Due to variability
   in field sampling conditions, some reproducibility issues were
   observed. Correlation heatmaps (Fig. [111]S2) showed that the
   correlation coefficient between samples within each altitude group was
   greater than 0.9, indicating high reproducibility and suitability for
   subsequent analysis. Additionally, PCA revealed potential differences
   between the three varieties at different altitudes (Fig. S3).

   DEGs analysis was performed to investigate the molecular mechanisms
   underlying the variation in phenotypes and the overall expression
   changes of three materials under different altitudes. To provide a
   visual representation of all DEGs across the groups, a volcano plot was
   constructed (Fig. S4). The number of DEGs across various groups and
   altitudes was summarized (Fig. [112]2B). Interestingly, the number of
   DEGs was in general higher between different high altitude across the
   three materials. For the parent K14, H_K14 showed 5,486 DEGs (2,245
   upregulated and 3,241 downregulated) compared to L_K14; H_E38 showing
   8,651 DEGs (4,406 upregulated and 4,245 downregulated) compared to
   L_E38; and H_D76 showing 7,756 DEGs (3,716 upregulated and 4,040
   downregulated) compared to L_D76. These data illustrated that that
   altitude significantly affects gene expression levels in highland
   barley. In addition, the E8-38 variant showed a relatively higher
   number of DEGs in comparison to D7-67. At highland, when comparing the
   two mutants to their parent, H_E38 showed upregulation of 1,123 genes
   and downregulation of 990 genes compared to H_K14, while H_D76
   exhibited 834 genes upregulated and 591 genes downregulated compared to
   H_K14. At lowland, L_E38 showed upregulation of 2,418 genes and
   downregulation of 2,696 genes compared to L_K14, and L_D76 exhibited
   1,071 genes upregulated and 994 genes downregulated compared to L_K14.
   Conclusively, two Venn diagrams were used to identify the relationships
   between the changing genes. Among the four/two comparisons between the
   high and low altitude and various plant lines, 167 genes were found to
   be commonly involved in trait changes (Fig. [113]2C). In the
   comparisons of altitude effects, 3,518 genes were found to have altered
   expression levels (Fig. [114]2D). These findings demonstrate that
   altitude significantly influences gene expression as well as the growth
   in highland barley.

Genetic basis of phenotypic variations in highland barley introduced by
nitrogen ion implantation

   RNA-seq data from the D7-67 and E8-38, obtained from high and low
   altitude locations, were compared pairwise with the parent “K14” to
   analyze the genetic underpinnings of their phenotypic variations. The
   D7-67 and E8-38 mutants were analyzed from different perspectives due
   to their distinct traits. The E8-38 mutant, characterized by its higher
   1000-grain weight, was analyzed in terms of metabolic pathways; while
   the D7-67 mutant, with its unique two-row phenotype, was examined
   through SNP variations and genes involved in row type regulation.

   First, a comprehensive GO annotation was performed, categorizing DEGs
   into three main categories: molecular function, biological process, and
   cellular component. The top 10 GO terms with the smallest p-values for
   each group were displayed (Fig. [115]3A). In the biological process
   category, the most enriched and largest number of DEGs across all four
   groups were related to “response to abiotic stimulus” and “response to
   oxygen-containing compound”. In the cellular component category, the
   most confidently enriched term in all four groups was “cell periphery”.
   Molecular function terms showed a lower level of enrichment. In
   comparison of high and low altitude E8-38 with the K14 parent, the GO
   terms commonly enriched included “response to abiotic stimulus”,
   “response to hormone”. Further analysis of DEGs within the “response to
   oxygen-containing compound” GO term revealed that high-altitude E8-38
   exhibited 69 upregulated and 39 downregulated genes, while low-altitude
   E8-38 showed 144 upregulated and 130 downregulated genes. Notably, 31
   genes were identified as responsive to oxygen-containing compounds in
   both high and low altitude E8-38 mutants, categorized into nine major
   functional groups (Fig. [116]3B and Table S4). Among these,
   Cytokinin-responsive GATA transcription factor 1-like, DNA
   (cytosine-5)-methyltransferase (DRM), Carbonic anhydrase, Chitin
   elicitor receptor kinase 1-like, Flavone O-methyltransferase 1-like,
   and iron ascorbate-dependent oxidoreductase (ANS) displayed significant
   upregulation. These genes are implicated in molecular processes such as
   cytokinin signaling, DNA modification and epigenetic regulation, grain
   development and cell expansion, immunity and defense, and metabolic
   regulation, all of which are highly active and critically influence
   seed development and weight, thereby contributing to the high
   1000-grain weight phenotype of E8-38. Additionally, GO enrichment
   analysis indicated a substantial number of genes involved in
   phytohormone response in E8-38. Given that plant hormones are known to
   regulate developmental processes and traits directly linked to
   productivity, and single mutations in phytohormone regulatory pathway
   genes can alter multiple phenotypes [[117]59], we further analyzed the
   expression patterns of hormone-responsive DEGs in E8-38 grown under
   high and low altitude conditions compared to K14 (Fig. [118]3C and
   Table S5). The heatmap revealed that cytokinin biosynthesis-related
   genes were the most upregulated, while genes associated with auxin and
   ethylene were downregulated. These differential expression patterns of
   phytohormones are likely key contributors to the elevated 1000-grain
   weight observed in E8-38.

Fig. 3.

   [119]Fig. 3
   [120]Open in a new tab

   GO enrichment analysis of DEGs associated with phenotypic variation of
   highland barley mutation materials (A). The expression levels of
   response to oxygen-containing compound related genes shared by the
   different altitude E8-38 (B). The expression levels of plant
   hormone-related genes shared by the different altitude E8-38 (C)

   KEGG pathway analysis was subsequently performed to identify enriched
   pathways (Fig. [121]4A). Common pathways across all comparison groups
   included “Plant hormone signal transduction”, “Starch and sucrose
   metabolism”, and “Phenylpropanoid biosynthesis”, suggesting their
   pivotal roles in high-altitude adaptation in these mutation materials.
   Notably, the “Starch and sucrose metabolism” and “Phenylpropanoid
   biosynthesis” pathways contained the largest number of DEGs. A heatmap
   of gene expression levels for shared genes in E8-38 and D7-67 under
   both high and low altitude conditions was generated. An expression
   value of 0 indicates the absence of the gene in that group
   (Fig. [122]4B). Overall, E8-38 and D7-67 exhibited widespread DEGs in
   both metabolic pathways. Intriguingly, in the H_D76vsH_K14 and
   L_D76vsL_K14 comparison groups, three genes in the “Starch and sucrose
   metabolism” pathway were significantly upregulated, indicating enhanced
   starch and sucrose metabolism in the D7-67. This observation implies
   that the two-rowed spike phenotype of D7-67 may be associated with
   enhanced starch and sucrose accumulation, although further experiments
   are needed to validate these findings to confirm the actual increase in
   nutrient content.

Fig. 4.

   [123]Fig. 4
   [124]Open in a new tab

   KEGG enrichment analysis of DEGs associated with phenotypic variation
   induced by highland barley mutations (A). Heatmap of DEGs expression in
   the pathway with the highest number of KEGG in E8-38 and D7-67 (B)

   Using SSR analysis, we confirmed genomic-level sequence differences
   between the mutation lines and the control parent, induced by nitrogen
   ion beam implantation. To further investigate the characteristics of
   gene mutations induced by nitrogen ion beams, we aligned the
   transcriptome data of the H_K14 parent, H_D76, and H_E38 to the
   reference genome, generating SAM files. SNP counts were then analyzed
   using GATK, SnpEff, and bcftools. Initially, the number of SNPs was
   statistically analyzed with H_K14 as the control (Table [125]3). A
   total of 45,707 SNP sites were identified in the comparison between
   H_E38 and H_K14, while 35,850 SNP sites were found in the comparison
   between H_D76 and H_K14. This difference may explain the higher number
   of DEGs identified in E8-38 compared to D7-67. Overall, these gene
   mutations were distributed across all chromosomes, although there were
   differences in mutation distribution between chromosomes. Specifically,
   the highest number of mutations were observed on the 4 H chromosome,
   while the fewest mutations occurred on the 2 H chromosome. It is
   evident that ionizing radiation easily induces nucleotide variations in
   chromosomes, which may also result from the DNA repair of damage caused
   by ionizing radiation [[126]60].

Table 3.

   Statistics on the numbers of SNPs in each sample
   Sample H_E38 vs. H_K14 H_D76 vs. H_K14
   No. of variants 54,207 42,724
   Total SNPs 45,707 35,850
   INS 4,445 3,670
   DEL 4,055 3,204
   Genetic mutation counts in each chromosome (Variants) Chromosome Size
   (bp) SNP density
   1 H 6,236 4,383 516,505,932 105,774
   2 H 9,635 9,168 665,585,731 70,795
   3 H 8,409 6,450 621,516,506 83,655
   4 H 5,689 4,750 610,333,535 116,933
   5 H 8,170 7,296 588,218,686 102,602
   6 H 8,180 6,321 561,794,515 77,483
   7 H 7,831 4,252 632,540,561 104,699
   Un 57 104 2,974,422 36,495
   [127]Open in a new tab

   Next, we extracted protein-coding genes with SNP sites and compared
   them to DEGs identified in the parent control group. In the comparison
   between H_E38 and H_K14, 2,113 DEGs were identified, with 704 genes
   containing SNP sites. In the comparison between H_D76 and H_K14, 643
   DEGs were associated with SNP sites. These SNP-containing DEGs were
   further annotated using GO terms. The results showed that, within the
   biological process category, genes in the E8-38 were predominantly
   enriched in photosynthesis-related processes (Fig. [128]5A). In
   contrast, genes in the D7-67 mutant were primarily enriched in lipid
   metabolism and fatty acid metabolism, both closely associated with
   nutrient metabolism (Fig. [129]5B). At the molecular function level,
   the E8-38 was mainly enriched in hydrolase activity and monooxygenase
   activity, while the D7-67 was enriched in water channel activity and
   passive transmembrane transporter activity.

Fig. 5.

   [130]Fig. 5
   [131]Open in a new tab

   Annotated categorization of DEGs identified as having SNP sites present
   at the level of cellular composition, molecular function, and
   biological function (A and B). Detailed genetic variations of the
   row-type controlling gene in H_D76 are depicted. The notation indicates
   the location of the SNP, while the box represents the full-length
   structure of the gene, with distinct regions of the sequence
   highlighted in different colors (C). The expression levels of the genes
   perhaps controlling row-type in H_D76 (D). The qRT-PCR validates the
   expression levels of genes related to control row types in H_D76 and
   L_D76, * indicate significant differences (E)

   When any one of the five genes associated with the six-row spike (VRS
   genes, VRS1-VRS5) undergoes a loss-of-function mutation, the fertility
   of the lateral spikelets is significantly enhanced. This leads to a
   complete or partial restoration of seed-setting ability, and
   consequently, the formation of multi-row spikes [[132]61]. Notably, 38
   genes have been identified in K14 that are implicated in the regulation
   of row-type formation, including 6 VRS1 (Homeobox-leucine zipper
   protein), 29 VRS4 (Lateral organ boundaries), and 3 VRS5 (Teosinte
   branched 1) genes. Through SNP analysis, we identified two VRS1 genes
   (HORVU.MOREX.r3.4HG0394590 and HORVU.MOREX.r3.2HG0198150) that exhibit
   genetic variations and differential expression in D7-67
   (Fig. [133]5C-D). Specifically, VRS1_1 harbors 5 SNP loci distributed
   across different regions of the gene in chromosome 4H, with upregulated
   expression observed in H_D76. Among these, the mutation at locus
   520509790 is an insertion variant (T→TCTCAA) located in the 5’ UTR
   region, which may influence transcriptional regulation. The loci
   520510084 and 520511314 are synonymous mutations (A→T and G→A) at
   positions 72 and 672 of the coding region, respectively, altering
   codons without changing the amino acid sequence. Most of these SNP loci
   have high-quality scores (QD > 28), indicating high reliability. In
   contrast, VRS1_2 on chromosome 2H contains 33 SNP loci, of which 24 are
   classified as low-impact variants, 1 is located in the 3’ UTR region,
   and 3 are situated in the 5’ UTR region, potentially affecting gene
   transcription, splicing, or protein function. It is noteworthy that we
   detected a proximal SNP at position 620318407 (663 bp upstream of 5’
   end (620319070) at VRS1_2). A distal 3’ SNP was found at position
   620379239 (53,961 bp downstream). While the 3’ SNP is too distant for
   functional impact, the 5’ proximal SNP (663 bp) may potentially
   regulate VRS1_2 expression. Notably, this 5’ SNP is located within an
   exon of HORVU.MOREX.r3.2HG0198140, though this gene did not show
   differential expression in the H_D76 vs H_K14 comparison. The 663 bp
   proximity of the VRS1_2-associated SNP warrants further investigation
   as a potential regulatory element. The functional implications would
   benefit from experimental validation in future studies. These
   variations may collectively contribute to the phenotypic alterations
   observed in the H_D76. The upregulated expression of VRS1_1 is likely a
   primary factor driving the two-rowed phenotype in H_D76, highlighting
   the critical role of VRS1 in row-type determination and its potential
   as a key genetic target for highland barley breeding programs.

Functional analysis of DEGs involved in response to high altitude adaptation

   Highland barley is a crop adapted to high-altitude environments.
   Transcriptome analysis revealed that the number of significantly DEGs
   was higher across three comparison groups from different altitudes for
   the same material. To further investigate the potential biological
   significance of these DEGs, KEGG pathway enrichment analysis identified
   several significantly enriched pathways across all three highland
   barley materials, indicating shared biological responses to altitude
   conditions (Fig. [134]6A). The most prominent pathways included
   phenylpropanoid biosynthesis, starch and sucrose metabolism, and lipid
   biosynthesis. Additionally, several pathways associated with
   high-altitude adaptation, such as carbon fixation in photosynthetic
   organisms and flavonoid biosynthesis, were also enriched. These
   pathways may be linked to adaptation to UV radiation and high-altitude
   conditions. Interestingly, genes associated with photosynthesis,
   particularly those involved in carbon fixation and porphyrin
   chlorophyll metabolism, were upregulated under high-altitude conditions
   (Fig. [135]6B), reflecting a metabolic adaptation to high-altitude
   environments. The phenylpropanoid metabolic pathway is central to the
   biosynthesis of polyphenols, flavonoids, and other compounds, which
   function as antioxidants to mitigate stress in plants [[136]59]. This
   metabolic pathway initiates with phenylalanine, which is catalyzed by
   phenylalanine ammonia-lyase (PAL) to yield trans-cinnamic acid.
   Subsequently, cinnamic acid 4-hydroxylase (C4H) catalyzes the
   conversion of trans-cinnamic acid into compounds such as coumaric acid,
   ferulic acid, sinapic acid, and caffeic acid [[137]62]. Finally,
   enzymes such as 4-coumarate-CoA ligase (4CL) are involved in
   synthesizing phenylpropanoid metabolites, including flavonoids and
   polyphenols. As demonstrated in Fig. [138]6C and Fig. S5, over half of
   the genes associated with the phenylpropanoid pathway were
   significantly upregulated, particularly, PAL and 4CL showing pronounced
   increases. Given that the initial three catalytic reactions are central
   to this pathway [[139]63], these results strongly suggest enhanced
   phenylpropanoid metabolism in highland barley under high-altitude
   conditions. This metabolic shift is likely able to enhances the plant
   to synthesize more polyphenolic compounds, reducing intracellular free
   oxygen levels and improving stress resistance. Further enrichment was
   observed in pathways related to lipid biosynthesis, carotenoid
   biosynthesis, and linoleic acid metabolism. These findings suggest that
   high-altitude conditions promote the accumulation of beneficial
   compounds, such as pigments and unsaturated fatty acids, in highland
   barley. These compounds are known to contribute to plant stress
   resistance, indicating an adaptive mechanism for survival at higher
   elevations [[140]64]. These insights provide a valuable theoretical
   foundation for developing stress-resistant plant varieties.

Fig. 6.

   [141]Fig. 6
   [142]Open in a new tab

   KEGG functional enrichment analysis of DEGs involved in the
   high-altitude response of highland barley. A Summary of KEGG enrichment
   for three sets of sample DEGs. B In the three comparison groups of
   different highland barley species from varying altitudes, gene
   expression levels related to two photosynthesis-related processes and
   phenylpropanoid biosynthesis were analyzed. C Biosynthetic pathway of
   phenylpropanoid. The colored bar shows a heat map of the expression of
   key enzymes for the synthesis of compounds between the three groups.
   Key enzyme abbreviations: PAL, phenylalanine ammonia-lyase; C4H,
   cinnamic acid 4-hydroxylase; 4CL, 4-coumarate CoA ligase

   GO enrichment analysis also revealed that altitude elicited similar
   expression trends across different plant materials (Fig. [143]7A). In
   the Molecular Function category, genes were predominantly associated
   with transcription regulator activity, while the Biological Process
   category exhibited significant enrichment in “response to abiotic
   stimulus” (lowest p-value). Next, DEGs involved in transcription
   regulation and response to abiotic stimulus, common to all the three
   comparison groups, were extracted for PPI network analysis
   (Fig. [144]7B and Table S6). As a result, 67 interactions with scores
   exceeding 700 were revealed. The top three enriched functions in the
   network were: disaccharide metabolic process, response to light
   stimulus, and response to radiation. We further analyzed the
   interaction results using k-means clustering. K-means clustering of the
   PPI network identified a subnetwork of nine genes collectively involved
   in response to light stimulus (Fig. [145]7C and Table S7). These
   included three BZIP transcription factors (HORVU.MOREX.r3.5HG0433250,
   HORVU.MOREX.r3.6HG0577940, HORVU.MOREX.r3.7HG0713820). CPRF2, a BZIP
   transcription factor, was designated as BZIP1. In addition, two other
   BZIP transcription factors, BZIP2 and BZIP3, belong to the ELONGATED
   HYPOCOTYL 5 type. Expression data indicated that these BZIP factors
   positively regulate light responses during high-altitude growth. The
   PPI network highlighted key interactions: BZIP1 strongly interacted
   with BZIP2 (score: 0.963), suggesting heterodimer formation for
   enhanced regulatory capacity. BZIP1 also interacted with
   serine/threonine-protein kinase STN7 (score: 0.854), potentially
   linking chloroplast light signaling to nuclear gene phosphorylation
   modification regulation [[146]65]. Additionally, BZIP2 and BZIP3
   interacted with E3 ubiquitin-protein ligase COP1 (score: 0.836),
   implicating ubiquitin-mediated regulation of BZIP stability. Notably,
   BZIP2 and BZIP3 also interacted with WD repeat-containing protein
   RUP2-like (score: 0.914 and 0.966), a known modulator of
   photomorphogenesis [[147]66], further supporting their role in light
   signaling. Interestingly, cyclin-dependent kinase 5 (Cdk5) and
   mitogen-activated protein kinase 9-like, both downregulated, showed
   indirect associations with STN7 and BZIP1. These kinases may modulate
   the phosphorylation status of BZIP transcription factors, thereby
   influencing the dynamics of light signal transduction. These findings
   highlight the central role of BZIP transcription factors in light
   signaling and abiotic stress responses, providing mechanistic insights
   into high-altitude adaptation.

Fig. 7.

   [148]Fig. 7
   [149]Open in a new tab

   GO functional enrichment analysis of DEGs involved in the high-altitude
   response of highland barley. A Summary of GO enrichment results for
   three sets of sample DEGs. B PPI network diagram of transcription
   factors and response to abiotic stimulus shared by the three mutation
   materials (score > 700). Lines represent interactions, with red and
   blue indicating upregulated and downregulated expression, respectively.
   Circular nodes represent proteins related to response to abiotic
   stimulus, while diamond-shaped nodes represent transcription factors. C
   K-means clustering of PPI networks reveals genes jointly involved in
   the response to light stimulus. Red indicates up-regulation, blue is
   down-regulation. The size of the dots corresponds to the number of
   interacting proteins, with larger dots representing a greater number of
   interactions. In addition, the values adjacent to each circle on the
   graph denote the log[2]FC for the various comparison groups

Discussion

   Highland barley is a staple crop of the Qinghai-Tibet Plateau in China,
   known for its strong adaptability, stable yield, and ease of
   cultivation. As a model organism for studying adaptive evolution in
   plants under extreme environmental conditions, highland barley serves
   as an exemplary species for understanding plant resilience in
   high-altitude ecosystems [[150]67]. Over the years, breeding
   advancements have led to the development of numerous high-quality,
   high-yielding highland barley varieties, contributing significantly to
   increase per-unit yields. In this study, we employed ion beam
   implantation technology to induce mutations in the highland barley
   variety K14. This mutagenesis approach aims to expand the genetic base
   of highland barley, generating novel genetic resources for future
   breeding efforts. We conducted transcriptome sequencing on selected
   mutation lines exhibiting distinct phenotypic variations, such as the
   high 1000-grain weight mutant E8-38 and the dual-row trait mutant
   D7-67. Functional annotations, including GO and KEGG pathway
   enrichment, SNP analysis, and PPI network analysis, were performed.
   These analyses aimed to elucidate the molecular mechanisms underlying
   the observed phenotypic variations and their potential roles in
   high-altitude adaptation.

Ion beam implantation mutagenesis and screening of highland barley mutations

   Ion beam implantation mutagenesis has been widely applied in plant
   breeding since its introduction in rice in 1986. Ion beams exhibit a
   broader mutation spectrum than gamma rays in practical mutation
   breeding [[151]68]. The energy and dosage of the ion beam are critical
   factors influencing mutagenic effects, with low-dose radiation
   generally promoting growth, while high doses inhibit plant development
   [[152]69]. Zhao et al. (2018) irradiated rice with carbon ion beams at
   different dosages and found that low-dose heavy ion beams primarily
   induced hypermethylation at CG sites, which may activate the metabolic
   and biosynthetic mechanisms of the organism. In contrast, high-dose
   heavy ion beams induced hypomethylation at CG sites, potentially
   enhancing the plant’s resistance to stress [[153]70]. In this study,
   highland barley seeds were treated with nitrogen ion beams at two
   energy levels (15 KeV and 30 KeV) and three doses (2 × 10^16 N^+/cm²,
   5 × 10^16 N^+/cm², and 1 × 10^17 N^+/cm²). As expected, it was observed
   that as the ion implantation dose increased, seed germination rate
   decreased (Table [154]1). In the experimental group treated with a 30
   KeV ion beam at a dose of 1 × 10^17 N^+/cm², the germination rate was
   lowest at 62%, which did not reach the half-lethal dose, and the damage
   remained relatively low. Subsequently, six groups of mutagenized
   materials were subjected to multi-generation selection and SSR
   analysis, yielding 71 stable mutation lines with widespread genomic
   mutations. The experimental group treated with a 30 KeV ion beam at a
   dose of 1 × 10^17 N^+/cm² yielded the highest number of mutation
   materials and the most extensive range of variations.

   Low-energy nitrogen ion implantation has been shown to induce
   significant morphological variations in wheat, such as reduced plant
   height, awnless traits, and increased 1000-grain weight [[155]71].
   Statistical analysis of phenotypic variations and yield-related traits,
   such as spike length and 1000-grain weight, in the 71 mutation lines
   revealed broad variation in these traits. Compared to the parent, 37
   materials exhibited significant differences in 1000-grain weight, while
   20 and 26 materials showed significant differences in the number of
   spikes and spike length, respectively. Among these, the mutation
   material E8-38 exhibited a significant increase in 1000-grain weight
   and a 40% increase in single-spike yield compared to the parent,
   showing potential for further selection as a high-yield variety. UV-B
   radiation significantly inhibits highland barley growth, resulting in
   shorter plants, reduced green leaf area, and lower dry matter content.
   The inhibitory effect was more pronounced with stronger radiation, and
   the degree of inhibition varied with the growth stage [[156]72]. A
   comparative analysis of highland barley mutation populations grown at
   different altitudes revealed that those at higher altitudes generally
   exhibited poorer performance in spike length, spike number, and
   1000-grain weight compared to those at lower altitudes. This indicates
   that altitude is indeed an important environmental factor influencing
   highland barley growth and yield. Among these traits, the most
   significant differences were found in the number of spikes, suggesting
   that altitude has the most pronounced effect on this trial.

   SSR markers are reliable molecular tool for variety identification and
   genetic analysis [[157]73]. Therefore, we employed SSR analysis to
   examine 71 mutation lines, confirming genomic-level sequence
   differences between the mutations and the control parent induced by ion
   implantation. Secondly, by matching the experimental grouping,
   phenotypic variation, and yield traits of the materials, no clear
   association was found between the phenotypic variation and yield traits
   of the variant material samples and the clusters in the phylogenetic
   tree. However, when matching with the grouping of the variant materials
   during the ion beam implantation experiment, it was found that
   materials selected from four experimental groups (T1, T2, T4, T5) all
   clustered within the same major group I. In group VI, the variant
   materials were all derived from experimental group T3, while materials
   obtained from experimental group T6 showed the widest variation, being
   distributed across seven groups other than group VI. This phenomenon
   may be attributed to the fact that all mutation lines originated from
   the same parent, resulting in limited genomic differences.
   Additionally, the highland barley genome is extensive, and the 50 SSR
   primer pairs used cover only a small portion of it. Furthermore, by
   matching the phenotypic variations and yield-related traits of the
   mutation materials in group I, no significant associations were found.
   This suggests that higher energy and dose conditions of nitrogen ion
   implantation (experimental group T6) may be more conducive to highland
   barley mutagenesis, resulting in a broad range of mutations.

Influence of phytohormone signaling and metabolic regulation on the
1000-grain weight mutation in E8-38

   The mutation material E8-38 exhibits a significantly higher 1000-grain
   weight. Understanding the genetic and molecular mechanisms controlling
   grain weight is a prerequisite for improving grain yield potential
   [[158]74]. Grain weight is influenced by complex molecular and genetic
   factors that regulate processes such as cell division, proliferation,
   and differentiation. The remarkable increase in 1000-grain weight
   observed in the E8-38 is likely driven by a complex interplay between
   phytohormone signaling and metabolic regulation, as revealed by our
   transcriptomic and functional enrichment analyses. The significant
   upregulation of CKs-related genes (cytokinin-responsive gata
   transcription factor 1-like、B-ARR、AHP), coupled with the downregulation
   of auxin (auxin-responsive protein) and ethylene (ethylene receptor)
   pathways (Fig. [159]3B-C), suggests a finely tuned hormonal balance
   that promotes cell division, delays senescence, and enhances nutrient
   allocation to developing grains. CKs, known for their role in
   regulating sink strength and grain filling, were particularly prominent
   in E8-38, aligning with previous studies that highlight the importance
   of cytokinin signaling in improving grain yield [[160]59]. Moreover, it
   was shown that the cytokinin-responsive gata transcription
   factor1/gnc-like (CGA1/GNL) is involved in photosynthesis as well as
   amino acid and starch biosynthesis in rice [[161]75]. It is assumed
   that these changes also exist in E8-38. Ethylene receptor encodes a
   serine/threonine kinase whose overexpression results in reduced grain
   weight [[162]25]. The ethylene sensitivity is positively correlated
   with grain weight and grain size. It is also consistent with the
   findings of this paper. Auxin has also been reported to be involved in
   the control of seed development and grain yield [[163]76]. The
   suppression of auxin and ethylene pathways may reduce inhibitory
   effects on grain development, further contributing to the increased
   grain weight. The upregulation of cytokinin signaling and suppression
   of auxin and ethylene pathways in E8-38 likely results in pleiotropic
   effects that influence various growth and stress response traits.
   Conversely, the suppression of auxin signaling could alter root
   architecture and impact the plant’s ability to acquire water and
   nutrients, especially under stressful conditions such as high-altitude
   environments [[164]77]. Additionally, ethylene suppression may result
   in delayed senescence, improving photosynthetic duration and overall
   growth under shorter growing seasons or high stress [[165]78]. These
   changes in hormonal balance are likely to provide E8-38 with advantages
   in stress resilience, especially under harsh environments like high
   altitudes, where UV radiation, low oxygen levels, and limited resources
   prevail. However, under optimal growing conditions, these same hormonal
   modifications could result in unbalanced growth, possibly leading to
   excessive vegetative growth or water/nutrient limitations. Thus, the
   balance of these hormonal pathways also plays a crucial role in the
   overall fitness of E8-38 under different environmental conditions.

   The enrichment of genes associated with “response to oxygen-containing
   compounds” highlights the critical role of metabolic regulation in the
   high 1000-grain weight phenotype of E8-38 (Fig. [166]3B). Key genes
   such as carbonic anhydrase and ANS were significantly upregulated,
   suggesting enhanced photosynthetic efficiency and redox homeostasis in
   E8-38. Carbonic anhydrase, in particular, plays a critical role in
   CO[2]hydration and photosynthetic carbon assimilation [[167]79],
   processes essential for sustaining high 1000-grain weight. The
   upregulation of flavone O-methyltransferase 1-like further highlights
   the importance of secondary metabolite biosynthesis in supporting grain
   development and stress resilience. E8-38 may exhibit enhanced
   resistance to pathogens and pests, reducing the impact of biotic stress
   on grain development [[168]80]. These metabolic adaptations likely
   enable E8-38 to maintain high productivity under varying environmental
   conditions, particularly in high-altitude regions. The significant
   upregulation of DRM in E8-38 suggests a critical role for epigenetic
   modifications in shaping the high 1000-grain weight phenotype.
   DRM-mediated DNA methylation may regulate the expression of genes
   involved in cell division, grain development, and stress responses. For
   instance, DRM could promote the methylation and silencing of negative
   regulators of grain size, while activating genes that enhance cell
   expansion and nutrient accumulation in developing grains. Furthermore,
   DRM may modulate the expression of hormone-responsive genes to optimize
   hormonal balance during grain filling. This epigenetic regulation,
   combined with the mutant’s enhanced metabolic efficiency, likely
   underpins the robust grain development observed in E8-38.

   In conclusion, the high 1000-grain weight phenotype of E8-38 likely
   results from a synergistic interplay between phytohormone signaling,
   metabolic regulation, and epigenetic modifications. The role of DRM
   underscores the importance of epigenetic regulation in optimizing gene
   expression networks for grain development. These findings elucidate the
   molecular mechanisms underlying the high-yield phenotype of E8-38 and
   identify valuable targets for future crop improvement programs aimed at
   enhancing grain yield and stress tolerance in highland barley and
   related cereals.

SNP variations and the molecular basis of the two-row phenotype in D7-67

   Most of the reported mutants in the literatures exhibit changes in
   quality or traits, such as the single-tiller mutant and the
   non-flowering mutant in rice, with these studies primarily focusing on
   understanding the genes controlling target traits [[169]81]. In this
   study, we identified 45,707 SNP loci in the ion beam-induced H_E38
   mutant and 35,850 SNP loci in the H_D76 mutant. The distribution of
   SNPs was not uniform across chromosomes or throughout the genome, which
   is inconsistent with the characterization of gamma-irradiation-induced
   genetic mutations in rice [[170]60]. The reason for this may be that
   our study only focused on specific two typical mutants. The
   identification of two VRS1 genes (HORVU.MOREX.r3.4HG0394590 and
   HORVU.MOREX.r3.2HG0198150) with significant genetic variations and
   differential expression in D7-67 highlights their pivotal role in
   regulating spike architecture. VRS1_1, located on chromosome 4H,
   harbors five SNP loci, including an insertion variant in the 5’ UTR
   region that may influence transcriptional regulation. Although two
   mutations are synonymous and do not alter the amino acid sequence,
   their presence in the coding region may still affect mRNA stability or
   translation efficiency (Fig. [171]4C). The high-quality scores
   (QD > 28) of these SNP loci further confirm their reliability and
   potential functional significance. Furthermore, VRS1_2, located on
   chromosome 2H, contains 33 SNP loci, including several in the 5’ and 3’
   UTR regions that may impact gene transcription, splicing, or protein
   function. The upregulation of VRS1_1 in D7-67 is particularly
   noteworthy, as it aligns with the known role of VRS1 in suppressing
   lateral spikelet fertility, thereby promoting the two-row phenotype
   [[172]12, [173]82]. The functional enrichment analysis of
   SNP-containing DEGs in D7-67 further supports the importance of lipid
   metabolism and fatty acid metabolism in nutrient allocation and spike
   development. These metabolic processes are closely linked to the
   regulation of spikelet fertility and grain filling, which are critical
   for determining row type. The enrichment of water channel activity and
   passive transmembrane transporter activity at the molecular function
   level suggests enhanced nutrient and water transport capabilities in
   D7-67. These capabilities may contribute to its distinct phenotypic
   characteristics.

   Overall, the nitrogen ion beam-induced mutation materials E8-38 and
   D7-67 exhibit substantial genomic alterations, affecting multiple genes
   that contribute to their phenotypic divergence. While whole-genome
   resequencing is a preferred method for mutation identification and gene
   isolation in species with small or simple genomes [[174]83, [175]84],
   its application in large and complex genomes remains cost-prohibitive.
   In this study, we employed transcriptome sequencing to simultaneously
   detect genomic variations and differential gene expression between
   mutants and wild-type plants, providing critical insights into the
   genetic and biological underpinnings of the observed phenotypes.
   Although transcriptome sequencing may overlook intronic and non-coding
   regulatory variants (with relatively low probability), it serves as a
   cost-effective alternative to whole-genome sequencing for large-scale
   mutation screening. The two-rowed phenotype of D7-67 appears to result
   from the synergistic effects of SNP-induced allelic variation and VRS1
   upregulation, which together disrupt the regulatory mechanisms
   controlling lateral spikelet fertility. These findings not only
   elucidate the molecular basis of the two-rowed phenotype but also
   highlight VRS1 as a promising target for genetic manipulation in
   highland barley breeding programs. Future studies should prioritize
   functional validation of the identified SNP loci to define their
   specific roles in inflorescence architecture regulation. Additionally,
   the broader implications of these mutations for enhancing yield and
   environmental adaptability in highland barley should be explored.

Molecular basis of high-altitude adaptation in highland barley

   Highland barley is a crucial crop in the high-altitude regions of
   China, known for its rich nutritional profile and strong resilience,
   enabling normal growth under adverse conditions [[176]67]. However,
   current research on highland barley stress tolerance has primarily
   focused on metabolite content analysis [[177]85], gene mining and
   validation through genomics [[178]3], and simulating adverse
   environments via salt stress [[179]59]. Studies analyzing highland
   barley growth under real high-altitude conditions are still lacking.
   Our study involved cultivating three highland barley varieties at both
   high and low altitudes. The transcriptomic and network analyses
   presented here reveal a multifaceted molecular strategy underpinning
   highland barley’s adaptation to high-altitude environments.
   Polyphenolic compounds, including flavonoids and lignin, have been
   recognized as critical secondary metabolites in plants, playing pivotal
   roles in stress tolerance by scavenging ROS to enhance UV resistance
   and mitigate oxidative damage under adverse environmental conditions
   [[180]86]. Based on our results, genes associated with photosynthesis
   and phenylpropanoid biosynthesis were significantly co-enriched and
   predominantly upregulated in three highland barley varieties under
   high-altitude conditions (Fig. [181]6A-B). This coordinated
   upregulation likely enhances highland barley’s ability to survive
   extreme abiotic stresses, including intense UV radiation, fluctuating
   temperatures, and oxidative stress. Notably, the phenylpropanoid
   biosynthesis pathway serves as a cornerstone of this adaptive response
   [[182]59]. The first three catalytic reactions in the phenylpropanoid
   metabolic pathway have been identified as the core steps driving this
   process [[183]63]. Among these, PAL, the initial enzyme in the pathway,
   catalyzes the rate-limiting conversion of phenylalanine to cinnamic
   acid, a precursor for downstream secondary metabolite synthesis
   [[184]87]. Under high-altitude conditions, three PAL genes exhibited
   significant transcriptional upregulation across the three barley
   mutation materials, with HORVU.MOREX.r3.2HG0182090 showing the most
   pronounced increase in K14 (Log2FoldChange = 3.303) (Fig. [185]5C).
   4CL, the final catalytic enzyme in the three core steps of the
   phenylpropanoid pathway, also serves as a major branch point enzyme
   that generates activated thioesters [[186]88]. Moreover, 4CL is
   valuable in metabolic pathway engineering for curcuminoids,
   resveratrol, biofuel production and nutritional improvement. Our study
   revealed significant upregulation of two 4CL genes in high-altitude
   barley (Fig. [187]6C). The marked upregulation of key enzymes such as
   PAL and 4CL suggests an increased synthesis of polyphenols and
   flavonoids, compounds widely reported for their roles in UV shielding
   and antioxidant defense [[188]89, [189]90]. These findings align with
   metabolite profiling studies in highland barley, which demonstrate a
   synergistic selection for constitutive and inducible phenylpropanoids
   in UV-B protection [[190]3]. Furthermore, the results indicate that
   high-altitude conditions induce the production of greater quantities of
   medicinal compounds in barley [[191]91], highlighting its potential for
   applications in the production of bioactive phytochemicals. The
   co-enrichment of carotenoid and linoleic acid biosynthesis further
   supports a synergistic mechanism for photoprotection, as carotenoids
   quench singlet oxygen and dissipate excess light energy [[192]92],
   while unsaturated fatty acids enhance membrane fluidity to withstand
   temperature fluctuations [[193]93]. Equally notable is the observed
   upregulation of photosynthetic pathways. Elevated expression of genes
   involved in porphyrin-chlorophyll metabolism and Calvin cycle enzymes
   implies an adaptive enhancement of photosynthetic efficiency,
   potentially counterbalancing the reduced partial pressure and oxygen
   radical burden at high altitudes [[194]94]. This metabolic
   prioritization may reflect a compensatory strategy to maintain energy
   production under hypoxic conditions while managing photoinhibition
   risks.

   Research on high-altitude highland barley revealed that many genes
   involved in abiotic stress responses are transcription factors. This
   indicates that highland barley’s stress tolerance is robust and closely
   linked to transcriptional regulation. PPI and K-means clustering
   analyses identified three BZIP transcription factors co-participating
   in highland barley’s light response (Fig. [195]7B-C). Studies have
   shown that BZIP transcription factors are essential regulators of
   various biological processes in plants, including light signaling,
   stress responses, disease resistance, seed maturation, and flower
   development [[196]95]. In Arabidopsis, the activity of BZIP
   transcription factors has been implicated in responses to UV-B
   radiation [[197]9]. This study also confirms that the stability and
   activity of BZIP in highland barley are closely associated with light
   response mechanisms, particularly in high-altitude environments. The
   central role of BZIP transcription factors in the regulatory network
   highlights their importance as integrators of light and stress
   signaling. The interaction between STN7 kinase and BZIP transcription
   factors (CPRF2) in the PPI network further suggests a
   chloroplast-nuclear signaling axis that dynamically adjusts
   photosynthetic gene expression in response to light shifts, a critical
   adaptation given the altered spectral composition at high elevations
   [[198]96]. BZIP transcription factors (ELONGATED HYPOCOTYL 5) coupled
   with their interactions with COP1 and RUP2, reveals a sophisticated
   regulatory module balancing photomorphogenesis and stress acclimation.
   COP1-mediated ubiquitination typically destabilizes
   photomorphogenesis-promoting factors in darkness [[199]97], but its
   association with BZIP(ELONGATED HYPOCOTYL 5) in this context may
   represent a mechanism to accurately sense light intensity
   hierarchically regulating light-responsive gene expression under
   intense high-altitude irradiation. This interaction mechanism was
   predicted for the first time in highland barley. The concurrent
   downregulation of Cdk5 and mitogen-activated protein kinase 9-like
   suggests a phosphorylation-dependent modulation of BZIP activity. Such
   kinase-transcription factor interactions have been implicated in stress
   responses [[200]98], and here we propose their role in linking
   circadian or stress signals to adaptive transcriptional rewiring under
   high-altitude conditions.

   Interesting, from an evolutionary perspective, the observed metabolic
   shifts mirror convergent strategies seen in alpine plants, where
   phenylpropanoid overproduction and photosynthetic optimization are
   recurrent themes [[201]99]. However, the unique prominence of lipid
   remodeling pathways in highland barley suggests taxon-specific
   adaptations, possibly related to its domestication history in the
   Qinghai-Tibetan Plateau. The integration of these molecular
   adaptations—enhanced antioxidant capacity, optimized light harvesting,
   and dynamic transcriptional regulation—enables highland barley to
   thrive in extreme environments characterized by high UV radiation, low
   oxygen, and low temperatures. These findings advance our understanding
   of plant altitude adaptation by elucidating how metabolic plasticity
   interfaces with regulatory networks. Future studies should validate the
   functional roles of key candidates like BZIP heterodimers through
   mutagenesis or overexpression, particularly their influence on
   flavonoid biosynthesis and photoprotection. The insights gained not
   only deepen fundamental knowledge of plant-environment interactions but
   also provide actionable targets for engineering crops resilient to
   climate change-induced stresses.

Conclusion

   This study employed nitrogen ion beam mutagenesis on the highland
   barley cultivar ‘Kunlun 14’ (K14), generating 71 phenotypically
   distinct mutants. Notably, E8-38 exhibited a high 1000-grain weight,
   while D7-67 transitioned from a six-rowed to a two-rowed spike
   phenotype, highlighting their potential for high-yield breeding.
   Transcriptomic analysis revealed that the high 1000-grain weight of
   E8-38 was driven by synergistic regulation of phytohormone signaling,
   metabolic pathways, and epigenetic modifications, particularly the
   upregulation of the CKs signaling pathway and starch metabolism genes.
   In D7-67, the two-rowed spike phenotype was underpinned by the
   upregulation of VRS1 genes. Adaptive mechanisms to high-altitude
   environments included the upregulation of PAL and 4CL genes in
   phenylpropanoid biosynthesis, enhancing UV resistance and antioxidant
   capacity, as well as optimization of the photosynthetic pathway. BZIP
   transcription factors facilitated adaptation of highland barley to
   light changes and oxidative stress by regulating downstream gene
   expression. The COPI-RUP2-BZIP interaction module emerged as a key
   regulator balancing photomorphogenesis and stress responses in highland
   barley under high-altitude conditions, with phosphorylation-dependent
   regulation of BZIP activity linking circadian rhythms and stress
   signaling. Future work should prioritize regional trials of E8-38 and
   D7-67, functional validation of VRS1 in D7-67, and investigation of the
   STN7-BZIP and COPI-RUP2-BZIP interaction modules to elucidate
   chloroplast-nuclear signaling under high-altitude stress. These efforts
   will facilitate the development of high-yield, stress-resilient barley
   varieties for extreme environments.

Supplementary Information

   [202]Supplementary Material 1.^ (792.1KB, docx)
   [203]Supplementary Material 2.^ (42.8KB, xlsx)

Acknowledgements