A schematic diagram illustrating the overall experimental design is shown in Figure 1. We used coix germplasm YY03-03 provided by the Hunan Provincial Research Center of Engineering Technology for Grass Crop Germplasm Innovation and Utilization. Our research group conducted extensive preliminary evaluations of coix germplasm resources, which revealed that the YY03–03 germplasm exhibits three key advantageous characteristics: high yield, strong cadmium tolerance, and low cadmium accumulation capacity. These findings formed the basis for its selection in this study. The soil was taken from the rice experimental field (0∼20 cm) of Hunan Agricultural University. The basic soil characteristics were determined by conventional methods after soil drying (Table 1). The soil sample with no added Cd was used as control (CK). The remaining soil samples were treated with either 15 or 30 mg kg⁻¹ Cd and expressed as T1 and T2, respectively. Briefly, CdCl₂·2.5H₂O was ground in a mortar, and its solution was prepared with a 1/10000 dilution. The solution was loaded into a spray pot and mixed with crushed soil by spraying and mixing in plastic pots (diameter × height = 50 cm × 35 cm, leaky bottom belt tray). A total of 25 kg of soil was filled in each basin, and 50 g of compound fertilizer (NPK; 15% each of nitrogen, phosphorus, and potassium) was applied to each basin. Three biological repeats were used for each treatment, amounting to three plants per pot, nine plants per biological repeat, and 27 plants per treatment. Healthy seeds with full grains and uniform size were disinfected with 0.1 mol/L nitric acid solution for 24 h. The seeds were sown in a nutrition bowl and raised into seedlings till the six-leaf stage. Seedlings with uniform growth were shifted into pots with three plants per pot. The plants were grown under natural light, and conventional water and fertilizer management practices were employed.

During the filling stage, transpiration rate (Tr), stomatal conductance (GS), intercellular carbon dioxide concentration (Ci), and photosynthetic rate (Pn) of the coix plants were measured using a photosynthesis meter (Li6400XT). The second fully expanded leaf of the main stem from the top was analyzed, and the analyses were done between 9:00 am h and 11:00 am h on a sunny day.

After 111 days of Cd stress treatment, plant height, tiller number, and leaf color, morphology, and growth were measured under varying Cd stresses.

Dry root, stem, leaf, and grain samples (about 0.5 g per replicate) were crushed and digested separately in 5 mL of a digestive solution containing HNO₃ and HClO₄ in a 4:1 ratio. Following digestion, the samples were diluted to 25 mL of deionized H₂O. Cd concentration was measured using a flame atomic absorption spectrophotometer (ZEEnit 700, Analytik, Germany). We measured Cd (root), Cd (soil), Cd (stem), Cd (leaf), and Cd (grain), which refer to the Cd concentrations in root, soil, stem, leaf, and grain, respectively. Next, we calculated root [REF; Cd (root)/Cd (soil)], stem [SEF; Cd (stem)/Cd (soil)], leaf [LEF; Cd (leaf)/Cd (soil)], and grain enrichment factor [GEF; Cd (grain)/Cd (soil)], stem [STF; Cd (stem)/Cd (root)], leaf [LTF; Cd (leaf)/Cd (root)], and grain translocation factor [GTF; Cd (grain)/Cd (root)].

Total RNA was extracted from 100 mg coix root, stem, leaf and grain samples using the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The extracted RNA samples were stored at -80°C for subsequent steps such as library construction and sequencing.

The extracted mRNA (500 ng) was enriched using mRNA Capture Beads. After purification with beads, the mRNA was fragmented using high temperatures. The fragmented mRNA was then used as a template to synthesize the first strand of cDNA in a reverse transcription enzyme mixture system. While synthesizing the second strand of cDNA, end repair and A-tailing were completed. Next, adapters were ligated, and Hieff NGS. DNA Selection Beads were used for purification to select target fragments. PCR library amplification was then performed, and finally, detection was carried out using the Illumina NovaSeq 6000 platform at Gene Denovo Biotechnology Co. (Guangzhou, China).

