Seeds of 83 worldwide natural A. thaliana accessions (Atwell et al., 2010; Cao et al., 2011; Horton et al., 2012) used in our previous GWAS (Sadhukhan et al., 2020; Nakano et al., 2020b) and T-DNA insertion lines were obtained from the Arabidopsis Biological Resource Center (ABRC, Columbus, OH, United States), the Nottingham Arabidopsis Stock Centre (NASC, Nottingham, United Kingdom), and the RIKEN BioResource Center (RIKEN BRC, Tsukuba, Japan). Prior to experimental use, the procured seeds were multiplied by a single-seed descent process. Homozygosity was confirmed in T-DNA insertion mutant line using primers from the SALK database, following their protocols.¹ The sequences of the primers are given in Supplementary Table S1. The T-DNA insertion mutants used in this study were pgip1 (SALK_001662), STOP1-KO (SALK_114108), at1g64105 (SAIL_235_D03), at2g21950 (GK-707C08), at3g24480 (GABI_017A08), at5g15300 (SALK_044494C), at5g38900 (SAIL_453_G03), at5g43460 (SALK_067877C), at5g58900 (SALK_084867C), at5g58910 (SALK_064093C), at1g51070 (SALK_104253C), at2g38090 (SALK_127250C), at2g04780 (SALK_113729C) and plc9 (SALK_021982C).

Seedlings were grown on nylon mesh floating on modified MGRL solution (Fujiwara et al., 1992; 2% solution with 200 μM CaCl₂; initial pH 5.6) for 10 days at 22°C using a 12-h photoperiod with 37 μMol m⁻² s⁻¹ photon flux density. The culture solution was renewed every 2 days. After 10 days, the seedlings were transferred to another modified MGRL solution (without P and pH 5.0) containing 25 μM AlCl₃·6H₂O (Sawaki et al., 2016). The shoots were harvested after 24 h of Al treatment and immediately frozen with liquid nitrogen for RNA extraction. The same Al toxic solution containing 50 μM 2-(4-carboxyphenyl)-4,4,5,5-tetramethyl-imidazoline-1-oxyl-3-oxide (cPTIO), an NO scavenger, was used to evaluate the effect of NO (D’Alessandro et al., 2013; Sun et al., 2016). The concentration of cPTIO was based on preliminary experiments from which maximum suppressed responses were obtained without affecting the plant root (Supplementary Figures S1, S2).

Soil culture was conducted using commercial acidic soil (PROTOLEAF, Tokyo, Japan; pH 4.2; 1:2.5 w/v soil/water solution; Agrahari et al., 2020). The acidic soil was neutralized by the addition of CaCO₃ (Koyama et al., 2000; Sawaki et al., 2016; Agrahari et al., 2020; 4.0 g kg⁻¹; pH 5.1; 1:2.5 w/v soil/water solution) and used as the control soil. Plants (100 seeds) were grown for 2 weeks at 22°C during a 12-h photoperiod with 37 μMol m⁻² s⁻¹. Throughout the experiment, the plants were irrigated daily with deionized water to maintain soil moisture. The soil pH (water) and exchangeable Al were determined using the method described by Koyama et al. (2000).

Total RNA was isolated from the shoots using Sepasol-RNA I Super G (Nacalai Tesque, Kyoto, Japan) according to the manufacturer’s instructions. The RNA quality was analyzed using the A260/A280 ratio on a NanoVue Plus spectrophotometer (Biochrom, Holliston, United States). Total RNA was reverse-transcribed using ReverTra Ace quantitative PCR master mix with genomic DNA remover (Toyobo, Osaka, Japan) following the manufacturer’s instructions. The gene expression levels were quantified using SYBR Premix Ex Taq II (Takara Bio, Otsu, Japan) with a Dice Real Time System II MRQ thermal cycler (Takara Bio, Otsu, Japan) according to the manufacturer’s instructions. Briefly, all quantifications were carried out based on the real-time quantitative reverse transcription PCR (qRT-PCR) standard curve method of Bustin et al. (2009), as described by Kobayashi et al. (2014). For all quantifications, a standard curve was constructed using a cDNA dilution series, and the transcript levels of selected genes were quantified relative to that of the stable internal reference gene, Ubiquitin 1 (UBQ1; AT3G52590; Kobayashi et al., 2007, 2014). We have checked the invariant expression of UBQ1 for all experimental condition in used lines in this study (Supplementary Figure S3). We included a control with no reverse transcriptase to assess genomic DNA contamination, and the amplification efficiency of all primers was confirmed. The primer sequences used to amplify the selected genes are shown in Supplementary Table S2.

