The pho1;2a'-m1.1 allele of Phosphate1 conditions misregulation of the phosphorus starvation response in maize (Zea mays ssp. mays L.).

Plant PHO1 proteins play a central role in the translocation and sensing of inorganic phosphate. The maize (Zea mays ssp. mays) genome encodes two co-orthologs of the Arabidopsis PHO1 gene, designated ZmPho1;2a and ZmPho1;2b. Here, we report the characterization of the transposon footprint allele Zmpho1;2a'-m1.1, which we refer to hereafter as pho1;2a. The pho1;2a allele is a stable derivative formed by excision of an Activator transposable element from the ZmPho1;2a gene. The pho1;2a allele contains an 8-bp insertion at the point of transposon excision that disrupts the reading frame and is predicted to generate a premature translational stop. We show that the pho1;2a allele is linked to a dosage-dependent reduction in Pho1;2a transcript accumulation and a mild reduction in seedling growth. Characterization of shoot and root transcriptomes under full nutrient, low nitrogen, low phosphorus, and combined low nitrogen and low phosphorus conditions identified 1100 differentially expressed genes between wild-type plants and plants carrying the pho1;2a mutation. Of these 1100 genes, 966 were upregulated in plants carrying pho1;2a, indicating the wild-type PHO1;2a to predominantly impact negative gene regulation. Gene set enrichment analysis of the pho1;2a-misregulated genes revealed associations with phytohormone signaling and the phosphate starvation response. In roots, differential expression was broadly consistent across all nutrient conditions. In leaves, differential expression was largely specific to low phosphorus and combined low nitrogen and low phosphorus conditions. Of 276 genes upregulated in the leaves of pho1;2a mutants in the low phosphorus condition, 153 were themselves induced in wild-type plants with respect to the full nutrient condition. Our observations suggest that Pho1;2a functions in the fine-tuning of the transcriptional response to phosphate starvation through maintenance and/or sensing of plant phosphate status.

INTRODUCTION

Phosphorus (P) is essential for crop growth and productivity. P is required for photosynthesis and respiration, and is a structural component of both nucleic acids and phospholipids (Chien et al., 2018; Péret et al., 2011; Plaxton & Tran, 2011; Veneklaas et al., 2012). Elemental P is highly reactive and readily oxidized to orthophosphate (PO4 3−; hereafter, inorganic phosphate, Pi) which occurs in neutral soil solution as a mixture of hydrogen phosphate (HPO4 2−) and dihydrogen phosphate (H2PO4 −). To access soil P, plants must acquire it in the form of free Pi (Schachtman et al., 1998). However, Pi has a strong tendency to react with other soil components (predominantly with iron and aluminum at acidic pH and with calcium and magnesium at alkaline pH), reducing both its mobility and its availability to plants (Hinsinger, 2001; Vance et al., 2003). The concentration of plant available Pi in the soil ranges between 2 and 10 μM (Raghothama, 1999).

It is estimated that Pi availability in ~70% of agricultural soils is below the level required to realize yield potential (López‐Arredondo et al., 2014). For optimal growth, plants must maintain an intracellular Pi concentration of 5–20 mM, far higher than the concentration available in the soil (Raghothama, 1999; Vance et al., 2003). To acquire Pi against this concentration gradient, plants actively transport Pi through high‐affinity transporters of the PHOSPHATE TRANSPORTER1 (PHT1) family (Misson et al., 2004; Rausch & Bucher, 2002). Once acquired, Pi is distributed around the plant through further action of the PHT1 proteins (e.g., Ai et al., 2009; Chang et al., 2019). Pi translocation within the plant is additionally supported by the action of PHOSPHATE1 (PHO1) proteins. In Arabidopsis, PHO1 is indispensable for the translocation of Pi from roots to shoots, acting as a Pi efflux transporter (Bulak Arpat et al., 2012; Ma et al., 2021; Poirier et al., 1991; Secco et al., 2010; Zhao et al., 2019). PHO1 is expressed in the pericycle and root xylem parenchyma cells where it loads Pi to the xylem for translocation (Bulak Arpat et al., 2012). PHO1 proteins contain a tripartite N‐terminal SPX domain and a hydrophobic C‐terminal EXS domain (Ried et al., 2021; Secco et al., 2012; Stefanovic et al., 2011). The EXS domain is required for the subcellular localization of PHO1 to the Golgi and trans‐Golgi networks and is necessary for Pi export (Wang et al., 2004; Stefanovic et al., 2007: Secco et al., 2012; Wege et al., 2016). The SPX domain is characteristic of a broader family of SPX domain proteins that act in Pi sensing and signaling through physical interaction with other proteins and inositol pyrophosphate molecules (Ried et al., 2021).

Plants have evolved developmental, molecular, physiological, and metabolic strategies to optimize P acquisition and internal P use efficiency. These strategies are modulated by plant P status and collectively referred to as the P starvation response (PSR). At the molecular level, P starvation promotes widespread changes in both transcriptional and post‐transcriptional gene regulation (Chang et al., 2019; Rausch & Bucher, 2002; Shin et al., 2004; Torres‐Rodríguez et al., 2021). Physiological experiments have demonstrated that plants respond both locally to environmental Pi concentration and systemically to internal Pi status (Lin, Huang, et al., 2014; Thibaud et al., 2010). The PSR is tightly regulated by a transcriptional network that converges on the MYB transcription factor PHOSPHATE STARVATION RESPONSE 1 (PHR1), a master regulator of P sensing and signaling (Barragán‐Rosillo et al., 2021; Bustos et al., 2010; Paz‐Ares et al., 2022; Puga et al., 2014; Ried et al., 2021; Wang et al., 2014). In Arabidopsis, the absence of full PHO1 function in under‐expressor lines, or pho1 mutants partially complemented by heterologous expression of rice PHO1, uncouples low internal Pi from the transcriptional and physiological symptoms of P starvation (Rouached et al., 2011). In wild type Arabidopsis plants a mild reduction in external Pi availability and shoot Pi acquisition results in reduced shoot growth, even though unused Pi remains in reserve in the vacuole. Although the shoots of plants with partial PHO1 function are Pi deficient at high external Pi availability, the plants grow normally by using the majority of Pi available and significantly reducing vacuolar reserves (Rouached et al., 2011). Similarly, PHO1 under‐expressors do not exhibit the lipid remodeling or complete transcriptional PSR seen in wild‐type plants grown under low external Pi to establish an equivalent rate of shoot Pi acquisition (Rouached et al., 2011). As such, PHO1 links plant Pi status, the PSR, and the development of P deficiency symptoms.