Reads obtained from the sequencing machines included raw reads containing adapters or low-quality bases which would affect the following assembly and analysis. Thus, to get high quality clean reads, reads were further filtered by fastp (version 0.18.0) with the following parameters:

The cleaning reads were then assembled de novo using the Trinity software (Grabherr et al., 2011) with the parameter –KMER_SIZE 31 –min_kmer_cov 12 –normalize_reads –normalize_max_read_cov 50.

For each transcription region, a FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify its expression abundance and variations using RSEM software.

The DEGs in the Cd-treated samples were compared with those in the control samples using the DESeq2 software with significance criteria set as p-value < 0.05 and |log₂FC| ≥ 1.

Functions of identified genes were annotated using BlastX with an E value of less than 10⁻⁵ against gene ontology (GO), Swiss Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups (COG), and NCBI non-redundant (Nr) protein databases. GO and KEGG pathway enrichment analyses were performed using Blast2GO and KOBAS2.0, respectively (Kanehisa and Goto, 2000; Conesa et al., 2005).

The expression of 12 Cd stress-related DEGs identified by transcriptomics was confirmed by qRT-PCR. Gene expression profiles were determined in triple biological replicates using the actin gene as the internal control (Wei et al., 2019). The sequences of the primers used for qRT-PCR are listed in Supplementary Table S1.

Compared with controls, under Cd stress of 15 and 30 mg kg^-¹, Pn, GS, and Tr (19.97 μmol m⁻² s⁻¹, 0.21 mol m⁻² s⁻¹, and 9.50 mmol m⁻² s⁻¹, respectively) decreased significantly by 27.7% and 37.6%, 40.4% and 47.9%, and 54.0% and 64.8%, respectively (P < 0.05), while Ci (119.00 μmol mol⁻¹) increased significantly by 82.1% and 111.5% (P < 0.05). (Figure 2). The photosynthesis of coix was significantly inhibited under 15 and 30 mg kg^-¹ Cd stress.

Cd-treated plants exhibited significantly reduced growth compared to the control, including reduced plant height and tiller number and yellowish and curled leaves. Compared with controls (with a plant height of 126.37 cm and 17.83 tillers per plant), under Cd stress of 15 and 30 mg kg^-¹, the plant height and tiller number decreased significantly by factors of 19.0% and 24.2%, and 57.0% and 57.8%, respectively (P < 0.05). The degree of growth inhibition increased with increasing Cd concentration (Figure 3). The growth of coix was significantly inhibited under Cd stress at concentrations of 15 and 30 mg kg^-¹.

The concentration of Cd in roots, stems, leaves, and grains was significantly decrease in an order, which was (19.5 and 32.1), (7.3 and 14.8), (3.8 and 6.9), and (0.5 and 0.8) mg kg⁻¹ under 15 and 30 mg kg⁻¹ Cd, respectively and increased as the Cd stress concentration increased, and the peak of the Cd content was observed under 30 mg kg⁻¹ CdCl₂ treatment (Figure 4).

Furthermore, all enrichment and translocation factors decreased with increasing Cd concentration. These findings showed that coix has low Cd absorption and transportation capacity. The Cd enrichment capacity (enrichment factor) was highest in the roots of coix, followed by the stems and leaves, while the grains showed the lowest accumulation ability. The Cd translocation capacity (translocation factor) was highest in the stems of coix, followed by the leaves, while the grains showed the weakest translocation ability.

Using the Illumina NovaSeq 6000 platform, a total of 6.28–8.79, 7.93–9.37, 8.51–9.79, 6.62–7.38, 7.75–8.47, 7.31–8.37, 6.92–7.75, and 6.40–8.51 billion bases were acquired from control roots (CK-G), stems (CK-J), leaves (CK-Y), and grains (CK-Z), and Cd-treated roots (C-G), stems (C-J), leaves (C-Y), grains (C-Z), respectively (Table 2).