The eGWAS analyses were carried out using TASSEL v3.0 software following a mixed linear model (MLM; Bradbury et al., 2007) using a total of 160,748 genome-wide single nucleotide polymorphism (SNP) information from public databases (Atwell et al., 2010; Cao et al., 2011; Horton et al., 2012),^{2 3} which excluded SNPs of missing data or those with less than 5% minor allele frequency, as described earlier (Nakano et al., 2020b). The heritability (h²) was estimated by the following formula h² = (the additive genetic variance)/(the additive genetic variance + the residual variance). The suggestive SNPs were determined by quantile–quantile (Q–Q) plot analysis (Zhang et al., 2018) using a free statistics software,⁴ and the genes closest to the SNPs (Table 1) were identified using the TAIR 10 database.⁵

Gene ontology (GO) analysis was performed using an online tool available in the TAIR database.⁶ Gene polymorphisms were mined from the 1,001 Genomes database⁷ and POLYMORPH database.⁸ The genes upregulating PGIP1 expression in T87-cultured cells of Arabidopsis were identified by the Regulatory-network Research (RnR) database (Sakurai et al., 2014).⁹ Co-expression network analysis was conducted on the eGWAS-detected and RnR database-listed genes using the ATTED-II database (Obayashi et al., 2018).¹⁰ Cis-elements and corresponding transcription factors (TFs) of the promoters were predicted using PlantPAN 3.0 (Chang et al., 2008).¹¹

The STOP1 complementation Arabidopsis transgenic plant was constructed as described by Ohyama et al. (2013). STOP1 genomic DNA containing the promoter (−2,848 from the first ATG) and downstream (+626 from the stop codon) regions was cloned into a binary vector (promoterless pBIG2113SF). This construct was introduced into Agrobacterium tumefaciens strain GV3101 and transformed into STOP1-KO plants by the floral dip method (Clough and Bent, 1998). A T3 homozygous line was used for the experiments.

The binding of STOP1 to double-stranded, synthetic promoter fragments was studied using an amplified luminescent proximity homogeneous assay (AlphaScreen^™, PerkinElmer, Waltham, MA, United States) and a 276 EnSpire Multimode Plate Reader (PerkinElmer) as described by Enomoto et al. (2019). Competition assays were performed with non-biotinylated probes (450 nM) according to the manufacturer’s instructions. The competitor probes, consisting of the promoter fragments −193 to −222 bp upstream from the PGIP1 and − 2,694 to −2,723 bp upstream from the AT5G38900 start codons, respectively, were designed around the STOP1-binding site according to the Plant Cistrome Database (Sadhukhan et al., 2019).¹² The mutated probes were designed following the method of Tokizawa et al. (2015). The forward and reverse probe sequences are listed in Supplementary Table S3.

Pectin was extracted from 250 mg powdered Arabidopsis shoot tissue (500 seedlings were grown for 10 days in the hydroponics system mentioned above) in buffer containing 50 mM Tris–HCl (pH 7.2) and 50 mM cyclohexane-trans-1, 2-diamine tetra-acetate (CDTA; Bethke and Glazebrook, 2014). The extraction was continued for 15 min at 95°C with intermittent vortexing, and the sample was then centrifuged at 10,000 × g for 10 min. The supernatant containing pectin was analyzed for Al content using inductively coupled plasma mass spectrometry as described by Watanabe et al. (2015). For nuclear magnetic resonance (NMR) analysis, the supernatant containing pectin was lyophilized and dissolved in D₂O. ¹H-¹³C-Heteronuclear Single Quantum Coherence (HSQC) NMR was performed at 600.17 MHz on a JEOL ECA 600 NMR spectrometer (JEOL, Tokyo, Japan) equipped with a 5-mm FG/TH tunable probe, using the pulse sequence ‘hsqc_dec_club_pn’. NMR measurements were recorded at 70°C (Siedlecka et al., 2008). Sweep widths of 15 and 170 ppm were used to acquire the ¹H and ¹³C spectra, respectively. For each NMR experiment, 88 scans were collected using a relaxation delay of 1.5 s.

