Genome-wide identification of RING finger genes in flax (Linum usitatissimum) and analyses of their evolution

Identification of RING finger proteins in L. usitatissimum

In order to identify all potential RING finger proteins in the flax genome, the 469, 508, and 509 RING proteins reported in Arabidopsis (Jimenez-Lopez et al., 2018; Pavicic et al., 2017; Stone et al., 2005) were retrieved from the TAIR database (Cheng et al., 2017), merged, and used as queries in BlastP and tBlastN (Camacho et al., 2009) against the latest L. usitatissimum genome database (De Santana Lopes et al., 2018; Wang et al., 2012b) in Phytozome v13 (Goodstein et al., 2012) with an E-value of < 10⁻⁵. The potential flax RING finger proteins determined from the Blast results were subsequently verified using the InterProScan5 program (Jones et al., 2014) to confirm the presence of the RING domain using the SMART (Letunic & Bork, 2017), PROSITE (Sigrist et al., 2013), and Pfam (El-Gebali et al., 2018) databases. Then, the RING domain sequences of these verified potential RING finger proteins were extracted using custom Perl scripts (File S7) and aligned using MUSCLE v3.8.31 with default parameter (Edgar, 2004) and MAFFT v7 with L-INS-i method (Rozewicki et al., 2019) to confirm the position and composition of the eight ml residues. Sequences with an incomplete RING domain or those lacking classical ml residues were excluded from the final sequence dataset.

Phylogenetic, gene structure, and domain organization analyses of RING finger proteins

To examine the evolutionary relationships of the flax RING finger genes, the retrieved represented RING domain sequences were aligned using the MEGAX software (Kumar et al., 2018) and the MUSCLE algorithm with the default parameters. The neighbor-joining (NJ) tree of the RING domain sequences was constructed using the MEGAX software with 1,000 bootstrap repetitions. The maximum-likelihood (ML) tree of the RING domain sequences was constructed using the MEGAX software based on JTT+G+I model with 100 bootstrap repetitions. Phylogenetic trees were annotated by the ITOL v6 Web server (Letunic & Bork, 2019) and EvolView 3.0 (Subramanian et al., 2019). The gene structures of the RING finger genes were prepared using the Gene Structure Display Server 2.0 (Hu et al., 2015) through comparisons of the CDS (coding sequences) with the corresponding genomic sequences. The domain structures of the putative flax RING proteins were analyzed using SMART annotation with Genomic mode combined with the aforementioned annotation of InterProScan5 results and manually curated.

Chromosomal locations, duplications, and microsynteny analyses of flax RING finger genes

The chromosomal locations of flax RING finger genes were obtained using BlastN with an E-value of < 10⁻¹⁰ against the representative genome (GCA_000224295.2) from NCBI. The MCScanX software (Wang et al., 2012a) was used to confirm duplicated and syntenic flax RING finger genes using the default settings. The chromosomal locations and microsyntenic relationships of the RING finger genes were described using the Circos−0.69-9 software (Krzywinski et al., 2009).

To show the syntenic relationships of the homologous RING finger genes in flax and other plants, syntenic analysis maps were constructed using the TBtools software (version 1.098) (Chen et al., 2020). Ks and Ka substitution rates were calculated as previously described (Khan et al., 2020). The divergence time (T, MYA) of the duplications for each paralogous gene pair was estimated by the averaged Ks values from T = Ks/2 λ, in which (λ) represents the averaged Ks substitution rate for flax, which is 1. 5 × 10⁻⁸ (Koch, Haubold & Mitchell-Olds, 2000; Wang et al., 2012b).

Expression analysis of flax RING finger genes during seed developmental stages and in the AR and BR of the shoot apex

The RNA-seq data (accession number GSE130378 and GSE80178) (Li et al., 2021; Zhang & Deyholos, 2016) were downloaded from the NCBI GEO dataset (Tables S5 and S6) and filtered to explore the expression patterns of RING finger genes in flax. These transcript data (Li et al., 2021; Zhang & Deyholos, 2016) were obtained for seeds at four developmental stages (DAP5, DAP10, DAP20, and DAP30) and for the AR and BR of the shoot apex. The expression data represented log2 based Fragments per Kilobase of Exon per Million Fragments Mapped FPKM) values, and the heatmap of the expression profiles of flax RING finger genes was drawn using the heatmap.2 function in R.