After filtration and trimming, 6.20–8.68, 7.63–8.34, 7.84–9.27, 7.22–8.28, 8.41–9.68, 6.84–7.66, 6.54–7.29, and 6.33–8.41 billion bases were acquired from CK-G, C-G, CK-J, C-J, CK-Y, C-Y, CK-Z, and C-Z, respectively. The Q20 and Q30 values were more than 97% and 93%, respectively, demonstrating the high quality of transcriptome sequencing data. The GC contents were about 52.54%–56.61% (Table 2).

In total, 1144, 2924, 3818, and 1702 genes were differentially expressed in the Cd-treated roots, stems, leaves, and grains compared to controls, with 682 (59.62%), 942 (32.22%), 1907 (49.95%), and 877 (51.53%) upregulated and 462 (40.38%), 1982 (67.78%), 1911 (50.05%), and 825 (48.47%) downregulated genes, respectively (Figure 5; Supplementary Tables S2–S5).

The number of downregulated genes in the Cd-treated stems was 2.10 times the number of upregulated genes, indicating that most genes were repressed under Cd stress. These findings indicated that Cd stress substantially influenced gene expression.

In root, stem, leaf, and grain, the DEGs were categorized into 51, 58, 58 and 51 functional groups, respectively, including “biological process” (BP), “cellular component” (CC), and “molecular function” (MF), with 22, 23, 24, and 21; 18, 21, 19, and 19; and 11, 14, 15, and 11 subcategories, respectively (Figure 6; Supplementary Tables S6–S9).

In BP, the majority of GO terms were clustered into metabolic process, biological regulation, localization, single organism process, regulation of biological process, cellular process, cellular component organization or biogenesis subcategories, and response to stimulus. In CC, the top seven subcategories were organelle part, membrane part, macromolecular complex, organelle, membrane, cell part, and cell. In MF, the transporter activity, catalytic activity, and binding subcategories were dominant. In particular, “antioxidant activity” in MF and “response to stimulus,” “signaling,” and “detoxification” in BP were functionally annotated with five, 83, 24, and six DEGs in root, respectively. These DEGs play vital roles in plant stress tolerance. A total of 35, 274, 62, and 16; 37, 280, 58, and 23; and ten, 114, 20, and nine DEGs were identified in stems, leaves, and grains, respectively.

The first four GO terms significantly enriched for DEGs between control and Cd-treated roots were transposition, DNA-mediated (GO: 0006313), transposition (GO: 0032196), transposase activity (GO: 0004803), and transferase activity (GO: 0016740) (Figure 7). In stems, the significantly enriched GO terms were cell wall (GO: 0005618), heme binding (GO: 0020037), iron ion binding (GO: 0005506), and tetrapyrrole binding (GO: 0046906). In leaves, the significantly enriched GO terms were thylakoid part (GO: 0044436), thylakoid (GO: 0009579), photosynthetic membrane (GO: 0034357), and photosynthesis (GO: 0015979). In grains, the significantly enriched GO terms were chloroplast part (GO: 0044434), thylakoid (GO: 0009579), photosynthesis (GO: 0015979), and thylakoid part (GO: 0044436). These data showed variable responses of roots, stems, leaves and grains to Cd stress. The responses of both leaves and grains collectively revealed the core harm of Cd stress—a severe inhibition of photosynthesis, thereby directly affecting plant growth.

All DEGs in roots, stems, leaves, and grains were mapped to 87, 120, 116, and 96 KEGG pathways, respectively. The pathways are divided into five categories, which were further divided into 18, 19, 19, and 19 subcategories for roots, stems, leaves, and grains, respectively (Figure 8; Supplementary Tables S10–S13). Cd stress affected subcategories “carbohydrate metabolism,” “translation” (root, stem, and leaf) or “replication and repair” (grain), “signal transduction,” “transport and catabolism,” and “environmental adaptation.” In addition, the DEGs were concentrated in the “energy metabolism,” “lipid metabolism,” and “amino acids metabolism” subcategories. These subcategories are related to plant responses to plant growth and development and environmental stimuli.

