A total of 38 DIR and DIR-like members in S. italica were confirmed by HMMER and BLASTP methods. According to their chromosomal location, the 38 DIR genes were named SiDIR1 to SiDIR38. The protein sequence length of all the SiDIRs was less than 350 amino acids, ranging from 158 (SiDIR10) to 340 (SiDIR38); the molecular weight ranged from 16.8 (SiDIR8) to 34.9 (SiDIR7) kDa; the theoretical PI ranged from 4.79 (SiDIR38) to 9.82 (SiDIR9) ( Supplemental Data 1 ). The predicted subcellular localization of most SiDIR proteins was in the cell membrane, while the remaining family members were located in the cell wall, chloroplast, mitochondrion, nucleus, etc. ( Supplemental Data 1 ). The N-glycosylation sites of SiDIR protein sequences were also analyzed. It is interesting to find that 60% of proteins have the distribution of N-glycosylation sites ( Supplemental Data 1 ). These results showed the divergence of millet DIR genes and hinted at their distinct functional potential.

To analyze the homology of SiDIRs, we used ClustalW and found that the protein sequences within the same subfamily have high similarity ( Figure 1A ). For example, SiDIR7, SiDIR8, and SiDIR9 have high sequence similarity (78%–86%); SiDIR19, SiDIR20, SiDIR21, SiDIR23, SiDIR25, and SiDIR26 have high sequence similarity (72%–89%). Also, we found that individual amino acids within the conserved motifs have undergone specific mutations. The amino acid of SiDIR35 changed from “FG” to “FS”, and SiDIR2, SiDIR10, and SiDIR11 changed from “FG” to “LG”. In addition, the conservative amino acid of SiDIR2, SiDIR5, SiDIR8, SiDIR9, SiDIR10, and SiDIR11 changed from “VGRAQG” to “VARAQG” ( Supplemental Data 2 ). Due to the differences in amino acid sequences, different DIR proteins may have diverse functions.

To further elucidate the evolutionary relationship of SiDIRs, a phylogenetic tree containing 38 SiDIRs, 26 AtDIRs, 55 OsDIRs, 13 ZmDIRs, 11 GmDIRs, 4 GhDIRs, 3 LuDIRs, 2 TpDIRs, 1 FiDIR, and 1 ScDIR were constructed ( Figure 1B ). The DIR members could be grouped into DIR-a, DIR-b/d, DIR-c, DIR-e, DIR-f, and DIR-g. DIR-e (blue in Figure 1B ) contained the largest number at 11 SiDIRs, while DIR-a (pink purple in Figure 1B ) contained the smallest number at three SiDIRs. DIR-b/d (green in Figure 1B ), DIR-c (purple in Figure 1B ), DIR-f (pale blue in Figure 1B ), and DIR-g (red in Figure 1B ) contained four, nine, six, and five SiDIRs, respectively. Interestingly, no AtDIR member was found in the subgroup of DIR-c.

We used NCBI-CDD and MEME databases to analyze the gene structure and the distribution of conserved motifs, respectively. We found that the gene structure and conserved motifs of SiDIRs were similar within the same subgroups ( Figures 2A–D ). For example, SiDIR4, SiDIR7, SiDIR8, SiDIR9, and SiDIR38 contained five to six identical conserved motifs ( Figure 2B ). Except for the dirigent domains, the DIR family also contains three other specific domains including jacatin, dirigent superfamily domain, and tudor_AtPTM-like domain ( Figure 2C ). Some DIR-c subfamily members contain an N-terminal DIR domain and C-terminal end of a jacalin-related lectin (JRL) domain ( Figure 2C ), which show high specificity to bind mono- or oligo-saccharides (Kittur et al., 2007; Huwa et al., 2021; Huwa et al., 2022). Among the 38 SiDIRs, four members contained the jacatin domain, two members contained the dirigent superfamily domain, and no member contained the tudor_AtPTM-like domain ( Figure 2C ). Additionally, nine members (SiDIR17, SiDIR32, SiDIR31, SiDIR37, SiDIR10, SiDIR28, SiDIR5, SiDIR2, and SiDIR11) contained introns ranging from one to three, which is higher than that of the other DIR members. Most of the DIR and DIR-like family members contained exons ranging from two to three ( Figure 2D ). The differences in gene structure or conserved motifs might be due to various biological functions of SiDIRs.