Identification of RING finger proteins in L. usitatissimum

In Arabidopsis, 469, 508, and 509 RING finger proteins have been identified or curated (Jimenez-Lopez et al., 2018; Pavicic et al., 2017; Stone et al., 2005). In the present study, using the protein sequences of merged Arabidopsis RING finger proteins as queries, we identified 587 RING domains (File S2) in 574 flax RING finger proteins (File S1) via Blast (Camacho et al., 2009), InterProScan (Jones et al., 2014), and domain alignment analysis (Table S1). These predicted flax RING finger proteins included 563 proteins containing only a single RING finger domain, nine containing two RING finger domains, and two containing three RING finger domains. The protein lengths in the identified flax RING finger proteins ranged from 108 AA for Lus10029456 to 4872 AA for Lus10001530. The 587 putative RING domains can be classified into seven RING types according to the amino acid residues at the eight ml locations and the distances between them: RING-H2 (292), RING-HC (204) (RING-HCa (181), RING-HCb (23)), RING-v (53), RING-C2 (31), RING-D (2), RING-S/T (3), and RING-G (2).

The 587 RING domains comprised 292 (49.74%) of the RING-H2 type, the largest domain, followed by 204 domains (34.75%) of the RING-HC type (Table 1). Based on the spacing patterns between ml7 and ml8, RING-HC was classified into two subgroups: RING-HCa (181) and RING-HCb (23) (Table 1). In addition to the two canonical RING domain types, the five modified RING domain types (RING-v, RING-C2, RING-D, RING-S/T, and RING-G) represented only 15.50% of the total identified RING domain types. RING-v represented 9.03% (53) of the total RING domains and was featured by a cysteine (Cys, C) residue at the ml4 site and a histidine (His, H) residue at the ml5 site. This pattern was reversed in RING-HC. The fourth RING type, RING-C2, represented 5.28% of the RING domains and was featured by Cys residues at the ml4 and ml5 sites, rather than the His residue observed in the RING-H2 type. In the flax genome, the three RING-S/T-type protein sequences differed from the RING-HC proteins, with a serine (Ser, S) residue at the ml2 site or ml6 site instead of a Cys residue. The two RING-D-type and the two RING-G-type proteins differed from the RING-HC proteins by an aspartic acid (Asp, D) or a glycine (Gly, G) residue at the ml5 site instead of a Cys residue.

Conserved spaces between ml residues in flax RING domains

The eight ml residues of the RING domain coordinate two zinc ions in a unique cross-brace structure; ml1-ml2 and ml5-ml6 bind with one zinc ion, and ml3-ml4 and ml7-ml8 bind with the other (Freemont, 2000). This unique structure requires highly conserved spaces of the ml pairs ml1–ml2, ml3–ml4, ml4–ml5, ml5–ml6, and ml7–ml8, while relatively variable spaces of the ml pairs ml2-ml3 and ml6-ml7. The spaces patterns between the distinct ml residues of each of the eight putative RING types in flax are shown in Fig. 1. All 587 flax identified RING domains showed two amino acid residues (100%) between ml1 and ml2, whereas 99.66% (585/587) demonstrated two amino acid residues between ml5 and ml6, 93.19% (547/587) showed one amino acid between ml3 and ml4, 95.51% (561/587) contained two amino acids between ml7 and ml8, and 77.00% (452/587) showed two amino acid residues between ml4 and ml5 (Fig. 1A). The space between ml2 and ml3 comprised 8 to 33 residues, with 11 (140/587), 14 (118/587), and 15 (168/587) residues being most commonly observed in this space. The space between ml6 and ml7 comprised up to 64 residues, with 10 (272/587) being the most commonly observed number of residues (Fig. 1B).

Analysis of the variations in the spaces indicated patterns in the same RING types. For example, RING-H2 showed the highest number of amino acid residues between ml2 and ml3 (14 [102/292] or 15 [133/292]) sites and between ml6 and ml7 (10 [208/292]) sites. However, most of RING-HCa had only 11 (126/181) amino acid residues between ml2 and ml3 sites and 6 (24/181), 10 (60/181), 11 (20/181), or 12 (17/181) amino acids between ml6 and ml7 sites. Among RING-v, 87% showed seven amino acid residues between ml4 and ml5 sites, while 65% of RING-C2 had four amino acid residues between these sites (Fig. 1).