In Cd-treated roots, the DEGs primarily enriched in biosynthesis of secondary metabolites (51 DEGs), phenylpropanoid biosynthesis (ten DEGs), MAPK signaling pathway–plant (11 DEGs), plant-pathogen interaction (nine DEGs), plant hormone signal transduction (seven DEGs), flavonoid biosynthesis (three DEGs), and others (Figure 9; Supplementary Table S14).

In Cd-treated stems, the main altered pathways included biosynthesis of secondary metabolites (192 DEGs), metabolic pathways (387 DEGs), phenylpropanoid biosynthesis (44 DEGs), flavonoid biosynthesis (ten DEGs), glutathione metabolism (28 DEGs), photosynthesis (eight DEGs), photosynthesis-antenna proteins (four DEGs), tryptophan metabolism (22 DEGs), and others (Figure 9; Supplementary Table S15). A total of 22 DEGs (5 up-regulated and 17 down-regulated) were identified in the tryptophan metabolic pathway. Tryptophan serves as a precursor for the synthesis of key hormones such as auxin (IAA). Auxin is responsible for regulating plant growth, development, and stress responses. Alterations in this pathway indicate that Cd stress interferes with auxin synthesis, thereby inhibiting the growth of coix seeds.

In Cd-treated leaves, the primarily altered pathways included photosynthesis (24 DEGs), photosynthesis-antenna proteins (ten DEGs), phenylpropanoid biosynthesis (30 DEGs), biosynthesis of secondary metabolites (193 DEGs), metabolic pathways (334 DEGs), plant hormone signal transduction (24 DEGs), glutathione metabolism (24 DEGs), flavonoid biosynthesis (ten DEGs), and others (Figure 9; Supplementary Table S16).

In Cd-treated grains, the main altered pathways included biosynthesis of secondary metabolites (78 DEGs), photosynthesis (14 DEGs), photosynthesis-antenna proteins (four DEGs), plant-pathogen interaction (11 DEGs), metabolic pathways (128 DEGs), and others (Figure 9; Supplementary Table S17).

Under abiotic stress, transporters play vital roles in absorbing and transporting heavy metals and plant metabolism. The present study identified 30 transporter genes, including ABC (ATP-binding cassette transporter), in the root. Among them are four heavy metal transporters (HMTs). We observed a significant upregulation of all these HMTs under Cd stress, including Unigene0015490, Unigene0048533, Unigene0063149, and Unigene0101091 (Log₂FoldChange = 1.79, 2.18, 1.96, and 2.19, respectively; Supplementary Tables S18, S19). Furthermore, six genes were identified to be associated with detoxification, including five upregulated genes and one downregulated gene.

In total, 165 transporter genes, including MTP (metal tolerance protein), ABC, ZIP (zinc and iron regulated transporter proteins), and CAX (cation exchanger), were identified in the stem. This subset included 18 HMTs. Among them, six genes were significantly upregulated under Cd stress, including Unigene0004542, Unigene0063912, Unigene0071261, Unigene0092412, Unigene0095544, and Unigene0095545 (Log₂FoldChange = 1.25, 1.67, 1.83, 3.35, 1.93, and 1.55, respectively; Supplementary Tables S20 and S21). Furthermore, 16 genes were identified to be associated with detoxification, including six upregulated and ten downregulated genes.

Furthermore, 159 transporter genes, such as ABC, ZIP, YSL (yellow stripe-like transporter), and CAX, were identified in the leaf. This subset included 13 HMTs. Among them, seven genes were significantly upregulated under Cd stress, including Unigene0002853, Unigene0044711, Unigene0051971, Unigene0090192, Unigene0092438, Unigene0099561, and Unigene0099564 (Log₂FoldChange = 1.46, 2.16, 1.31, 1.97, 7.37, 1.42, 2.04, respectively; Supplementary Tables S22, S23). Furthermore, 23 genes were identified to be associated with detoxification, including nine upregulated and 14 downregulated genes.

Finally, 53 transporter genes, such as ABC, were identified in the grain. This subset contained nine HMTs. Of them, four genes were significantly upregulated under Cd stress, including Unigene0006961, Unigene0016620, Unigene0051971, Unigene0094764 (Log₂FoldChange = 2.49, 2.15, 1.37, and 1.34, respectively; Supplementary Tables S24, S25). Furthermore, nine genes were identified to be associated with detoxification, including five upregulated and four downregulated genes.