FiDIR1 and DRR206 participate in the formation of (+)-pinoresinol (Pickel et al., 2010; Seneviratne et al., 2015). AtDIR6 and LuDIR5/LuDIR6 are able to form (−)-pinoresinol (Dalisay et al., 2015; Gasper et al., 2016). Based on the results, SiDIR7, SiDIR8, SiDIR9, TpDIR5/8, LuDIR1/5/6, AtDIR5/6, FiDIR, and ScDIR belonged to the same DIR-a subfamily ( Figure 1 ). We aligned the amino acids and found that some sites were highly conserved, such as alanine (A) located in β3 ( Figure 3A ). Phenylalanine (F) sites located in the β4 structure were similar, except for their mutation to isoleucine (I) in SiDIR9 ( Figure 3A ). Also, the hydrophilic amino acid “Y” located in the β4 structure mutated into unhydrophobic amino acid “F” (phenylalanine), while the hydrophobic amino acid “I” located in β5 structure mutated to hydrophobic amino acid “L” (leucine) ( Figure 3A ).

It has been known that the tertiary structure of the AtDIR6 trimer is an eight-stranded antiparallel β-barrel with spatially well-separated cavities for substrate binding. The binding cavity is composed of two lobes, and each of the two lobes is lined with a set of hydrophilic and potentially catalytic residues that are conserved in (+)- and (−)-pinoresinol-forming DIRs (Gasper et al., 2016). A further forecast showed that the tertiary structure of SiDIR7, SiDIR8, and SiDIR9 consisted of eight-stranded antiparallel β-barrels ( Figures 3B, D, F ). SiDIR7/8/9 were able to form trimers by themselves through homologous modeling prediction ( Figures 3C, E, G ). These suggested that SiDIR7, SiDIR8, and SiDIR9 may be able to direct the formation of (−)-pinoresinol. Meanwhile, the different locations of hydrophobic cavities and diverse substrate binding sites ( Figures 3B–G ) also potentially revealed their different functions.

According to genome annotation, we analyzed the distribution of 38 SiDIR genes. The 38 SiDIR genes were randomly distributed on nine chromosomes (Chr) ( Figure 4 ). Both Chr3 and Chr4 contain three members (~7.89%), while Chr9 contains two members (~5.26%). Compared to Chr5 (1 gene, ~2.63%), Chr8 contained the largest number of SiDIR family (15 genes, ~39.47%), which appears in the form of gene clusters ( Figure 4 ). Although there was no distribution on the first and sixth chromosomes, most SiDIR genes were distributed on the ends of the chromosomes.

Tandem duplication (TD) and whole-genome duplication (WGD)/segmental duplication (SD) of genes drive the evolution and expansion of gene family (Weidenbach et al., 2016; Zhang et al., 2022; Guo et al., 2023). We analyzed the duplication events of SiDIR genes and found 10 tandem duplication events involving 17 SiDIR genes on Chr2, 3, 4, 7, and 8 ( Figure 4 ; Supplemental Data 3 ). The genes of SiDIR7/8/9, SiDIR10/11, SiDIR13/15, SiDIR20/21, and SiDIR25/26 in tandem duplication events were from the same subgroup ( Figure 4 ), indicating the accuracy of the group division of the phylogenetic tree. However, no genome duplication (WGD)/SD events were found. These results revealed that TD events contribute largely to expanding SiDIR gene family.

Meanwhile, we counted the Ka, Ks, and Ka/Ks ratios of the duplication gene pairs using DNASP to analyze the evolutionary selection of duplication pairs in SiDIR gene family. We found that the Ka/Ks ratios of most gene pairs were less than 1 ( Supplemental Data 3 ), implying that these SiDIRs have undergone negative selection. Only two gene pairs, SiDIR13 and SiDIR14, and SiDIR23 and SiDIR24 ( Supplemental Data 3 ) had Ka/Ks ratios greater than 1, indicating that SiDIR13 and SiDIR14, and SiDIR23 and SiDIR24 may undergo positive selection and that they are important for the evolution of millets.