To examine whether any amino acid residues other than those involved in ml pairing are highly conserved in the flax RING domain sequences, a multiple sequence alignment of all RING domain sequences was conducted, and domain sequence logos of the distinct representative RING types were performed (Fig. 2). Ile (I) or Val (V) was the most frequently observed amino acid residue ahead of ml2 in the distinct RING types. Another highly conserved amino acid was Pro (P), which was always observed at the second site after ml7 site in over 94% of the putative RING domains, except RING-v, which showed an E amino acid instead. In the RING-H2 type, a phenylalanine (Phe, F) most often preceded the ml5 site, a leucine (Leu, L) was often present next to the ml2 site, and a D was typically present at the second site after the ml6 site. A G regularly existed in front of the ml4 site in RING-HC. More than 88% of RING-H2 and RING-v included a tryptophan (Trp, W) at the fourth site after the ml6 site. The residue after the ml1 site in RING-v was almost often an arginine (Arg, R).

Domain organization in flax RING finger proteins

To better classify the flax RING finger proteins, full sequences of the 574 identified RING finger proteins were examined using the SMART database and InterProScan to verify the conserved domain organization. Besides the RING domain, 68 additional domains were identified in these proteins, and based on their organization, they were classified into 45 groups and corresponding subgroups (Table S2). Group 1 (without additional domains, included 188 members) and group 2 (one or more transmembrane domains, included 157 members) were the two largest groups, accounting for 60% of the total proteins. The other groups or subgroups usually contained fewer members, and most of them had only several members. These additional domains in the putative flax RING finger proteins were predicted to play essential functional roles. The predicted protein-binding domains might be related to substrate recognition, including ANK, BRCT, SPRY, Vwaint, and WD40. The domains participating in ubiquitination involved CUE, GIDE, and RWD. The predicted nucleic acid-binding domains involved DEXDc, HIRAN, RRM, ZnF-NFX, ZnF-C2H2, and ZnF-C3H1. In this research, several domains, for examples Zn ²⁺ binding domains (ZnF-CHY, ZnF_NFX, ZnF_RBZ, Zinc_ribbon_6, and Zinc_ribbon_9) and heavy metal ion binding domains (HMA), were predicted to be involved in binding metal ions. The previously verified domains such as IBR, PA, Pep3_Vps18, SPRY, Zinc_Ribbon_6, Zinc_Ribbon_9, DUF1117, DUF1232, and DUF3675, were found to be associated with RING domains. Most of these additional domains were observed in diverse plants such as A. thaliana, O. sativa, V. vinifera, and P. trichocarpa, indicating that their functions could be conserved, and they may perform the same or similar functions among these species.

Phylogenetic and structural analyses of flax RING finger genes

To examine the evolutionary relationships among the RING genes in flax, the 587 RING domain sequences were aligned using the MEGAX software employing the MUSCLE algorithm with default parameters to construct a neighbor-joining (NJ) phylogenetic tree with 1,000 bootstrap repetitions (Fig. 3, Figs. S1, S2, and File S2) and a ML tree based on JTT+G+I model with 100 bootstrap repetitions (Fig. S4 and File S6). Based on the evolutionary analysis, the flax RING genes could be classified into three clades: Clade I (RING-H2 Clade), Clade II (RING-HC Clade), and Clade III (RING-temp Clade). Clade I comprised the highest number of RING-H2, RING-v, and RING-D genes and several RING-HC genes. Clade II included the highest number of RING-HC, RING-C2, RING-G, and RING-ST genes and several RING-H2 genes. Clade III consisted of RING-H2 and RING-HC genes.

Chromosome locations and duplication of flax RING finger genes