OsPHO1;2, the functional rice ortholog of AtPHO1, is abundantly expressed in roots and is required for Pi translocation (Secco et al., 2010). Interestingly, OsPHO1;2 produces both sense transcripts and a cis‐acting natural antisense transcripts (cis‐NAT ). The cis‐NAT is upregulated under P starvation to promote translation of sense OsPHO1;2 transcripts through the formation of double‐stranded RNA and ribosome association (Jabnoune et al., 2013; Reis et al., 2021). Maize has two co‐orthologs of OsPHO1;2, designated Pho1;2a and Pho1;2b (Salazar‐Vidal et al., 2016). The coding sequences of Pho1;2a and Pho1;2b are highly similar, reflecting their presumed origin during an ancient whole‐genome duplication (Salazar‐Vidal et al., 2016). Pho1;2a and Pho1;2b transcripts accumulate preferentially in the roots (Salazar‐Vidal et al., 2016), although their absolute expression level differs, with Pho1;2b transcript accumulation ~30 fold greater than that of Pho1;2a (Salazar‐Vidal et al., 2016; Woodhouse et al., 2021). Reverse‐transcription PCR has indicated the production of a cis‐NAT associated with Pho1;2a homologus to the rice product, although this remains to be definitively demonstrated (Salazar‐Vidal et al., 2016). To date, there is no evidence of cis‐NAT transcripts associated with the maize Pho1;2b gene. CRISPR‐Cas9 editing has suggested maize Pho1;2a and Pho1;2b to play non‐redundant roles in grain filling, with edited lines showing reduced starch synthesis and small, shrunken kernels (Ma et al., 2021). However, there are no reports concerning the role of maize PHO1 genes in root‐to‐shoot Pi translocation or the PSR.

To study PHO1 function in maize, we previously identified insertion of an Activator (Ac) transposon into the sixth exon of the maize Pho1;2a gene (pho1;2a‐m1::Ac; Salazar‐Vidal et al., 2016). The original pho1;2a‐m1::Ac allele showed a high level of somatic instability, prompting us to generate the stable footprint derivative pho1;2a′‐m1.1 (Bai et al., 2007; Salazar‐Vidal et al., 2016). Here, we report the characterization of pho1;2a′‐m1.1 mutants. We found the pho1;2a′‐m1.1 mutation to have only subtle effects on plant growth and no statistically significant impact on total P concentration in the leaf. We characterized root and leaf transcriptomes in a stock segregating pho1;2a′‐m1.1 under different nutrient regimes, observing misregulation of several hundred transcripts in both heterozygous and homozygous mutant plants. Many of the transcripts upregulated in the leaves of plants carrying pho1;2a′‐m1.1 were also induced by low P availability in wild‐type plants, suggesting a role for Pho1;2a in feedback regulation of the transcriptional PSR, potentially through sensing of Pi status.

MATERIALS AND METHODS

Plant material and genotyping

Analysis of the wild‐type siblings of the mutant segregants described here was reported previously in Torres‐Rodríguez et al. (2021). Three self‐pollinated sibling families segregating for the pho1;2a′‐m1.1 mutation in the Zea mays ssp. mays var. W22 background were used for the experiments (Salazar‐Vidal et al., 2016). Plants were genotyped as previously reported (Salazar‐Vidal et al., 2016) with small modifications. Briefly, DNA was extracted from a punch of cotyledon tissue and a fragment was amplified spanning the position of the pho1;2a′‐m1.1 8‐bp insertion using the primers MS002 and MS126 (Table S1). PCR was performed using Kapa Taq DNA polymerase (Kapa Biosystems), under the following cycling conditions: denaturation at 95°C for 5 min; 35 cycles of 95°C for 30 s, 60°C for 40 s, 72°C for 30 s; final extension at 72°C for 5 min. PCR products were digested with BseYI enzyme (New England BioLabs) according to the manufacturer's protocol and analyzed by gel electrophoresis to determine the genotype.

Growth conditions

Plants were grown under a greenhouse system as previously described (Torres‐Rodríguez et al., 2021). Sand was used as substrate with nutrient conditions maintained by use of a combination of solid‐phase P buffer (loaded with 209‐μM KH2PO4 for high P treatments and 11‐μM KH2PO4 for low P treatments; Lynch et al., 1990) and fertilization with modified Hoagland solution (5‐mM KNO3, .25‐mM Ca (NO3)2, 2‐mM MgSO4, 1‐mM KH2PO4, 20‐μM FeC6H6O7, 9‐μM MnSO4, 1.2‐μM ZnSO4, .5‐μM CuSO4, 10‐μM Na2B4O7, .008‐μM (NH4)6Mo7O24; Hoagland & Broyer, 1936). For low N treatments, Hoagland was adjusted by substitution of KNO3 and Ca (NO3)2 with KCl and CaCl2, respectively (Baxter et al., 2008; Escobar et al., 2006; Reddy et al., 1996). For low P treatments, Hoagland P concentration was adjusted by substitution of KH2PO4 with KCl (Diepenbrock, 1991). Hoagland solution was applied at 1/3 strength with final N/P concentrations in the different treatments as follows: Full treatment—1750‐μM NO3 −, 333‐μM PO₄3−; LowN—157.5‐μM NO3 −, 333‐μM PO₄3−; LowP—1750‐μM NO3 −, 10‐μM PO₄3−; LowNP—157.5‐μM NO3 −, 10‐μM PO₄3−.

For growth evaluation, plants were grown in PVC tubes (15‐cm diameter; 1 m tall) until 40 days after emergence (DAE). Tubes were filled with ~17 L of washed sand. In the upper third of the tube, sand was amended with 1.5% of the full or low P loaded solid‐phase P buffer. Four seeds were selected randomly from segregating seed stock, imbibed, and planted at 4‐cm depth per tube. The experiment was established over four planting dates spaced at 4‐day intervals. On each day, nine tubes were established per nutrient treatment, arranged as a 3 × 3 group. The resulting 16 groups were arranged in a Latin square with respect to nutrient treatment, although row number was confounded with planting date. Plants were thinned to a single plant a week after emergence, and the retained plants were genotyped for the pho1;2a′‐m1.1 mutation. Plants were irrigated with distilled water up until 10 DAE after which Hoagland treatments were applied, at a rate of 200 ml every third day. Plants were evaluated every 5 days by non‐destructive measurement of stem width, leaf number, leaf length, and leaf width for each fully expanded leaf. Leaf area was calculated as green leaf length × leaf width × .75 (Francis et al., 1969). Total leaf area was the sum of all leaves per plant. At 40 DAE, plants were carefully removed from the tubes, minimizing damage to the root system, washed in distilled water, and dried with paper towels. The cleaned root system was placed in a water‐filled tub and photographed using a digital Nikon camera D3000. Raw images were individually processed using Adobe Photoshop CC (Version 14.0) to remove the background and maximize the contrast between foreground and background non‐root pixels. Processed images were scaled and analyzed using GiA Roots software (Galkovskyi et al., 2012) as previously described (Torres‐Rodríguez et al., 2021). After photography, the root system was divided into segments corresponding to increments of 15‐cm depth (numbered 1 to 6) and each segment weighed individually. Shoot and root tissue were dried for 1 week at 42°C using a drying oven before taking dry weights and collecting tissue for element quantification. Trait description is provided in Data Set S1A.

