- 1Grupo de Biología Molecular y Biotecnología de Plantas, Departamento de Biología Molecular y Bioquímica, Universidad de Málaga, Málaga, Spain
- 2Departamento de Ecología y Geología, Universidad de Málaga, Málaga, Spain
- 3Integrative Molecular Biology Lab, Universidad de Málaga, Málaga, Spain
Epitranscriptome constitutes a gene expression checkpoint in all living organisms. Nitrogen is an essential element for plant growth and development that influences gene expression at different levels such as epigenome, transcriptome, proteome, and metabolome. Therefore, our hypothesis is that changes in the epitranscriptome may regulate nitrogen metabolism. In this study, epitranscriptomic modifications caused by ammonium nutrition were monitored in maritime pine roots using Oxford Nanopore Technology. Transcriptomic responses mainly affected transcripts involved in nitrogen and carbon metabolism, defense, hormone synthesis/signaling, and translation. Global detection of epitranscriptomic marks was performed to evaluate this posttranscriptional mechanism in un/treated seedlings. Increased N6-methyladenosine (m6A) deposition in the 3’-UTR was observed in response to ammonium, which seems to be correlated with poly(A) lengths and changes in the relative abundance of the corresponding proteins. The results showed that m6A deposition and its dynamics seem to be important regulators of translation under ammonium nutrition. These findings suggest that protein translation is finely regulated through epitranscriptomic marks likely by changes in mRNA poly(A) length, transcript abundance and ribosome protein composition. An integration of multiomics data suggests that the epitranscriptome modulates responses to nutritional, developmental and environmental changes through buffering, filtering, and focusing the final products of gene expression.
Introduction
Discoveries in the nascent field of molecular biology culminated in the central dogma of molecular biology during the 1970s (Crick, 1970). Currently, it is known that transcription and translation processes are not always directly linked. Multiple factors intervene in the gene response and its regulation. Two good examples of this are long noncoding RNAs (lncRNAs) (Liu et al., 2015) and microRNAs (miRNAs) (Paul et al., 2015). These types of RNA regulate important aspects of both development and response to external stimuli in plants (Liu et al., 2015; Paul et al., 2015). However, these studies represent only some aspects of overall RNA metabolism. In recent years, an increasing interest has focused on determining the biological role of RNA chemical modifications emerging as a new field of study under the term epitranscriptomics. These modifications are found in all RNA types, such as transfer RNAs (tRNAs), ribosomal RNAs (rRNAs), messenger RNAs (mRNAs) and small RNAs (RNAs) (Xiong et al., 2017). To date, more than 160 different modifications have been identified in RNA (Shen et al., 2019). In Arabidopsis thaliana, m7G, m6A, m1A, m5C, hm5C, and uridylation have been identified as modifications in mRNA (Shen et al., 2019). The marriage between classical detection techniques and high-throughput sequencing has allowed to determine N6-methyladenosine (m6A) as the most prevalent chemical modification in mRNAs, both in animals and plants (Fray and Simpson, 2015). Transcriptome-wide analysis revealed that the m6A mark in transcripts is predominantly located near the stop codon and throughout the 3′ untranslated region (UTR) (Shen et al., 2019). An m6A methylation motif (RRACH [R = A/G; H = A/U/C; A = m6A]) that is conserved between plants and other eukaryotic organisms has been described (Shen et al., 2019). m6A deposition, recognition and elimination are carried out by different proteins commonly called writers, readers, and erasers, respectively (Shen et al., 2019). Several cellular functions have been observed to be affected by m6A modification, such as mRNA stability (Wei et al., 2018) or translational efficiency (Luo et al., 2014). In addition, proper m6A deposition has been reported to be essential during Arabidopsis embryo development (Zhong et al., 2008) and to take part in biotic and abiotic plant stress responses (Martínez-Pérez et al., 2017; Anderson et al., 2018), fruit ripening (Zhou et al., 2019), flowering transition (Duan et al., 2017), leaf morphogenesis (Arribas-Hernández et al., 2018), trichome development (Bodi et al., 2012) and apical shoot meristem development (Shen et al., 2016).
Nitrogen (N) is an essential element for plant growth and development, and a key component of cellular constituents such as nucleic acids, proteins, and chlorophylls (Hawkesford et al., 2012). However, little is known about the potential role of the epitranscriptome in the regulation of N nutrition, with only a recent study on the involvement of m6A in nitrate assimilation and signaling in Arabidopsis (Hou et al., 2021). In soils, plants can take up N inorganic forms such as nitrate () or ammonium () and organic forms such as urea, peptides, and amino acids (Jia and von Wirén, 2020). Plants such as rice, tea and maritime pine prefer over as the main N source (Sasakawa and Yamamoto, 1978; Ruan et al., 2016; Ortigosa et al., 2020). In many plants, important changes in the transcriptome, proteome and metabolome have been described in relation to nitrogen nutrition and have mainly focused on the supply of and (Patterson et al., 2010; Yang et al., 2018; Ravazzolo et al., 2020). In this way, it has been observed that these N forms trigger both shared and differential responses involving different pathways and many result in phenotypic differences such as specific changes in the root system architecture (RSA) and growth (Jia and von Wirén, 2020). Therefore, it is reasonable to hypothesize that some of these described responses to N nutrition may be influenced by epitranscriptome regulatory processes.
Oxford Nanopore Technology (ONT) is a third-generation sequencing platform that is currently the only option for direct sequencing of RNA samples without the requirement of reverse transcription and amplification steps (Parker et al., 2020). These features are of great relevance in transcriptomics as they reduce sequencing biases and maintain nucleoside modifications that enable epitranscriptomic studies (Gao et al., 2021).
Research efforts of this work were focused on maritime pine (Pinus pinaster Aiton), which is a conifer typically found in the western Mediterranean region. In these areas maritime pine constitutes extensive forests being mainly located in Portugal, Spain, and France where it has been used for raw material obtention such as timber, pulp, and resin. This tree is a model species for research on functional genomics, drought resistance or nitrogen nutrition and metabolism in conifers (Sterck et al., 2022; Ávila et al., 2022). Maritime pine is a plant with a preference for ammonium over nitrate nutrition, which promotes an increase in biomass accumulation, especially in the roots where a higher number of lateral roots has been observed (Ortigosa et al., 2020; Ortigosa et al., 2022). In addition, ammonium supply promotes transcriptomic changes in several phytohormone-related transcripts such as ACC oxidase, as well as localization of phytohormones in root tips (Ortigosa et al., 2022). One of its main attractions is that it can provide an evolutionary insight into different processes studied in other model organisms, since the maritime pine is included in the group of gymnosperms whose appearance on Earth is estimated to be about 100 million years before the appearance of angiosperms (Clarke et al., 2011).
The aim of the present work is to shed light on the short-term response of maritime pine roots to nutrition, elucidating what kind of regulatory relationship exists between transcriptomics, epitranscriptomics and proteomics. For this purpose, cutting edge and commonly used omics approaches, such as comprehensive transcriptome analyses by direct RNA sequencing (DRS) using ONT, epitranscriptomic modification detection focused on m6A assisted by the ONT platform, and quantitative proteomic, were combined in the present study.
Material and methods
Plant material
Seeds from maritime pine (Pinus pinaster Aiton) from “Sierra Segura y Alcaraz” (Albacete, Spain) were provided by the Área de Recursos Genéticos Forestales of the Spanish Ministerio de Agricultura, Pesca y Alimentación. Maritime pine seed germination was carried out following the protocol described in (Cañas et al., 2006). Seedlings were grown in vermiculite in plant growth chambers under 16 h light photoperiod, a light intensity of 100 μmol m−2 s−1, constant temperature of 25 °C and watered twice a week with distilled water. One-month old maritime pine seedlings were used for the experiment. Pine seedlings were randomly subdivided into two different groups, relocated into forestall seedbeds and watered with 80 mL distilled water. After three days of acclimation, the control group was irrigated with 80 mL of water (C) and the experimental group with 80 mL of 3 mM NH4Cl. This ammonium concentration (3mM) is N-sufficient for the growth of maritime pine (Canales et al., 2010). Root samples were collected at 24 hours post-irrigation and immediately frozen in liquid N. This experiment was carried out three independent times. The adequate development of each experiment was verified through the gene expression analysis by RT-qPCR of two control genes, PpAMT1.3 and PpAMP1 following previous results (Canales et al., 2011; Castro-Rodríguez et al., 2016). The sequences of this genes can be found in Genbank under the following accession numbers: KC807909 (PpAMT1.3) and HM210085 (PpAMP1).
Metabolite profile
The metabolites for 1H-NMR analysis were extracted following the protocol previously described by Kruger et al. (2008) and according to Ortigosa et al. (2020). Extended method descriptions are in the Supplementary Methods. All data and results have been included in Supplementary Figure 1 and Supplementary Dataset 1.
Total RNA isolation
Total root RNA from maritime pine seedlings was isolated following the protocol described by Liao et al. (2004) and modified by Canales et al. (2012). The RNA concentration and purity were determined via spectrophotometry on a Nanodrop ND-1000 (Thermo Scientific, Waltham, MA, USA). Purity was determined through the 260/280 and 260/230 ratios. RNA quality was also determined in a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). The concentration was verified with a Qubit 4 Fluorometer (Invitrogen, Paisley, UK) and Qubit RNA BR, Broad-Range, Assay Kit (Cat. No. Q10210, Invitrogen, Paisley, UK).
mRNA isolation and preparation
Samples with a RIN value > 7 were selected to mRNA isolation. The poly(A)-RNA isolation was performed using Dynabeads™ mRNA Purification Kit (Cat. No. 61006, Invitrogen, Paisley, UK) following the manufacturer’s instructions. This process was carried out twice per sample to avoid rRNA contamination. poly(A)-RNA quality was determined in a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). The concentration was verified with a Qubit 4 Fluorometer (Invitrogen, Paisley, UK) and Qubit RNA HS, High Sensitivity, Assay Kit (Cat. No. Q32852, Invitrogen, Paisley, UK).
Direct RNA sequencing (DRS) and differential epitranscriptomic analysis
Nanopore libraries for DRS were prepared from 1.65 up to 2.18 µg of isolated poly(A)-RNA using the Nanopore Direct RNA Sequencing kit (Cat. No. SQK-RNA001, Oxford Nanopore Technologies, ONT, Oxford, UK) according to manufacturer’s instructions. The DRS libraries were loaded onto a R9.4 SpotON Flow Cells (Cat. No. FLO-MIN106D, Oxford Nanopore Technologies, Oxford, UK) and sequenced until complete depletion of the nanopores. Extended method descriptions are in the Supplementary Methods. The DRS data have been deposited in the NCBI's Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series with the accession number GSE174830 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE174830). Basecalling was carried out with ONT Guppy software (https://community.nanoporetech.com). The resultant reads were filtered by quality (Q>9). Read alignment was made with minimap2 software (Li, 2018) using root transcriptome of maritime pine as reference (Ortigosa et al., 2022). The alignment parameters were adjusted for DRS (-uf and -k14). Differentially expressed (DE) transcripts were identified using the edgeR package for R, the transcripts were normalised by count per million mapped reads (cpm) and filtered (2 cpm in at least 2 samples) (Robinson et al., 2010). The transcripts with False Discovery Rate < 0.05 (FDR < 0.05) were considered as differentially expressed.
ONT-DRS reads were used for de novo modification detection with the TOMBO software (Stoiber et al., 2017). The total mapped reads per base and the number of modified bases in each position were obtained using the text_output_browser_file method with the options coverage and fraction. Only the positions with at least 50 mapped reads were employed for subsequent analyses. A Fischer exact test was carried out for each transcript position to determine the differential expression of modified bases among control and supplied samples. The transcript positions with a P-value < 0.05 and an absolute logFC > 0.5 were considered as differentially modified.
Transdecoder software (https://github.com/TransDecoder/transdecoder.github.io) were used to determine modification positions and nucleobases in the transcripts of the reference transcriptome. Identification of m6A sites were carried out with the bioinformatic pipeline Nanom6A using default parameters (Gao et al., 2021). The length of poly(A) tails was determined using Nanopolish 0.11.1 software package (https://github.com/adbailey4/nanopolish) with the polya function.
Functional annotation and enrichment analyses
The transcriptome was functionally annotated with BLAST2GO (Götz et al., 2008) using DIAMOND software with blastx option (Buchfink et al., 2015) against the NCBI’s plants-nr database (NCBI Resource Coordinators, 2016). Blast results were considered valid with e-value < 1.0E-6. Singular enrichment analysis (SEA) of the GO terms was made in the AGRIGO v2.0 web tool under standard parameters using as GO term reference the whole assembled transcriptome annotation (Tian et al., 2017). Representative enriched GO was determined using REVIGO (Supek et al., 2011).
RT-qPCR
The cDNA synthesis was performed using 1 μg of total RNA and iScript™ cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA) following manufacturer’s instructions. The qPCR primers were designed following the MIQE guidelines (Bustin et al., 2009). The primers are listed in the Supplementary Table 1. qPCRs were carried out using 10 ng of cDNA and 0.4 mM of primers and 2X SsoFast™ EvaGreen® Supermix (Cat. No. 1725204, Bio-Rad, Hercules, CA, USA) in a total volume of 10 μL. Relative quantification of gene expression was performed using thermocycler CFX 384™ Real-Time System, Bio-Rad, Hercules, CA, USA). The qPCR program was: 3 min at 95°C (1 cycle), 1 s at 95°C and 5 s at 60°C (50 cycles) and a melting curve from 60 to 95°C, to generate the dissociation curve in order to confirm the specific amplification of each individual reaction. The analyses were carried out as described by Cañas et al. (2014) using the MAK3 model in the R package qpcR (Ritz and Spiess, 2008). Expression data were normalized to two reference genes, SKP1/ASK1 and SLAP that were previously tested for RT-qPCR experiments in maritime pine (Granados et al., 2016). The qPCR analyses were made with three biological replicates and three technical replicates per sample.
Validation of differential deposition of m6A by RT-qPCR
The validation of the differential deposition of m6A in the transcripts were made using the SELECT method (Xiao et al., 2018). A differential cDNA synthesis was made per each transcript with 30 ng of total RNA. qPCR determinations were made with 2 μL of cDNA. The expression level of each transcript was determined in parallel by RT-qPCR and their result were used to normalize the SELECT results. The primers are listed in the Supplementary Table 1. Extended method description can be found in the Supplementary Methods.
Differential proteomics analysis
The proteins were extracted following the protocol described by González Fernández et al. (2014). The extractions were carried out with 200 mg of sample. Protein content was determined using a commercially kit (Cat. No. 5000006, Protein Assay Dye Reagent; Bio-Rad, CA, USA) and bovine serum albumin as a standard (Bradford, 1976). Protein extracts were cleaned-up in 1D SDS-PAGE at 10% polyacrilamyde as described in Valledor and Weckwerth, 2014. Protein bands were cut off, diced, and kept in water at 4°C until digestion.
Protein digestion and nLC-MS2 analysis were carried out in the Proteomics Facility at Research Support Central Service, University of Cordoba. Nano-LC was performed in a Dionex Ultimate 3000 nano UPLC (Thermo Scientific, Waltham, MA, USA) with a C18 75 μm x 50 Acclaim Pepmam column (Thermo Scientific, Waltham, MA, USA). Eluting peptide cations were converted to gas-phase ions by nano electrospray ionization and analyzed on a Thermo Orbitrap Fusion (Q-OT-qIT, Thermo Scientific) mass spectrometer operated in positive mode. Extended method descriptions are in the Supplementary Methods.
Root transcriptome from Pinus pinaster was translated into the six open reading frames with transeq tool (Madeira et al., 2019). The output peptides chains were filtered by length, deleting those less than 50 amino acids (Romero-Rodríguez et al., 2014). To reduce the redundancy of proteins in the database, CD-HIT-EST with a 99% identity filter was used (Fu et al., 2012). The raw data were processed using Proteome Discoverer (version 2.3, Thermo Scientific). MS2 spectra were searched with SEQUEST engine against the reference proteome database. In silico peptide lists were created using the followings settings: trypsin digestion, a maximum of two missed internal cleavage sites per peptide, precursor mass tolerance of 10 ppm and fragment mass tolerance of 0.75 Da per fragment ions. Only peptides with a high confidence (FDR ≤ 0.01) and minimum XCorr of 2 were selected. The identified proteins were filtered by a minimum of two different peptides. Minora algorithm was used to determine relative quantification. The protein identification with redundancy is considered by Proteome Discoverer and SEQUEST software. Proteins sharing peptides were grouped and all those groups without a unique peptide were removed. Only proteins detected in 5 of the 6 samples were considered for subsequent analyses. The resultant proteins were annotated using BLAST with blastp option against NCBI’s non-redundant database (Camacho et al., 2008) and BLAST2GO software (Götz et al., 2008). Proteins with the same sequence annotation and quantification profile across the samples was considered the same protein selecting the protein with the longest sequence. Normalized data from Proteome Discoverer software were used for a differential expression analysis using the edgeR package for R (Robinson et al., 2010). Only the proteins with P-value < 0.05 were considered as differentially expressed. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository (Vizcaino et al., 2013) with the Data identifier PXD025331 and 10.6019/PXD025331.
Results
Direct RNA sequencing (DRS)
The global sequencing results are shown in Table 1. The mean read sizes were between 908 bp to 1059 bp (Figure 1A). The longest reads ranged from 10298 bp to 14299 bp. Differential expression analyses identified 350 differentially expressed (DE) transcripts. From the total DE transcripts obtained, 106 were upregulated and 244 were downregulated (Figure 1B; Supplementary Dataset 2). Singular enrichment analysis (SEA) was performed individually for each gene expression regulation (up- and downregulated) to classify the biological functions under nutrition. The SEA global results are shown in Supplementary Dataset 3. Upregulated transcripts were significantly enriched with GO terms for Biological Process (BP) such as ammonia assimilation cycle, protein glutathionylation, and developmental growth; for Cellular Component (CC) such as chloroplast stroma and cytosolic ribosome; and for Molecular Function (MF) such as glutamate synthase (NADH) activity, and phospholipase activity (Figure 1C; Supplementary Figure 2). The downregulated transcripts were enriched in BP terms such as protein folding and ethylene-activated signaling pathway and MF terms such as transcription coactivator activity, mRNA binding, and chaperone binding. A more detailed study of the upregulated transcriptomic response revealed that PpGS1b (pp_68481) and PpNADH-GOGAT (pp_238920) were upregulated. Interestingly, transcripts for genes involved in defense-related response were upregulated: PpAMP1 (pp_58005, pp_58008), different class IV chitinases (pp_239593, pp_239598, pp_239600, pp_117809), different splicing coding forms of patatin-like protein 2 (pp_71017, pp_71018, pp_71019, pp_71020, pp_71022), a PR-1 pathogenesis-related protein (pp_87427), a defensin coding transcript (pp_92119) and an RPW8 domain-containing protein (pp_142311) (Supplementary Dataset 2). Cell wall-related transcripts were also upregulated, such as those encoding expansin-A18 (pp_134987), nonclassical arabinogalactan protein 30 (pp_66323), probable prolyl 4-hydroxylase 4 (pp_235715), xyloglucan:xyloglucosyl transferase (pp_68519) and different forms of xyloglucan endotransglucosylase/hydrolase (pp_66707, pp_66708, pp_68517), was also observed. In contrast, the downregulation of different transcription factors (TFs) was observed (Supplementary Dataset 2), such as ethylene response factors (ERFs) (pp_58625, pp_58626, pp_86737, pp_96228, pp_96234), a trihelix transcription factor (pp_59947), and a MYB coding transcript (pp_202778). The repression of different splicing forms of an auxin-repressed protein/dormancy-auxin associated protein coding transcript (pp_58457, pp_58458, pp_58459, pp_58461, pp_58462, pp_58463), the SnRK1 regulator FCS-like zinc finger transcript (pp_202096), and transcripts encoding carbohydrate metabolism enzymes such as pyruvate decarboxylase (pp_78343, pp_123611, pp_126347) and sucrose synthase (pp_144843) was also observed. The differential expressions observed in DRS was confirmed by RT-qPCR analyses of seven transcripts including SAM synthase, MYB5, PFK, AMP1, NADH-GOGAT, GS1b, and ASPG. The results showed the same trend between the DRS and RT-qPCR results for all DE transcripts, with some differences in the logFC values (Figure 1D).
Figure 1 Main DRS transcriptomics results. (A) Size distribution of reads obtained through DRS-ONT sequencing. (B) MA plot containing all detected transcripts. Blue points correspond to significantly downregulated transcripts in maritime pine roots at 24 h after the treatment with 3 mM ammonium. Red points correspond to significantly upregulated transcripts in maritime pine roots at 24 h after the treatment with 3 mM ammonium. (C) Significant GO terms from significant DE transcripts after a SEA analysis. (D) Comparison between DRS and RT-qPCR results from different transcripts with significant differential expression in the DRS analysis.
Differential DRS epitranscriptomics
To determine the differential epitranscriptomic marks, DRS results were explored to identify the chemically modified nucleosides in the mRNAs using Tombo software. The number of putative modified nucleosides in the transcripts ranged from 722,655 in sample C1 to 4,063,005 in sample N3 (Supplementary Datasets 4–9). After differential deposition analysis, 513 significant differentially modified nucleosides were obtained in 283 transcripts (Figure 2A; Supplementary Dataset 10). There were 221 undermodified positions in 161 transcripts and 292 overmodified positions in 184 transcripts. Among them, 58 transcripts had significant over- and undermodified positions, including PpGS1b (pp_68474) and translationally-controlled tumor protein (TCTP, pp_72505) (Supplementary Dataset 11). The percentage of modifications for each kind of nucleoside was similar for adenosine, guanosine, and uridine (28%, 29% and 25%) but lower for cytidine (17%), without any obvious difference between sample conditions (Figure 2B). When the global set of modification ratios and transcript amounts were compared (Figure 2C), significant and negative correlations were found to be stronger under supply (-0.36) than under the control (-0.30) (Figure 2C). Modification ratios of individual nucleosides compared to transcript amounts were also similar (Supplementary Figure 3).
Figure 2 Epitranscriptomics modification results from Tombo software analysis. (A) Volcano plot of the detected epitranscriptomics modifications. Yellow points correspond to modifications with P-values lower than 0.05 but with an absolute logFC value lower than 0.5. Blue points correspond to modifications with an absolute logFC value higher than 0.5 but with P-values higher than 0.05. Red points correspond to significant DE modifications with an absolute logFC value higher than 0.5 and P-values lower than 0.05. (B) Percentage of detected modifications per nucleobase. Yellow columns correspond to control samples. Purple columns correspond to 3 mM ammonium treated samples. (C) Scatter plot and correlation between detected modifications using Tombo software and transcript amounts determined by DRS sequencing. (D) Significant GO terms from significant DE modifications after a SEA analysis. (E) Comparison between DRS and SELECT results from different transcript modified positions with significant differential modification in the Tombo analysis.
The functions of the transcripts with DE modifications were analyzed using SEA (Figure 2D; Supplementary Figure 4 and Supplementary Dataset 12). A total of 116 significant GO terms were obtained. At the BP level, the main functions were related to the terms ribosome biogenesis, translation, protein ubiquitination, and glutamine metabolic process. Similarly, the ribosome term was the main function at the CC level. Finally, the terms translation factor activity, ubiquitin protein ligase binding, and endopeptidase inhibitor activity were the main functions at the MF level.
The differential modifications were verified using SELECT. This qPCR-based technique was initially designed to determine differential deposition of m6A in total RNA mixtures. However, nucleoside modifications putatively detected in cytidine and uridine by Tombo were also detected with SELECT (Figure 2E). The obtained results for each position showed a similar trend between the differential epitranscriptomic results from SELECT and Tombo (Figure 2E).
m6A identification
The bioinformatic pipeline Nanom6A was used to specifically identify m6A modifications in the RRACH sites from DRS data. Statistical analysis identified 176 nucleosides with significant (P-value < 0.05) differential deposition of m6A, but only 29 were considered as having a |logFC| > 0.5 (Supplementary Dataset 13). The distribution of m6A sites along the full-length transcripts showed a higher accumulation of marks over two-thirds of the relative length in control samples while a higher accumulation of marks in the final portion of the transcripts was observed in -treated samples (Figure 3A). A more detailed distribution study showed that m6A sites were more abundant in the 5’-UTR and coding (CDS) regions under the control condition, while m6A sites tended to be more abundant in the 3’-UTR regions of transcripts isolated from -treated seedlings (Figure 3B). The m6A frequency was higher in the CDS, mainly in the central portion, and at the beginning of the 3’-UTR regions, but it was lower at the transcript ends (Figures 3A, B).
Figure 3 m6A identification results using the bioinformatics pipeline Nanom6A. (A) Distribution of the identified m6A all along the transcripts. (B) Distribution of the identified m6A along each part of the transcripts: 5’UTR, CDS and 3’UTR respectively. (C) Percentage of identified m6A for each different variant (kmer) of the consensus sequence RRACH for adenine methylation. (D) Scatter plot and correlations between m6A modification ratio and transcript amount from DRS sequencing. (E) Poly(A) length of the transcripts identified in the DRS sequencing. (F) Scatter plot and correlations between poly(A) length and m6A modification ratio. (G) Scatter plot and correlations between poly(A) length and transcript amount from DRS sequencing.
The most predominant RRACH sequence was AAACA (>15%), while GGACC was the least abundant (< 5%) (Figure 3C). The comparison between m6A modification ratios and transcript amounts showed significant negative correlations for the control (-0.42) and conditions (-0.46) (Figure 3D). The lengths of the poly(A) tails were determined from DRS data (Figure 3E). As expected, most of the poly(A) tails had a size between 40-250 nt with similar means in both conditions, 91.01 nt and 91.8 nt. The poly(A) tail lengths had significant positive correlations with m6A ratios in both conditions (0.2) (Figure 3F). Significant but lower positive correlations were found when poly(A) tail lengths were compared with global and nucleoside modification ratios obtained with Tombo (Supplementary Figure 5). Finally, the poly(A) tail lengths and transcript amounts exhibited significant negative correlations for the control (-0.17) and conditions (-0.18) (Figure 3G).
Differential proteomics
A total of 2,385 proteins were identified in the shotgun proteomics analysis (Supplementary Dataset 14). Among the identified proteins, 114 were differentially regulated by : 38 were more abundant, while 76 were less represented (Figure 4A; Supplementary Dataset 14). To elucidate the biological roles of the identified proteins, SEA analyses were performed (Figure 4B; Supplementary Figure 6 and Supplementary Dataset 15). The upregulated proteins showed as representative BP terms cell redox homeostasis, protein complex assembly, and translation. At CC level, ribosome, nucleolus and chloroplast were the enriched terms. Finally, structural constituent of ribosome and enzyme regulator activity were the significant MF terms. Among the downregulated proteins, the representative enriched functions were ribosome assembly and translation among the BP terms. At the CC level, the significant terms most interesting was ribosome. Structural constituent of ribosome and GTPase activity were the enriched MF terms.
Figure 4 Differential expression proteomics results. (A) Volcano plot of the identified proteins. Yellow points correspond to modifications with P-values lower than 0.05 but with an absolute logFC value lower than 0.5. Red points correspond to significant DE modifications with an absolute logFC value higher than 0.5 and P-values lower than 0.05. (B) Significant GO terms from significant DE proteins after a SEA analysis. (C) Scatter plot and correlations between m6A modification ratios and protein amounts. (D) Scatter plot and correlations between transcript amounts from DRS sequencing and protein amounts.
The putative relationship between the m6A modification ratio and protein abundance was determined through Pearson correlation analysis (Figure 4C). In -treated seedlings, a significant positive correlation was observed (0.18), while there was no correlation in control seedlings. The same effect was observed between nucleoside modification ratios from Tombo and protein amounts (Supplementary Figure 7). In addition, similar correlation results were obtained when poly(A) tail lengths and protein amounts were compared, and only -treated roots had a significant positive correlation (0.12) (Figure 4D).
Integration of data from omics approaches
The results of the different omics approaches employed in the present work were integrated to explore possible regulatory steps in response to supply in maritime pine. The global data showed 30 different elements/genes with significant results based on at least two approaches; 14 of them were significant in DRS and epitranscriptomics, 5 in DRS and proteomics, 9 in epitranscriptomics and proteomics, and 2 in all approaches (Figure 5A; Supplementary Dataset 16). Altogether, the genes/proteins identified were involved in N metabolism (PpASPG, PpGS1b, alanine-glyoxylate aminotransferase and isocitrate dehydrogenase), defense (PpAMP1, a ginkbilobin and a chitinase), oxidative stress response (alcohol dehydrogenase, aldehyde dehydrogenase, peroxidase and glutaredoxin), translation (ribosomal proteins and elongation factors) and RNA binding (cold shock proteins and a BURP domain protein RD22). Among them, the presence of the 1-aminocyclopropane-1-carboxylate (ACC) oxidase, which had a putative modification at the position 462 on the contig, must be highlighted. This epitranscriptomic mark was overexpressed (logFC 0.84) under nutrition, while the ACC oxidase protein was underexpressed (logFC -2.69). Similarly, changes in transcript level, protein level and transcript modification ratio can show opposite trends with variations between genes. For the aldehyde dehydrogenase, the transcript expression and modification were overexpressed (0.58 and 0.73, respectively), but protein expression was underexpressed (logFC -1.20). However, for the CSP/GRP (pp_211512 and pp_211516), transcript accumulation was underexpressed (-0.43 and -0.7) and the epitranscriptomic modifications and protein expression were overexpressed (0.88-0.73 and 2.11 respectively).
Figure 5 Common results among different omics approaches. (A) Genes with significant expressions in more than one omics approach. Red boxes correspond to upregulated transcripts and proteins. Blue boxes correspond to downregulated transcripts and proteins. Yellow boxes correspond to significant differentially modified positions in the transcripts. (B) Venn diagram of significant GO terms in the different omics approaches. (C) Heatmap of significant GO terms shared by at least two omics approaches. (D) Heatmap of genes that produce proteins involved in ribosome and translation. Ribosomal protein equivalences obtained from Martinez-Seidel et al. (2020). Red boxes correspond to upregulated transcripts and proteins. Blue boxes correspond to downregulated transcripts and proteins. Yellow boxes correspond to significant differentially modified positions in the transcripts.
The comparison between the significant GO terms in the omics approaches revealed that 94 of the 287 GO terms were shared (Figures 5B, C; Supplementary Dataset 17). Interestingly, 20 (7%) of them were common to the three sets of results. The epitranscriptomics and proteomics comparison included the highest number of GO terms (59, 20.6%), and the transcriptomics and proteomics comparison included the lowest number, with only 5 (1.7%). The most representative GO terms common to the three omics data were ribosome, structural constituent of ribosome and translation. Between the transcriptomics and epitranscriptomics data, mRNA binding and glutamine metabolic process were the main GO terms. In the transcriptomics and proteomics comparison, oxoacid metabolic process was the main GO term. Finally, between epitranscriptomics and proteomics the main GO terms were small ribosomal subunit, GTPase activity, ribosome assembly, defense response and response to external stimulus. According to these functional results, several transcripts coding for eukaryotic ribosomal proteins (38) and translation factors (8) had differential epitranscriptomic marks based on the Tombo results (Figure 5D).
Discussion
The response of maritime pine roots to nutrition has been studied from a multiomics perspective that includes direct RNA sequencing using the ONT platform, which has allowed a global epitranscriptomics analysis. Although N is an essential nutrient for proper plant growth and development (Xu et al., 2012), little is known about the role of epitranscriptomic marks in the regulation of gene expression in response to N nutrition. The results in the present work highlight the importance of epitranscriptomic marks in the regulation of gene expression.
Epitranscriptome changes in response to NH4+nutrition
Correlation of global epitranscriptomics results with transcript abundance (Figure 2C), especially m6A modifications (Figure 3D), are consistent with those from previous works in Populus (Gao et al., 2021) and Arabidopsis (Luo et al., 2014; Wan et al., 2015) and support the role of m6A in mRNA turnover, as previously described in mammals (Wang et al., 2014; Li et al., 2019). Additionally, m6A identification revealed that in the roots of maritime pine, AAACA and AAACU were the most abundant RRACH 5-mer motifs (Figure 3C), as previously reported in Arabidopsis, maize, and poplar (Wan et al., 2015; Miao et al., 2020; Parker et al., 2020; Gao et al., 2021). These findings strongly suggest conservation of the RNA m6A methylation machinery in plants. However, pine m6A distribution, with similar levels of enrichment in the middle of the CDS and at beginning of the 3’-UTR (Figures 3A, B), slightly differs from the obtained results in angiosperms, where a greater enrichment of m6A was observed in the 3’-UTR (Miao et al., 2020; Parker et al., 2020; Gao et al., 2021). This discrepancy may be attributed to the lack of a well-curated reference genome for maritime pine. Interestingly, the distribution of m6A sites in pine transcripts is dynamically regulated by increasing the m6A deposition in the 3’-UTR in response to , which seems to be correlated with poly(A) length and changes in protein abundance (Figures 5C, D), as well as the identification of transcripts with both up- and down-regulated epitranscriptomic marks (Supplementary Dataset 11). This observation agrees with previous results in other eukaryotic organisms (Meyer et al., 2012; Anderson et al., 2018).
Regarding poly(A) tail length, a previous work in Caenorhabditis elegans reported that highly expressed transcripts contained a relatively short and well-defined poly(A) tail (Lima et al., 2017), which has been related to translational efficiency (Wu et al., 2020). In yeast, it was demonstrated that poly(A) tail lengths also have negative correlations with transcript abundance (Tudek et al., 2021). Our results reveal that increased poly(A) tail lengths correlated with higher m6A and lower transcript abundances (Figures 3F, G). These findings are consistent with the effect of m6A modifications on transcript levels (Figure 3D), which follow the same trend as previously published results (Luo et al., 2014; Wan et al., 2015; Parker et al., 2020; Gao et al., 2021). The integration of all these results with the positive correlations of m6A ratios and poly(A) tail length with protein amounts (Figures 4C, D) suggests that m6A could affect the translation efficiency of the differentially modified transcripts, as previously described in mammals (Meyer, 2019). This evidence could also explain, at least in part, the lack of a relationship between the transcriptomic and proteomic data.
Epitranscriptomic marks, mainly m6A, seem to modulate protein synthesis through mRNA stability and modify translation processes. This is in line with a previous hypothesis considering that initial responses caused by environmental factors are modulated at the level of final products of gene expression to maintain cellular homeostasis (Cañas et al., 2015). Thus, the variable transcriptomic response can be controlled through intermediate steps, such as epitranscriptomic regulation, to generate an appropriate, quick, and stable cellular response.
Carbon and nitrogen metabolism
As expected, N assimilation was induced by , as shown by the significant increase in GS1b and NADH-GOGAT transcripts (Figure 6; Supplementary Dataset 3), consistent with previous transcriptomic reports (Li et al., 2017; Sun et al., 2017). However, no overrepresentation of GS1b and NADH-GOGAT was observed at the proteomic level (Supplementary Dataset 14), even though a trend towards higher GS activity was observed in the presence of (Supplementary Figure 8). Similar proteomics results for nutrition have been described before (Marino et al., 2016; Coleto et al., 2019), as well as the lack of correlation between relative expression levels of GS transcripts and GS activity and Gln/Glu levels. To explain these discrepancies, it was hypothesized that the response to N nutrition could be regulated through a posttranscriptional (translational or posttranslational) mechanism (Ortigosa et al., 2020). In this way, several nucleosides with differential epitranscriptomic marks, including m6A, have been identified in GS1b transcripts, which could explain the differences between transcriptomics and proteomics data (Supplementary Datasets 10 and 13). Accordingly, in Arabidopsis, m6A can regulate the alternative polyadenylation and transcript abundance of the GLN1;1 and GLN1;3 genes during nitrate assimilation and signaling (Hou et al., 2021).
Figure 6 Schematic representation of functional results obtained through the omics integrative approach. triggers carbon and nitrogen metabolism, ethylene biosynthesis and signaling and translation responses. Geometric red and blue forms indicate upregulation and downregulation respectively of transcripts, RNA nucleoside modifications, proteins, and metabolites.
Additionally, the results of this study suggest the existence of a carbon flux through glycolysis and the TCA cycle to provide carbon skeletons for N assimilation (de la Peña et al., 2019). Interestingly, the fermentation pathway seems to be repressed, as suggested by the downregulation of pyruvate decarboxylase (PDC) transcripts and the decrease in alcohol dehydrogenase (ADH) and aldehyde dehydrogenase (ALDH) protein levels (Supplementary Datasets 2 and 14). It is well known that plants under low oxygen conditions regulate their metabolism, inducing the fermentation pathway in which PDC and ADH activities are essential, which results in pyruvate consumption, the production of ethanol and the concomitant oxidation of NADH to NAD+ (Mithran et al., 2014). The consumption of soluble sugars (sucrose and D-glucose) (Supplementary Figure 1 and Supplementary Dataset 1), the repression of the fermentative pathway and the overrepresentation of TIM proteins (inner membrane translocase proteins) (Supplementary Dataset 14) highlight the need for organic acids from TCA cycle to assimilate and produce energy for plant growth. In addition, our results suggest that the biosynthesis of organic acids could be regulated through the epitranscriptome according to the differentially overmodified sequence found in the coding transcript for the pyruvate kinase enzyme (pp_236058), which could imply an additional regulation mechanism together with those previously described (Wulfert et al., 2020).
Ethylene
Transcriptomic studies carried out in rice described that under excess , ethylene (ET) could be one of the major regulatory molecules in roots (Sun et al., 2017). Additionally, Li et al. (2013) described in Arabidopsis that shoot-supplied promoted ET biosynthesis only in shoots, resulting in a reduction in the lateral root formation process due to auxin transporter 1 (AUX1) repression. These data might indicate that ET biosynthesis is involved in detrimental -related phenotypes, such as a reduction in the number of lateral roots (LRs). Some reports have described that the application of inhibitors of the ET biosynthesis relieved symptoms of toxicity affecting RSA (Li et al., 2013). However, previous results in maritime pine seedlings indicated that promotes root growth (Ortigosa et al., 2020) and a wide repression of the ACC oxidase in the root apex during early transcriptome response (Ortigosa et al., 2022).
A comparison of the root transcriptomic response and differential proteomics results revealed that promotes a decrease in ET-related transcripts and proteins. Several ET-responsive TFs were downregulated at the transcriptomic level, while proteomic results revealed that ACC oxidase was downregulated (Supplementary Datasets 2 and 14). ACC oxidase is the enzyme responsible for the final step in the biosynthesis of ET in plants (John et al., 1999). Furthermore, ACC oxidase transcripts exhibited a nucleoside putative modification in the 3’-UTR (position 462) that was differentially increased under treatment (Supplementary Dataset 16). Since the levels of ACC oxidase transcripts had no significant changes, it is possible that this epitranscriptomic mark could be involved in the translational regulation of these transcripts, therefore affecting the ET levels of the organ. Interestingly, it is reasonable to think that this kind of response could be related to conifer preference for . Additional technical approaches are required to identify the exact modification and validate its biological relevance in the regulation of ACC oxidase and ET biosynthesis. Finally, the relationship between ET and N nutrition is remarked by the differential expression (negative) of a trihelix family TF that has a very high similarity degree to AT3G54390 in Arabidopsis (Figure 6). This TF in Arabidopsis has been shown to interact with the promoters of genes involved in N metabolism and ET biosynthesis including glutamine synthetase (AT1G66200), nitrite reductase (AT2G15620), ET-responsive transcription factor (ERF012; AT1G21910) and 1-aminocyclopropane-1-carboxylate synthase (AT5G65800) (Gaudinier et al., 2018).
Translation and growth
The results obtained also indicate a reconfiguration of the ribosomal proteins and elongation factors at all biological levels, although this was more visible in epitranscriptomics and proteomics results (Figures 5D, 6; Supplementary Datasets 2, 10, 14 and 17). This process is clearly related to the correlation between m6A ratios and protein amounts (Figure 4C), suggesting that nutrition promotes general translation activation to support the root growth of maritime pine seedlings. This kind of effect on the ribosomal protein composition and proteins involved in translation has been previously observed under different conditions, including plant mineral nutrition (Wang et al., 2013; Prinsi and Espen, 2018; Ferretti and Karbstein, 2019). The data suggest that changes in the profiles of ribosomal proteins can be significantly mediated by epitranscriptomic modifications of their transcripts. These posttranscriptional changes can influence ribosome function and performance. One example is the RACK1 proteins, which are components of the small subunit of the ribosome and are involved in translation quality control promoting the endonucleolytic cleavage of nonstop mRNA (Ikeuchi and Inada, 2016). These proteins in plants are involved in the response to phytohormones, including ET, and in the regulation of growth and development (Wang et al., 2019). Interestingly, the mRNA of a RACK1 gene in maritime pine roots has several differential epitranscriptomic marks (Supplementary Datasets 10 and 11). Although it was not significant, RACK1 protein expression was repressed under conditions, suggesting a change in the quality control mechanisms and even in the integration of the plant signaling pathways. This result agrees with the negative regulation of ACC oxidase and repression of several ERFs (Supplementary Datasets 2 and 14), since ET signaling affects gene translation in different ways (Merchante et al., 2015).
Additionally, related to the regulation of translation and root growth, among the transcriptomic results, the repression of an FCS-like zinc finger was observed (Supplementary Dataset 2). FCS proteins interact with SnRK1 to mediate the interaction of SnRK1 with other proteins (Jamsheer et al., 2018). SnRK1 is a kinase that phosphorylates RAPTOR, a regulatory element of the TOR complex, repressing TOR activity (Fu et al., 2020). Considering the ribosomal protein changes and the increase in root growth observed under nutrition in maritime pine seedlings (Ortigosa et al., 2020), it can be hypothesized a decrease of SnRK1 activity over the TOR complex in these conditions. The TOR complex integrates different developmental and environmental signals, including nutritional status, promoting ribosome biogenesis, translation, and plant growth (Henriques et al., 2022). Finally, a transcript coding for a TCTP protein had several differential modifications (Supplementary Dataset 11). It is well established in plants that TCTP mRNA is synthetized in shoots and transported to roots, where it is translated and promotes lateral root formation (Branco and Masle, 2019). Interestingly, transport of this mRNA is mediated through epitranscriptomic marks (Yang et al., 2019). This suggests additional mechanisms controlling growth and development in response to N nutrition in maritime pine roots.
Therefore, the results presented in this work suggest that protein translation and growth are finely regulated through epitranscriptomic marks, including m6A, to acquire an optimum response to N supply (Figure 6). More research efforts are required to corroborate this hypothesis and to investigate whether ET could act as a modulator of the integrated response observed due to its effects on the translation process.
Conclusions
Finally, and according to our results in maritime pine, triggers a root systemic response at short-term that mainly involved changes in key pathways such as carbon and nitrogen metabolism, ET signaling pathway, translation and root growth (Figure 6). Interestingly, ET-related response observed was different from that of previously reported in other tolerant plants such as rice (Sun et al., 2017) and supports previous findings in maritime pine (Ortigosa et al., 2022). Additionally, the obtained results strongly suggest that protein translation is finely regulated through the epitranscriptome affecting mRNA turnover and probably the ribosome performance. Although the gene transcription is very reactive to development and environmental changes the processes that regulate mRNA metabolism and translation, such as epitranscriptomics marks, can modulate the response buffering, filtering, and focusing the final products of the gene expression. In this case, the epitranscriptomic regulation seems directed to acquire a proper response to promote root growth under supplementation.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
FO have performed the experiments; CL-F has performed the bioinformatic analyses; JP-C and FRC have performed the statistical data analysis; FO and RC have wrote the manuscript; CA and FMC made additional contributions and edited the manuscript. FO and RC have planned and designed the research. RC, CA, and FMC were the responsible of the funding acquisition. All authors contributed to the article and approved the submitted version.
Funding
This work was funded by Spanish Ministerio de Economía y Competitividad, grants numbers BIO2015-73512-JIN MINECO/AEI/FEDER, UE, BIO2015-69285-R and RTI2018-094041-B-I00, and by Spanish Ministerio de Ciencia e Innovación, grant numbers PID2021-122641NB-C21 and PID2021-125040OB-I00. FO was supported by a grant from the Universidad de Málaga (Programa Operativo de Empleo Juvenil vía SNJG, UMAJI11, FEDER, FSE, Junta de Andalucía) and funds from the research group BIO-114.
Acknowledgments
We thank Dr. Bertrand Hirel for his helpful comments on the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.1102044/full#supplementary-material
References
Anderson, S. J., Kramer, M. C., Gosai, S. J., Yu, X., Vandivier, L. E., Nelson., A. D., et al. (2018). N6-methyladenosine inhibits local ribonucleolytic cleavage to stabilize mRNAs in Arabidopsis. Cell Rep. 25, 1146–1157. doi: 10.1016/j.celrep.2018.10.020
Arribas-Hernández, L., Bressendorff, S., Hansen, M. H., Poulsen, C., Erdmann, S., Brodersen, P. (2018). An m6A-YTH module controls developmental timing and morphogenesis in Arabidopsis. Plant Cell 30, 952–967. doi: 10.1105/tpc.17.00833
Ávila, C., Cañas, R. A., de la Torre, F. N., Pascual, M. B., Castro-Rodríguez, V., Cantón, F. R., et al. (2022). “Functional genomics of Mediterranean pines,” in The pine genomes. compendium of plant genomes. Ed. de la Torre, A. R. (Cham: Springer). doi: 10.1007/978-3-030-93390-6_9
Bodi, Z., Zhong, S., Mehra, S., Song, J., Graham, N., Li, H., et al. (2012). Adenosine methylation in Arabidopsis mRNA is associated with the 3’ end and reduced levels cause developmental defects. Front. Plant Sci. 3. doi: 10.3389/fpls.2012.00048
Bradford, M. M. (1976). A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal. Biochem. 72, 248–254. doi: 10.1006/abio.1976.9999
Branco, R., Masle, J. (2019). Systemic signalling through translationally controlled tumour protein controls lateral root formation in Arabidopsis. J. Exp. Bot. 70, 3927–3940. doi: 10.1093/jxb/erz204
Buchfink, B., Xie, C., Huson, D. H. (2015). Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60. doi: 10.1038/nmeth.3176
Bustin, S. A., Benes, V., Garson, J. A., Hellemans, J., Huggett, J., Kubista, M., et al. (2009). The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55, 611–622. doi: 10.1373/clinchem.2008.112797
Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. (2008). BLAST+: architecture and applications. BMC Bioinf. 10, 421. doi: 10.1186/1471-2105-10-421
Canales, J., Avila, C., Canovas, F. M. (2011). A maritime pine antimicrobial peptide involved in ammonium nutrition. Plant Cell Environ. 34, 1443–1453. doi: 10.1111/j.1365-3040.2011.02343.x
Canales, J., Flores-Monterrosso, A., Rueda, M., Avila, C., Cánovas, F. M. (2010). Identification of genes regulated by ammonium availability in the roots of maritime pine trees. Amino Acids 39, 991–1001. doi: 10.1007/s00726-010-0483-9
Canales, J., Rueda-López, M., Craven-Bartle, B., Avila, C., Cánovas, F. M. (2012). Novel insights into regulation of asparagine synthetase in conifers. Front. Plant Sci. 3. doi: 10.3389/fpls.2012.00100
Cañas, R. A., Canales, J., Gómez-Maldonado, J., Ávila, C., Cánovas, F. M. (2014). Transcriptome analysis in maritime pine using laser capture microdissection and 454 pyrosequencing. Tree Physiol. 34, 1278–1288. doi: 10.1093/treephys/tpt113
Cañas, R. A., Canales, J., Muñoz-Hernández, C., Granados, J. M., Ávila, C., García-Martín, M. L., et al. (2015). Understanding developmental and adaptive cues in pine through metabolite profiling and co-expression network analysis. J. Exp. Bot. 66, 3113–3127. doi: 10.1093/jxb/erv118
Cañas, R. A., de la Torre, F., Cánovas, F. M., Cantón, F. R. (2006). High levels of asparagine synthetase in hypocotyls of pine seedlings suggest a role of the enzyme in re-allocation of seed-stored nitrogen. Planta 224, 83–95. doi: 10.1007/s00425-005-0196-6
Castro-Rodríguez, V., Assaf-Casals, I., Pérez-Tienda, J., Fan, X., Avila, C., Miller, A., et al. (2016). Deciphering the molecular basis of ammonium uptake and transport in maritime pine. Plant Cell Environ. 39, 1669–1682. doi: 10.1111/pce.12692
Clarke, J. T., Warnock, R. C., Donoghue, P. C. (2011). Establishing a time-scale for plant evolution. New Phytol. 192, 266–301. doi: 10.1111/j.1469-8137.2011.03794.x
Coleto, I., Vega-Mas, I., Glauser, G., González-Moro, M. B., Marino, D., Ariz, I. (2019). New insights on Arabidopsis thaliana root adaption to ammonium nutrition by the use of a quantitative proteomic approach. Int. J. Mol. Sci. 20, 814. doi: 10.3390/ijms20040814
de la Peña, M., González-Moro, M. B., Marino, D. (2019). Providing carbon skeletons to sustain amide synthesis in roots underlines the suitability of Brachypodium distachyon for the study of ammonium stress in cereals. AoB Plants 11, plz029. doi: 10.1093/aobpla/plz029
Duan, H. C., Wei, L. H., Zhang, C., Wang, Y., Chen, L., Lu, Z., et al. (2017). ALKBH10B is an RNA N6-methyladenosine demethylase affecting Arabidopsis floral transition. Plant Cell 29, 2995–3011. doi: 10.1105/tpc.16.00912
Edgar, R., Domrachev, M., Lash, A. E. (2002). Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 30, 207–210. doi: 10.1093/nar/30.1.207
Ferretti, M. B., Karbstein, K. (2019). Does functional specialization of ribosomes really exist? RNA 25, 521–538. doi: 10.1261/rna.069823.118
Fray, R. G., Simpson, G. G. (2015). The arabidopsis epitranscriptome. Curr. Opin. Plant Biol. 27, 17–21. doi: 10.1016/j.pbi.2015.05.015
Fu, L., Niu, B., Zhu, Z., Wu, S., Li, W. (2012). CD-HIT: accelerated for clustering the next generation sequencing data. Bioinformatics 28, 3150–3152. doi: 10.1093/bioinformatics/bts565
Fu, L., Wang, P., Xiong, Y. (2020). Target of rapamycin signaling in plant stress responses. Plant Physiol. 182, 1613–1623. doi: 10.1104/pp.19.01214
Gao, Y., Liu, X., Wu, B., Wang, H., Xi, F., Kohnen, M. V., et al. (2021). Quantitative profiling of N6-methyladenosine at single-base resolution in stem-differentiating xylem of Populus trichocarpa using nanopore direct RNA sequencing. Genome Biol. 22, 22. doi: 10.1186/s13059-020-02241-7
Gaudinier, A., Rodriguez-Medina, J., Zhang, L., Olson, A., Liseron-Monfils, C., Bågman, A. M., et al. (2018). Transcriptional regulation of nitrogen-associated metabolism and growth. Nature 563, 259–264. doi: 10.1038/s41586-018-0656-3
González Fernández, R., Redondo, I., Jorrín-Novo, J. V. (2014). “Making a protein extract from plant pathogenic fungi for gel- and LC-based proteomics,” in Plant proteomics: Methods and protocols, 2nd edition, (Totowa, NJ: Humana) vol. 1072. , 93–109. doi: 10.1007/978-1-62703-631-3_8
Götz, S., García-Gómez, J. M., Terol, J., Williams, T. D., Nagaraj, S. H., Nueda, M. J., et al. (2008). High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 36, 3420–3435. doi: 10.1093/nar/gkn176
Granados, J. M., Avila, C., Cánovas, F. M., Cañas, R. A. (2016). Selection and testing of reference genes for accurate RT-qPCR in adult needles and seedlings of maritime pine. Tree Genet. Genomes 12, 60. doi: 10.1007/s11295-016-1018-7
Hawkesford, M., Horst, W., Kichey, T., Lambers, H., Schjoerring, J., Møller, I. S., et al. (2012). “Functions of macronutrients,” in Marschner's mineral nutrition of higher plants (Waltham, MA, USA: Academic Press), 135–189. doi: 10.1016/B978-0-12-384905-2.00006-6
Henriques, R., Calderan-Rodrigues, M. J., Crespo, J. L., Baena-González, E., Caldana, C. (2022). Growing of the TOR world. J. Exp. Bot. 73, 6987–6992. doi: 10.1093/jxb/erac401
Hou, Y., Sun, J., Wu, B., Gao, Y., Nie, H., Nie, Z., et al. (2021). CPSF30-l-mediated recognition of mRNA m6A modification controls alternative polyadenylation of nitrate signaling-related gene transcripts in Arabidopsis. Mol. Plant 14, 688–699. doi: 10.1016/j.molp.2021.01.013
Ikeuchi, K., Inada, T. (2016). Ribosome-associated Asc1/RACK1 is required for endonucleolytic cleavage induced by stalled ribosome at the 3' end of nonstop mRNA. Sci. Rep. 6, 28234. doi: 10.1038/srep28234
Jamsheer, K. M., Sharma, M., Singh, D., Mannully, C. T., Jindal, S., Shukla, B. N., et al. (2018). FCS-like zinc finger 6 and 10 repress SnRK1 signalling in Arabidopsis. Plant J. 94, 232–245. doi: 10.1111/tpj.13854
Jia, Z., von Wirén, N. (2020). Signaling pathways underlying nitrogen-dependent changes in root system architecture: from model to crop species. J. Exp. Bot. 71, 4393–4404. doi: 10.1093/jxb/eraa033
John, P., Reynolds, E. A., Prescott, A. G., Bauchot, A. D. (1999). “ACC oxidase in the biosynthesis of ethylene,” in Biology and biotechnology of the plant hormone ethylene II. Eds. Kanellis, A. K., Chang, C., Klee, H., Bleecker, A. B., Pech, J. C., Grierson, D. (Dordrecht: Springer). doi: 10.1007/978-94-011-4453-7_1
Kruger, N. J., Troncoso-Ponce, M. A., Ratcliffe, R. G. (2008). 1H NMR metabolite fingerprinting and metabolomic analysis of perchloric acid extracts from plant tissues. Nat. Protoc. 3, 1001–1012. doi: 10.1038/nprot.2008.64
Li, H. (2018). ). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094–3100. doi: 10.1093/bioinformatics/bty191
Liao, Z., Chen, M., Guo, L., Gong, Y., Tang, F., Sun, X., et al. (2004). Rapid isolation of high-quality total RNA from taxus and ginkgo. Prep. Biochem. Biotech. 34, 209–214. doi: 10.1081/PB-200026790
Li, G., Li, B., Dong, G., Feng, X., Shi, W. (2013). Ammonium-induced shoot ethylene production is associated with the inhibition of lateral root formation in Arabidopsis. J. Exp. Bot. 64, 1413–1425. doi: 10.1093/jxb/ert019
Lima, S. A., Chipman, L. B., Nicholson, A. L., Chen, Y. H., Yee, B. A., Yeo, G. W., et al. (2017). Short poly(A) tails are a conserved feature of highly expressed genes. Nat. Struc. Mol. Biol. 24, 1057–1063. doi: 10.1038/nsmb.3499
Liu, J., Wang, H., Chua, N. H. (2015). Long noncoding RNA transcriptome of plants. Plant Biotech. J. 13, 319–328. doi: 10.1111/pbi.12336
Li, W., Xiang, F., Zhong, M., Zhou, L., Liu, H., Li, S., et al. (2017). Transcriptome and metabolite analysis identifies nitrogen utilization genes in tea plant (Camellia sinensis). Sci. Rep. 7, 1–12. doi: 10.1038/s41598-017-01949-0
Li, F., Yi, Y., Miao, Y., Long, W., Long, T., Chen, S., et al. (2019). N6-methyladenosine modulates nonsense-mediated mRNA decay in human glioblastoma. Cancer Res. 79, 5785–5798. doi: 10.1158/0008-5472.CAN-18-2868
Luo, G. Z., MacQueen, A., Zheng, G., Duan, H., Dore, L. C., Lu, Z., et al. (2014). Unique features of the m6A methylome in Arabidopsis thaliana. Nat. Commun. 5, 1–8. doi: 10.1038/ncomms6630
Madeira, F., Park, Y. M., Lee, J., Buso, N., Gur, T., Madhusoodanan, N., et al. (2019). The EMBL-EBI search and sequence analysis tools APIs in 2019. Nucleic Acids Res. 47, 636–641. doi: 10.1093/nar/gkz268
Marino, D., Ariz, I., Lasa, B., Santamaría, E., Fernández-Irigoyen, J., González-Murua, C., et al. (2016). Quantitative proteomics reveals the importance of nitrogen source to control glucosinolate metabolism in Arabidopsis thaliana and Brassica oleracea. J. Exp. Bot. 67, 3313–3323. doi: 10.1093/jxb/erw147
Martínez-Pérez, M., Aparicio, F., López-Gresa, M. P., Bellés, J. M., Sánchez-Navarro, J. A., Pallás, V. (2017). Arabidopsis m6A demethylase activity modulates viral infection of a plant virus and the m6A abundance in its genomic RNAs. P. Natl. Acad. Sci. U.S.A. 114, 10755–10760. doi: 10.1073/pnas.1703139114
Martinez-Seidel, F., Beine-Golovchuk, O., Hsieh, Y. C., Kopka, J. (2020). Systematic review of plant ribosome heterogeneity and specialization. Front. Plant Sci. 11. doi: 10.3389/fpls.2020.00948
Merchante, C., Brumos, J., Yun, J., Hu, Q., Spencer, K. R., Enríquez, P., et al. (2015). Gene-specific translation regulation mediated by the hormone-signaling molecule EIN2. Cell 163, 684–697. doi: 10.1016/j.cell.2015.09.036
Meyer, K. D. (2019). m6A-mediated translation regulation. BBA Gene Regul. Mech. 1862, 301–309. doi: 10.1016/j.bbagrm.2018.10.006
Meyer, K. D., Saletore, Y., Zumbo, P., Elemento, O., Mason, C. E., Jaffrey., S. R. (2012). Comprehensive analysis of mRNA methylation reveals enrichment in 3′ UTRs and near stop codons. Cell 149, 1635–1646. doi: 10.1016/j.cell.2012.05.003
Miao, Z., Zhang, T., Qi, Y., Song, J., Han, Z., Ma, C. (2020). Evolution of the RNA N6-methyladenosine methylome mediated by genomic duplication. Plant Physiol. 182, 345–360. doi: 10.1104/pp.19.00323
Mithran, M., Paparelli, E., Novi, G., Perata, P., Loreti, E. (2014). Analysis of the role of the pyruvate decarboxylase gene family in Arabidopsis thaliana under low-oxygen conditions. Plant Biol. 16, 28–34. doi: 10.1111/plb.12005
NCBI Resource Coordinators (2016). Database resources of the national center for biotechnology information. Nucleic Acids Res. 44, 7–19. doi: 10.1093/nar/gkv1290
Ortigosa, F., Lobato-Fernández, C., Shikano, H., Ávila, C., Taira, S., Cánovas, F. M., et al. (2022). Ammonium regulates the development of pine roots through hormonal crosstalk and differential expression of transcription factors in the apex. Plant Cell Environ. 45, 915–935. doi: 10.1111/pce.14214
Ortigosa, F., Valderrama-Martín, J. M., Urbano-Gámez, J. A., García-Martín, M. L., Ávila, C., Cánovas, F. M., et al. (2020). Inorganic nitrogen form determines nutrient allocation and metabolic responses in maritime pine seedlings. Plants 9, 481. doi: 10.3390/plants9040481
Parker, M. T., Knop, K., Sherwood, A. V., Schurch, N. J., Mackinnon, K., Gould, P. D., et al. (2020). Nanopore direct RNA sequencing maps the complexity of Arabidopsis mRNA processing and m6A modification. eLife 9, e49658. doi: 10.7554/eLife.49658
Patterson, K., Cakmak, T., Cooper, A., Lager, I., Rasmusson, A. G., Escobar, M. A. (2010). Distinct signalling pathways and transcriptome response signatures differentiate ammonium- and nitrate-supplied plants. Plant Cell Environ. 33, 1486–1501. doi: 10.1111/j.1365-3040.2010.02158.x
Paul, S., Datta, S. K., Datta, K. (2015). miRNA regulation of nutrient homeostasis in plants. Front. Plant Sci. 6. doi: 10.3389/fpls.2015.00232
Prinsi, B., Espen, L. (2018). Time-course of metabolic and proteomic responses to different nitrate/ammonium availabilities in roots and leaves of maize. Int. J. Mol. Sci. 19, 2202. doi: 10.3390/ijms19082202
Ravazzolo, L., Trevisan, S., Forestan, C., Varotto, S., Sut, S., Dall’Acqua, S., et al. (2020). Nitrate and ammonium affect the overall maize response to nitrogen availability by triggering specific and common transcriptional signatures in roots. Int. J. Mol. Sci. 21, 686. doi: 10.3390/ijms21020686
Ritz, C., Spiess, A. N. (2008). qpcR: an r package for sigmoidal model selection in quantitative real-time polymerase chain reaction analysis. Bioinformatics 24, 1549–1551. doi: 10.1093/bioinformatics/btn227
Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Romero-Rodríguez, M. C., Pascual, J., Valledor, L., Jorrín-Novo, J. V. (2014). Improving the quality of protein identification in non-model species. characterization of Quercus ilex seed and Pinus radiata needle proteomes by using SEQUEST and custom databases. J. Proteomics 105, 85–91. doi: 10.1016/j.jprot.2014.01.027
Ruan, L., Wei, K., Wang, L., Cheng, H., Zhang, F., Wu, L., et al. (2016). Characteristics of and fluxes in tea (Camellia sinensis) roots measured by scanning ion-selective electrode technique. Sci. Rep. 6, 1–8. doi: 10.1038/srep38370
Sasakawa, H., Yamamoto, Y. (1978). Comparison of the uptake of nitrate and ammonium by rice seedlings: influences of light, temperature, oxygen concentration, exogenous sucrose, and metabolic inhibitors. Plant Physiol. 62, 665–669. doi: 10.1104/pp.62.4.665
Shen, L., Liang, Z., Gu, X., Chen, Y., Teo, Z. W., Hou, X., et al. (2016). N6-methyladenosine RNA modification regulates shoot stem cell fate in Arabidopsis. Dev. Cell 38, 186–200. doi: 10.1016/j.devcel.2016.06.008
Shen, L., Liang, Z., Wong, C. E., Yu, H. (2019). Messenger RNA modifications in plants. Trends Plant Sci. 24, 328–341. doi: 10.1016/j.tplants.2019.01.005
Sterck, L., de María, N., Cañas, R. A., de Miguel, M., Perdiguero, P., Raffin, A., et al. (2022). “Maritime pine genomics in focus,” in The pine genomes. compendium of plant genomes. Ed. de la Torre, A. R. (Cham: Springer). doi: 10.1007/978-3-030-93390-6_5
Stoiber, M. H., Quick, J., Egan, R., Lee, J. E., Celniker, S. E., Neely, R., et al. (2017). De novo identification of DNA modifications enabled by genome-guided nanopore signal processing. BioRxiv. doi: 10.1101/094672
Sun, L., Di, D., Li, G., Kronzucker, H. J., Shi, W. (2017). ). spatio-temporal dynamics in global rice gene expression (Oryza sativa l.) in response to high ammonium stress. J. Plant Physiol. 212, 94–104. doi: 10.1016/j.jplph.2017.02.006
Supek, F., Bošnjak, M., Škunca, N., Šmuc, T. (2011). REVIGO summarizes and visualizes long lists of gene ontology terms. PloS One 6, e21800. doi: 10.1371/journal.pone.0021800
Tian, T., Yue, L., Hengyu, Y., Qi, Y., Xin, Y., Zhou, D., et al. (2017). agriGO v2.0: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 45, 122–129. doi: 10.1093/nar/gkx382
Tudek, A., Krawczyk, P. S., Mroczek, S., Tomecki, R., Turtola, M., Matylla-Kulińska, K., et al. (2021). Global view on the metabolism of RNA poly(A) tails in yeast saccharomyces cerevisiae. Nat. Commun. 12, 4951. doi: 10.1038/s41467-021-25251-w
Valledor, L., Weckwerth, W. (2014). An improved detergent-compatible gel-fractionation LC-LTQ-Orbitrap-MS workflow for plant and microbial proteomics. Methods Mol. Biol. (Clifton N.J.) 1072, 347–358. doi: 10.1007/978-1-62703-631-3_25
Vizcaino, J. A., Cote, R. G., Csordas, A., Dianes, J. A., Fabregat, A., Foster, J. M., et al. (2013). The proteomics identifications (PRIDE) database and associated tools: status in 2013. Nucleic Acids Res. 41, 1063–1069. doi: 10.1093/nar/gks1262
Wang, J., Lan, P., Gao, H., Zheng, L., Li, W., Schmidt, W. (2013). Expression changes of ribosomal proteins in phosphate- and iron-deficient Arabidopsis roots predict stress-specific alterations in ribosome composition. BMC Genomics 14, 783. doi: 10.1186/1471-2164-14-783
Wang, Y., Li, Y., Toth, J. I., Petroski, M. D., Zhang, Z., Zhao, J. C. (2014). N6-methyladenosine modification destabilizes developmental regulators in embryonic stem cells. Nat. Cell Biol. 16, 191–198. doi: 10.1038/ncb2902
Wang, W., Wang, X., Wang, X., Ahmed, S., Hussain, S., Zhang, N., et al. (2019). Integration of RACK1 and ethylene signaling regulates plant growth and development in Arabidopsis. Plant Sci. 280, 31–40. doi: 10.1016/j.plantsci.2018.11.009
Wan, Y., Tang, K., Zhang, D., Xie, S., Zhu, X., Wang, Z., et al. (2015). Transcriptome-wide high-throughput deep m6A-seq reveals unique differential m6A methylation patterns between three organs in Arabidopsis thaliana. Genome Biol. 16, 272. doi: 10.1186/s13059-015-0839-2
Wei, L. H., Song, P., Wang, Y., Lu, Z., Tang, Q., Yu, Q., et al. (2018). The m6A reader ECT2 controls trichome morphology by affecting mRNA stability in Arabidopsis. Plant Cell 30, 968–985. doi: 10.1105/tpc.17.00934
Wulfert, S., Schilasky, S., Krueger, S. (2020). Transcriptional and biochemical characterization of cytosolic pyruvate kinases in Arabidopsis thaliana. Plants 9, 353. doi: 10.3390/plants9030353
Wu, X., Wang, J., Wu, X., Hong, Y., Li, Q. Q. (2020). Heat shock responsive gene expression modulated by mRNA poly(A) tail length. Front. Plant Sci. 11. doi: 10.3389/fpls.2020.01255
Xiao, Y., Wang, Y., Tang, Q., Wei, L., Zhang, X., Jia, G. (2018). An elongation-and ligation-based qPCR amplification method for the radiolabeling-free detection of locus-specific N6-methyladenosine modification. Angew. Chem. Int. Ed. 57, 15995–16000. doi: 10.1002/anie.201807942
Xiong, X., Yi, C., Peng, J. (2017). Epitranscriptomics: toward a better understanding of RNA modifications. Genomics Proteomics Bioinf. 15, 147. doi: 10.1016/j.gpb.2017.03.003
Xu, G., Fan, X., Miller, A. J. (2012). Plant nitrogen assimilation and use efficiency. Annu. Rev. Plant Biol. 63, 153–182. doi: 10.1146/annurev-arplant-042811-105532
Yang, L., Perrera, V., Saplaoura, E., Apelt, F., Bahin, M., Kramdi, A., et al. (2019). m5C methylation guides systemic transport of messenger RNA over graft junctions in plants. Curr. Biol. 29, 2465–2476. doi: 10.1016/j.cub.2019.06.042
Yang, Y., Wang, F., Wan, Q., Ruan, J. (2018). Transcriptome analysis using RNA-seq revealed the effects of nitrogen form on major secondary metabolite biosynthesis in tea (Camellia sinensis) plants. Acta Physiol. Plant 40, 127. doi: 10.1007/s11738-018-2701-0
Zhong, S., Li, H., Bodi, Z., Button, J., Vespa, L., Herzog, M., et al. (2008). MTA is an Arabidopsis messenger RNA adenosine methylase and interacts with a homolog of a sex-specific splicing factor. Plant Cell 20, 1278–1288. doi: 10.1105/tpc.108.058883
Keywords: nitrogen nutrition, ammonium, epitranscriptomics, ONT sequencing, Pinus pinaster, translation
Citation: Ortigosa F, Lobato-Fernández C, Pérez-Claros JA, Cantón FR, Ávila C, Cánovas FM and Cañas RA (2022) Epitranscriptome changes triggered by ammonium nutrition regulate the proteome response of maritime pine roots. Front. Plant Sci. 13:1102044. doi: 10.3389/fpls.2022.1102044
Received: 18 November 2022; Accepted: 08 December 2022;
Published: 22 December 2022.
Edited by:
Guangjie Li, Institute of Soil Science, Chinese Academy of Sciences, ChinaReviewed by:
Yufang Lu, Institute of Soil Science (CAS), Nanjing, ChinaRunlai Hang, University of California, Riverside, Riverside, United States
Copyright © 2022 Ortigosa, Lobato-Fernández, Pérez-Claros, Cantón, Ávila, Cánovas and Cañas. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Rafael A. Cañas, rcanas@uma.es
†These authors have contributed equally to this work and share first authorship