Previous studies reported several DETFs under Cd stress, such as MYB, WRKY, bHLH, NAC, and C2H2 family. In the present study, 1549 TFs from 54 families were identified in Cd-treated samples. Among them, C2H2, bHLH, bZIP, and ERF were identified as the four largest TF families with 153, 123, 119, and 119 TFs, respectively. The HSF, WRKY, GRAS, MYB, and NAC TFs, being responsive to a variety of stresses, were also characterized as large families, with 25, 87, 70, 68 and 90 members, respectively (Table 3).

Under Cd stress, expression analysis revealed that 39 TFs were significantly regulated in root. Among them, ERF, WRKY, GRAS, MYB, and Trihelix were ranked in the top five with three, four, three, three, and five DETFs, respectively. Interestingly, under Cd stress, most genes were found to be upregulated (682 out of 1144, accounting for 59.62%) and fewer genes were downregulated (Figure 5). Similarly, for most TF-encoding DEGs (27 out of 39, accounting for 81.69%) were upregulated (Table 3).

Under Cd stress, 94 TFs were significantly regulated in the stem. Among them, C2H2, bHLH, bZIP, WRKY, and NAC were ranked in the top five, with 14, 14, 11, 11, and six DETFs, respectively. Interestingly, under Cd stress, most genes were found to be downregulated and fewer genes (942 out of 2924, accounting for 32.22%) were upregulated (Figure 5). Similarly, fewer TF-related DEGs (37 out 94, accounting for 39.36%) were upregulated (Table 3).

Under Cd stress, 104 TFs were significantly regulated in the leaf. Among them, WRKY, bHLH, NAC, ERF, bZIP, and MYB were ranked in the top six with 11, ten, ten, eight, seven, and seven DETFs, respectively. Under Cd stress, about half of the genes (1907 out of 3818, accounting for 49.95%) were upregulated (Figure 5). Similarly, nearly half of the TF-encoding DEGs (44 out of 104, or 42.31%) were upregulated (Table 3).

Finally, 43 TFs were significantly regulated in grain. Among them, WRKY, bHLH, bZIP, ERF, NAC, and FAR1 were ranked in the top six with six, five, four, three, three, and three DETFs, respectively. Under Cd stress, around half of the genes (877 out of 1702, accounting for 51.53%) were upregulated (Figure 5). Similarly, nearly half of the TF-encoding DEGs (24 out of 43, or 55.81%) were upregulated (Table 3).

To experimentally validate the reliability of our transcriptome profiling data, we performed qRT-PCR analysis on 12 key DEGs implicated in Cd response and flavonoid biosynthesis. The selection criteria for these genes were based on their significant fold-changes and pivotal roles in the relevant biological pathways. The results demonstrated a strong concordance between the qRT-PCR data and the RNA-seq findings (Figure 10). This high degree of consistency confirms the accuracy and reproducibility of our transcriptomic data. Specifically, the validated genes include: Eight transporter-related genes: such as ABCG28 (Unigene0048533) and ABCC3 (Unigene0101091) in roots, MTP4 (Unigene0004542) in stems, and IRT1 (Unigene0090192), ABCA7 (Unigene0092438), YSL13 (Unigene0099561) and YSL12 (Unigene0099564) in leaves, ABCB11 (Unigene0006961) in grains. Four flavonoid biosynthetic genes: including HCT4 (Unigene0093309) in roots and CYP93G1 (Unigene0035743), CYP75B3 (Unigene0056120) and CYP73A12 (Unigene0090121) in stems. The successful validation of these key DEGs reinforces the credibility of our subsequent bioinformatic analyses and functional inferences regarding Cd detoxification and flavonoid accumulation in coix. These genes respond to Cd stress and play crucial roles in Cd uptake, efflux, sequestration, and antioxidation, thereby regulating Cd content in the roots, stems, leaves, and grains of coix.