For transcriptome analysis, plants were grown in shorter PVC tubes (15‐cm diameter; 50 cm tall). Substrate and nutrient treatments were the same as for the large tube experiment with the top 30 cm of the tube amended with 1.5% solid‐phase P buffer. At 25 DAE, the whole plant was harvested, separating the total root system, stem, and leaves. Leaf and root tissue for gene expression were immediately frozen in liquid nitrogen and stored at −80°C. Samples were homogenized with cooled mortar and pestle and aliquoted under liquid nitrogen for RNA extraction and transcriptome analysis.

For quantification of Pho1;2a and Pho1;2b transcripts, a third experiment was carried out with the same growth set‐up as the RNA‐Seq experiment. In this experiment, two treatments were applied, Full nutrients and LowP without the use of solid‐phase P buffer (alumina‐P). Plants were harvested at 25 DAE as described above.

Quantification of Pho1;2a and Pho1;2b transcripts

Total RNA was extracted using Trizol protocol (Invitrogen), and cDNA was synthesized using SuperScript II reverse transcriptase (Invitrogen) after DNase I treatment (Invitrogen) following the manufacturer's protocol. qRT‐PCR was performed using SYBR Green I‐based real‐time PCR reagent with the LightCycler 480 Instrument (Roche) using the following program: 95°C for 5 min, followed by 40 cycles of 95°C for 15 s; 60°C for 20 s; 72°C for 20 s. The ZmCDK gene (GRMZM2G149286) was used as a constitutive reference gene (Lin, Jiang, et al., 2014). Relative expression was calculated as 2ΔCt, where ΔCt = (Ct of constitutive gene − Ct of target gene). Three biological replicates with three technical replicates were analyzed. Statistical analysis was performed in R (R Core Team, 2022) using R/stats::lm in the model gene expression ~ genotype * treatment followed by Tukey's HSD test using R/stats::TukeyHSD. Primers used are indicated in Table S1.

Determination of elemental concentration by inductively coupled plasma mass spectrometry

Ion concentration was determined as described previously (Ramírez‐Flores et al., 2017). Briefly, dry leaf samples were analyzed by inductively coupled plasma mass spectrometry (ICP‐MS) to determine the concentration of 20 metal ions. Tissue samples were digested in 2.5‐ml concentrated nitric acid (AR Select Grade, VWR) with an added internal standard (20 ppb In, BDH Aristar Plus). The concentration of the elements Al, As, B, Ca, Cd, Co, Cu, Fe, K, Mg, Mn, Mo, Na, Ni, P, Rb, S, Se, Sr, and Zn were measured using an Elan 6000 DRC‐e mass spectrometer (Perkin‐Elmer SCIEX) connected to a PFA microflow nebulizer (Elemental Scientific) and Apex HF desolvator (Elemental Scientific). A control solution was run for every tenth sample to correct for machine drift both during a single run and between runs.

Statistical analysis of plant growth and ionomic data

For plants grown to 40 DAE, traits were obtained from 133 genotyped individuals across the wild‐type (wt), heterozygous (ht), and homozygous pho1;2a′‐m1.1 (mu/mt) genotypes, four nutrient treatments (Full, LowN, LowP, and LowNP), and four planting dates (Data Set S1B). Traits included direct measurements and derived values (e.g., total leaf surface area or biomass totals). Non‐destructive measurements were taken from 10 DAE and repeated at 5‐day intervals during the experiment. Destructive measurements were made for all 133 individuals at harvest. The data set included element concentrations determined by ICP‐MS and root architectural traits extracted by image analysis, as described above.

Statistical analysis was performed in R (R Core Team, 2022). Our focus was on differences between genotypes, and we performed separate analyses for each trait/treatment combination to limit model complexity. Traits were classified as dynamic (repeated measures) or endpoint (traits collected at harvest) and analyzed separately, adjusting for multiple testing within each group. Ionomic element concentrations were also analyzed separately. For visualization, we used R/stats::lm to fit the model trait value ~ 0 + genotype + planting date + error for each trait/treatment combination, extracting model coefficients and standard errors for plotting. For dynamic traits, we performed pairwise comparisons of wt, ht, and mu curves (trait over time) for each trait/treatment combination using R/statmod::compareGrowthCurves (Baldwin et al., 2007) under the meanT function, using 10,000 permutations for p value estimation; p values were adjusted using Holm's method with R/stats::p.adjust, both within each trait/treatment combination and globally. For endpoint analysis, trait values were square root transformed to improve normality and modeled as trait value  ~ planting date + genotype + error for each trait/treatment combination, extracting the genotype p value from the associated ANOVA table. Endpoint p values were adjusted for multiple testing using Holm's method with R/stats::p.adjust, applied separately to growth, GiA Roots, and ionomic data sets. Where the treatment effect was significant (adjusted p < .05), we used R/agricolae::HSD.test (de Mendiburu, 2020) to apply pairwise Tukey HSD tests at α = .1 to identify differences between genotypes.

RNA extraction and transcriptome sequencing

Transcriptome analysis was performed for each genotype (wt, ht, mt) in the four treatments (Full, LowN, LowP, and LowNP) from roots and leaf tissues represented by two biological replicates. RNA extraction and library generation were performed by Labsergen (Laboratorio de Servicios Genómicos, Langebio). Leaf libraries were sequenced using Illumina HiSeq4000 high‐throughput sequencing technology at the Vincent J. Coates Genomics Sequencing Laboratory (UC Berkeley). Root libraries were sequenced by Labsergen (Laboratorio de Servicios Genómicos, Langebio) using Illumina NextSeq 550 technology. Transcriptome data are available in the NCBI Sequence Read Archive under study SRP287300 at https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP287300. A separate analysis of the libraries generated from wild‐type segregants was published previously (Torres‐Rodríguez et al., 2021).

Analysis of differential gene expression