To investigate more deeply the evolution mechanisms of SiDIR genes, 12, 16, 14, and 8 ortholog DIR gene pairs were identified when compared millets with Arabidopsis thaliana, Oryza sativa, Zea mays, and G. max, respectively ( Figure 5 ; Supplemental Data 4 ). Interestingly, we found that some SiDIR genes were identified to be associated with at least three homologous gene pairs, such as SiDIR6 to AT4G11180.1, AT4G23690.1 and AT5G42500.1; SiDIR37 to AT1G65870.1, AT2G28670.1, AT2G39430.1 and AT3G55230.1; SiDIR3 to Glyma.07G157100.1, Glyma.08G258300.1, Glyma.18G282600.1 and Glyma.18G207700.1 ( Figure 5 ; Supplemental Data 4 ). Strikingly, some collinear gene pairs identified between Si and Os, and Zm were not found between Si and Os, and Zm were not found between Si and At, and Gm, such as SiDIR7/LOC_Os07g44250.1, and SiDIR7/Zm00001d022270 ( Figure 5 ; Supplemental Data 4 ). Some collinear gene pairs identified between Si and At were not found between Si and Os, Gm and Zm, such as SiDIR10/AT3G16450.1. We speculated that these orthologous genes may exert vital roles during the evolutionary process of DIRs.

To predict the functions and regulatory mechanisms of SiDIR genes, the cis-elements within their promoters were analyzed. A total of 23 different cis-elements in the promoter of SiDIRs were detected. These include cis-elements involved in phytohormone response (ABRE: abscisic acid, AuxRR-core: auxin, GARE-motif: gibberellin, CGTCA-motif: methyl jasmonate, and TCA-element: salicylic acid), light response, developmental regulation (circadian control, endosperm, flavonoid biosynthesis, and meristem), and environmental stress (defense and stress, drought, and low temperature) ( Figure 6 ). In plants, these cis-elements regulate various signaling pathways, hinting at the complicated regulatory function of SiDIRs.

In order to further clarify the functions and regulatory pathways among DIR family, a protein–protein interaction network of AtDIR family was analyzed and predicted by STRING software ( Figure S1 ). Nearly all of the AtDIR proteins could interact with other members. AtDIR9 (ESB1), At4g13580, and At2g39430 were central to the interaction network. At4g13580 might interact with AtDIR9 (ESB1), At2g39430, At5g42500, At3g24020, At4g11190, and AtDIR6. Similarly, At4g11190 may also interact with At3g24020 and AtDIR6. The interaction protein numbers of At3g13650, At1g22900, and At5g42510 might be the least. Interestingly, the interaction between At5g42510 and At2g39430 might be strong ( Figure S1 ). According to the protein similarity, we inferred that SiDIR family members might have a similar protein–protein interaction network to that of AtDIRs.

To gain insights into the functions of SiDIRs, the expression patterns of SiDIRs were analyzed based on RNA‐seq data. A total of 16 SiDIR genes including SiDIR1, SiDIR10, SiDIR11, SiDIR17, SiDIR19, SiDIR20, SiDIR21, SiDIR22, SiDIR23, SiDIR24, SiDIR26, SiDIR27, SiDIR31, SiDIR36, SiDIR37, and SiDIR38 expressed higher in root than other tissues ( Figure 7 ; Supplemental Data 5 ). Unlike the expression of SiDIR12, which was detected mainly in the stem and panicle, that of SiDIR3 was expressed mainly in the stem and root ( Figure 7 ; Supplemental Data 5 ). SiDIR10, SiDIR11, SiDIR19, and SiDIR31 are expressed only in root. Additionally, there were no transcripts of SiDIR5, SiDIR7, SiDIR9, SiDIR18, SiDIR28, SiDIR29, and SiDIR32 being detected in the organs of flag leaf, stem, root, panicle, and mesophyll ( Figure 7 ; Supplemental Data 5 ). These results hinted at the essential regulatory function of SiDIRs in plant development.

