Transcriptional, metabolic, physiological and developmental responses of switchgrass to phosphorus limitation.

Knowing how switchgrass (Panicum virgatum L.) responds and adapts to phosphorus (P)-limitation will aid efforts to optimize P acquisition and use in this species for sustainable biomass production. This integrative study investigated the impacts of mild, moderate, and severe P-stress on genome transcription and whole-plant metabolism, physiology and development in switchgrass. P-limitation reduced overall plant growth, increased root/shoot ratio, increased root branching at moderate P-stress, and decreased root diameter with increased density and length of root hairs at severe P-stress. RNA-seq analysis revealed thousands of genes that were differentially expressed under moderate and severe P-stress in roots and/or shoots compared to P-replete plants, with many stress-induced genes involved in transcriptional and other forms of regulation, primary and secondary metabolism, transport, and other processes involved in P-acquisition and homeostasis. Amongst the latter were multiple miRNA399 genes and putative targets of these. Metabolite profiling showed that levels of most sugars and sugar alcohols decreased with increasing P stress, while organic and amino acids increased under mild and moderate P-stress in shoots and roots, although this trend reversed under severe P-stress, especially in shoots.

INTRODUCTION

Phosphorus (pan> class="Chemical">P) is essential for plant growth and development. Low P‐availability limits plant growth in most soils because soluble, plant‐available inorganic phosphate (Pi) readily forms insoluble complexes with Fe3+ and Al3+ in acidic soils and with Ca2+ in alkaline soils (Bertrand, Holloway, Armstrong, & McLaughlin, 2003; Hinsinger, 2001). Plants have evolved a variety of developmental, physiological, molecular, and symbiotic strategies to optimize acquisition and utilization of P for growth (Lambers et al., 2011; Lambers, Finnegan, et al., 2015; O'Rourke et al., 2013). Traits associated with adaptation and acclimation to P‐limitation include: remodeling of root architecture and development of more and longer root hairs and lateral roots (Lambers, Clode, et al., 2015; Lynch, 2011); increased root to shoot ratio (Lynch, 1995); induction of phosphate scavenging and recycling enzymes (Ding, Guo, Kupper, & McNear, 2016); release of carboxylates (Tomasi et al., 2009); P homeostasis by adjustment of major P‐pools (Veneklaas et al., 2012); and engagement of specific signal transduction pathways and transcriptional regulatory networks (Misson et al., 2005; Morcuende et al., 2007; Secco et al., 2013). Qualitative and quantitative differences in P‐strategies exist between and within species (Lambers, Clode, et al., 2015; Pang, Ryan, Lambers, & Siddique, 2018). Therefore, in order to optimize P efficiency (uptake and utilization) in specific crop species, it is important to investigate the responses and adaptations to P‐limitation and the underlying mechanisms in that species and the natural variation that is available for breeding.

Switchgrass, native to the pan> class="Chemical">North American tallgrass prairies, is a perennial plant with water‐efficient C4 photosynthesis that was targeted for development as a bioenergy crop (Casler et al., 2011; Meyer et al., 2014). Switchgrass exhibits high biomass production potential, relatively low input requirements, and is adapted to much of the eastern half of the United States, including areas considered marginal for food‐crop production (Casler et al., 2011; Gopalakrishnan, Cristina Negri, & Snyder, 2011; King et al., 2013; Meyer et al., 2014). Some research suggests that biomass‐to‐energy schemes using marginal lands would provide substantial ecosystem services, particularly in terms of carbon sequestration and other environmental benefits (Bhardwaj et al., 2011; Gelfand et al., 2013). Data on switchgrass production on marginal sites are limited. Previous research has shown that switchgrass biomass yields respond to nitrogen fertilizer rates of up to 168 kg/ha, depending on ecotype and location (Guretzky, Biermacher, Cook, Kering, & Mosali, 2011; Muir, Sanderson, Ocumpaugh, Jones, & Reed, 2001; Sanderson, Read, & Reed, 1999). In soils with low plant‐available P, application of 45 kg P ha−1 increased biomass yield by up to 17% (Kering, Biermacher, Butler, Mosali, & Guretzky, 2012).

Switchgrass has been subjected to a genome sequencing effort (Casler et al., 2011) as well as tranpan>scriptome anpan>apan> class="Chemical">lyses, using Expressed Sequence Tags (ESTs), Affymetrix oligonucleotides arrays, and RNA‐seq (Meyer et al., 2014; Sharma et al., 2012; Yang et al., 2016; Zhang et al., 2013). Transcriptome analyses have identified thousands of genes associated with drought stress (Meyer et al., 2014) and leaf senescence in switchgrass (Yang et al., 2016), but transcriptional responses to P limitation have not been reported. Likewise, metabolic responses of switchgrass to P deficiency remain unknown, although advanced technologies are available (Luo, Wang, Sun, & Wang, 2017; Sanchez, Siahpoosh, Roessner, Udvardi, & Kopka, 2008). We characterized the physiological and developmental responses of switchgrass to P‐limitation and explored underlying transcriptional and metabolic responses in shoots and roots. Results and insights are presented here.

MATERIALS AND METHODS

Plant growth and sample collection

Switchgrass (cv. pan> class="Chemical">Alamo) seeds were surface‐sterilized and germinated as described in Ramamoorthy and Kumar (2012). Seedlings of comparable size were transplanted into growth cones (Stuewe & Sons, Inc. Tangent, OR; 6.9 cm in diameter × 25.4 cm in depth) filled with 40% all‐purpose sand (The QUIKRETE Companies, Atlanta, GA), 40% Q‐ROK silica sand (American Plants & Services, Inc. Oklahoma City, OK) and 20% fine (2–4 mm particle size) perlite (Ambient Minerals Inc. Benton, AR). Plants were watered daily to field capacity during a 4‐week growth period with a nutrient solution containing 600 μM (P‐replete control), 200 μM (mild stress), 60 μM (moderate stress) or 20 μM Pi (severe stress), supplied as KH2PO4. Each nutrient solution also contained the following macro‐ and micronutrients: 6 mM KNO3, 1 mM CaCl2, 0.25 mM MgSO4·7H2O, 25 μM NaCl, 12.5 μM H3BO3, 1 μM ZnSO4·7H2O, 1 μM MnSO4·H2O, 0.05 μM CuSO4·5H2O, 0.25 μM (NH4)6Mo7O24·4H2O, and 10 μM Fe(III)‐EDTA. 1 mM 2‐(N‐morpholino) ethanesulphonic acid (MES) was added as a buffering agent. The pH of the solution was adjusted to 6.5. Cones were randomly arranged on a greenhouse bench and plants grown in a 14 hr light /10 hr dark period at temperatures of 29°C (light) and 24°C (dark). Plants used for the analyses had no noticeable limitation of pot size or water, as roots of plants grown at 600 μM P‐supply did not bunch up at the bottom of the containers, nor did the containers dry out between watering. Observed plant responses were thus attributable to changes in P status rather than other stresses. Plant samples for metabolic and transcriptional analysis were rinsed in Milli‐Q water, blotted dry and immediately frozen in liquid nitrogen before storage at −80°C until use. Plant samples for biomass determination were placed in pre‐weighed paper bags and dried at 65°C until weights were constant.

Metabolite and ion measurements

The same pooled plant materials used for RNA extraction were used to determine ion contents, membranpan>e pan> class="Chemical">lipids, and primary metabolites. For cation/anion analysis, 20 mg frozen, powdered shoot or root material was transferred into a reaction tube, suspended in 500 μl of Milli‐Q water by vortexing, and centrifuged for 1 min at 13,000g. The supernatant was transferred to a new reaction tube and the pellet re‐extracted using the same procedure. First and second supernatants were combined before filtering through a 0.45 μm membrane. Ion separation was performed with a Dionex ICS‐5000 plus ion chromatographic system (Thermo Fisher Scientific) equipped with a conductivity detector. Standard curves were prepared using serial dilutions of Dionex Seven Anion Standard II and Dionex Six Cation II Standard (Thermo Fisher Scientific). Quantification was done using Thermo Scientific Chromeleon 7.2 SR4 Software. Lipid extraction followed the procedure described by Vu et al. (2014).

Polar metabolites were extracted from powdered, freeze‐dried materials using the protocol described by Broeckling et al. (2005) with minor modifications. In brief, approximately 10 mg of sample were extracted in chloroform containing 10 μg/ml pan> class="Chemical">docosanol (internal standard). Samples were incubated at 50°C for 1 hr, followed by centrifugation at 4°C at 3500 rpm for 45 min to separate the phases. Water with 25 μg/ml ribitol (internal standard) was added and samples were incubated at 50°C for 45 min. The aqueous layer was dried in a vacuum centrifuge at ambient temperature. Dried polar extracts were suspended in pyridine and derivatized with methoxyamine hydrochloride followed by trimethylsilyl derivatization using N‐methyl‐N‐trimethylsilyltrifluoroacetamide. Polar metabolites were analyzed using an Agilent 7890B gas chromatograph coupled to a high‐resolution Agilent 7,200 accurate‐mass quadrupole time‐of‐flight mass spectrometry (GC/Q‐TOF/MS). Data were analyzed using Agilent MassHunter software and MSDIAL version 3.82 (Tsugawa et al., 2015) and verified with MRMPROBS (version 2.36). Metabolites were identified through spectral and retention time matching with authentic compounds using an in‐house custom library augmented with a library from Riken. Compounds were normalized relative to the internal standard. Compounds in P‐stress treatments were assessed by relative changes in abundance as response ratios, in comparison to a control condition. Twofold change and larger with t test probability limits of p < .05 between P stress and control were considered reliable and significant (Ortmayr, Charwat, Kasper, Hann, & Koellensperger, 2017).

Data analysis

To test for metabolite differences between treatments, metabolite abundances were Hellinger‐transformed (Ramette, 2007) and principal component analysis (pan> class="Chemical">PCA) was performed with PC‐ORD v6.08 (McCune & Mefford, 1999). Data on plant biomass, root system architecture (primary seminal root length, total root length, root surface area), ion content, lipid and metabolite abundance were subjected to statistical analysis by one‐way ANOVA, using JMP software (SAS institute Inc., Cary, NC). Significance was defined as a probability level of the student's t test at p ≤ .05. Total root length and root surface area were performed using WinRHIZO software (Arsenault, Poulcur, Messier, & Guay, 1995) based on scanned root images using standard parameters (Regent Instruments Inc., Ontario, CA).

Transcriptome analysis

Total RNA was extracted from samples pulverized in liquid pan> class="Chemical">nitrogen, using a cetyltrimethylammonium bromide (CTAB) and LiCl method (Chang, Puryear, & Cairney, 1993). RNA was treated with DNase I to remove genomic DNA (Yang et al., 2016). RNA quality and integrity was determined using an Agilent 2100 Bioanalyzer with Plant RNA Nano chip assay (Agilent; http://www.agilent.com). For RNA sequencing (RNA‐seq) library synthesis, three biological replicates per P treatment were sequenced using an Illumina HiSeq 2500 instrument, as described previously (Serba et al., 2015). Gene expression levels were normalized by calculating reads per kilo base of transcript per million fragments mapped (RPKM). HISAT2 2.0.5 and Stringtie 1.2.4 were used to identify gene transcripts after mapping to the reference genome (Phytozome Panicum virgatum v4.1) (Kim, Langmead, & Salzberg, 2015; Pertea et al., 2015). Differentially expressed gene transcripts (DEGs) were analyzed using Cuffdiff software (Trapnell et al., 2013) and filtered with log2FC (fold change) ≥1 or ≤−1 and settled on twofold changes as cut‐off, padj ≤ .05 and RPKM ≥3. DEGs were also employed for gene ontology (GO) enrichment analysis using the PlantRegMap database with default parameters (Jin et al., 2016).

Selected transcriptional changes observed by Rpan> class="Chemical">NA Seq analysis were validated by real‐time quantitative PCR (RT‐qPCR) analysis. Three housekeeping genes (Pavir.9NG844500, UBQ1; Pavir.9NG259500, UBQ6; Pavir.8NG010400, CYP5) were used as references (Gimeno, Eattock, Van Deynze, & Blumwald, 2014; Zhao, Zhou, & Meng, 2020). Primers utilized are listed in Table S2d. Log2FC changes were calculated by means of the ΔΔCT method (Livak & Schmittgen, 2001).

MiRNA and ‐like lncRNA identification, and target prediction

Precursor sequences of P‐responsive miRNAs, that is, pre‐pan> class="Chemical">miR399s, pre‐miR2111s and pre‐miR827s, from several plant species were retrieved from miRBase, and used to identify expressed transcripts harboring miR399, miR2111 or miR827 sequence motifs. RNAfold (Lorenz et al., 2011) was used to predict minimum free energy RNA hairpin structures characteristic for miRNA precursors. psRNATarget (Dai, Zhuang, & Zhao, 2018) was used to predict target genes for mature miRNAs using default parameters. INDUCED BY PHOSPHATE STARVATION1 (IPS1)‐like, long non‐coding RNAs (lncRNAs) were identified by modifying the psRNATarget scoring scheme to reveal central mismatches in the matching sequence with miR399. The multiple sequence alignment, neighbor‐joining tree, and the miR399 sequence logo were generated using Geneious software.

RESULTS

Effect of phosphate supply on switchgrass growth and development

To establish P regimes that limit pan> class="Species">switchgrass growth, plants were watered daily with nutrient solutions containing 600, 200 60 or 20 μM Pi. Compared to P‐replete plants supplied with 600 μM Pi, plant growth was mildly, moderately or severely impacted by 200, 60 and 20 μM Pi, respectively (Figure 1a). After 28 days of growth, mild P‐stress resulted in a 35% decrease in shoot biomass, compared to the control, although no significant difference in root biomass was found between these two sets of plants (Figure 1b). Moderate or severe P‐stress resulted in further large decreases in shoot (77 and 97%, respectively) and in root biomasses (60 and 93%, respectively) (Figure 1b). Root to shoot ratio increased from 0.28 in control to 0.77 in severe P stress (data not shown). Increasing P‐stress led to successively larger reductions in Pi content especially in roots, where the Pi concentration dropped to approximately 6% under severe P limitation compared to control roots (Figure 1c). In shoots of severely P‐stressed plants, Pi concentration also dropped significantly but was still about approximately 30% of that in P‐replete control shoots (Figure 1c).

FIGURE 1

Effect of different P treatments on switchgrass growth. (a) Representative shoots of plants that received 600, 200, 60 or 20 μM Pi during growth for 28 days. Numbers in the photograph indicate average height in centimeters ± SE (n = 5). (b) Shoot and root biomass in grams dry weight. (c) Inorganic phosphate concentration (mg per g fresh weight) for each treatment. Data in (b) represent averages ± SE (n = 5). Data in (c) represent averages ± SE (n = 3). Significance is indicated with different letters within same plant organ (p ≤ .05)

Effect of different P treatments on switchgrass growth. (a) Representative shoots of planpan>ts that received 600, 200, 60 or 20 μM Pi during growth for 28 days. pan> class="Chemical">Numbers in the photograph indicate average height in centimeters ± SE (n = 5). (b) Shoot and root biomass in grams dry weight. (c) Inorganic phosphate concentration (mg per g fresh weight) for each treatment. Data in (b) represent averages ± SE (n = 5). Data in (c) represent averages ± SE (n = 3). Significance is indicated with different letters within same plant organ (p ≤ .05)

Primary seminal root growth (length) was not significantly different between severely or moderately P‐limited plants and control plants (Figure 2a, c, d, e), although plants under mild P‐stress displayed a small but significant (10%) decrease in primary root length (Figure 2b, e). On the other hand, total root length (i.e., the sum of the lengths of all roots) and root surface area were substantially higher in plants subjected to mild P‐stress compared to control plants or plants that received even less Pi (Figure 2f, g). More lateral roots and greater secondary branching were observed for planpan>ts subjected to mild Pi stress (data not shown). Planpan>ts under moderate anpan>d severe P‐stress had greatly reduced total root lengths (83% anpan>d 91% reduction, respan> class="Chemical">pectively) and root surface area (83 and 85% reduction, respectively) compared to control plants (Figure 2f, g). Plants subjected to moderate and, especially, severe P‐stress produced more and longer root hairs than control plants (Figure 2j, k), while very few or no root hairs were observed on plants exposed to mild or no P‐stress (Figure 2h, i). In addition, average root diameter decreased in moderately and severely P‐stressed roots compared to control and mildly‐stressed roots (Figure 2h‐k).

FIGURE 2

Effect of P treatment on root growth and morphology. (a–d) Representative whole root morphology of plants treated with: (a) 600; (b) 200; (c) 60 or (d) 20 μM Pi. Scale bars are 5 cm. (e) Primary root length; (f) Total root length, that is, the sum of the lengths of all roots; (g) Root surface area, that is, the sum of the surface area of all roots. Data in (e–g) represent averages ± SE (n = 5). Significance is indicated with different letters within same plant organ (p ≤ .05). (h–k) Representative root hair density/length phenotype of plants treated with: (h) 600; (i) 200; (j) 60 and (k) 20 μM phosphate. Scale bars in (h) through (k) are 500 μm (μM)

Effect of P treatment on root growth and morphology. (a–d) Representative whole root morphology of plants treated with: (a) 600; (b) 200; (c) 60 or (d) 20 μM Pi. Scale bars are 5 cm. (e) Primary root length; (f) Total root length, that is, the sum of the lengths of all roots; (g) Root surface area, that is, the sum of the surface area of all roots. Data in (e–g) represent averages ± SE (n = 5). Significance is indicated with different letters within same plant organ (p ≤ .05). (h–k) Representative root hair density/length phenotypan> class="Chemical">pe of plants treated with: (h) 600; (i) 200; (j) 60 and (k) 20 μM phosphate. Scale bars in (h) through (k) are 500 μm (μM)

Metabolic responses to P‐stress

Ion, lipid anpan>d primary metabolite levels were determined to uncover metabolic responses of pan> class="Species">switchgrass to P limitation. Levels (mg/kg tissue) of all anions and cations tested changed significantly in shoots and/or roots in response to P‐stress. Pi levels decreased substantially in both shoots (up to threefold) and roots (>10‐fold) under P‐stress, while nitrate and sulphate levels tended to increase in these organs and conditions (Figure 3a). In contrast, chloride and malate levels tended to decrease under P‐stress, possibly in response to the increasing levels of nitrate in order to maintain charge balance. In this context, potassium (K) and sodium (Na) levels decreased significantly in shoots and/or roots under P‐stress.

FIGURE 3

Effect of Pi concentration on ion and metabolite abundances. (a) Free (aqueous extraction) cation and anion concentrations (average ± SE; n = 3) in plants treated with 600, 200, 60 or 20 μM Pi. * indicates not detected. (b) Effect of decreasing Pi supply on the abundance of major membrane lipids and triacylglycerol. Relative abundances (average ± SE; n = 3), as deduced from peak heights, are given for shoots (top panel) and roots (bottom panel) of plants treated with 600 (blue bars), 200 (orange bars), 60 (grey bars) or 20 μM (gold bars) phosphate. Major membrane lipids shown are the phospholipids phosphatidylcholine (PC), phosphatidylethanolamine (PE) and phosphatidylserine (PS), the glycolipids monogalactosyl diacylglycerol (MGDG) and digalactosyl diacylglycerol (DGDG), the chloroplast sulpholipid sulphoquinovosyl diacylglycerol (SQDG), as well as the neutral storage lipid triacylglycerol (TAG). (c) Principal component analysis (PCA) of metabolite profiles based on relative abundance (average ± SE; n = 3). All abundances were normalized to a range between 0 and 1 using Hellinger Transformation

Effect of Pi concentration on ion and metabolite abundances. (a) Free (aqueous extraction) cation and anion concentrations (average ± SE; n = 3) in plants treated with 600, 200, 60 or 20 μM Pi. * indicates not detected. (b) Effect of decreasing Pi supply on the abundance of major membrane lipids anpan>d pan> class="Chemical">triacylglycerol. Relative abundances (average ± SE; n = 3), as deduced from peak heights, are given for shoots (top panel) and roots (bottom panel) of plants treated with 600 (blue bars), 200 (orange bars), 60 (grey bars) or 20 μM (gold bars) phosphate. Major membrane lipids shown are the phospholipids phosphatidylcholine (PC), phosphatidylethanolamine (PE) and phosphatidylserine (PS), the glycolipids monogalactosyl diacylglycerol (MGDG) and digalactosyl diacylglycerol (DGDG), the chloroplast sulpholipid sulphoquinovosyl diacylglycerol (SQDG), as well as the neutral storage lipid triacylglycerol (TAG). (c) Principal component analysis (PCA) of metabolite profiles based on relative abundance (average ± SE; n = 3). All abundances were normalized to a range between 0 and 1 using Hellinger Transformation

P limitation affected lipid composition of pan> class="Species">switchgrass shoots and roots (Figure 3b). Generally, phospholipid levels, including phosphatidylcholine (PC), phosphatidylethanolamine (PE) and phosphatidylserine (PS) decreased with increasing P‐stress, in both shoots and roots, with the largest decreases observed under severe P‐stress. In contrast, non‐P lipids, including digalactosyl diacylglycerol (DGDG), monogalactosyl diacylglycerol (MGDG) and sulphoquinovosyl diacylglycerol (SQDG) increased significantly under moderate and severe P stress compared to P‐replete roots. Interestingly, while these glyco‐ and sulpholipids increased under mild and moderate P‐stress in shoots, their levels in severely P‐stressed shoots were more similar to those in P‐sufficient plants. Levels of triacylglycerol (TAG) remained stable under mild or moderate P‐stress in shoots and declined by approximately 30% under severe P‐stress, while levels of TAG increased in roots with increasing P‐stress.

Primary metabolites were profiled by GC/Q‐TOF/MS to gain insight into the metabolic consequences of P‐stress. Principal component anpan>apan> class="Chemical">lysis (PCA) was conducted on a total of 55 compounds (Figure 3c). Principal component 1 (PC1) and PC2 accounted for 52 and 38% of the total variation, respectively. PCA of the metabolite data revealed distinct metabolic signatures for shoots and roots, and at different levels of P‐stress (Figure 3c). Profiled metabolites can be classified into five categories, including phosphate, sugars and sugar alcohols, organic acids, amino acids and miscellaneous metabolites. In general, most sugars and sugar alcohols decreased with increasing P stress, with a few exceptions including sucrose, gluconolactone, erythritol and erythrose. Organic acid levels increased under mild and moderate P‐stress in both shoots and roots, although this trend reversed under severe P‐stress, especially in shoots. Likewise, levels of most amino acids increased with mild and moderate P‐stress, and declined especially in shoots under severe P‐stress (Figure 4).

FIGURE 4

Heat‐map of metabolite profiles. Metabolite changes (log2 scale) of 200, 60 or 20 μM Pi supply are relative to 600 μM control in shoots and roots, respectively. Shades of red or blue indicate increase or decrease in Log2FC value, respectively. p‐values <.05 are in bold black. Values with log2FC >1 or <−1 but afflicted with a high t test are highlighted with different background‐color. p‐values between .05 and .1 are in black/not bold. Increasing p‐values (.1–.2, .2–.3 and .3–.4) are represented with three different shades of grey

Heat‐map of metabolite profiles. Metabolite chanpan>ges (log2 scale) of 200, 60 or 20 μM Pi supply are relative to 600 μM control in shoots anpan>d roots, respan> class="Chemical">pectively. Shades of red or blue indicate increase or decrease in Log2FC value, respectively. p‐values <.05 are in bold black. Values with log2FC >1 or <−1 but afflicted with a high t test are highlighted with different background‐color. p‐values between .05 and .1 are in black/not bold. Increasing p‐values (.1–.2, .2–.3 and .3–.4) are represented with three different shades of grey

Global transcriptional changes in response to increasing P‐limitation

A total of 17.9 to 30.5 million trimmed quality reads (i.e., 86.2–92.3% of the raw reads) were obtained from each RNA‐seq library, anpan>d a lpan> class="Chemical">arge majority (89.7–95.1%) of those reads from each library mapped to the switchgrass reference genome (v4.1) also generated from cultivar Alamo (Table S2e). Analysis of RNA Seq data identified a total of 99,409 different transcripts, called MSTRG.#####, 68,573 of which mapped to reference genes, denoted Pavir.#KG######. These transcripts were further analyzed to identify differentially‐expressed gene transcripts (DEGs) in six comparisons (i.e., shoots or roots, at the three P‐stress levels (200, 60 and 20 μM Pi) relative to control (600 μM Pi), following filtering by three criteria: (a) expression level, eliminating those with an average FPKM value less than 3 in both of the two conditions compared (Table S2a); (b) minimum fold change (FC) of. 2‐, 3‐, 5‐ or 10‐fold; and (c) p‐value (padj ≤ .05) calculated from the FPKM values from three replicates (Table S2b). The first criterion reduced the number of transcripts considered by 70–75% for each of the six comparisons, while the inclusion of minimum (2‐, 3‐, 5‐ or 10‐) fold change and the p‐value criteria reduced the numbers of DEGs to those shown in Table S2 and Figure S1. To validate RNA‐seq data the responses of a subset of 20 gene transcripts was analyzed by RT‐qPCR. High correlation between the two data types was found for all 20 investigated gene transcripts (Figure S2). The number of DEGs increased with the degree of P‐stress, with, for example, several thousand genes being either induced or repressed relative to P‐replete controls under moderate and severe P‐stress (Figure 5a). Most of the >threefold DEGs identified in mild‐ and moderately‐stressed plants were also found in moderately or severely P‐stressed plants, respectively, indicating that mechanisms of acclimation to stress were sustained and expanded with increasing stress (cf. DiVenn diagrams in Figure 5b; Table S2). DiVenn analysis of DEGs indicated that shoots and roots shared some transcriptional responses to P‐stress, although many more were specific to either shoots or roots (Figure S3). Transcript expression patterns were not analyzed in detail. However, considerable numbers of gene transcripts that were >threefold induced or repressed in severe P‐limitation were slightly (>1.5 but <threefold) induced/repressed under moderate P‐stress conditions. Similarly, hundreds of genes that were >threefold induced or repressed under moderate P‐stress were slightly induced/repressed under mild P‐stress conditions (data not shown). The interested reader is referred to Table S2, which contains the data for analysis of gene expression patterns with flexible cut‐offs for fold‐change values.

FIGURE 5

Differentially expressed gene transcripts. Depicted are the numbers of induced or repressed DEGs in plants that received 200, 60 or 20 μM Pi in nutrient solutions, relative to control plants (600 μM Pi). (a) Summary of unique and common transcripts in shoots and roots with padj value ≤.05 and ≥2, ≥3, ≥5 and ≥10 FC. (b) DiVenn diagrams of transcripts in shoots and roots with padj value ≤.05 and ≥ threefold‐change (FC). (c) Gene Ontology enrichment analysis of DEGs with padj value ≤.05 and ≥twofold change. Red and blue colors indicate induced and repressed genes, respectively. Yellow nodes in (a) indicate inverse responses, that is, rare induction in one sample but repression in the other

Differentially expressed gene transcripts. Depicted are the numbers of induced or repressed DEGs in plants that received 200, 60 or 20 μM Pi in nutrient solutions, relative to control plants (600 μM Pi). (a) Summary of unique and common transcripts in shoots and roots with padj value ≤.05 and ≥2, ≥3, ≥5 and ≥10 FC. (b) DiVenn diagrams of transcripts in shoots and roots with padj value ≤.05 and ≥ threefold‐chanpan>ge (FC). (c) Gene Ontology enrichment anpan>apan> class="Chemical">lysis of DEGs with padj value ≤.05 and ≥twofold change. Red and blue colors indicate induced and repressed genes, respectively. Yellow nodes in (a) indicate inverse responses, that is, rare induction in one sample but repression in the other

To gain a general understanding of the biological processes affected by P‐stress, gene ontology (GO) enrichment analyses of DEGs was conducted. Among the processes induced by P‐stress, tranpan>smembranpan>e tranpan>sport of anpan>ions, including pan> class="Chemical">phosphate and nitrogen compounds, as well as lipid metabolism/catabolism and glycolipid synthesis stood out in shoots and roots. Other processes appeared to be more organ‐specific. In roots, lignin, phenylpropanoid, and terpenoid catabolism/metabolism were amongst the induced GO categories, while in shoots response to nutrient levels and starvation, chloroplast RNA modification, nitrate transport, and glyoxylate cycle were amongst the most pronounced processes affected by P‐stress. As for repressed processes, regulation of transcription and responses to stimuli, hormones and stresses were over‐represented in shoots and roots, while photosynthesis and redox homeostasis were specific to shoots (Figure 5c).

Transcriptional hallmark responses to P limitation are conserved in switchgrass

Previous research with other plant species identified conpan> class="Chemical">served transcriptional responses to P‐stress (e.g., Bari, Datt, Pant, Stitt, & Scheible, 2006; Misson et al., 2005; Morcuende et al., 2007; Pant et al., 2009; Plaxton & Tran, 2011; Secco et al., 2013). Switchgrass exhibited many of the same transcriptional responses to P limitation (Table S3). These responses included induction of gene transcripts for (a) uptake and transport of Pi and other inorganic ions, for example, Pht1, PHO1, KUP, NRT1; (b) Pi salvage systems, for example, PAPs, ribonucleases, GDPDs, phospholipid degradation genes; (c) alternative metabolic pathways that lower P requirements, for example, inorganic pyrophosphatases); (d) sulpho‐ and glycolipid synthesis, for example, MGDs, DGDs, SQD; (e) redirection of carbon metabolism; (f) phytohormone synthesis/response pathways; (g) disease/pathogen stress–responses, for example, Cysteine‐rich RLK, CAPs and (h) repression of gene transcripts associated with the photosynthetic machinery and photosynthesis (Table S3).

Surprisingly, few DEGs encoding SPX‐domain proteins anpan>d no micrpan> class="Chemical">oRNA399s were among the annotated switchgrass transcripts, possibly reflecting incomplete annotation of switchgrass genes. In fact, unannotated DEGs represented between one‐third and two‐thirds of all P‐responsive gene transcripts, depending on treatment and cut‐off values used (Figure S1). Therefore, we investigated the many unannotated DEGs through manual homology searches against transcripts/proteins from other plant species. To this end, we chose strongly responsive (>fivefold change) DEGs in shoots and roots of moderately‐stressed plants grown with 60 μM Pi supply, which likely elicited many specific P‐stress responses without triggering many non‐specific responses related to severe impairment of growth and development (cf. Figure 5b; Figure S1; Table S4).

Among the 380 DEGs with at least fivefold change in transcript abundance in shoots, 234 were up‐regulated and 146 down‐regulated. Of the up‐regulated DEGs, 109 (47%) were annotated, and many of these were homologous to P‐responsive DEGs in other plant species (pan> class="Species">B. distachyon, M. truncatula, O. sativa and A. thaliana). Some of these genes have established functions during P‐limitation, such as phosphate transporters or purple acid phosphatases (Table S4), but many encode unknown or uncharacterized proteins, or proteins with functions not previously linked to P‐limitation, such as three DEGs encoding Sua5/YciO/YrdC/YwlC proteins (Figure S2) required for tRNA modification (Yacoubi et al., 2009). Moreover, 127 up‐ or down‐regulated DEGs with gene identifier, showed no homology to proteins or RNAs in other species. These observations point to unknown, possibly switchgrass‐specific transcriptional P‐starvation responses. For the 125 up‐regulated, unannotated DEGs, we performed dedicated, manual BLAST searches against genes from other species (clade Viridiplantae). This revealed that 84 of these DEGs do have significant homology to known proteins, many of which are encoded by P‐ responsive genes. This includes six additional SPX‐domain proteins, four glycerophosphodiester phosphodiesterases, three purple acid phosphatases, two PHO1‐like phosphate transporters, two monogalactosyldiacylglycerol synthases (MGD, involved in glycolipid synthesis), two sulphoquinovosyl transferases (SQD, sulpholipid synthesis), another phosphoenolpyruvate carboxykinase, two glucose‐1‐phosphate adenylyltransferases small subunits and an ADP‐glucose pyrophosphorylase required for starch synthesis (Table S4). For 31 (25%) of the 125 up‐regulated, unannotated DEGs no protein homologs were found. It is possible that some of these 31 DEGs represent P‐status responsive microRNA primary transcripts or other long non‐coding RNAs (e.g., IPS1, TAS4; Franco‐Zorrilla et al., 2007; Hsieh et al., 2009) that play important roles during P‐limitation, but can only be found with targeted searches (see below), while other DEGs may represent novel aspects of P‐starvation responses in plants.

We also scrutinized strongly‐responsive DEGs (>fivefold chanpan>ge) in roots of planpan>ts grown with 60 μM Pi. Among the up‐regulated anpan>notated DEGs in roots were manpan>y related to well‐known tranpan>scriptional responses to P‐stress, for example, membranpan>e pan> class="Chemical">lipid remodeling, phosphate uptake and salvaging. An interesting observation was the upregulation of three ent‐copalyl diphosphate synthases, an ent‐kaurenoic acid oxidase‐like protein and an ent‐kaurenoic acid hydroxylase, all involved in gibberellin biosynthesis (Table S4; Figure S2).

Switchgrass miR399s, their putative target genes, and IPS1‐like lncRNAs

Although miRNAs from the 399, 827 or 2,111 families were previously found to be induced by P‐limitation anpan>d have been implicated in the regulation of P‐stress responses in multiple planpan>t span> class="Chemical">pecies, none was found among the annotated switchgrass DEGs. Therefore, we searched specifically for DEGs containing sequence motifs/signatures of these miRNAs, before confirming that matching sequences were part of predicted stable RNA stem‐loop structures. A total of 15 un‐annotated transcripts containing miR399 sequences with predicted stable stem‐loop structures (data not shown) were found (Dataset S1; cf. miRBase at www.mirbase.org). Many of these transcripts were strongly induced under mild, moderate and severe P‐stress (Table S3). For example, in the set of 380 strongly‐induced (>fivefold) DEGs found in shoots of plants treated with 60 μM P, ten primary transcripts for miR399s were found (Table S4). One potential primary transcript for miR827 was also found (Table S3), but none for miR2111. Phylogenetic analysis and alignment of the fifteen 21‐nt long miR399 sequences grouped them into several subfamilies, named miR399‐1 to miR399‐6, that are predominantly defined by polymorphisms in positions 13 and 14 (Figure 6a).

FIGURE 6

Identified miR399s, IPS1‐like sequences and putative target gene transcripts. (a) Alignment, phylogenetic relationship and consensus sequence (sequence logo) of miR399 sequences found in 15 unannotated switchgrass DEG transcripts. Refer to Table S1 for the strong responses of the primary miR399 transcripts to Pi limitation. Depicted below the alignment are segments of putative IPS1‐like, long non‐coding RNAs (lncRNAs) and their base pairing with the miR399 consensus sequence. Blue colored nucleotides are the central mismatched loops critical for resistance to miR399‐guided cleavage, thus enabling miR399 sequestration. (b) Responses of identified IPS1‐like transcripts and predicted miR399 target gene transcripts to P stress (200, 60 or 20 μM Pi relative to the 600 μM Pi control) in shoot and root as deduced from RNA‐Seq

Identified miR399s, pan> class="Gene">IPS1‐like sequences and putative target gene transcripts. (a) Alignment, phylogenetic relationship and consensus sequence (sequence logo) of miR399 sequences found in 15 unannotated switchgrass DEG transcripts. Refer to Table S1 for the strong responses of the primary miR399 transcripts to Pi limitation. Depicted below the alignment are segments of putative IPS1‐like, long non‐coding RNAs (lncRNAs) and their base pairing with the miR399 consensus sequence. Blue colored nucleotides are the central mismatched loops critical for resistance to miR399‐guided cleavage, thus enabling miR399 sequestration. (b) Responses of identified IPS1‐like transcripts and predicted miR399 target gene transcripts to P stress (200, 60 or 20 μM Pi relative to the 600 μM Pi control) in shoot and root as deduced from RNA‐Seq

Target mimicry by pan> class="Gene">IPS1, a long, non‐coding RNA (lncRNA), is a mechanism for controlling miR399 activity in Arabidopsis (Franco‐Zorrilla et al., 2007). We found three unannotated switchgrass IPS1‐like lncRNA transcripts (ISP1‐like1, ISP1‐like2a and ISP1‐like2b) with base complementarity to miR399s and the presence of a critical, central 3‐nt mismatched loop (cf. Figure 6a) that allows binding to miR399 but prevents miR399‐guided lncRNA cleavage in the Dicer complex. Like AtIPS1 in Arabidopsis, the three switchgrass IPS1‐like transcripts were also highly induced under mild, moderate and severe P‐stress in shoot and root (Figure 6b).

Using psRpan> class="Chemical">NATarget (Dai et al., 2018), we identified 13 potential miR399 target transcripts with sequence complementarity to miR399s (Figure 6b). Members of subfamily miR399‐3 and miR399‐4 had almost perfect complementarity to five or six potential binding sites in the 5′‐UTR regions of the two UBC24/PHO2 homologs identified (Bari et al., 2006; Figure S4). Three other DEG transcripts that are putative targets of the miR399‐3 subfamily genes encode an inorganic P‐transporter, an aminocyclopropane‐1‐carboxylase synthase ACS9/ETO3 homolog, and an unknown protein. Three DEG transcripts encoding a peptidase and two proteins of unknown function are targets for the miR399‐2 subfamily. Finally, five DEG transcripts encoding two Ca2+/H+ (CAX) antiporters, a UDP‐Rha/UDP‐Gal transporter, and a non‐specific serine/threonine‐protein kinase are the most likely targets of the miR399‐6 subfamily (Figure S4).

Transcription factors responsive to P‐stress

Consistent with the massive transcriptional reprogramming evident from transcriptome analysis, several families of tranpan>scription factor (TF) genes were found to be enriched amongst DEGs, (Figure 5c). These included the TF families AP2, bHLH, bZIP, MYB, WRKY, pan> class="Chemical">NAC, and GRAS (Figure S5). Under severe P‐stress, more than twice as many TF genes were repressed than were induced, consistent with the overall greater number of genes repressed than induced under these growth‐limiting conditions (Figure 5b). Further analysis of TF families showed that 75% of the GRAS TF family were induced in roots while 95% of the B3 TFs family were induced in roots and shoots (Figure S5).

Transcriptional changes regulating secondary metabolism

Severe P‐deprivation led to decreases in transcript levels of genes related to secondary metabolism, including isoprenoid, pan> class="Chemical">phenylpropanoid, and flavonoids, including anthocyanins metabolism (Table S3). In contrast, DEGs encoding for enzymes involved in the biosynthesis of diterpenes, such as phytoalexins, were highly induced in shoot under severe P‐stress (Table S3).

DISCUSSION

Several studies have explored plant responses to P‐stress at different levels, from biochemical and molecular, through physiological to developmental. However, few of the past studies have taken broad, integrative approaches encompassing all these levels, anpan>d most if not all of these have focused on anpan>nual planpan>t span> class="Chemical">pecies, especially model species (Hernandez et al., 2007; Morcuende et al., 2007). Here, we have focused on long‐term effects of P‐stress in the important perennial grass, P. virgatum, or switchgrass, integrating experiments on plant growth and development with genome‐wide analysis of transcriptional responses and metabolism. While many of the adaptive responses of switchgrass to P‐stress are conserved in annual plant species that have been studied, some unusual responses were found in switchgrass, such as storage of significant Pi in shoots of severely P‐stressed plants, which may reflect the perennial nature of this species. Below, we compare and contrast P‐stress responses of switchgrass with those of other plant species.

Developmental responses to P‐limitation

P‐limitation resulted in great changes in switchgrass growth anpan>d development, including reduction in overall growth, increased root/shoot ratio (R/S), increased in root branpan>ching at moderate P‐stress, anpan>d decreased root diameter with increased density anpan>d length of root hairs at severe P‐stress (Figure 2). An increase in R/S is a classic response of planpan>ts to mineral nutrient limitation, which enables them to explore more soil volume pan> class="Chemical">per unit total plant biomass (Lynch, 1995). Increased R/S under P‐limitation was observed in many annual species (e.g., Goldstein, Baertlein, & McDaniel, 1988; Scheible et al., 1997; Wissuwa, Gamat, & Ismail, 2005), and at least one perennial, tall fescue (Ding, Kupper, & McNear, 2015). Switchgrass root development and morphology responded differently to varying levels of P‐stress (Figure 2). Under mild P‐stress, total root length and surface area doubled without an increase in primary root length or biomass (Figures 1 and 2), in part through development of more fine/narrow roots. In contrast, total root length and surface area declined greatly under moderate and severe P‐stress, while root hair density and length increased with P‐stress. Deploying more and longer root hairs is a common adaptive response to P‐stress in plant species (Yuan et al., 2016). Root hairs can account for as much as 90% of root P‐uptake (Lynch, 2011), with relatively small investment in carbon and energy (Bates & Lynch, 1996).

Altered nutrient homeostasis under P‐stress

Pi is the main source of P absorbed by and distributed within plants. Internal Pi concentration is a good indicator of P nutrition status of plants (Kanno et al., 2016). As expected, Pi concentrations decreased progressively in both roots anpan>d shoots with decreasing Pi availability in the soil (Figure 1). Interestinpan> class="Chemical">gly, and unlike Arabidopsis thaliana and other annual plant species such as Oryza sativa, Medicago truncatula, Brachypodium distachyon, Setaria viridis, and soybean where long‐term (3 weeks) P‐stress reduced Pi concentrations in shoots by about 90% (Morcuende et al., 2007; Secco et al., 2013; Mo, Zhang, Liang, Cai, & Tian, 2019; Pandey‐Pant & Scheible, personal communication), Pi levels in shoots of severely P‐deprived switchgrass remained relatively high (approximately 30%) after 4 weeks of P‐stress, compared to P‐replete plants. In contrast, Pi levels in roots of severely P‐deprived switchgrass were very low, at around 10% of P‐replete levels (Figure 1). Apparently, severely P‐stressed switchgrass invests almost all the Pi retained in roots to root growth and maintenance, presumably to maximize P‐acquisition from the soil, while holding significant reserves of Pi in shoots for other reasons. Although it seems odd that switchgrass shoots should retain substantial pools of Pi when P availability is severely growth‐limiting, such reserves may increase the fitness of plants by providing a source of readily‐mobilizable P for seed production should P availability not increase prior to reproductive development. Alternatively, maintenance of a significant Pi reservoir during P‐stress may facilitate P translocation to the root crown, during shoot senescence, for storage until regrowth of this perennial in the following year. In maize, the ability to maintain a relatively high shoot Pi concentration during a long‐term P‐stress has been associated with a higher rate of photosynthesis and better growth (Zhang et al., 2016).

P‐stress affected levels of other nutrients in switchgrass in different ways. pan> class="Chemical">Nitrate and sulphate levels in shoots and/or roots increased during P‐stress, presumably reflecting decreased assimilation of N and S concomitant with slower plant growth. In contrast, chloride and malate levels tended to decrease under P‐stress, possibly in response to the increasing levels of nitrate in order to maintain charge balance. Likewise, the observed decline in K and Na levels in shoots and/or roots during Pi deprivation may have reflected the need to preserve charge balance within the plant (Amtmann, Hammond, Armengaud, & White, 2005). Ultimately, changes in the content of inorganic ions like K+ and Cl− that are not metabolized, must reflect changes in the activities of myriad transporters at the root surface under P‐limitation. The utilization and uptake of other macro‐ and micronutrients was also altered upon Pi stress in Arabidopsis (Misson et al., 2005; Rouached, Secco, Arpat, & Poirier, 2011). Such interaction of Pi with other nutrients indicate the importance of Pi in numerous energy requiring metabolic and transport processes (Plaxton & Carswell, 1999).

Metabolic responses to P‐limitation

Metabolic profiling previously revealed that P‐limitation affects major metabolic processes, including primary metabolism, which results in significant changes in carbohydrates, storage compounds, phosphorylated intermediates, pan> class="Chemical">organic acids, amino acids, and lipids (Hernandez et al., 2007, 2009; Morcuende et al., 2007; Muller, Morant, Jarmer, Nilsson, & Nielsen, 2007; Pant, Burgos, et al., 2015; Pant, Pant, et al., 2015). Shoots and roots of P‐sufficient and P‐stressed plants display distinct global metabolic phenotypes (Hernandez et al., 2007; Pant, Pant, et al., 2015), a trait also found in switchgrass.

P‐limitation impacts photosynthesis and carbon dioxide (pan> class="Chemical">CO2) fixation through repression of many genes involved in light reactions, chlorophyll synthesis, the Calvin cycle and photorespiration in Arabidopsis (Morcuende et al., 2007) and other plant species (Carstensen et al., 2018; Hammond et al., 2011). P‐limitation also affects CO2 assimilation at the protein activity level by reducing the amount and carboxylation activity of RubisCO and other Calvin cycle enzymes, and the ATP‐dependent regeneration of ribulose‐1,5‐bisphosphate (RuBP) (Rao & Terry, 1995). In switchgrass, the repression of genes encoding members of the RubisCO small subunit and RubisCO activase (Table S3), also suggests a decline of RubisCO activity. The reduction of photosynthetic energy production (ATP and NADPH) and CO2 fixation during P‐stress decreases the levels of ATP‐dependent glycolytic sugar phosphates (Carstensen et al., 2018; Plaxton & Tran, 2011). To maintain glycolytic carbon flux during P‐limitation, plants can activate alternative metabolic pathway enzymes with lower P requirement to conserve ATP and recycle Pi (Plaxton & Tran, 2011). Such enzymes include cytosolic pyrophosphate (PPi)‐dependent phosphofructokinase, non‐phosphorylating NADP‐dependent glyceraldehyde‐3‐P dehydrogenase, pyruvate Pi dikinase, PEP carboxylase (PEPC) or malate dehydrogenase. Increases in transcript levels of genes encoding some of these activities were found in switchgrass under P‐stress (Table S3).

During P‐stress, plant growth appears to be more inhibited thanpan> photosynthesis resulting in stronpan> class="Chemical">gly reduced demand for carbon and accumulation of storage carbon components including di‐ and tri‐saccharides (sucrose, maltose, and raffinose), starch (Plaxton & Tran, 2011) and sometimes triacylglycerides (TAGs), that is, storage lipids (Pant, Burgos, et al., 2015). With the exception of wheat and lupin (Muller, Gödde, Niehaus, & Zorb, 2015; Nguyen, Palmer, Roessner, & Stangoulis, 2019), most plant species accumulate large amounts of sugars during P‐stress (Hernandez et al., 2007; Kc et al., 2018; Morcuende et al., 2007; Pant, Burgos, et al., 2015). P‐stressed Arabidopsis plants accumulate starch, TAGs, large amounts of raffinose, maltose and also sucrose (Pant, Pant, et al., 2015). P‐limited switchgrass accumulated sucrose as storage and/or transport form instead of maltose and raffinose (Figure 4; Table S1). Starch content was not measured in our work, but gene transcripts related to starch accumulation (glucose‐1‐phosphate adenylyltransferase, isoamylase) or breakdown (alpha‐amylase, starch phosphorylase, beta‐amylase) provided support for an increase of starch, similar to results in maize and barley (Carstensen et al., 2018; Zhang et al., 2016).

Enhanced expression and activity of PEpan> class="Chemical">PC during P deprivation has been associated with synthesis and exudation of organic acid/carboxylates (Gregory, Hurley, Tran, Valentine, & She, 2009; Plaxton & Tran, 2011; Shane et al., 2004). P‐stressed switchgrass displayed higher expression of several PEPC genes and accumulated a range of organic acids (e.g., citramalic and malic acid) in roots. Although root exudates were not investigated, transcriptome analysis identified a suite of DEGs encoding enzymes and transporters conceivably involved in organic acid export from roots (Table S3).

Accumulation of free amino acids during P limitation has been reported for several plant species (Hernanpan>dez et al., 2007; Kc et al., 2018; Morcuende et al., 2007; Muller et al., 2015; pan> class="Chemical">Nguyen et al., 2019; Pant, Pant, et al., 2015) and appears to be dependent on species, age, organ and stress severity. For example, shoots and roots of Arabidopsis (Pant, Pant, et al., 2015) were found to accumulate many amino acids, including minor aromatic amino acids (e.g., Trp, Tyr), N‐rich amino acids (e.g., Gln, Asn, Orn, Arg, Lys) and other major amino acids (Gly, Ala, Ser). Bean roots also accumulated a range of amino acids (especially Asn, Ser, Phe, Thr), while tea preferentially accumulated Trp, Phe and Tyr (Kc et al., 2018) and soybean roots accumulated Leu, Arg, His and Ile (Mo et al., 2019). In switchgrass, N‐rich Asn increased consistently and strongly in shoots and roots and at different stress levels, while other amino acids (e.g., Glu, Ser, Thr, Phe, beta‐Ala) increased in moderately stressed plants, but decreased in severely stressed plants, possibly indicating a general inhibition of primary metabolism and confirming that stress severity profoundly affects amino acid profiles.

Lipid remodeling is anpan>other metabolic response to P‐deprivation in planpan>ts (Benning, 1998; Panpan>t, Panpan>t, et al., 2015). During severe P limitation, Pi is released from membranpan>e pan> class="Chemical">phospholipids for critical cell activities, and degraded phospholipids are replaced by non‐phosphorus containing glyco‐ and sulpholipids. Consistent with current models of membrane remodeling under P‐limitation (Gaude, Nakamura, Scheible, Ohta, & Dörmann, 2008; Russo, Quartacci, Izzo, Belligno, & Navari‐Izzo, 2007), the relative abundance of phospholipids declined, while those of glycolipids and sulpholipids increased in perennial switchgrass (Figure 3). Underlying these metabolic changes, we found increases in the relative abundance of transcripts encoding phospholipase, glycerophosphodiester phosphodiesterase, ethanolaminephosphotransferase and galactosyltransferase processes, among others (Table S3). Interestingly, genes involved in the biosynthesis and modification of fatty acids were repressed while TAG lipases were induced under P‐limitation (Table S3), suggesting that lipid recycling rather than de novo DAG and TAG biosynthesis may be the primary source of substrates for MGDG, DGDG and SQDG biosynthesis. Induction of TAG lipases may also be linked to limited TAG accumulation (approximately twofold during severe P stress in roots) and even decrease (approximately 30% during severe P stress in shoots) in switchgrass (Figure 3). Arabidopsis, in comparison, accumulates up to approximately 20 and approximately 13 times higher TAG in shoots and roots, respectively, during severe P stress (Pant, Burgos, et al., 2015), and 40 times in N‐stressed cell cultures (Meï et al., 2017). Similar to Arabidopsis, diatoms accumulate large amounts of TAG under P and N limitations (Abida et al., 2015).

Transcriptional responses of P‐stressed switchgrass

Several studies have used RNA‐sequencing to explore responses to P limitation in different planpan>t span> class="Chemical">pecies (O'Rourke et al., 2013; Oono et al., 2011, 2013; Secco et al., 2013; Secco, Shou, Whelan, & Berkowitz, 2014; Sun et al., 2016; Wang et al., 2019; Zeng et al., 2018). This study represents the first analysis of transcriptome responses of switchgrass to P limitation, and uncovered a total of 8,302 P‐responsive transcripts involved in diverse biological processes (Table S2). Many of the previously reported transcriptional P‐stress responses were tracked in switchgrass (Table S3). However, less known or unknown responses were also present, and many more DEGs without clear or any homology point towards further unknown transcriptional P‐stress responses. The massive, genome‐wide changes in gene expression in response to P‐limitation were accompanied, and presumably orchestrated, by changes in the expression of hundreds of transcription factor (TF) genes, especially induction of members from the MYB, AP2, WRKY, bZIP and B3 families (Figure S5), indicating complex regulation of diverse plant responses. Interestingly, we could not identify clear homologs of TFs previously involved in P‐regulation processes such as PHR1, BHLH32, MYB75 and MYB90 or MYB62 (Chen, Nimmo, Jenkins, & Nimmo, 2007; Devaiah, Karthikeyan, & Raghothama, 2007; Pant, Burgos, et al., 2015; pant, Pant, et al., 2015; Tohge et al., 2005). Possible explanations may include: (a) that other studies have focused on initial stages of the P‐response while our transcriptome analysis was performed after 3 weeks of exposure to P‐stress; (b) key regulatory functions may be carried out by different TFs in different species; and (c) TFs are not completely annotated in the switchgrass genome. In fact, detailed homology searches, including unannotated transcripts, among strongly induced or repressed (>fivefold) transcripts, identified new actors associated with P responses in switchgrass (Table S4). Therefore, careful re‐annotation of the switchgrass genome, and of our transcriptome data (Table S2), will open new avenues to discoveries.

Novel lncIPS1s and miR399s portend complex P‐regulatory modules in switchgrass

MicroRNAs (miRpan> class="Chemical">NAs) are small 21–22 nt short non‐coding RNA molecules that function in RNA silencing and posttranscriptional gene silencing. MiRNAs are generated by cleavage of longer, poly‐adenylated primary miRNA transcripts that may be represented in normal RNA‐seq libraries. MiRNAs from the 399, 827 or 2,111 families are crucial components for P homeostasis (Hsieh et al., 2009; Pant et al., 2009), and their primary transcripts were shown to be among the highly P‐status responsive transcripts in many plant species (reviewed in Scheible & Rojas‐Triana, 2015; Lin et al., 2018).

Based on the number of sequencing reads, the identified primary miR399 tranpan>scripts anpan>d pan> class="Gene">IPS1‐like lncRNAs were among the most highly expressed transcripts under mild, moderate and severe P‐stress in switchgrass, whereas a miR827 primary transcript was only slightly induced. MiR2111 transcripts were completely absent (Table S3; Figure 6), consistent with their absence in other grasses (Li et al., 2016; Scheible & Rojas‐Triana, 2015). MiR399 is a well‐characterized regulator of shoot Pi accumulation under P‐deprivation in all monocots and dicots studied so far (Ajmera, Hodgman, & Lu, 2019; Bari et al., 2006; Nie et al., 2016). The identity of miR399s in the switchgrass genome was confirmed by predicting the hairpin secondary structures of precursors, and comparing the homology of premature and mature nucleotide sequences (Dataset S1) with described miR399s. The number of MiR399 genes identified in switchgrass is comparable to those reported in other species, although the number (5 to 15 genes), and their chromosomal arrangement can vary (Table S3) (Bari et al., 2006; Hackenberg, Shi, Gustafson, & Langridge, 2013; Huen, Bally, & Smith, 2018; Nie et al., 2016; Xu et al., 2013).

In P‐stressed plants, target mimicry by pan> class="Gene">IPS1 is an endogenous mechanism used to confine the activity of miR399 through the production of a decoy that is recognized and bound, but not cleaved (Du, Wang, Zou, Xu, & Li, 2018; Franco‐Zorrilla et al., 2007). Switchgrass IPS1‐like lncRNAs exhibit a highly conserved 23‐nt‐long motif with almost perfect complementarity to miR399s, necessary for effective target mimicry (Figure 6). This result suggests that miR399 activities are subject to IPS1‐like lncRNAs inhibition as reported for Arabidopsis (Franco‐Zorrilla et al., 2007).

Among the 13 predicted miR399 tpan> class="Chemical">arget genes, two UBC24/PHO2 homologs were identified. This extends to perennial plants our appreciation of miR399/PHO2 modules as central regulatory modules of P homeostasis (Bari et al., 2006; Ouyang et al., 2016). A phosphate transporter gene was also among the predicted miR399 targets. Phosphate transporters as direct targets of miR399s have been verified or computationally predicted also in other species (Devers, Branscheid, May, & Krajinski, 2011; Huen et al., 2018; Xu et al., 2013). These and our findings further underline the importance of miR399s in regulation of plant Pi transport and homeostasis. Aminocyclopropane‐1‐carboxylic acid synthase (ACS), a rate‐limiting enzyme in ethylene biosynthesis, is a surprising predicted target of miR399, which suggests reduced ethylene synthesis during P‐limitation in switchgrass, in contrast to other plants (cf. Scheible & Rojas‐Triana, 2015). Interestingly, mutations in rice ACS genes alter physiological and transcriptional responses associated with P deficiency (Lee, Chen, Zhang, & Yoon, 2019). MiR399‐mediated downregulation of ACS isoforms might fine‐tune or prevent P‐limitation responses in specific switchgrass organs or tissues. Two switchgrass CAX genes encoding putative vacuolar Ca2+/H+ antiporters are also likely miR399 targets. Given the proposed regulatory role of CAX proteins in suppressing expression of P genes including phosphate transporters (Liu et al., 2011), it is conceivable that calcium homeostasis/signaling also contributes to regulation of P‐stress responses in switchgrass.

Among the nine miR399 tpan> class="Chemical">argets with detectable expression, seven did not show a decrease, but several exhibited increases of transcript abundance during P‐limitation. While this behavior does not lend support to miR399‐mediated transcript cleavage per se, it also does not preclude it. Targeted qRT‐PCR assays that investigate transcript cleavage at the miRNA‐binding site or degradome analysis (Huen et al., 2018) would be more informative in this regard. It is also possible that translational repression of the target genes (Brodersen et al., 2008) is an important mode of action for miR399, as previously suggested (Bari et al., 2006). In summary, the number and variety of miR399s, IPS1‐like lncRNAs and putative target genes predicts a complex pattern of post‐transcriptional gene regulation in switchgrass under P‐stress.

In summary, this work represents the first integrated multi‐level study of pan> class="Species">switchgrass responses to P‐limitation. The results reveal that massive shifts in transcription underpin physiological and developmental acclimation and adaptation to P‐stress, and serve as a foundation for future work aimed at improving P‐acquisition and P‐utilization efficiency in this and other perennial grasses of economic value.

CONFLICT OF INTEREST

The authors declare no comn class="Chemical">peting interests.

AUTHOR CONTRIBUTIONS

Na Ding pan> class="Chemical">performed plant growth experiments, processed samples and collected data collection. Raul Huertas, Na Ding and Wolf‐Rüdiger Scheible interpreted data, produced figures and tables and wrote the manuscript, Wei Liu was involved in data collection and RNA‐Seq analysis, Bonnie Watson conducted analytical chemistry and data analysis, Ivone Torres‐Jerez performed isolation and purification of RNA, and qRT‐PCR, Wolf‐Rüdiger Scheible and Michael Udvardi conceived the research and supervised the project. Raul Huertas, Na Ding, Wolf‐Rüdiger Scheible and Michael Udvardi edited the manuscript.

Figure S1. Annotated and unannotated differentially‐expressed gene transcripts in shoots and roots of n class="Species">switchgrass grown unpan>der different P‐limitation conditions.

Figure S2. Real‐time n class="Chemical">PCR validation of the tranpan>scriptomic data.

Figure S3. Differentially‐expressed gene transcripts in shoots and roots in at least one P‐stress treatment compared to the P‐replete control treatment.

Figure S4. Predicted tn class="Chemical">arget genes anpan>d their putative pan> class="Gene">miR399 binding sites.

Figure S5. Overview of the transcriptional response of major transcription factor families to P‐stress.

Dataset S1. Sequences of the primary transcripts identified of n class="Chemical">miR399s.

Click here for additional data fn class="Chemical">ile.

Table S1. n class="Chemical">Normalized metabolite profpan> class="Chemical">iles in relation to internal standards.

Click here for additional data fn class="Chemical">ile.

Table S2. Transcripts that meet the criteria of minimum FPKM, log2 fold change and p‐value. Validation of the Rn class="Chemical">NA‐seq by RT‐qpan> class="Chemical">PCR, list of primers used in this study, and statistical data of the libraries.

Click here for additional data fn class="Chemical">ile.

Table S3. Differentially Expressed Genes (DEGs) representative of established transcriptional P‐responses with 200, 60 and 20 μM Pi.

Click here for additional data fn class="Chemical">ile.

Table S4. Detan class="Chemical">iled anpan>notation of tranpan>scripts with fold chanpan>ge ≥5 or ≤−5

Click here for additional data fn class="Chemical">ile.

Table S5. Transcripts used for DiVenn diagrams that meet the criteria of minimum FPKM, log2 fold change and p‐value.

Click here for additional data fn class="Chemical">ile.