We analyzed the expression of PGIP1 in the shoots of wild-type (WT) Arabidopsis, Columbia (Col-0) at different time points after root exposure to 25 μM Al. The expression of PGIP1 was induced after 12 h and was markedly induced (five-fold on average) after 24 h of treatment (Figure 1A). It is also a condition of Al accumulation in the shoots (Supplementary Figure S4). Hence, in this study, we chose 24 h as the time point for evaluating PGIP1 gene expression under Al stress to conduct the eGWAS of PGIP1 in the shoots. Next, we analyzed the expression of PGIP1 in 83 Arabidopsis accessions (Supplementary Figure S5; Supplementary Table S4) that had been treated with treated with 25 μM Al for 24 h and found that the expression range between log₂ RFC −3.3 and log₂ RFC 0.1 (RFC: relative fold change; compared to Col-0; Figure 1B; h² = 88.3%). We performed a GWAS following MLM using the PGIP1 expression data, but after Bonferroni correction for multiple testing, we could not detect any significant SNPs at the genome-wide significance level. This could be due to the dependence of the statistical power of GWAS on the population size and allele frequency. Although the Q–Q plot showed only a slight deviated plot, in this study we set a suggestive threshold based on this result (p < 10^–3.5) and selected the top-ranked 17 SNPs as suggestive SNPs associated with PGIP1 expression level variation (Figures 1C,D; Table 1). These potentially associated SNPs were selected for further analysis.

We first characterized 17 protein-coding genes carrying the 17 SNPs directly in their exons, introns, and promoters (Table 1) and then considered a priori candidate genes related to the regulation of PGIP1 expression out of these 17 genes. Out of these, seven genes belonged to the GO term of “regulation of gene expression,” “DNA-binding transcription factor activity,” “intracellular signal transduction,” and “hormone-mediated signaling pathway” and may be functionally associated with PGIP1 expression (Table 1; Supplementary Figure S6). Out of the a priori candidate genes, AT2G21950 (SKP1 interacting partner 6: SKIP6; Farrás et al., 2001) and AT3G24480 (leucine-rich repeat extension 4: LRX4; Zhao et al., 2018) are involved in hormonal signaling and cell-wall integrity, respectively. AT5G38900 (Thioredoxin superfamily protein: TRX SF) is involved in regulation of gene expression via redox signaling (Sevilla et al., 2015). On the other hand, five among the 17 genes were induced more than 1.2-fold by Al in the shoot, according to our earlier transcriptome analysis, carried out under the same experimental conditions (Sawaki et al., 2016; Table 1). Therefore, we finally selected 12 genes with functions related to transcriptional regulation and Al-responsive expression, as a priori candidate genes associated with PGIP1 expression, for subsequent analyses.

We examined the 12 candidate genes for expression level polymorphisms (ELPs) and amino acid polymorphisms associated with the PGIP1 expression level. The expression levels of these genes were compared between representative accession groups that carried different detected SNP alleles (Figure 2, Supplementary Figure S7). Five genes, AT1G64105 (NAC027), AT5G38900 (TRX SF) AT5G43460 (HR-like lesion-inducing protein-like protein: HR-like protein), AT5G58900 (Homeodomain-like transcriptional regulator: R-R MYB), and AT5G58910 (laccase 16: LAC16) exhibited significant differences in the level of expression between accessions carrying different SNP alleles (Figure 2). When comparing the different alleles, the minor allele group with elevated PGIP1 expression showed higher expression levels of each gene than the major allele group with low PGIP1 expression (Table 1; Figure 2). In the case of genes with SNP alleles located directly in exons, we examined the amino acid polymorphisms caused by the detected SNPs using reliable DNA sequences from the 1,001 genome database. The SNPs detected in the exons of five genes caused amino acid polymorphisms, viz. Ser29Cys in NAC027, Ala187Ser in SKIP6, Asn60Lys in LRX4, Phe134Ile in AT5G15300 (pentatricopeptide repeat superfamily protein) and His246Gln in R-R MYB. Among them, NAC027 and R-R MYB showed both ELP and amino acid polymorphism. These polymorphisms might affect the expression level variation of PGIP1 as cis-factors. In this way, eight genes were selected as the first group of possible candidate genes, which were further studied using reverse genetics.