To further characterize the SiDIR genes in response to abiotic stresses, the millets were treated with CaCl₂, NaCl, CdCl, and PEG6000, and the transcriptional analysis of six SiDIRs in roots was carried out. After treatment with 20 mM of CaCl₂, the transcriptional level of SiDIR19 and SiDIR36 reached the maximum at 24 h, while SiDIR10 and SiDIR20 reached the maximum at 48 h ( Figure 8A ). Specifically, the expression level of SiDIR10 (at 48 h) and SiDIR36 (at 24 h) was approximately eight times higher than that of control. Unlike this, the relative expression of SiDIR22 and SiDIR27 was downregulated ( Figure 8A ). For 1 mM CdCl treatment, the transcription level of SiDIR10/19/20 was upregulated, while SiDIR22/27 was significantly downregulated. However, the expression induction of SiDIR36 was not obvious under CdCl treatment ( Figure 8B ). Upon treatment with 150 mM of NaCl, the expression of SiDIR10/19/22/27/36 was upregulated at 24 h and 48 h, and SiDIR20 was upregulated only at the time point of 48 h ( Figure 8C ). In addition, the expression levels of SiDIR19/20/22/27/36 were upregulated after being treated with 10% PEG6000 for 24 h, and SiDIR10 was not upregulated until being treated for 48 h ( Figure 8D ). These results showed that different SiDIR genes with different expression patterns may function differently during millet growth.

Co-expression analysis can help find genes that were closely co-regulated during the physiological process. Based on the MDSI database (Li et al., 2023), we constructed co-expression networks centered on the SiDIR10, SiDIR19, SiDIR20, SiDIR22, SiDIR27, and SiDIR36. As shown in Figure 9 , we obtained a total of six co-expression networks. Among them, the network centered on SiDIR19 is the largest (including 21 genes). In contrast, the network centered on SiDIR20 and SiDIR27 is the smallest (including one gene).

As shown in Supplemental Data 6 , the co-expressed gene network centered on SiDIR10 is significantly enriched in the carbohydrate metabolic process, oxidation-reduction process, ion binding, response to cadmium ion, and oxidation-reduction process. The network centered around SiDIR19 showed significant enrichment in disease resistance, sterol biosynthetic process, malate transport, metal ion binding, response to oxidative stress, cell wall biogenesis, cell wall organization process, response to salt stress, response to abscisic acid, phosphate starvation, and response to stimulus process. Meanwhile, the network centered on SiDIR36 showed significant enrichment in response to water deprivation, ion binding, lignin catabolic process, cell wall biogenesis, and carbohydrate metabolic process. Overall, these findings present an interesting phenomenon that warrants further investigations.

To deeply analyze the function of SiDIRs, the fusion expression vectors of 35S-SiDIR7-YFP, 35S-SiDIR19-YFP, and 35S-SiDIR22-YFP were constructed, with an empty vector of 35S-YFP used as a negative control and DAPI fluorescence signal used to indicate the nucleus. These vectors were transferred into the leaves of Nicotiana tabacum L., and the fluorescence was observed by ×20 laser confocal microscopy. The negative YFP signal was expressed in the membrane and nucleus, while the fluorescence signals of 35S-SiDIR7-YFP, 35S-SiDIR19-YFP, and 35S-SiDIR22-YFP were mainly expressed in the cell membrane ( Figure 10 ). This was coherent with the predicted analysis ( Supplemental Data 1 ).

Yugu 1 was used as the experimental material in this study. Millets were grown in a greenhouse in Wuhan, Hubei Province, China. Millets were grown in Hoagland nutrient solution (Li et al., 2022). For CaCl₂, NaCl, CdCl, and PEG6000 treatments, millet seedlings were grown in Hoagland solution for 10 days and then treated with 20 mM of CaCl₂, 150 mM of NaCl, 1 mM of CdCl, and 10% PEG6000, respectively. Plant roots were collected after treatment for 0 h, 24 h, and 48 h. Three biological replicates were carried out for each treatment.