RNA sequencing reads were aligned to cDNA sequences from the AGPV3.30 version of the maize genome available at Ensembl Plants (ftp://ftp.ensemblgenomes.org/pub/plants/release-30/) using kallisto (version .43.1; Bray et al., 2016). Transcript abundance data were pre‐processed using R/tximport (Soneson et al., 2015) with gene‐level summarization. Transcript counts were analyzed with a GLM approach using R/edgeR (Robinson et al., 2010; McCarthy et al., 2012) and R/limma (Ritchie et al., 2015). To explore the impact of the pho1;2a′‐m1.1 mutation on gene regulation, we compared genotypes in separate analyses for each tissue (root or leaf) and treatment (Full, LowN, LowP, LowNP) combination (condition). We then pooled information across all analyses using R/mashr (Urbut et al., 2019) to generate the final list of differentially expressed genes (DEGs). Initial analyses were performed by generating a series of DGEList objects, one for each condition, each list containing data for wild‐type, heterozygous and pho1;2a′‐m1.1 plants. Normalization matrices and offsets were calculated for each DGEList as described in the R/edgeR documentation. Genes with fewer than 10 counts per million (CPM) in all, or all but one, of the libraries in a given DGEList were removed before normalization with R/edgeR::calcNormFactors. Data were prepared for linear modeling using R/limma::voom with no additional normalization. Global patterns of gene expression within each condition were visualized with R/edgeR::plotMDS. We estimated the effect of the genotype (GEN) at Pho1;2a for each condition using the linear model: where y is the normalized log2CPM value for a given gene in a given library, μ is the intercept for that gene in a particular condition, β is the effect of the genotype, and e is the model residual, which is assumed to be independent for each y and to have a variance proportional to Φ , the empirical weight factor calculated by voom. We fitted models (Model1) using R/limma::lmfit, sent the output to R/limma::eBayes, and extracted the estimates of β for both heterozygous and homozygous effects, along with their corresponding standard errors, calculated as the product of the standard deviation and the residual for each gene. We generated separate objects for heterozygous and homozygous effects, each containing effect estimates and standard errors for each gene in each of the eight conditions. Where we had no estimate for a given gene in a given condition (i.e., because of filtering for low expression), we set the estimate to 0 and the standard error to 1000. Estimates and standard errors were passed to R/mashr::mash using a canonical covariance matrix. For each gene, we extracted posterior estimates for heterozygous and homozygous mutant effects in each condition, along with standard deviations and local false sign rates (lfsr). For any given gene, we reset the posterior estimate to 0 for any condition for which the value had originally been set to NA because of low expression. We combined heterozygous and homozygous mutant effects into a single object and selected genes with posterior effects >|1|and lfsr < .01 as DEGs. The sign of the posterior effect was used to assign genes as up‐ or down‐regulated in any given comparison. DEGs identified in more than one comparison were called globally up‐ or down‐regulated based on the most common call, although only seven genes in our final set had both up‐ and down‐calls in different comparisons. Genes called as DEG in at least one condition for either heterozygous or homozygous effect were retained in our set of 1100 Pho1;2a‐regulated genes (Data Set S2A). Gene annotation was assigned as previously described (Gonzalez‐Segovia et al., 2019). Gene ontology (GO) analysis was performed with BiNGO (v. 3..3; (Maere et al., 2005) in Cytoscape (v. 3.7.2) using the set of 1100 Pho1;2a‐regulated genes in the context of 39,468 maize B73 RefGen_v3 gene models. Significant GO terms were obtained using the hypergeometric test with the Benjamini and Hochberg FDR correction and a significance level of .05 (Data Set S2B). We used posterior effects to visualize differential expression between genotypes. To visualize expression differences between nutrient treatments, we ran a second linear model (Model2) using all samples, considering each combination of nutrient treatment and genotype as a different level of a single grouping variable and extracting coefficients as a measure of normalized expression (Data Set S2A). Subsets of the 1100 Pho1;2a‐regulated genes were selected based on tissue/nutrient specific calls. Assignment of these genes as responsive to nutrient treatments was made using comparison of Model2 coefficients and taking an absolute difference >1 as evidence of differential expression.

RESULTS

The pho1;2a′‐m1.1 mutation introduces a premature stop codon

The maize pho1;2a′‐m1.1 allele (hereafter, pho1;2a) contains an 8‐bp (CTGCCCAG) insertion at the point of Activator (Ac) excision from the progenitor allele pho1;2a‐m1::Ac (Figure 1a; Salazar‐Vidal et al., 2016). This 8‐bp insertion shifts the reading frame relative to the wild‐type gene and generates a premature stop codon after residue 419 (based on the W22 gene model Zm00004b022763_T001), resulting in a mutant protein that is predicted to be truncated between the SPX and EXS domains (Figure 1a). In Arabidopsis, the EXS domain is crucial for subcellular localization of PHO1, Pi export, and root‐shoot signaling (Wege et al., 2016). We predict that the pho1;2a mutation knocks‐out native protein function, although we cannot exclude the possibility that a novel truncated SPX domain containing product is produced (Figure 1a). Our pho1;2a allele is similar to a recently reported CRISPR‐Cas9 editing event that deleted 2 bp in the sixth exon of Pho1;2a, shifting the reading frame after residue 391 and generating a premature stop codon at residue 627 (based on the W22 gene model Zm00004b022763_T001; Ma et al., 2021).

FIGURE 1

The pho1;2 am1.1 mutation introduces a shift in reading frame and a premature stop codon. (a) Structure of the ZmPho1;2a gene showing the site of transposon insertion (*) in the allele pho1;2a‐m1::Ac. The DNA sequence surrounding the insertion site in the wild type and the pho1;2a′‐m1.1 footprint allele is shown below, along with predicted protein sequence. Sequence adjacent to the insertion site is shown in blue. Eight base pairs of additional sequence generated at Ac insertion and retained following Ac excision are shown in red. The DNA and predicted protein sequences are based on the v2.0 of the W22 genome (Zm00004b022763). Black rectangles represent exons, black lines represent introns, and white rectangles represent UTRs. The black arrow indicates the position of a putative cis‐Natural Antisense Transcript. The positions of regions encoding the conserved SPX and EXS protein domains are marked. P1, P2 indicate the positions of primers used to quantify Pho1;2a sense transcripts in RT‐qPCR analysis. P3, P4 indicate primers used to quantify putative Pho1;2a antisense transcripts. (b) Quantification of Pho1;2a transcripts in the roots of wild type (WT, green), heterozygous (HT, yellow), and homozygous mutant (MT, brown) individuals of a family segregating the pho1;2′am1.1 mutation, grown under sufficient (Full) and or low phosphate (low P) treatment. RT‐qPCR was used to quantify relative expression in three individuals per condition. Transcript accumulation is shown relative to the maize CDK gene (GRMZM2G149286). (c and d) As (b), showing quantification of Pho1;2a cis‐NAT and Pho1;2b transcripts, respectively. In (b)–(d), significant differences among genotypes (G) or treatments (T) are indicated as *, p < .05; **, p < .01; ***, p < .01.

To evaluate the effect of the pho1;2a mutation on the Pho1;2a gene, we used quantitative PCR to assay Pho1;2a transcripts in a stock segregating pho1;2a. We included a primer set designed to amplify a previously described putative anti‐sense transcript (Salazar‐Vidal et al., 2016). We grew plants under both Full and LowP treatments (see Section 2) and assessed transcript levels at 25 DAE (Figure 1b–d). Although we could detect Pho1;2a sense transcripts in all genotypes under both P conditions, there was a dosage‐dependent reduction in accumulation associated with the pho1;2a allele (Figure 1b). Our quantitative PCR primers were located in the 3′‐UTR, indicating that full length transcripts are produced from pho1;2a. We did not detect a statistically significant effect of pho1;2a on the putative Pho1;2a antisense transcript (Salazar‐Vidal et al., 2016), although there was some evidence of increased accumulation in homozygous pho1;2a mutant plants under low P (Figure 1c). Consistent with our previous report (Salazar‐Vidal et al., 2016), the level of accumulation of Pho1;2b sense transcripts in the roots was greater than that of Pho1;2a (Figure 1d). We did not detect any effect of the pho1;2a mutation on Pho1;2b transcript accumulation (Figure 1d).

The pho1;2a mutation is linked to a reduction in seedling growth that is intensified under low nitrogen availability

There was no obvious whole‐plant phenotype associated with pho1;2a during generation and propagation of homozygous stocks in the field (cf. pho1 mutants in Arabidopsis, rice, or tomato; Poirier et al., 1991; Secco et al., 2010; Zhao et al., 2019). Despite the molecular similarity of our pho1;2a allele to previously reported CRISPR‐Cas9 events, we did not see a clear shrunken kernel phenotype (Ma et al., 2021). We weighed individual kernels from pho1;2a segregating stocks prior to planting for seedling characterization. Based on the seedling genotype data, there was no effect of the embryo/endosperm genotype on kernel weight (n = 133; mean kernel weight: wild‐type = 197 mg, heterozygote = 194 mg, mutant = 189 mg; Kruskal–Wallis p = .15). We also compared 100 kernel‐weight of a homozygous pho1;2a stock to that of a wild‐type sibling family. Here, genotype was significant, although potentially acting as a maternal effect (n = 5; mean 100 kernel weight: wild‐type = 22.5 g, mutant = 20.2 g; Kruskal–Wallis p = .01).

The SPX domain protein family is important in P nutrition and in crosstalk between P and nitrogen (N) signaling pathways (Hu et al., 2019; Ueda et al., 2020). We speculated that any pho1;2a phenotype might be enhanced under nutrient deficiency. We used controlled greenhouse conditions to characterize young (up to 40 DAE) plants from a stock segregating the pho1;2a mutation under different levels of N and P availability (Full, LowN, LowP, and combined LowNP; see Section 2). We followed plant growth by manual measurement of green leaf area (LA) every 5 days, starting at 10 DAE to 40 DAE, followed by endpoint measurements at harvest at 41 DAE. The first two leaves were fully expanded when measurements began at 10 DAE.

We observed a subtle reduction in seedling growth in pho1;2a mutants (Figures 2 and S1 and Data Set S1). Under Full and LowN conditions, LA of leaf (L)1 was smaller in homozygous mutants and heterozygotes than in wild‐type siblings at 10 DAE, with evidence of more rapid senescence in homozygous mutants and heterozygotes, especially under Full conditions (Figures 2a and S1 and Data Set 1C). We did not observe differences in L1 among genotypes under LowP or LowNP. L3 was initiated during the experiment, with plants reaching L8 or L9 by 40 DAE. Leaf growth was reduced in homozygous mutants and heterozygotes compared with wild‐type from L4 onwards under LowN, and from L6 under LowP (Figures 2a and S1 and Data Set S1C). Statistical support for growth effects was limited, with only effects on L1 and L2 under Full and L1‐5 under LowN being significant at 5%. Significant differences in Total LA (TLA) were observed under LowN and robust to adjustment for multiple testing across traits and treatments (Figure 2b). In all significant cases, homozygous mutants and heterozygotes were equivalent, and different to wild‐type, with the exception of L1 and L2 under LowN where heterozygotes followed wild‐type growth. Overall, our analysis indicated a subtle dominant effect of the pho1;2a mutation on vegetative growth that was enhanced under LowN availability. LA measurements were mirrored by endpoint shoot biomass (Figure S2—SFW LowN; Data Set S1D).

FIGURE 2

Seedling growth is reduced in pho1;2′am1.1 mutants. (a) Fully expanded green area (cm2) of first (LA1) and fifth (LA5) leaves from 10 days after emergence (DAE) until day 40. Data collected every 5 days from 10 DAE. Points show the coefficient estimated for each treatment, with bars extending ±1 standard error (SE). Colored polygons follow SE bars. (b) As (a) for total leaf area (TLA). Endpoint coefficients (±1 SE) for root occupancy (root fresh weight/total plant fresh weight as a percentage) in (c) the upper 45 cm (R0–45) of the soil and (d) below 45 cm (R45+) in the soil, for LowN or LowP treatments. Significant pairwise differences in Tukey HSD test (α = 0.1) are indicated by lowercase letters

To assess the importance of Pho1;2a in Pi translocation and more broadly in mineral nutrition, we used inductively coupled plasma mass spectrometry to quantify the concentration of total P and 19 further elements in leaf tissue harvested from the segregating mutant stock (Data Set S1E). There was no significant effect of pho1;2a on total P concentration in the leaf. A previously reported increase in leaf P concentration under lowN (Schlüter et al., 2013; Torres‐Rodríguez et al., 2021) was observed in all genotypes, further indicating that pho1;2a mutants are not compromised in the translocation of Pi from roots to leaves. Although we saw treatment effects on several elements, differences between genotypes were subtle, with limited statistical support that was not robust to correction for multiple testing (Figure S3 and Data Set S1E).

Root distribution is altered in plants carrying the pho1;2a mutation

The plants in our greenhouse evaluation were grown in 1‐m‐tall tubes, allowing us to characterize both shoot and root phenotypes. Total root biomass showed a mild reduction in homozygous pho1;2a mutants and heterozygotes compared with wild‐type plants, with the greatest statistical support under LowN (Figure S2—RFW LowN; Data Set S1D). Crown root number was reduced in plants carrying pho1;2a specifically under LowP (wild‐type: 15.02 ± 1.15, heterozygote: 14.21 ± .93, mutant: 12.14 ± 1.28; here and below, we give model coefficients and associated standard errors, Figure S2—CN LowP), although without strong statistical support. When we examined the distribution of roots by depth, mutants showed a significant reduction in root growth in the upper 45 cm of the soil column under LowN, in terms of both absolute biomass and proportion of total biomass (Figures 2c and S2 and Data Set S1D). A similar, although non‐significant, trend was also seen under LowP (Figure 2c). In contrast to the shoot growth traits, the action of the pho1;2a mutation on root biomass appeared additive, with heterozygous plants intermediate to homozygous mutants and wild‐type. Interestingly, in LowN and LowP, homozygous mutant and heterozygous plants showed greater root occupancy in the lower 45 cm of the soil column than wild‐type, both in terms of absolute mean biomass and mean proportion of total biomass (Figure 2d). The trend toward deeper rooting in mutants was reversed under LowNP, with homozygous mutant and heterozygous plants showing greater biomass than wild‐type in the top 45 cm and less than wild‐type in the lower 45 cm, the latter difference being detected as significant for both absolute and proportional biomass.

Root and leaf transcriptomes are modified in plants carrying pho1;2a

The Arabidopsis PHO1 gene is involved not only in Pi translocation but also in nutrient signaling and regulation of the PSR (Rouached et al., 2011). Although maize pho1;2a mutants lacked a strong Pi translocation phenotype, we hypothesized that subtle growth differences in the mutant might result from a signaling effect. To evaluate nutritional signaling in pho1;2a mutants, we grew a second batch of segregating plants under the same treatments used for growth evaluation. At 25 DAE (the point at which the nutrient treatments first influenced growth; Figure S1 and Data Set S1), we harvested leaf and root tissue from homozygous mutants, heterozygotes, and wild‐type plants and performed transcriptome analysis by RNA sequencing. To identify Pho1;2a regulated genes in our transcriptome data, we first compared genotypes within any given tissue/treatment combination and then pooled statistical support to identify genotype differences. Across all tissue/treatment combinations and considering both wild‐type to heterozygote (wt‐ht) and wild‐type to mutant (wt‐mt) comparisons, we identified a set of 1100 Pho1;2a‐regulated genes (Figure 3, Data Set S2, and Figure S5). Of these 1100 genes, 966 were up‐regulated in plants carrying pho1;2a with respect to wild‐type, 126 were down‐regulated, and 7 were observed to be up‐ or down‐regulated depending on the treatment/tissue combination. We saw a greater number of expression differences in roots than in leaves (Figure 3). The number of expression differences in roots was similar under all four nutrient treatments (Figure 3a) with a high proportion shared across treatments (Figure 3b). In leaves, differential expression was largely specific to LowP and LowNP treatments (Figures 3a and 3c). Patterns of differential expression were similar in wt‐ht and wt‐mt comparisons, indicating a dominant effect of pho1;2a on global transcript accumulation (Figures 3 and Figure S5).

FIGURE 3

The pho1;2′am1.1 mutation is linked to transcriptional misregulation. (a) Effect size estimates (log2 fold change; homozygous–wild‐type) for genes differentially expressed in plants homozygous for pho1;2a compared with segregating wild‐type siblings. Points corresponding to significant genes in any given tissue‐treatment combination are filled and outlined in bold. A single gene with an effect size >5 in the roots is not shown. (b) Venn diagram showing the overlap (count) of genes differentially expressed in the root (mutant–wild‐type) among the four nutrient treatments. The count of genes from the 1100 list not called in this genotype‐tissue combination is given outside the circles. (c) As (b), for genes differentially expressed in the leaf

We performed a GO analysis on the Pho1;2a regulated gene set and found enrichment for 620 GO terms, many implicated in stress responses including terms related to phytohormone (e.g., GO:9737 response to abscisic acid, GO: 9753 response to jasmonic acid) and nutrient (e.g., GO: 10167 response to nitrate, GO16036: cellular response to P starvation) signaling (Figure 4a–d and Data Set S2B). PHO1 has previously been linked to abscisic acid and jasmonic acid signaling in Arabidopsis (Huang et al., 2017; Khan et al., 2016; Zimmerli et al., 2012). In our data, genes associated with response to abscisic acid and response to jasmonic acid GO terms followed the global trend of upregulation in the roots of plants carrying pho1;2a across all treatments, and upregulation specifically under LowP and LowNP in the leaves (Figure 4a,b). For the terms response to nitrate and cellular response to P starvation, a number of the associated genes were upregulated in plants carrying pho1;2a specifically under LowP in both roots and leaves (Figure 4c,d).

FIGURE 4

The phosphate starvation response is enhanced in plants carrying pho1;2a. (a) Effect size estimates (standardized by row; z score) of the pho1;2a mutation on expression of called genes associated with abscisic acid (ABA) response (GO:9737). Columns represent conditions (root or leaf as stated; treatments indicated as green: Full, red: LowN, yellow: LowP, blue: LowNP; genotype as orange: heterozygous–wild type, brown: homozygous–wild type). Each row is a different gene. (b–d) As (a), for genes associated with jasmonic acid (JA; GO: 9753), nitrate (NO3 −; GO: 10167), and the cellular response to phosphorus starvation (PS; GO: 16036) responses, respectively. (e) Scatter plot of root expression coefficients comparing control (root full) and LowP (root low P) for a set of 581 genes identified as Pho1;2a‐regulated in LowP roots. Expression shown for homozygous (MT), heterozygous (HT), or wild type (WT) plants. A linear fit for each genotype is shown by the colored solid lines. The dashed red line marks equal expression in the two treatments. Marginal box plots indicate the distribution of values with respect to the facing axis. Where there was a significant genotype effect (Kruskal–Wallis test, p < .001), lower case letters indicate groups defined by Dunns test (p < .001). (f) As (e), showing a set of 366 genes identified as Pho1;2a‐regulated in LowP leaves. (g) Expression across treatments/genotypes (expression coefficient) of 276 genes identified as Pho1;2a‐upregulated in LowP leaves. The shape of the violin plot represents probability density and the bar the median values. The red dashed line shows the median for wild‐type plants under Full conditions. Lower case letters indicate groups defined by Tukey HSD at α = .05.

To investigate more closely the effect of pho1;2a on phosphate signaling and the PSR, we focused on genes that were Pho1;2a regulated under LowP and compared their expression across nutrient treatments in wild‐type and mutant plants (Data Set S2). Of the 1100 gene set, 581 genes were identified as Pho1;2a‐regulated in LowP roots. Overall, these 581 genes were expressed similarly under both LowP and Full conditions (Figure 4e). In contrast, the 366 genes of the 1100 gene set that were identified as Pho1;2a‐regulated in LowP leaves were largely not differentially expressed among genotypes under Full conditions—that is, starting from a common baseline in Full nutrient conditions, differences among genotypes were the result of a variable response to the LowP treatment (Figure 4f). Of 276 genes upregulated in the leaves of plants carrying pho1;2a under LowP, 153 (55%) were themselves part of the wild‐type leaf PSR (Log2 fold change >1; Figure S6 and Data Set S2). In heterozygous and homozygous mutant plants, 270 and 251, respectively, of the 276 Pho1;2a‐upregulated gene set were also induced by LowP treatment. The 276 gene set includes genes that, in addition to our transcriptome data, are implicated in the PSR by annotation, including the single SPX domain genes Spx1 (GRMZM2G083655), Spx3 (GRMZM2G024705), Spx5 (GRMZM5G828488), and Spx6 (GRMZM2G065989; Data Set S2; Torres‐Rodríguez et al., 2021; Xiao et al., 2021). Collectively, we interpret our data as evidence of an exaggerated transcriptional PSR in which the level of upregulation of responsive genes in plants carrying pho1;2a is higher than in wild‐type.

We have previously observed antagonism in wild‐type plants between transcriptional responses to P and N limitation, along with suppression of the PSR under a combined low NP treatment (Torres‐Rodríguez et al., 2021). This broader trend was conserved in the 276 gene set for wild‐type plants, with 56 genes down‐regulated under LowN and only 36 induced under LowNP, in contrast to the 153 induced under LowP (Data Set S2). In plants carrying pho1;2a, however, a high proportion of those genes induced under LowP were also induced under LowNP (208 in heterozygotes, and 174 in homozygotes). Consideration of the 276 gene set as a whole under the different nutrient treatments supported this general trend (Figure 4g). The genes observed to be upregulated under LowP in plants carrying pho1;2a were themselves LowP responsive in all genotypes, the level of LowP induction being higher in mutant plants. In wild‐type plants, the combined LowNP treatment suppressed upregulation, while in mutant plants, median expression remained higher than in the Full condition. There was evidence of downregulation of the 276 gene set under LowN. Interestingly, this signal was slightly stronger in plants carrying pho1;2a, suggesting that both the positive PSR and the antagonistic LowN repression of these genes was enhanced in the mutants.

DISCUSSION

We observed that an insertion of 8 bp in the 6th exon of the maize Pho1;2a gene was correlated with mild growth‐reduction and transcriptional misregulation, notably the upregulation of genes associated with the PSR. In Arabidopsis, the PHO1 protein functions in the translocation of Pi from roots to shoots and in P signaling (Hamburger et al., 2002; Rouached et al., 2011; Secco et al., 2010). Consequently, loss‐of‐function pho1 mutants in Arabidopsis show a dramatic reduction in shoot Pi concentration and growth (Poirier et al., 1991). A similar, non‐redundant role in Pi translocation is played by the orthologous proteins PHO1;2 and PHO1;1 in rice and tomato (Secco et al., 2010; Zhao et al., 2019), respectively, although shoot growth can recover by maturity in rice (Ma et al., 2021). The single copy rice PHO1;2 is represented by the paralogous gene pair Pho1;2a and Pho1;2b in maize (Salazar‐Vidal et al., 2016). Characterization of the pho1;2a allele presented here indicates that the maize PHO1;2a protein does not contribute significantly to Pi translocation. There was no significant reduction in leaf total P concentration in pho1;2a mutants nor any severe impact on growth, such as has been seen in other plant species. We observed Pho1;2b transcripts in the root to accumulate to levels ~40‐fold greater than those of Pho1;2a, consistent with Pho1;2b playing the primary role in Pi translocation. We saw no evidence of transcriptional up‐regulation of Pho1;2b in pho1;2a mutants, suggesting that there was no masking of a loss‐of‐function phenotype by the compensatory action of the paralogous gene. In mature plants, Pho1;2a and Pho1;2b have been shown to play non‐redundant roles in grain filling (Ma et al., 2021), and their expression levels are more nearly equivalent in the endosperm (Zhan et al., 2015). Interestingly, the previous study did not report any dramatic growth phenotype associated with either pho1;2a or pho1;2b single mutants, suggesting that double or higher‐order mutations in the maize pho1 family might be required to produce the Pi translocation phenotype seen in other plants.

We have interpreted our observations under the assumption that the pho1;2a′‐m1.1 allele conditions a complete loss of PHO1;2a protein function. Of course, this would be better demonstrated by isolation of additional mutant alleles. Similarly, the level of redundancy between the Pho1;2a and Pho1;2b paralogs cannot be fully addressed without generation of pho1;2b mutants. The CRISPR/Cas9 edit reported by Ma et al. (2021) is clearly an important additional reference, although the molecular nature of the Ma et al. pho1;2a allele—a small deletion in Exon 6—is very similar to the pho1;2a′‐m1.1 footprint event. PHO1 proteins contain both an N‐terminal SPX domain and a C‐terminal EXS domain. The molecular lesions in both pho1;2a′‐m1.1 and the Ma et al. pho1;2a allele are located downstream of the SPX domain‐encoding portion of the gene with the potential to generate truncated SPX containing proteins. In the absence of the EXS domain, the function of a truncated maize pho1;2a product would clearly be greatly compromised. Nevertheless, there is an entire class of single SPX domain proteins that play an important role in P signaling mediated through protein–protein interaction (Poirier, 2019; Torres‐Rodríguez et al., 2021; Xiao et al., 2021). Experimental manipulation has demonstrated the stability and nuclear localization of an isolated PHO1 SPX domain in Arabidopsis (Wege et al., 2016). Furthermore, overexpression of the SPX domain of Arabidopsis PHO1;H4 is sufficient to promote constitutive hypocotyl elongation, mimicking the dominant Pho1;H4 shb1‐D mutant (Zhou & Ni, 2010). The seedling growth and transcriptional effects we observed to be associated with pho1;2a were predominantly dominant. As for redundancy, isolation of additional mutant alleles will clarify the interpretation of dominance and dosage.

The presence of both Pho1;2a and Pho1;2b in the maize genome, coupled with their syntenic relationship to the single Pho1;2 gene of the closely related plant sorghum, suggests that the two maize genes are a paralog pair retained following a lineage specific whole‐genome duplication (Salazar‐Vidal et al., 2016). For the majority of genes, fractionation following whole‐genome duplication restores gene copy number to pre‐duplication levels (Schnable et al., 2011). The retention of Pho1;2a, despite a low level of expression relative to Pho1;2b in most tissues analyzed, suggests the gene to have functional significance, although the fractionation process may be ongoing. In addition to the reported role in grain filling (Ma et al., 2021), it is interesting to note that Pho1;2a has been identified as a candidate in several quantitative trait loci mapping and genome wide association studies, including studies of variation in P efficiency (Zhang et al., 2014) and root system architecture (Wang et al., 2019). Furthermore, Pho1;2a is located within a large‐scale inversion polymorphism that is characteristic of the native maize varieties and maize wild‐relatives of the central Mexican highlands (Crow et al., 2020) and has been associated with low temperature and low soil P availability (Aguirre‐Liguori et al., 2019). To our knowledge, Pho1;2b has not been associated with quantitative trait variation in maize. It may be significant that Pho1;2a has no major role in Pi translocation, freeing it up for possible sub‐ or neo‐functionalization or a role in adaptation to different environments.

In Arabidopsis, PHO1 plays a role in systemic P signaling (Rouached et al., 2011; Wege et al., 2016). PHO1 is not itself a transcriptional regulator and the implication is that PHO1 acts through maintenance and/or sensing of plant Pi status. We found evidence for transcriptional misregulation in plants carrying pho1;2a, predominantly examples of upregulation. Given that many of the Pho1;2a upregulated genes were associated with phytohormone or nutrient signaling GO terms, we interpret this result as evidence of a role for PHO1;2a in the repression of transcriptional stress responses. In the leaves, we saw little differential expression under full nutrient conditions, but genes that were induced by LowP in wild‐type were induced to greater levels in plants carrying pho1;2a, suggesting a disruption of feedback regulation. Arabidopsis plants with a partial loss of PHO1 function do not fully induce the PSR even though their rate of shoot Pi acquisition is at a level that would be low enough to induce the PSR in wild type plants (Rouached et al., 2011). Although, in our case, we observed an increased transcriptional response under low P in pho1;2a mutants, the data would be consistent with a general role for PHO1 in maintaining and sensing Pi concentration. As mentioned above, we cannot rule out the production of a truncated PHO1;2a product in mutant plants. Such a product would be missing the C‐terminal EXS domain, which may be significant given that we see leaf specific effects on transcription and that the EXS domain has been shown to play a role in root‐to‐shoot signaling in Arabidopsis (Wege et al., 2016). It was also interesting to observe that in roots, the presumed primary site of PHO1;2a accumulation in vegetative stage plants, there was significant transcriptional misregulation in pho1;2a mutants under all four nutrient treatments while, in the leaves, the effects were specific to LowP and combined LowNP treatments. It has previously been hypothesized that Arabidopsis PHO1 might transport a systemic signal indicating root Pi status (Rouached et al., 2011). The conditional effect of pho1;2a on the leaf transcriptome would be broadly consistent with this model.

In rice, PSR induction is contingent on sufficient nitrate availability to promote turnover of the negative regulator SPX4 by the NRT1.1b nitrate transporter (Hu et al., 2019). In our data, genes induced in wild‐type plants under LowP were not responsive under LowNP, suggesting an analogous mechanism to be acting in maize. We have previously observed that even a mild reduction in nitrate availability is sufficient to prevent transcriptional upregulation of PSR genes under P limitation (Torres‐Rodríguez et al., 2021). In plants carrying pho1;2a, the majority of genes that were induced in LowP remained induced in LowNP, suggesting that aspects of NP crosstalk were modified.

In summary, we show that full‐length PHO1;2a does not play a major role in root to shoot Pi translocation in maize, although the pho1;2a allele is associated with transcriptional misregulation indicating a role in nutrient signaling. In the future, it will be informative to evaluate additional alleles of pho1;2a and the role of the paralogous Pho1;2b gene. Characterization of pho1;2a; pho1;2b double mutants will provide a better understanding of functional divergence and redundancy between the two maize Pho1;2 genes.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Data Set S1A. Description of the traits collected to characterize the growth development of the pho1;2a´‐m1.1 segregating plants. Plants with a segregating pho1;2a´‐m1.1 mutation were evaluated for 40 DAE under four nutrient conditions to follow seedling growth by manual measurements. The data represents direct and non‐destructive measurements during the experiment. Also, derived values and destructive measurements calculated at the endpoint of the experiment. Traits abbreviations, name, units of measure and descriptions are displayed.

Click here for additional data file.

Data Set S2A. Set of Pho1;2a‐regulated genes. Transcriptome analysis was performed in roots (R) and leaves (L) of wild type (WT), heterozygous (HT) and homozygous pho1;2a´‐m1.1 mutant (MT) plants. Plants were grown under four nutrient conditions (Full, LowN, LowP, LowNP) until 25 DAE. Mixture model output (Model1) shows posterior estimates (pm) for HT and MT effects (Log2 fold change) in each condition, along with standard deviations (sd) and local false sign rates (lfsr). Differentially expressed genes (DEGs) was selected according to|pm| > 1 and lfsr < .01. Genes called as DEG in at least one condition for either HT or MT effect are shown. Calls flag DEGs +/− 1 for up‐ or down‐ regulated, respectively, based on the sign on the pm. Gene annotation was assigned as previously described (Gonzalez‐Segovia et al., 2019). A second linear model (Model2) was run to visualize expression differences between conditions, considering each combination of condition (Full, LowN, LowP, and LowNP) and genotype (WT, HT and MT) as a different level of a single grouping variable and extracting coefficients. Gene sets for Figure 4 and associated discussion were selected based on DEG Calls. Expression values were based on Model2 output. Nurtient treatment effects (up‐ or down‐regulation) discussed in gene sets were called based on an adsolute difference of Model2 coefficients >1.

Click here for additional data file.

Figure S1. Seedling growth is reduced in A) Fully expanded green area (LA; cm2; y‐axis) of leaves 1 to 7 from 10 days‐after‐emergence (DAE; days; x‐axis) until day 40, under the four nutrient treatments. Data collected every 5 days from 10 DAE. Points show the coefficient estimated for each treatment, with bars extending +/− 1 standard error (SE). Colored polygons follow SE bars. Total leaf area (TLA) calculated as the sum of individual LA measurements. Green shading: wild‐type; Yellow shading: heterozygous pho1;2a; Brown shading: homozygous pho1;2a.

