Genome-wide identification and characterization of m6A regulatory genes in Soybean: Insights into evolution, miRNA interactions, and stress responses
Sabrina Bintay Sayed, Md. Afser Rabbi, Joy Prokash Debnath, Kabir Hossen, Ajit Ghosh, Kai Huang, Kai Huang, Kai Huang
TL;DR
This study identifies and characterizes m6A regulatory genes in soybean, revealing their roles in stress responses and evolution.
Contribution
The first genome-wide analysis of m6A regulatory genes in soybean, including their interactions and stress-related expression patterns.
Findings
42 m6A regulatory genes were identified in soybean, with 18 gene duplication events observed.
miRNA regulation was stronger for writer components than for erasers and readers.
Specific genes showed distinct expression patterns under stress conditions and mosaic virus infection.
Abstract
N6-methyladenosine (m6A) is one of the most prevalent mRNA modifications in eukaryotes, playing a crucial role in plant development and stress responses. The m6A modification is regulated by three key components: writers (which add methyl groups), erasers (which remove them), and readers (which interpret the modification). Despite its significance, the role of m6A regulatory genes in plants, particularly in soybean, remains largely unexplored. This study identified 42 m6A regulatory genes in soybean through a comprehensive genome-wide analysis. Structural analysis revealed diverse gene architectures and functional variations across different subgroups. A total of 18 gene duplication events were identified, predominantly evolving under purifying selection. Network analysis and hub gene identification suggested weaker interactions among eraser proteins. Moreover, interaction analysis…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Fig 1
Fig 2
Fig 3
Fig 4
Fig 5
Fig 6
Fig 7
Fig 8Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsRNA modifications and cancer · HVDC Systems and Fault Protection · Cancer-related gene regulation
Introduction
Plants have continuously suffered from various biotic and abiotic stressors such as pathogens, nutritional imbalance, salt, light intensity, and drought. To cope with these environmental variations, plants must evolve multiple mechanisms [1]. Epigenetic modification is one of them, which is associated with gene expression regulation under stress [2]. Epigenetics refers to stimuli-triggered changes in gene expression without altering the underlying DNA sequence that can be passed to offspring [3]. This process is mainly involved in DNA/RNA methylation, histone modifications, chromatin remodeling, and noncoding RNAs.
To date, RNA modification is an emerging major epigenetic modification in the plant life cycle and plays a critical role in the regulation of gene expression, both at the transcriptional and post-transcriptional levels [4]. Like DNA and protein, RNA goes through 160 types of post-translational modifications, where mRNA modification is found to be less common in comparison to tRNA and rRNA. Among mRNA modifications, N6-methyl-adenosine (m6A) is one of the most abundant, inner, dynamic, and reversible posttranscriptional modifications in eukaryotes and serves as a novel epigenetic marker that is involved in various biological processes [5]. Furthermore, three enzyme complexes, namely methyltransferases (writers), demethylases (erasers), and m6A-binding proteins (readers), are found to be involved in introducing, deleting, and interpreting specific methylation marks on mRNAs, respectively [6].
The writer complex was first discovered in mammals that included methyltransferase-like 3 (METTL3) [7], methyltransferase-like 14 (METTL14) [8], wilms tumor 1-associated protein (WTAP) [9], RNA binding motif protein 15/15B (RBM15/RBM15B) [10], vir-like m6A methyltransferase associated (VIRMA or KIAA1429) [11], zinc finger CCCH-type containing 13 (ZC3H13) [12], methyltransferase-like 16 (METTL16) [13], methyltransferase-like 4 (METTL4) [14], methyltransferase-like 5 (METTL5) [15], and zinc finger CCHC-type containing 4 (ZCCHC4) [16]. The orthologs of these complexes have also been identified in Arabidopsis, such as MTA (ortholog of METTL3), MTB (ortholog of METTL14), VIR (ortholog of VIRMA), and FIP37 (ortholog of WTAP) [17–20]. Moreover, E3 ubiquitin ligase HAKAI, an additional m6A writer component, has been identified in Arabidopsis [21]. METTL3 (S-adenosyl-methionine-binding protein) is the key catalytic component of the m6A methyltransferase complex. METTL14 and WTAP, making complex with METTL3, are found to be involved in M6A installation and modification, respectively [8]. Besides this, VIRMA is a regulatory subunit of m6A methyltransferase, recruiting the m6A complex to the special RNA site and facilitating m6A installation [11].
Conversely, the m6A demethylase erases the methyl group from the mRNA, indicating the reversibility of m6A modification. FTO (fat mass and obesity-associated protein) is the first eraser complex discovered in mammals that facilitates in restoration of the methylated base to the adenine base [22]. Subsequently, ALKBH5 (alkB homolog 5), the second mammalian m6A demethylase facilitating the removal of m6A, potentially affects different subsets of target mRNAs [23]. Arabidopsis contains several putative m6A eraser ALKBH family proteins [20]. ALKBH9B and ALKBH10B are the only two eraser proteins found in Arabidopsis that have been functionally investigated concerning viral infection and floral transition [24].
While writers and erasers are actively engaged in their dynamic action, the reader family regulates mRNA nuclear export, splicing, degradation, translation, and stability. The YT521-B homology (YTH) domain family is the first identified mammalian reader complex, which includes YTHDF1, YTHDF2, YTHDF3, YTHDC1, and YTHDC2 [25]. Different reader proteins have different functions on m6A-modified RNAs. For example, YTHDF2 regulates the degradation [26], YTHDF1 modulates translation [27], FMRP enhances nuclear export and stability [28], and YTHDC1 modulates nuclear export and splicing [29]. Furthermore, YTHDC2 and Eukaryotic initiation factor 3 (eIF3) are also crucial in mRNA translation [27,30]. Notably, three evolutionarily conserved C-terminal regions (ECT) family proteins have recently been functionally characterized in Arabidopsis as YTHD homologs [31,32].
The m6A modification is a crucial regulator performing three main functions: mRNA processing, plant growth and development, and stress response [33]. Researchers suggest that m6A regulatory writer proteins are essential for shoot and root growth, leaf, cotyledon, and floral development [18]. Partial loss of FIP37 causes huge over-proliferation of shoot meristems. Moreover, depletion of FIP37 results in embryo lethality [17]. The m6A modification-associated genes are expressed differentially in response to biotic and abiotic stresses, including drought [34], salt [35], cold [36], and bacterial infections [37]. Nevertheless, m6A modifications have been studied in limited plant species, highlighting the need for future investigation in other plants.
Soybean (Glycine max), a Leguminosae family member, has become one of the main economic oilseed beans [38]. Currently, researchers are shedding light on pinpointing stress-responsive writer as well as eraser genes in soybean [39,40]. Nevertheless, comprehensive or systematic analysis of all the m6A gene families remains scarce. Therefore, we were encouraged to conduct genome-wide identification of writers, erasers, and readers genes of the m6A pathway in soybean. A total of 12 writers, 11 erasers, and 19 readers were identified. A comprehensive analysis of the m6A regulatory genes in soybeans, including classification, evolution, gene structure, and potential interaction network, was performed. Expression patterns in response to a range of biotic and abiotic stresses were also studied.
Materials and methods
Genome-wide identification of the m6A gene family in Soybean
To identify m6A components and their protein families in soybeans, the amino acid sequences of m6A-related proteins reported in Arabidopsis thaliana [41], including writers, erasers, and readers, were used as queries to perform BLASTP against the soybean genomic sequences in the plant genomics resource database phytozome [42]. Putative genes and corresponding protein sequences were retrieved from the same database, and the presence of conserved domains was confirmed by NCBI Conserved Domain (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and Pfam (http://pfam.xfam.org/). The identity of the “writer” family was confirmed by the presence of the MT-A70 domain (PF05063), Wtap domain (PF17098), and Vir-N domain (PF15912). The presence of 2-OG Fe (II) oxygenase superfamily (PF13532) and YTH domain (PF04146) confirmed the identity of the “eraser” and “reader” families, respectively. Candidates without the m6A conserved and typical domains were not further considered. The prefix “Gm” for Glycine max was added to all members, followed by subclass IDs (MT, FIP37, VIRILIZER, ECT, CPSF, and ALKBH). Chromosomal location, strand position, CDS coordinate (5’ to 3’), gene length, CDS length, and protein length were retrieved from the Phytozome database [42]. Furthermore, the physicochemical properties of the identified proteins, such as molecular weight and theoretical isoelectric point, were calculated using the ExPASy-ProtParam tool [43]. Subcellular localization of proteins was predicted using BUSCA [44].
Chromosomal localization and gene duplication events
The GFF (General Feature Format) file of soybean was retrieved from the phytozome database [42]. Chromosomal locations were visualized by TBtool [45]. Using the Plant Genome Duplication Database (PGDD, http://pdgd.njau.edu.cn:8080/) [46], paralogous genes were found by calculating the number of duplication events among soybean genes. From the same database, the synonymous substitution rate (Ks) and non-synonymous substitution rate (Ka) were obtained. Using the Ka/Ks ratio, the selection pressure of duplicated genes was computed. Ka/Ks < 1 denotes purifying selection, Ka/Ks = 1 implies neutral drifting, and Ka/Ks > 1 signifies positive or Darwinian selection for a pair of genes [47]. Using T = Ks/ (2 × 6.1 × 10^−9^) × 10^−6^ Mya formula, the approximate date (Mya, million years ago) of each duplication event was estimated [48].
Gene structure, conserved domain, and conserved motif analysis
Genomic and CDS sequences were retrieved from the phytozome database [42]. Subsequently, the physical mapping of the exon-intron sequence was plotted using GSDS (Gene Structure Display Server 2.0, http://gsds.gao-lab.org/). The conserved domains of m6A regulatory protein sequences were analyzed in the Batch CD-Search program [49], and then the output file was submitted to TBtools for visualizing analyses. The conserved motifs were analyzed using TBtools with a maximum number of motifs of 10 [45].
Phylogenetic analysis
Protein sequences of writer, eraser, and reader of 13 species, namely Linum usitatissimum, Acorus americanus, Anacardium occidentale, Aquilegia coerulea, Brachypodium hybridum, Coffea arabica, Gossypium barbadense, Oryza sativa, Panicum hallii, Solanum lycopersicum, Triticum aestivum, Zea mays, and Porphyra umbilicalis were retrieved from the phytozome database. Multiple sequence alignment was performed using the MUSCLE algorithm [50], and phylogenetic trees for the protein families were constructed using Molecular Evolutionary Genetics Analysis (MEGA) 11 software [51]. The maximum likelihood method was employed with node reliability assessed through bootstrap analysis based on 1000 replicates. The resulting phylogenetic trees were further refined and graphically represented using the Interactive Tree of Life (iTOL) platform (https://itol.embl.de).
Promoter analysis in Soybean m6A regulatory genes
The 1000 bp upstream sequences from the transcription start site were collected and analyzed through the PlantCare database [52] for cis-acting element analysis. The presence of identified cis-acting elements on the putative promoter of every gene had been depicted using TBtool [45].
Functional enrichment visualization of m6A regulatory genes
Gene ontology enrichment in three categories, such as cellular component (CC), biological process (BP), and molecular function (MF), was executed using the STRING v12 database (https://string-db.org/). FDR (False discovery rate) =<0.05 was set. The number of terms was shown as 10. The results were shown as a bar chart and a bubble chart.
Protein-protein interactions, cluster prediction, and hub gene identifications
The STRING v12 database (https://string-db.org/) was used to construct the protein-protein interaction network among m6A regulatory genes in soybeans. The output file from STRING was visualized by Cytoscape software [53]. The network’s co-expressing network clusters (highly linked areas) were examined using the Molecular Complex Detection (MCODE) (v2.0.2) plugin [54]. Hub genes were recognized using the cytoHubba plugin [55] in Cytoscape [53] with the shortest path and ranked by degree.
Structure construction by homology modeling
The homology modeling approach was used to predict the three-dimensional (3D) structures of the m6A regulatory proteins. The amino acid sequence was queried against the SWISS-MODEL server (https://swissmodel.expasy.org/) [56] to search for templates, and the best templates with a similar amino acid sequence and known three-dimensional structures were used to build the models. These proteins’ structures were examined with the UCSF ChimeraX [57] visualization tools.
microRNA targets prediction
Mature miRNAs were downloaded from the miRBase database [58]. The psRNAT-target database [59] was used to search the regulatory relationship between miRNAs and soybean m6A regulatory genes, with an Expectation threshold of < 5 and other parameters at their default values. Cytoscape [53] was used to visualize the interaction network.
Expression analysis of m6A regulatory genes
RNAseq datasets (GSE137263 [60], GSE129509 [61], and GSE186317 [62]) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) of the NCBI (https://www.ncbi.nlm.nih.gov/) to observe the expression pattern of m6A genes in different organs, biotic, and abiotic stress. Normalized FPKM (fragments per kilobase of exon per million) data were used to observe expression levels in emerging nodules, mature nodules, emerging lateral roots, and mature lateral roots. To identify the expression patterns in abiotic stress, including salt stress, heat stress, dehydration, water deficit, and combined water deficit and heat stress, Log Fold Change (LFC) values were extracted using the R package DESeq2 [63]. The expression pattern of m6A regulatory genes in soybean mosaic virus (SMV) infection was also observed. The heatmaps were plotted using TBtool with normalized FPKM (roots and nodules) and LFC values (abiotic and biotic stress).
Results
Genome‑wide identification and characterization of m6A regulatory genes in Soybean
After a comprehensive screening of the soybean genome, 12 m6A writers, 11 m6A erasers, and 19 m6A readers were identified. All m6A regulatory genes were named according to their subclass identifier. The amino acid sequence length, relative molecular weights (MWs), isoelectric points (pI), subcellular localization, and gravy were calculated (Table 1). It was found that the soybean writer’s CDS length ranged from 624 bp to 6561 bp. GmVIRILIZER2 is the biggest member of the writer family and has a polypeptide length of 2187 aa and a molecular weight of 240.88 kDa. Meanwhile, GmFIP37c is the smallest, having a polypeptide length of 208 aa and a molecular weight of 23.13 kDa. The pI values of the writer family differed widely, ranging from 5.04 (GmFIP37a) to 8.81 (GmMTC2). Among the 12 writer members, ten exhibit acidic pI, while the remaining two show basic features. Likewise, the CDS length of the m6A eraser varies from 957 bp to 2052 bp. GmALKBH9B3 is the longest eraser protein with a 684 aa long and a molecular weight of 73.46 kDa. On the other hand, the smallest one, GmALKBH10B4*,* is 957 aa long with a molecular weight of 35.39 kDa. In terms of pI, eight of eleven erasers have acidic pI, whereas three others possess basic pI. The predicted CDS length of the m6A reader also varies from 387 bp (GmECT15) to 2127 bp (GmECT9). The largest one (GmECT9) is 77.17 kDa in size and 709 aa long. The smallest one (GmECT15) is 129 aa long with a molecular weight of 14.59 kDa. In terms of localization, most m6A regulatory proteins were in the nucleus, and only several proteins were in the plasma membrane (GmECT13, GmALKBH10B5, GmECT9) and chloroplast (GmMTC2, GmECT6). The Grand Average of Hydropathicity (GRAVY) values for all m6A proteins are less than zero, indicating that all soybean m6A proteins are hydrophilic.
Table 1: List of identified RNA m6A genes in Glycine max along with their detailed information and subcellular localization.
Chromosomal localization and gene duplication
The identified putative 42 genes have been distributed throughout the 17 chromosomes of soybeans (Table 1 and S1 Fig). The distribution of genes across the chromosomes can be summarized as follows: Chromosomes 17 and 8 each have the highest count, with four genes. Following closely, chromosomes 16, 10, 14, 5, and 2 each have three genes. Chromosomes 20 and 19 have two genes each. Finally, chromosomes 4, 6, 3, 1, 12, 9, and 15 each have a single gene in their respective location. Different duplication events of m6A regulatory genes in soybeans were identified through PGDD [46], as gene duplication is an important mechanism for gene family expansion. A total of 18 duplication events were observed (Table 2 and S1 Fig). Six writer’s sister gene pairs (GmCPSF30a-GmCPSF30b, GmMTA1-GmMTA2, GmMTB1-GmMTB2, GmFIP37a-GmFIP37c, GmFIP37b-GmFIP37d, GmVIRILIZER1-GmVIRILIZER2) were observed to be duplicated. Alongside, five erasers were discovered to be duplicated: GmALKBH10B7-GmALKBH9B1, GmALKBH10B1-GmALKBH10B3, GmALKBH10B4-GmALKBH10B5, GmALKBH10B2-GmALKBH10B6, and GmALKBH9B2-GmALKBH9B3. In addition, sister gene pairs were shown to be duplicated in seven readers (GmECT4-GmECT11, GmECT15-GmECT16, GmECT5-GmECT12, GmECT9-GmECT13, GmECT8-GmECT17, GmECT1-GmECT3, and GmECT10-GmECT14). All the duplication events of m6A regulatory genes were segmental duplications (Table 2). To explore the evolutionary history and direction of selection, non-synonymous to synonymous substitution (Ka/Ks) ratios were calculated. Purifying selection, neutral selection, and positive selection were indicated by Ka/Ks ratios below 1, equal to 1, and above 1, respectively. In this study, all analyzed duplication events yielded Ka/Ks ratios that were less than 1, indicating the influence of purifying selection on these genes throughout their evolutionary history. The time of divergence of these gene pairs varied from 5.71 to 15.50 Mya (Table 2).
Table 2: List of duplicated RNA m6A genes identified in soybean.
Structural analysis of writer, eraser, and reader members
Using amino acid sequences, maximum likelihood phylogenetic trees were built to investigate the structure and sequencing properties of the writer, eraser, and reader genes in soybeans (Fig 1A). All the writer genes, such as GmMTA, GmMTB, GmMTC, GmVIRILIZER, and GmFIP37 clustered according to their group. Likewise, all the GmALKBH10B genes clustered together except GmALKBH10B7, which formed a cluster with GmALKBH9B1. In the reader group, most of the GmECT genes clustered in the same clade. However, GmECT15 and GmECT16 were placed away from those clusters. Besides, GmECT10 and GmECT14 were clustered with GmCPSF30 genes rather than GmECT genes. The exon-intron distribution indicates the diversity of gene functions. According to the gene structure, the exon numbers of soybean m6A writers, erasers, and readers varied from 6–27, 6–8, and 2–9, respectively (Fig 1B). Amongst the m6A writers, the highest number of exons was observed in GmVIRILIZER1 (27), and the lowest number was seen in GmMTB1 (6) and GmMTB2 (6). In terms of eraser, all the genes contain only 6 to 8 exons, indicating the lowest number of exons among the three m6A regulatory proteins. Within readers, GmECT15 contains only two exons.
The evolutionary relationship, gene structure, and the appropriate domain architecture of the m6A-regulating genes in soybean.A. Using MEGA11, three groups’ phylogenetic trees were constructed. B. The gene structures were plotted by using GSDS. The exon-intron length was proportionately shown and could be approximated using the scale below. Exons are shown as blue boxes, upstream/downstream sequences are represented as dark red boxes, and introns are shown as grey lines. C. Functional domains of m6A regulatory genes in soybean. The YTH reader domains were marked in lavender, while the Eraser domain 2OG-Fe (II) oxy and superfamily domains were indicated by sky blue and cyan color, respectively. The MTA-70, MTA-70 superfamily, WTAP, and VIR-N writer domains were visually distinguished by light grey, dark grey, red, and yellow, respectively.
Group-specific conserved domains analysis indicates MT-A70, Wtap, HAKAI, and VIR_N as m6A writers. MT-70 domain was found in MTA, while the MT-70 superfamily was found in both MTB and MTC proteins (Fig 1C). In addition, GmFIP37a, GmFIP37b, GmFIP37c, and GmFIP37d were characterized with the WTAP superfamily. The VIR_N superfamily was found in GmVIRILIZER1 and GmVIRILIZER2. We found GmHAKAI according to phytozome. However, Pfam analysis did not show a HAKAI-related domain for GmHAKAI. Therefore, this gene was excluded from downstream analysis. The eraser family is characterized by the 2OG-FeII_Oxy superfamily (Fe (II) dependent Oxygenase superfamily) and 2OG-FeII_Oxy, all the identified putative eraser genes containing these conserved domains (Fig 1C). YTH (YT-521-B-like domain) was used to identify reader proteins, and all the identified readers contained that domain (Fig 1C). In the motif analysis, MTC1 and MTC2 contain only one motif (S2A Fig). Other members of the writer contain 4 to 6 motifs. All the eraser families contain 5–9 conserved motifs, while 2–5 motifs were observed in the reader family (S2B and S2C Fig).
The phylogenetic classifications of m6A regulatory genes in soybean
Demonstrating phylogenetic relationships is crucial as it gives insight into deducing the gene origins, examining their molecular adaptation, learning about the evolution of morphological traits, and rebuilding demographic variations in recently differentiated species [64]. Phylogenetic trees for m6A writers, erasers, and readers were constructed to study the evolutionary relationships of m6A modifiers. All the writer family, MT, FIP37(WTAP), VIR_N, and HAKAI formed four groups named I.II.III and IV (Fig 2A). The phylogeny revealed that GmVIRILIZER, GmMTA, GmMATB, and GmMTC have three common closest orthologs of Gossypium barbadense, Linum usitatissimum, and Anacardium occidentale. Another writer family, GmFIP37, has the closest ortholog of Coffea arabica and Solanum lycopersicum, while Linum usitatissimum is the distantly related ortholog. Zea mays is the most distantly related ortholog of GmVIRILIZER, demonstrating the plant kingdom’s evolutionary divergence. On the other hand, the red algae species Porphyra umbilicalis has the most distant association with the orthologs of GmMTA, suggesting a wider evolutionary range. The orthologs of GmMTB and GmMTC are most distantly linked to Coffea arabica, highlighting the evolutionary connection among flower species.
Phylogenetic analysis of m6A regulatory genes from six dicotyledons, six monocotyledons, and one red alga.A. Phylogenetic tree of writer genes. B. Phylogenetic tree of m6A eraser genes. C. Phylogenetic tree of m6A reader genes.
The YTH domain-containing eraser protein orthologs also form 4 different groups named I, II, III, and IV (Fig 2B). The various species’ readers are categorized according to their class (Fig 2C). Most of these readers’ orthologs of Oryza sativa are out-grouped from their respective clades. While monocot-monocot and dicot-dicot couples are found in certain sister taxa, monocot-dicot pairs are found in others. The closely related species A. americanus and A. coerulea, for example, show strong evidence for their nodes and a recent shared ancestor. The implication here is that monocots and dicots may have split out from a common ancestor. These three evolutionary trees demonstrate how m6A dynamics are conserved across a variety of kingdoms. The presence and functional relevance of the m6A alteration in different species within these categories have been demonstrated by these phylogenies.
Identification of cis-regulatory elements in the putative promoter region
The spatiotemporal specific expression of protein-coding genes can be uncovered by analyzing the cis-acting regulatory elements (CREs) in the promoter region. Using the PlantCare database [52], CREs in the promoter regions (1000 bp upstream of the transcription start site) of m6A regulatory genes were predicted to better understand the transcriptional regulation and potential biological functions of soybean m6A regulatory genes (Fig 3). These CRE elements were divided into four groups, including light-responsive, phytohormone-responsive, stress-responsive, and development-related. Light-responsive elements were predominant in all the m6A regulatory genes, suggesting that these gene sets may be regulated by light signaling. ABRE and TGA are predominant in phytohormone elements. This substantial number of phytohormone-responsive elements indicated the possible activation or inhibition of the expression level of the m6A regulatory gene by the hormone signals. In terms of development, CAT-box, GCN4-motif, ARE, and O_2_-site were enriched in the promoter regions. These CREs were respectively associated with meristem expression, endosperm expression, anaerobic induction, and metabolism regulation. MBS (related to drought response), TC-rich repeat (related to defense and stress), and LTR (related to cold response) elements were also highly enriched.
The distribution of cis-acting elements in the promoter regions of m6A regulatory genes.The 1000 bp 5′ upstream region of all the identified m6A genes was retrieved and analyzed using the PlantCARE database to determine the presence and number of cis-acting regulatory elements. The presence of selected motifs (x-axis) in the promoter of corresponding genes (Y-axis) was represented by the number of motifs in the heatmap.
Functional analysis of m6A regulatory genes
The GO analysis of the m6A regulatory genes revealed that the biological process was mostly engaged in the regulation of gene expression (Fig 4A). Nucleic acid metabolic process, regulation of mRNA metabolic process, negative regulation of gene expression, and mRNA metabolic process were the dominant molecular processes. RNA binding, mRNA binding, Nucleic acid binding, and binding were observed in molecular function (Fig 4B). Intracellular anatomical structure was seen in the cellular component (Fig 4C). These results highlight the critical function of m6A regulatory genes in coordinating cellular and molecular processes that influence the dynamics of genes and RNA.
Functional analysis of m6A regulatory genes in Soybean.A. Bubble diagram illustrating the enrichment of the gene ontology for the category of biological processes. B. A bubble chart for the Gene Ontology enrichment analysis term for the category of molecular function C. Bubble diagram for the term cellular component category from the analysis of Gene Ontology enrichment. Bubble size is proportional to the number of genes, and bubble color represents the value of False Discovery Rate (FDR). Lighter green bubbles indicate lower FDR values, meaning the results are statistically significant.
Protein-protein interaction and cluster identification
A multi-protein complex performs together in m6A modification, highlighting the prediction of a potential protein-protein interaction network in this regard. There were 192 edges and 42 nodes in the protein-protein interaction (PPI) network (Fig 5). Two clusters were identified by using MCODE [54]. Cluster 1 contains writer proteins, including MTA, MTB, and VIRILLIZER (S3A Fig). Cluster 2 contains reader and writer proteins named GmFIP37a, GmFIP37b, GmFIP37c, GmFIP37d, MTC1, MTC 2, CPSF30a and ECT14 (S3B Fig). The cluster formation gives valuable insight into the biological mechanism of m6A regulatory proteins. These proteins are vital for plant growth and development, as well as for how well plants adapt to their surroundings. These findings shed light on the molecular mechanism and possible biological purposes of soybean m6A regulatory genes.
PPI network of m6A regulatory genes.PPI networks of writer eraser and reader genes were generated using the STRING database. Nodes represent the proteins, and edges represent the protein-protein associations.
Identification of hub genes and their 3D structure modeling
The top ten hub genes were extracted by the cytohubba plugin [55] in Cytoscape [53]. Notably, there was no eraser member in this hub gene (S3C Fig). GmVIRILIZER1 may interact with GmVIRILIZER2, GmECT10, GmECT14, GmMTB1, and GmMTB2, according to the network. GmECT10 and GmECT14 interact with GmMTB1. Furthermore, GmMTA1 engages in interactions with every hub gene in this network (S3C Fig). However, GmMTA2 has no connection with GmMTA1 or GmCPSF30a. Furthermore, GmMTB2 engages in interactions with GmMTB1, GmECT14, and GmECT10. No interaction was found between GmCPSF30a and reader hub genes. However, GmFIP37d interacts with two reader families, including GmECT10 and GmECT14. GmVIRILIZER2 interacts with GmMTB1, GmECT10, GmECT14 and GmMTB2. The protein structure of 10 hub members was predicted using SWISS-MODEL [56]. Protein structure has a polypeptide chain of amino acids specific to each protein type. The chain is folded into a secondary structure characterized by alpha-helices and beta-sheets [65] (Fig 6).
3D structure prediction of 10 hub genes using SWISS-MODEL.Blue and green colors represent the alpha helix and beta sheets, respectively.
MicroRNA target prediction
Plant microRNAs (miRNAs), small non-coding RNAs which are usually form duplexes with their targets. They play a crucial role in regulating a wide range of biological metabolic activities, including hormone signaling, plant growth, flowering time, leaf morphogenesis, and reactions to environmental stressors. Therefore, potential miRNA binding sites of m6A regulatory genes were examined in this study by using psRNATarget [59]. After filtering interactions with an Expectation penalty score below 5, it was found that 10 writers had 81 miRNAs, 10 erasers had 31 miRNAs, and 18 readers had 60 miRNAs (S4 Fig). It is noteworthy that soybean m6A writers (n = 81) seemed to be more regulated by miRNA than m6A readers (n = 60) and erasers (n = 31). Interestingly, our research did not identify any miRNAs that regulate the two writers (GmMTC1 and GmMTC2), one eraser (ALKBH9B2), and one reader (ECT16) in soybeans. However, more research is required to fully understand the complex interactions between m6A regulatory genes and miRNA.
Expression analysis of m6A regulatory genes in roots and nodules of Soybean
Symbiotic legume nodules and lateral roots arise away from the root meristem via dedifferentiation events in soybean [61]. These organs share some morphological and developmental similarities. We explored the expression level of m6A in emerging nodules (EN), mature nodules (MN), emerging lateral roots (ELR), and young lateral roots (YLR) of soybeans. Most of the genes displayed a broad expression range across all the samples, indicating that they were extensively involved in the growth and development of roots and nodules (Fig 7). GmMTA1 displayed vast expression levels in mature nodules. Notably, GmMTB1, GmFIP37a, GmFIP37d, GmALKBH10B6, GmECT1, and GmECT12 showed a medium level of expression in all the selected samples (EN, MN, ELR, and YLR). GmALKBH10B5 showed higher expression levels in mature nodules compared to EN, ELR, and YLR. Additionally, GmECT9 was expressed tremendously in all the selected tissues, while GmECT13 exhibited immense expression in EN and ELR compared to MN and YLR. Importantly, GmECT17 showed higher expression in ELR but displayed medium expression in the other three tissues.
Expression profiling of m6A regulatory genes in roots and nodules of soybean.Soybean RNA writer, eraser, and reader gene’s genome-wide RNA-Seq data were obtained from the GEO of NCBI. Expression levels in emerging nodules, mature nodules, emerging lateral roots, and mature lateral roots were observed by using normalized FPKM (fragments per kilobase of exon per million). The color scale provided at the right of the figure represents the level of expression.
Expression patterns of m6A regulatory genes under various abiotic and biotic stress
The expression pattern of m6A regulatory genes was observed under various abiotic and biotic stresses (Figs 8A and 8B). The majority of m6A regulatory genes remained unresponsive when subjected to salt stress and dehydration. GmMTB2 and GmECT9 exhibited significant upregulations in heat stress and combined water deficit and heat stress (Fig 8A). On the contrary, GmALKBH9B2 showed downregulation in combined water deficit and heat stress. The remaining soybean’s m6A genes displayed a medium level of expression except GmECT4 and GmECT2. These two genes did not show any changes across all the abiotic stresses (Fig 8A).
Expression profiling of m6A regulatory genes under abiotic and biotic stress.A. The expression pattern of 42 genes was analyzed in response to five abiotic stress conditions, including salt, heat, dehydration, water deficit, and combined water deficit and heat stress. B. The expression profile of writer, eraser, and reader was analyzed in response to SMV (Soybean Mosaic Virus). Heatmaps were generated using TBtool with the LFC (Log Fold Change), which delineates the expression pattern of 42 m6A regulatory genes. The color scale provided at the right of the figure represents the level of expression. The stress-induced upregulation and downregulation of all the genes are indicated by the red and blue colors, respectively.
The soybean mosaic virus (SMV) treatment resulted in the significantly higher up-regulation of multiple genes, such as GmMTB1, GmMTB2, and GmALKBH9B1 (Fig 8B). In the writer’s family GmMTA2, GmCPSF30a, GmVIRILIZER1, GmMTA1, GmVIRILIZER2, and GmCPSF30b showed upregulation compared to control. In the ALKBH group, three genes, such as GmALKBH10B3, GmALKBH10B7, and GmALKBH10B6 upregulated in SMV treatment. Additionally, in the readers group, GmECT9, GmECT2, GmECT17, GmECT15, and GmECT7 displayed upregulation. GmALKBH9B4 was the only member of the m6A regulatory gene that displayed significant downregulation in SVM treatment. Overall, the expression patterns of m6A regulatory genes in soybeans under different stress conditions indicated that they might play critical biological roles in various stress responses.
Discussion
RNA N6-methyladenosine plays a crucial regulatory role in plant growth and development. The methylation levels of target transcripts are dynamically regulated by three types of m6A regulatory proteins. A total of 42 candidate m6A regulatory genes were identified in soybean, including 12 writers, 11 erasers, and 19 readers (Table 1). The gene sets were confirmed by the conserved domains of writer, reader, and eraser. Soybean (42) has more m6A regulatory genes than others, such as Arabidopsis (33), tomato (25), grape (40), rice (33), S. moellendorffii (22), M. polymorpha (16), P. patens (18), and Chinese pine (36) [66]. On the other hand, soybean falls short of tobacco (52) and upland cotton (75), common wheat (85), and maize (55) [66]. According to Su et al [66], a low number of m6A regulatory genes were found in Chinese pine (36) compared to species such as common wheat and upland cotton, though Chinese pine has the largest genome size (25.4 Gb). The plants that have experienced polyploidization events, such as common wheat (2n = 6x = 42), upland cotton (2n = 4x = 52), and tobacco (2n = 4x = 48), exhibit a higher abundance of m6A regulatory genes [67–69]. Likewise, soybean (2n = 4x = 40) [70], a polyploid species, also possesses a higher number of m6A regulatory genes. Additionally, maize, a segmental allopolyploid, also shows a higher number of m6A regulatory genes [71]. Using PGDD, a total of 18 duplication events were identified in this study. All the duplication events were segmental duplication (Table 2, and S1 Fig**).** From the Ka and Ks analysis, duplicated gene pairs went through purifying selection, a finding consistent with studies on tomatoes [72], tea plants [73], and tobacco [66].
To explore the evolutionary patterns of m6A writer, eraser, and reader genes, phylogenetic trees were constructed across 13 plant species, which included Rhodophyta (red algae) and Angiosperms (comprising eudicots and monocots) (Fig 2). Monocot-dicot, monocot-monocot, and dicot-dicot pairs were spotted in this study. For instance, Acorus americanus and Aquilegia coerulea share a recent ancestor, implying shared monocot-dicot divergence. In the previous study on m6A genes of Oryza sativa, a similar evolutionary pattern was observed [74]. Cis-regulatory elements in the promoter region were identified. According to this analysis, CREs were associated with light response, hormone response, plant growth and development, and stress response. Similar findings were also observed in other plant species such as tobacco [66], tea [73], tomato [72], and poplar [75]. Light-responsive CREs were predominant in soybean m6A regulatory genes (Fig 3).
GO terms in three categories, including biological process, cellular components, and molecular function, were documented. These results indicate the involvement in RNA modification of writer, eraser, and reader genes of soybean (Fig 4). The PPI network revealed that MTA1, MTA2, MTB1, and MTB2 are highly interconnected genes, indicating their importance as the catalytic core (Fig 5). Most of the erasers did not interact with the writer and the reader. On the other side, CPSF30a and CPSF30b, two readers, had stronger interaction than any other reader components (Fig 5). Cluster analysis showed two cluster formations of m6A regulatory genes of soybeans (S3 Fig). Six writer components formed cluster 1, highlighting their dynamic action in RNA modification (S3A Fig). Writer components also predominated in cluster 2; nevertheless, no erasers were found in any cluster (S3A and S3B Fig). Previous studies also found similar results that there were no interactions among eraser proteins in tobacco [66]. Furthermore, hub gene identification exhibited that there were no eraser genes in the top 10 hub genes (S3C Fig). Using homology modeling, the 3D structure of 10 hub genes was predicted as these genes might be the core catalytic component in m6A methylation (Fig 6). The secondary structures, such as alpha helices and beta sheets, are important to maintain protein stability and functionality [76]. To explore the mutual interaction of miRNA and m6A soybean, miRNA target prediction was performed. Interestingly, the writer component had greater miRNA regulation compared to readers and erasers (S4 Fig). These findings indicate that soybean writers may prefer to be regulated by miRNA.
Researchers have demonstrated that altering the m6A system may cause irregular growth and development. Many previous studies have mentioned the association of m6A modification with reproductive cell development [77] and fruit ripening [78]. In this study, soybean m6A expression was observed in roots and nodules. Readers and erasers have a higher expression than writers regarding roots and nodules (Fig 7). Two reader genes named GmECT9 and GmECT13 are highly expressed in the respective samples. According to prior studies, CPSF30L in Arabidopsis might regulate APA by binding to m6A-modified RNAs that are involved in nitrate signaling, particularly WRKY1 and NRT1.1 [79]. As GmECT9 and GmECT13 were highly expressed in roots and nodules, it can be hypothesized that these genes have a role in nitrogen absorption. However, further investigation is needed in this regard.
CRE elements such as MBS (related to drought response), TC-rich repeat (related to defense and stress), and LTR (related to cold response) were highly enriched in soybean m6A genes according to Promoter analysis (Fig 3). This result indicates the involvement of soybeans’ m6A genes in environmental stimuli. Salt stress has been shown to have a major impact on the m6A methylation levels of mRNA in Arabidopsis [35]. The expression patterns of GmMTAs and GmMTBs under various abiotic stressors have been documented in earlier research, suggesting their possible role in stress tolerance, particularly in the reaction to darkness or alkalinity [39]. According to another research, GmALKB10B genes were probably triggered by alkalinity, cold, and dryness [40]. In this study, five abiotic stress responses, such as salt stress, dehydration, heat stress, water deficit, and combined water and heat have been reported. Most of the soybean’s m6A genes did not show any changes in salt and dehydration stress, indicating their stress tolerance possibility (Fig 8A). Notably, two genes named GmMTB2 and GmALKBH9B2 showed upregulation and downregulation, respectively, in combined water and heat stresses. Interestingly, GmECT4 and GmECT2 did not show any changes across all the abiotic stresses. Therefore, we assumed that these two genes might be stress tolerance m6A genes in soybeans. Validating these findings requires a wet lab inquiry. Distinct biotic stresses also have a great impact on the plant life cycle. NtALKBH10 and NtVIR1 of tobacco exhibited robust responses to specific biotic stresses (R. solanacearum and black shank) but showed insensitivity to the Cucumber Mosaic Virus (CMV) [66]. In this study, biotic stress responses (soybean mosaic virus) have been documented. Most of the soybean’s m6A genes displayed upregulation in SMV treatment (Fig 8B). Nevertheless, one member called GmALKBH9B4 was downregulated in this regard. In conclusion, the m6A regulatory genes of soybeans showed a variety of expression patterns under stress, nodules, and roots, suggesting that these genes have a multitude of roles. The precise processes of the putative m6A regulatory genes in soybeans require further investigation.
Conclusion
This study presents a comprehensive and systematic investigation of m6A regulatory genes in soybean, and a total of 42 m6A regulatory genes were identified. We analyzed the features of these genes in terms of gene structure, conserved domains, and motifs. Phylogenetic study and duplication events were also examined. Eighteen duplication events that evolved through purifying selection were found. Functional analysis of m6A regulatory genes in cis-elements and interaction networks demonstrated their crucial roles in soybeans. The altered expression pattern of these genes had been profiled in roots, nodules, and under distinct stresses. For instance, the three readers—GmECT9, GmECT13, and GmECT17—were highly expressed in roots and nodules. GmMTB2 exhibited the highest upregulation in response to combined water deficit and heat stress, while GmALKBH9B2 showed the greatest downregulation. Furthermore, GmALKBH9B4 was significantly downregulated, whereas GmMTB1, GmMTB2, and GmALKBH9B1 were upregulated under soybean mosaic virus treatment. This study gives a reference framework for investigating the functional variety of soybean m6A regulatory genes at the epigenetic level.
Supporting information
S1 FigChromosomal distribution of m6A regulatory genes.All m6A regulatory genes are found to be located on different chromosomes of soybeans. The relative size of the corresponding chromosomes and the position of the respective genes could be estimated by using the scale provided left side of the figure. The chromosome numbers are provided in the middle of each bar. The straight lines connect the duplicated gene pairs.(PDF)
S2 FigOrganization and distribution of the conserved motifs.The organization and distribution of the conserved motifs in the m6A writer (A), eraser (B), and reader (C) genes. The squares in the motif represent the positions of conserved domains.(PDF)
S3 FigCluster and hub genes identification in soybeans.Two clusters (A and B) were identified using MCODE. C. The top 10 hub genes were extracted by the Density of Maximum Neighborhood Component method in the cytohubba plugin in Cytoscape.(PDF)
S4 FigGenome-wide network of soybean miRNA.Genome-wide miRNA-regulated networks of the writer (A), eraser (B), and reader (C). Green nodes: miRNAs, Pink nodes: genes that may be miRNA targets, and Black edges: correlations.(PDF)
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Nawaz M, Sun J, Shabbir S, Khattak WA, Ren G, Nie X, et al. A review of plants strategies to resist biotic and abiotic environmental stressors. Sci Total Environ. 2023;900:165832. doi: 10.1016/j.scitotenv.2023.165832 37524179 · doi ↗ · pubmed ↗
- 2Thiebaut F, Hemerly AS, Ferreira PCG. A Role for Epigenetic Regulation in the Adaptation and Stress Responses of Non-model Plants. Front Plant Sci. 2019;10:246. doi: 10.3389/fpls.2019.00246 30881369 PMC 6405435 · doi ↗ · pubmed ↗
- 3Deans C, Maggert KA. What do you mean, “epigenetic”?. Genetics. 2015;199(4):887–96. doi: 10.1534/genetics.114.173492 25855649 PMC 4391566 · doi ↗ · pubmed ↗
- 4Tang J, Chen S, Jia G. Detection, regulation, and functions of RNA N 6-methyladenosine modification in plants. Plant Commun. 2023;4(3):100546. doi: 10.1016/j.xplc.2023.100546 36627844 PMC 10203383 · doi ↗ · pubmed ↗
- 5Wu R, Jiang D, Wang Y, Wang X. N (6)-Methyladenosine (m(6)A) Methylation in m RNA with A Dynamic and Reversible Epigenetic Modification. Mol Biotechnol. 2016;58(7):450–9. doi: 10.1007/s 12033-016-9947-9 27179969 · doi ↗ · pubmed ↗
- 6Fang X, Li M, Yu T, Liu G, Wang J. Reversible N 6-methyladenosine of RNA: The regulatory mechanisms on gene expression and implications in physiology and pathology. Genes Dis. 2020;7(4):585–97. doi: 10.1016/j.gendis.2020.06.011 33335958 PMC 7729110 · doi ↗ · pubmed ↗
- 7Bokar JA, Shambaugh ME, Polayes D, Matera AG, Rottman FM. Purification and c DNA cloning of the Ado Met-binding subunit of the human m RNA (N 6-adenosine)-methyltransferase. RNA. 1997;3(11):1233–47. 9409616 PMC 1369564 · pubmed ↗
- 8Liu J, Yue Y, Han D, Wang X, Fu Y, Zhang L, et al. A METTL 3-METTL 14 complex mediates mammalian nuclear RNA N 6-adenosine methylation. Nat Chem Biol. 2014;10(2):93–5. doi: 10.1038/nchembio.1432 24316715 PMC 3911877 · doi ↗ · pubmed ↗