To examine the effect of the candidate genes on PGIP1 expression, its expression level was quantified in the T-DNA insertion knockout (KO) or knockdown mutants (KD; Supplementary Figure S8) of the eight candidate genes. PGIP1 expression was significantly lower in the three T-DNA insertion mutants of nac027, trx sf and r-r myb than in the WT under Al treatment (Figure 3A). In particular, the reduction of PGIP1 expression in the r-r myb was as high as approximately 40%. In contrast, other mutants (at2g21950, at3g24480, at5g15300, at5g43460, at5g58910) showed no significant difference of PGIP1 expression compared to the WT under Al treatment (Figure 3A). From this analysis, we found that NAC027, TRX SF and R-R MYB are involved in the regulation of PGIP1 expression.

We also analyzed the patterns of expression of NAC027, TRX SF, and R-R MYB under Al stress in the shoots of WT A. thaliana after 24 h of treatment. We observed a significant induction of these genes relative to the control (Figure 3B). The Al-induced responses of NAC027 were weak, but their gene expression showed Al responses similar to PGIP1 expression. These results suggest that NAC027, TRX SF, and R-R MYB are related to the regulation of Al-induced PGIP1 expression. In addition, Al-responsive PGIP1 and its regulatory system were involved in Al tolerance. In acidic soils containing higher exchangeable Al, pgip1, like the Al-hypersensitive stop1, was much more inhibited in growth than the WT (Supplementary Figure S9). On the other hand, the growth was recovered in the neutralized soil (Supplementary Figure S9). Similarly, the growth of trx sf and r-r myb was inhibited than the WT in acidic soil, although the growth of nac027 was not severely inhibited (Supplementary Figure S9). This may be consistent with the lower degree of repression of PGIP1 in nac027 and the weaker induction of Al on the expression level of NAC027 compared to the other two genes (Figure 3).

We investigated the relationship between STOP1 regulation and the eGWAS-detected genes because it has been reported that STOP1 regulates PGIP1 expression (Sawaki et al., 2009). There was no difference in expression levels of STOP1 in the nac027, trx sf, and r-r myb (Supplementary Figure S10A). In contrast, we found a significant reduction in the gene expression level of TRX SF in the STOP1-KO, whereas the other two genes showed similar expression levels compared to WT (Figure 4A). In addition, in the STOP1-complemented line, the expression of TRX SF was fully recovered similar to Col-0 (Figure 4B). These results suggest that TRX SF is regulated by the STOP1 transcription factor, while none of the three genes affect the expression level of STOP1.

In our previous eGWAS of ALS3, we found the involvement of phosphoinositide (PI)-dependent phospholipase C9 (PLC9) signaling upstream of the STOP1 regulation system, which regulates the expression of ALS3 and PGIP1 in the shoots of A. thaliana (Sadhukhan et al., 2020). Therefore, we examined the expression of TRX SF, NAC027, and R-R MYB in the plc9. We found that the expression of TRX SF was significantly suppressed in plc9 (Figure 4A), similar to the downregulation of PGIP1 (Supplementary Figure S10B). These findings suggest that TRX SF is regulated by a PI signaling-mediated STOP1-dependent pathway. In contrast, we found that the expression of ALS3 remained unchanged in the trx sf, nac027, and r-r myb mutants (Supplementary Figure S10A). These results suggest that TRX SF differentially regulates transcription of PGIP1 and ALS3, which are co-regulated by STOP1 in Arabidopsis shoots.