The genome data of A. thaliana, S. italica, Z. mays, G. max, and O. sativa spp. japonica were downloaded from Ensembl Plants (http://plants.ensembl.org/index.html/). The AtDIR protein sequence was downloaded from TAIR (https://www.Arabidopsis.org/). The hidden Markov Model (HMM) file of the dirigent domain (PF03018) was downloaded as reported previously (Duan et al., 2023). HMMER 3.0 (E-value ≤ 1e⁻⁵, similarity > 50%) was used to search the DIR protein from the S. italica protein database. Further, based on the BLASTP method, we searched SiDIR protein sequences using AtDIR protein sequences (E-value ≤ 1e⁻⁵, similarity > 50%). All candidate DIR protein sequences were used to verify the DIR domain as analyzed previously (Duan et al., 2023). The longest transcript was obtained using the R package seqfinder (https://github.com/yueliu1115/seqfinder).

The phylogeny tree of identified SiDIRs and DIRs of rice, Arabidopsis, maize, soybean, cotton, etc., were constructed using the neighbor-joining (NJ) method of MEGA7.0 (bootstrap: 1,000 replications) (Kumar et al., 2016). The website of iTOL (Interactive Tree of Life, https://itol.embl.de/) was used to enhance the evolutionary tree.

The sequence alignment of DIR proteins was carried out by ClustalW, and ESPript 3.0 (https://espript.ibcp.fr/ESPript/ESPript/) was used to illustrate the conserved residues. For the prediction of tertiary structure, the protein sequences were input and analyzed through homologous modeling in SWISS-MODEL (https://swissmodel.expasy.org/).

The conserved motifs of S. italica DIR proteins were determined by MEME (http://meme-suite.org/) with a conserved motif number of 10. The gene structure information was acquired from GFF data. The conserved domains were obtained from NCBI-CDD (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and were subsequently visualized using TBtools software (Chen et al., 2020).

Segmental and tandem duplications were detected by MCScanX (Wang et al., 2012). The non-synonymous (Ka)/synonymous (Ks) ratios of duplication gene pairs were calculated using TBtools software. TBtools software was used to visualize the duplication events (Chen et al., 2020). The divergence time of all duplicate gene pairs was estimated as previously (Deng et al., 2019).

Gene expression level data of different tissues were downloaded from MDSi: Multi-omics Database for S. italica (http://foxtail-millet.biocloud.net/page/tools/expressionVisualization) (Yang et al., 2020; Li et al., 2023). The heatmap was generated by TBtools software (Chen et al., 2020).

The primer 5.0 software was used to design specific primers for SiDIR genes in this study ( Supplemental Data 7 ). Total RNA was extracted using the KKFast Plant RNApure Kit (ZOMANBIO, ZP405K-2). The cDNA was synthesized by PrimerScript™ IV 1st strand cDNA Synthesis Mix (TaKaRa, Mountain View, CA, USA; 6215A). The quantitative real-time PCR (RT-qPCR) system program was performed according to the previous research (Duan et al., 2023). The gene expression was analyzed by the 2^−ΔΔCT method as used in a previous study (Gong et al., 2022).

The CDS sequence of SiDIR7/19/22 was cloned from the cDNA of S. italica by the primers of SiDIR7-F/R, SiDIR19-F/R, and SiDIR22-F/R, respectively. The vector of 35S-YFP was digested with BamHI. The amplified products were then inserted into the linearized carrier of 35S-YFP by the kit of ClonExpress^® MultiS One Step Cloning (Vazyme, Nanjing, China; C113) and verified by DNA sequencing. These four vectors were transformed into Agrobacterium tumefaciens strain GV3101. Agrobacterium cultures harboring each construct were resuspended and mixed before being infiltrated into Nicotiana benthamiana leaves as described previously (Gong et al., 2022). After 48 h, the fluorescence signals were detected using a Leica TCS SP8 confocal microscope, and images were captured by LAS‐X software (Leica, Wetzlar, Germany).