Click here for additional data file.

Figure S2. Seedling growth is reduced in Endpoint coefficients (+/− 1 SE) of traits at harvest, under the four nutrient treatments, for wild‐type (WT), heterozygous (HT) and homozygous pho1;2a (MT) plants. Trait codes and measurement units are given in Supplemental Data Set 1. Significant pairwise differences in Tukey HSD test (α = .1) are indicated by lowercase letters.

Click here for additional data file.

Figure S3. Leaf ionome nutrient responses are modulated in Total concentration of selected elements in the leaf displayed as heatmaps of standardized values (z value, standardized across genotypes for each element/treatment combination) and as plots (ppm). The genotype term was significant in the element/treatment combinations shown as *** p < .001, ** p < .01, * p < .05, p < .1 (value before adjustment for multiple testing). Significant pairwise differences (p < .05) recovered from Tukey analysis are indicated with lowercase letters. Plots show estimated coefficients with bars extending +/− 1 standard error (SE).

Click here for additional data file.

Figure S4. Giaroots traits at harvest Coefficients (+/− 1 SE) for Giaroots traits at harvest, under the four nutrient treatments, for wild‐type (WT), heterozygous (HT) and homozygous pho1;2a (MT) plants. Trait codes and measurement units are given in Supplemental Data Set 1. Significant pairwise differences in Tukey HSD test (α = .1) are indicated by lowercase letters.

Click here for additional data file.

Figure S5. Plants heterozygous for the A) Effect size estimates (log2 fold change; heterozygous ‐ wild‐type) for genes differentially expressed in plants heterozygous for pho1;2a compared with segregating wild‐type siblings. Only significant genes are shown in any given tissue‐treatment combination. Numbers at the base of the plot show the count of differentially expressed genes. A gene with an effect size >5 in the roots is not shown. B) Venn diagrams showing the overlap (count) of genes differentially expressed in the root (heterozygous pho1;2a compared with segregating wild‐type siblings) among the four nutrient treatments. C) as B, for genes differentially expressed in the leaf.

Click here for additional data file.

Figure S6. The majority of genes mis‐regulated in Counts of genes called as upregulated in pho1;2a leaves with respect to wild type under LowP; upregulated in WT leaves in Low P with respect to Full; upregulated in the leaves of pho1;2a homozygotes with respect to Full.

Click here for additional data file.

Table S1. Primers used to quantify Pho1;2a and Pho1;2b transcripts

Click here for additional data file.