Next, we performed a promoter binding analysis to reveal whether STOP1 directly regulates TRX SF and PGIP1 in the STOP1-dependent pathway. The STOP1-binding positions were predicted by searching enriched sequences (octamer units) in the stress-inducible promoters.¹³ They were identical to the binding sites provided by the Plant Cistrome Database. Therefore, we searched for putative STOP1-binding sites in the promoters of PGIP1 and TRX SF using the Plant Cistrome Database. Based on DNA affinity purification sequencing (DAP-seq), the database identified the “GGNVS” consensus sequence in the PGIP1 and TRX SF promoters for binding STOP1-like proteins (Figure 5A), as previously identified in rice by Tsutsui et al. (2011). The binding capacity of STOP1 to the sequences available from the Plant Cistrome Database was validated by an in vitro competitive binding assay using the AlphaScreen^™ system (Figure 5B). A 30-bp synthetic double-stranded DNA probe, designed around the binding site in the PGIP1 (−193 to −222 bp) and TRX SF (−2,694 to −2,723) promoters, could compete for STOP1 protein binding with a known STOP1-binding site in the AtALMT1 promoter (Tokizawa et al., 2015). On the other hand, replacing the “GGNVS” consensus sequences within the PGIP1 and TRX SF probe with A/T stretches abolished STOP1 binding (Figure 5B). These results indicate that STOP1 binds directly to the PGIP1 and TRX SF promoters.

NO generation is positively correlated with cell wall pectin demethylation and alteration of cell wall metabolism under Al stress (Zhou et al., 2012; Sun et al., 2016). To establish whether Al accumulation induces NO signaling, we examined NO-inducible marker gene (AT2G06050, AT3G45140, and AT5G42650; Huang et al., 2004) expression in shoots after exposing the roots to Al or Al plus cPTIO (NO scavenger; Shi et al., 2015) for 24 h. We found that the expression of the NO marker genes was substantially induced in the Al-treated samples, whereas their expression was suppressed in the samples treated with Al plus cPTIO (Figure 6A). This suggests that Al induces NO signaling in the shoots. Next, we examined whether PGIP1 and ALS3, which are regulated by STOP1, function downstream of NO signaling by quantifying their expression levels in the shoots after exposing the roots to Al or Al plus cPTIO for 24 h. The transcript levels of PGIP1 were suppressed significantly in the Al plus cPTIO samples, whereas ALS3 expression was unchanged (Figure 6B). STOP1 expression was neither Al-induced nor affected by cPTIO.

We also assessed whether the eGWAS-identified genes that regulate PGIP1 expression function downstream of NO signaling by quantifying their expression levels in the shoots after exposing the roots to Al or Al plus cPTIO for 24 h. We found that only NAC027 expression and R-R MYB expression were significantly suppressed in the Al plus cPTIO samples, while TRX SF remained unchanged (Figure 6B). These results clearly indicate that Al-inducible PGIP1 expression is regulated by the NO signaling pathway through NAC027 and R-R MYB, and this is not regulated by STOP1 (Figure 4A).

The expression of NAC027 and R-R MYB was measured in each KO line. A significant reduction in the expression of NAC027 was observed in the r-r myb compared with the WT expression (Figure 7A). In contrast, the expression of R-R MYB in the nac027 showed no difference compared with wild type (Figure 7B). In addition, the PlantPAN3.0 database (provides TF binding sites in genome-wide promoters based on DAP-seq analysis of various transcription factors) identified that R-R MYB directly binds to the promoter of NAC027 (Figure 7C).

Interestingly, we found regulators of PGIP1 using the RnR database that is different from the eGWAS-identified factors. The RnR database indicated that overexpression of AT1G51070 [BASIC HELIX-LOOP-HELIX 115 (BHLH115)], AT2G38090 (duplicated homeodomain-like superfamily protein), and AT2G04780 [FASCICLIN-LIKE ARABINOGALACTAN 7 (FLA7)] upregulated PGIP1 expression, which was at the approximately 99th percentile of expression regulation and about 1.3–2.0-fold expression compared with the control. In fact, we confirmed the downregulation of PGIP1 expression in the KO or KD mutants of these three genes under Al stress (Figure 7D). Among them, AT2G38090 is a member of the R-R-type MYB family (Yanhui et al., 2006) and is a close homolog of R-R MYB that we found in the eGWAS. This TF also binds to the NAC027 promoter as assessed by a PlantPAN3.0 database search, similar to R-R MYB (Figure 7C). Co-expression gene network analysis using the six genes (TRX SF, NAC027, R-R MYB, BHLH115, AT2G38090, and FLA7), whose PGIP1 expression levels were decreased in these KO or KD lines, revealed that R-R MYB was co-expressed with AT2G38090 (R-R type MYB) and FLA7 (Supplementary Figure S11). These results suggest that these genes are involved in the direct/indirect regulation of PGIP1 expression.