To investigate the chromosome locations of the predicted 574 RING finger genes in flax, the DNA sequences of all flax genes in Phytozome v13 were queried using BlastN against the representative genome assemble (GCA_000224295.2, chromosome level) from NCBI (E-value <10⁻¹⁰). The chromosomal locations of the 574 RING genes (Files S4 and File S5) were extracted from the Blast results under strict control. All the RING finger genes were distributed among 15 flax chromosomes at different densities: Lu1-Lu15 comprised 54, 47, 56, 39, 34, 35, 27, 45, 23, 15, 32, 47, 49, 34, and 37 RING finger genes, respectively (mean: 38; max: Lu3, 56; min: Lu10, 15) (Fig. 4). To examine the gene duplication events, all the flax RING finger genes were analyzed using BlastP and the MCScanX software. A total of 312 segmental duplication events and 27 tandem duplication events were confirmed in the whole flax genome (Table S3). To examine the selection pressure on the flax RING finger genes, the synonymous (Ks) and nonsynonymous (Ka) substitution rates and the Ka/Ks ratio of the RING finger gene pairs were calculated (Table S3). The Ka/Ks ratio was estimated to be less than 1, indicating that the duplicated RING finger genes in flax were under strong negative selection. Gene duplication analysis indicated that both segmental and tandem duplication events contributed to the expansion of the flax RING finger gene family.

To investigate the phylogenetic mechanisms of the flax RING finger genes, a comparative syntenic analysis (Fig. 5) of flax and the other four representative plants, including three dicots (P. trichocarpa, A. thaliana, and V. vinifera) and one monocot (O. sativa), was performed. A total of 572, 304, 314, and 42 orthologous RING finger gene pairs were observed between flax and the other four species (P. trichocarpa, A. thaliana, V. vinifera, and O. sativa), respectively. Additionally, most of the orthologous RING finger genes demonstrated a Ka/Ks ratio of less than 1, indicating that the RING finger gene family was under a strong negative selection pressure during the evolution process (Table S4).

Expression of flax RING finger genes during seed developmental stages

To examine the potential functions of RING finger genes during the seed development process in flax, a published flax RNA sequencing (RNA-seq) dataset (NCBI GEO accession: GSE130378) (Li et al., 2021) was used to confirm the expression profiles of these genes. A total of 505 RING finger genes were expressed at four stages (days after pollination [DAP]5, DAP10, DAP20, and DAP30) during the seed development process in flax, determined based on the transcriptome data; the other 69 genes were not found and excluded from the analysis (Fig. 6, Table S5). According to the hierarchical clustering of the expression patterns, these 505 genes were classified into six clusters. The RING finger genes in Cluster I (including 171 members) exhibited relatively lower expression levels at most stages, whereas Lus10023055, Lus10041156, Lus10040617, and Lus10025145, and Lus10012819 showed relatively higher expression at a specific stage than other genes in Cluster I. In Cluster II (including 74 members), typically, genes were highly expressed at different stages. For instance, Lus10007538 and Lus10028597 showed stronger expression than other clusters. Lus10011416 and Lus10009695 showed less abundance in seeds at DAP5 or DAP10. In Cluster III, IV, and V (including 56, 88, and 50 members), with an increase in the number of days, gene expression showed a downward trend, whereas the expression of genes in Cluster III and IV seemed to decrease faster than that of those in Cluster V. Additionally, genes in Cluster V showed relatively higher expression than those in Cluster III and IV at the same stages. The expression pattern of genes in Cluster VI (including 66 members) was different from that of others. The expression levels remained high at DAP5, decreased at DAP10, and increased again at DAP20 or DAP30.

Differential expression patterns of flax RING finger genes in the apical region (AR) and basal region (BR) of the shoot apex

To investigate RING finger genes that may contribute to the specification of phloem fiber identity in flax, publicly available RNA-seq data (NCBI GEO accession: GSE80718) (Zhang & Deyholos, 2016) were used to determine the differential expression patterns of these genes. 556 RING finger genes were expressed in the differentially expressed transcriptome data of the flax shoot apex. A total of 216 genes showed differential expression patterns in the apical and basal tissues (Fig. 7, Table S6). A total of 49 genes, including Lus10040716 and Lus10002639, showed high expression in the AR and showed low expression in the BR, whereas the other 167 genes, including Lus10013210, Lus10012811, and Lus10030728, demonstrated low expression in the AR and high expression in the BR. Differential transcript expression data corresponding to the shoot apex in which fiber specification is observed would complement other approaches aimed at understanding primary phloem fiber differentiation.