ALUMINUM RESISTANCE TRANSCRIPTION FACTOR 1 (ART1) contributes to natural variation in aluminum resistance in diverse genetic backgrounds of rice (O. sativa)

Abstract Transcription factors (TFs) regulate the expression of other genes to indirectly mediate stress resistance mechanisms. Therefore, when studying TF‐mediated stress resistance, it is important to understand how TFs interact with genes in the genetic background. Here, we fine‐mapped the aluminum (Al) resistance QTL Alt12.1 to a 44‐kb region containing six genes. Among them is ART1, which encodes a C2H2‐type zinc finger TF required for Al resistance in rice. The mapping parents, Al‐resistant cv Azucena (tropical japonica) and Al‐sensitive cv IR64 (indica), have extensive sequence polymorphism within the ART1 coding region, but similar ART1 expression levels. Using reciprocal near‐isogenic lines (NILs) we examined how allele‐swapping the Alt12.1 locus would affect plant responses to Al. Analysis of global transcriptional responses to Al stress in roots of the NILs alongside their recurrent parents demonstrated that the presence of the Alt12.1 from Al‐resistant Azucena led to greater changes in gene expression in response to Al when compared to the Alt12.1 from IR64 in both genetic backgrounds. The presence of the ART1 allele from the opposite parent affected the expression of several genes not previously implicated in rice Al tolerance. We highlight examples where putatively functional variation in cis‐regulatory regions of ART1‐regulated genes interacts with ART1 to determine gene expression in response to Al. This ART1–promoter interaction may be associated with transgressive variation for Al resistance in the Azucena × IR64 population. These results illustrate how ART1 interacts with the genetic background to contribute to quantitative phenotypic variation in rice Al resistance.


| INTRODUCTION
Plants achieve abiotic stress resistance through diverse mechanisms that evolved through both natural and artificial selection. The characterization of quantitative trait loci (QTL) and their underlying genes has led to insights about the genetic architecture and the biological mechanisms of resistance to several abiotic stresses. Two major classes of genes predominate as determinants of stress resistance: membrane transporters and transcription factors (TFs) (Mickelbart, Hasegawa, & Bailey-Serres, 2015). The reprogramming of gene expression through transcriptional regulation, mediated by TFs, is a finely orchestrated and tightly regulated process, and is one of the hallmarks of plant response to stress (Vaahtera & Brosch e, 2011). To achieve specificity in the transcriptional response, a TF must activate only target genes involved in adaptation to a particular stress or combination of stresses (Vaahtera & Brosch e, 2011). TFs are activated through posttranslational modification and nuclear import and/ or may be transcriptionally activated by other TFs. Genetic variation in the TF loci, or in any element along the regulatory chain, can contribute to phenotypic variation in plant stress resistance. Understanding the quantitative nature of TF-mediated stress resistance in plants is essential when predicting the impact of introducing alleles into a new genetic background in the context of plant breeding and crop improvement.
Aluminum toxicity severely limits plant growth on acidic soils (pH <5). Under these conditions, the rhizotoxic Al species Al 3+ is solubilized, inhibiting root growth and function (Kochian, Piñeros, Liu, & Magalhaes, 2015), and leaving plants more vulnerable to drought and mineral nutrient deficiencies. Approximately 30% of the earth's land area consists of highly acid soils, and as much as 50% of all potentially arable lands are acidic (von Uexk€ ull & Mutert, 1995). As vast areas of acid soils in the tropics and subtropics are critical foodproducing regions, Al toxicity constitutes a food security threat exceeded only by drought among the abiotic limitations to crop production.
Rice (Oryza sativa L.) is the most Al-resistant species among the small grain cereals (Foy, 1988;Ma et al., 2002). A comparative cross-species study in hydroponics showed that rice is two-to sixfold more resistant than maize, wheat, and sorghum (Famoso et al., 2010). This high level of resistance is likely achieved by the pyramiding of multiple mechanisms conferred by multiple genes, a hypothesis supported by results from both genome-wide association (GWAS) and QTL studies (Famoso et al., 2011). Mapping of Al resistance QTL in a rice recombinant inbred line (RIL) population derived from a cross between Al-resistant Azucena (O. sativa L., tropical japonica) and Al-sensitive IR64 (O. sativa L., indica) cultivars identified four genomic regions associated with Al resistance, on chromosomes 1, 2, 9, and 12 (Famoso et al., 2011;Spindel et al., 2013). The QTL Alt12.1 on chromosome 12, which explains a large proportion (>19%) of the variation in Al resistance, encompasses a genomic region that includes the gene ALUMINUM RESISTANCE TRANSCRIP-TION FACTOR 1 (ART1), encoding a C2H2-type zinc finger transcription factor . The art1 mutant, producing a truncated version of the ART1 protein, is sensitive to Al stress.
A total of 32 genes were up-regulated in response to Al in the wild type but were not up-regulated in the art1 mutant . Some of these genes have been shown to play a role in rice Al resistance. These "ART1-regulated genes" have received this denomination because they are mis-regulated in the art1 mutant; in other words, their expression is up-regulated by Al stress in the wild type, but not in the mutant . It is worth noting that direct binding of the ART1 protein to upstream regulatory regions has so far been experimentally demonstrated only for STAR1 (Tsutsui, Yamaji, & Ma, 2011) and OsFRDL4 (Yokosho, Yamaji, Kashino-Fujii, & Ma, 2016).
The genetic architecture of rice Al resistance is complex and quantitative in nature (Famoso et al., 2011). In the present study, we describe how natural variation in the ART1 locus affects transcriptional responses to Al stress in two distinct genetic backgrounds of rice. We fine-mapped the Al resistance QTL Alt12.1 to a small region surrounding ART1 and used reciprocal near-isogenic lines (NILs) to examine the effect of two different ART1 alleles in japonica (Al-resis-  2.2 | Fine-mapping of the Alt12.1 QTL The Alt12.1 locus was fine-mapped using a substitution mapping approach (Appendix S2) to a 0.2 cM region (44.7 kb; chromosome 12: 3,578,363 -3,623,299 bp), between markers K_3.57 and in-del_3.62 (Figure 2a, Fig. S1, Table S1). According to the rice MSUv.7 genome assembly (Kawahara et al., 2013), this region encompasses six gene models: ART1 (LOC_Os12g07280), LOC_Os12g07290, LOC_Os12g07300, LOC_Os12g07310, LOC_Os12g07340, and LOC_Os12g07350.
To assess the levels of structural variation between japonica and indica across the 44.7-kb fine-mapped Alt12.1 locus, we aligned sequences from the japonica rice reference genome, Nipponbare (which served as a proxy for Azucena) to a long-read, de novo assembly of the IR64 genome (Schatz et al., in preparation; http://sc hatzlab.cshl.edu/data/ir64/). A single IR64 scaffold was identified that covered >36 kb of the Alt12.1 fine-mapped locus (44.7 kb). No major rearrangements were detected in the region, and both gene order and gene content were preserved (Fig. S2). Based on previous reports showing the major role played by ART1 in rice Al resistance, and the extreme Al-sensitive phenotype of the art1 mutant , we considered ART1 to be the primary candidate gene underlying the Al resistance QTL Alt12.1.

| Allelic variation in ART1
To analyze the extent of allelic variation in the ART1 gene, we examined both sequence and expression level polymorphisms.  (Garris, Tai, Coburn, Kresovich, & McCouch, 2005;Zhao et al., 2011). Azucena and Nipponbare, both Al-resistant japonica cultivars (Famoso et al., 2011), differ from each other by a single, nonsynonymous amino acid substitution at position 436 of the ART1 protein   2.6 | The effect of the reciprocal Alt12.1 QTL

introgressions on gene expression
We first examined the RNA-seq results for each NIL alongside its recurrent parent grown under control conditions (absence of Al).
Of these, 40 are located within the introgression from IR64 in Plasmids carrying the Azucena, Nipponbare, IR64, Kasalath, or art1 mutant allele of ART1 were individually infiltrated into tobacco leaves. The top row shows the YFP fluorescence, while the middle row shows the superimposed chloroplast fluorescence and bright field channels. The third row shows merged images from all three channels (YFP fluorescence, chloroplast fluorescence, and bright field). Scale bar = 20 lm. (b) Schematic diagram of the constructs used in yeast one-hybrid assays. For the reporter construct, a fragment of the STAR1 promoter was cloned in front of the lacZ gene and integrated into the yeast genome. Individual yeast colonies carrying the reporter construct were selected and transformed with one of five different transcription factor constructs. Each transcription factor construct contained the constitutive ADH1 promoter driving the expression of the GAL4 activation domain (GAL4 AD) fused to the ART1 CDS from Azucena (tropical japonica; Al-resistant), Nipponbare (temperate japonica; Al-resistant), IR64 (indica; Al-sensitive), Kasalath (aus Al-sensitive), and the art1 mutant allele. (c) Yeast one-hybrid assay comparing promoter activation by different ART1 alleles. A 164-bp fragment of the OsSTAR1 promoter was used. Constructs carrying an empty vector (control), or the Azucena, Nipponbare, IR64, Kasalath, or art1 mutant allele were assayed for beta-galactosidase reporter gene activity (n = 10). Error bars indicate standard deviation. All pairs of means in (c) were compared using a Tukey-Kramer HSD test; lines not connected by the same letter are significantly different (p < .05) AZU [IR6412.1] , and are interpreted to represent genotypic differences between indica and japonica. Three of them were also found to be differentially regulated by Al in Azucena: two were down-and one was up-regulated (File S2   In previous mutant studies, 32 genes were reported to be up-regulated by Al stress in wild type but not in the art1 mutant Yokosho, Yamaji, & Ma, 2011); these 32 genes have been designated ART1-regulated genes. Several of these were also identified in our study. The parental lines Azucena (Al-resistant) and IR64 (Al-sensitive) up-regulated 17 and 10 "ART1-regulated" genes, The genes that respond more strongly to the ART1 allele from than in Azucena (Fig. S7, File S4). Variation in OsFRDL4 expression was previously associated with the presence of a 1.2-kb insertion in the promoter (Yokosho et al., 2016). This insertion is present in Azucena and absent in IR64 (Fig. S7), suggesting that, in the genetic backgrounds used our study, variation in both the ART1 protein and the cis-regulatory region of OsFRDL4 may contribute to the regulation of OsFRDL4 expression in response to Al (Fig. S7).

| ART1 cis element analysis
The cis element binding site of ART1 has been experimentally defined as the short, degenerate sequence motif GGN(T/g/a/C)V(C/ A/g)S(C/G). This motif was identified based on protein gel-shift assays using the promoter of STAR1 (Tsutsui et al., 2011). We set out to determine whether the putative regulatory regions of the Al-  Because ART1 was shown to bind with lesser affinity when certain bases were found in positions 3 and 4 of the motif (Tsutsui et al., 2011), we also mapped the less degenerate motif GGYMS, in which the bases predicted to have less affinity to ART1 were excluded.
The motif GGYMS was found in 55,460 of 55,554 putative regulatory sequences, with an average of 10 hits per sequence (Fig. S8b).
We determined that enrichment of the 5-bp ART1-binding cis element in the regulatory region of putative Al response genes is insufficient criteria for identifying ART1-regulated genes, due to the widespread presence of the motif in the background of putative regulatory regions.  (Figure 7a). In other words, these two genes were upregulated only in the IR64 genetic background and not in the Azucena genetic background, and the up-regulation is triggered by either of the ART1 alleles. We hypothesized that IR64, the Al-sensitive parent, carries positive alleles conferring enhanced Al resistance at both of these loci (and possibly others). In line with this hypothesis, we previously reported transgressive segregation for rice Al resistance (Famoso et al., 2011).
To further explore this hypothesis, we investigated the expres-  The RNA-seq results presented here indicate that the ART1 allele from Azucena leads to a stronger transcriptional response to Al stress in rice roots. Using yeast one-hybrid assays, we demonstrate that the ART1 proteins encoded by the alleles found in indica/aus (IR64; Kasalath) and japonica (Azucena; Nipponbare) differ significantly in their DNA-binding ability. In the yeast one-hybrid assay, the ART1 proteins encoded by the IR64 and Kasalath alleles activated the reporter gene more strongly than the ones from Azucena and Nipponbare (Figure 3c). This is a nonintuitive result, given that the Azucena allele drives gene expression more strongly than the IR64 allele in planta. Yeast one-hybrid is a heterologous system, in which the TF of interest is fused with the strong GAL4 activation domain and the DNA sequence of interest is placed upstream of a reporter gene. As such, other elements that may affect TF activity in planta are missing. Further studies using promoter-reporter systems in plant protoplasts will shed more light into the transcriptional activation ability of the different ART1 proteins.

| Using reciprocal NILs to study transcriptomic responses to Al in IR64 and Azucena
We generated reciprocal NILs in which the Alt12.1 QTL region from Azucena was introgressed into IR64, and vice versa (Figure 1) in each of the NILs, suggesting that these genes represent "background" genotypic differences between indica and japonica. However, five genes showed differential expression in response to Al treatment (File S2). This included LOC_Os12g07340, which is one of six genes located within the fine-mapped region, only 28 kb away from ART1 (Figure 2a). The function of these genes is unknown, and although none are predicted to encode a regulatory protein or TF, we cannot rule out the possibility that one or more may contribute to Al resistance. Our study illustrates some of the complexities and limitations of using NILs to molecularly ARBELAEZ ET AL.  that the Azucena ART1 allele is more effective in conferring Al resistance than the IR64 allele in both genetic backgrounds. The small sample size (n = 4 genotypes) limits our ability to statistically analyze the differences between genetic backgrounds; nevertheless, our data suggest that the identity of the genes affected by the two ART1 alleles differed in Azucena and IR64.
3.5 | ART1 allele swapping between Azucena and IR64 affects the expression pattern of many genes not previously designated as ART1-regulated Among the genes up-regulated by Al stress in our RNA-seq study are many previously characterized Al resistance genes, including OsNrat1, OsFRDL4, OsFRDL2, OsALS1, STAR1, STAR2, and OsCDT3.
OsNrat1 encodes a plasma membrane-localized Al uptake transporter (Li et al., 2014;Xia, Yamaji, Kasai, & Ma, 2010). OsNrat1 function is likely coupled with a mechanism of internal detoxification, involving another ART1-regulated gene: OsALS1. This gene encodes a half-size ABC transporter localized to the tonoplast of root cells, where it is thought to sequester Al 3+ into the vacuole (Huang, Yamaji, Chen, & Ma, 2012). STAR1 and STAR2 were also shown to play a role in rice Al resistance: STAR1 encodes the nucleotide-binding domain of a bacterial-type ABC transporter that interacts with the transmembrane domain encoded by STAR2 ). Disruption of either gene results in higher Al sensitivity; however, the mechanism by which the STAR1-STAR2 complex confers Al resistance is unknown. When expressed in Xenopus oocytes, STAR1-STAR2 facilitates the export of UDP-glucose; this substrate is proposed to modify the cell wall in a way that reduces Al 3+ accumulation. OsCDT3 encodes a small cysteine-rich peptide that exhibits Al-binding properties (Xia, Yamaji, & Ma, 2013).
The genes that respond more strongly to the ART1 allele from Azucena than to the ART1 allele from IR64 are presumably the drivers of the phenotypic difference in Al resistance observed between NILs and their corresponding recurrent parents. Among the genes previously designated as ART1-regulated, only a few were among those that responded more strongly to the ART1 allele from Azucena: OsFRDL4, Al-activated release of organic acids from the root is a major physiological mechanism of plant Al resistance (Kochian et al., 2015). OsFRDL4 is the major candidate for the gene underlying an Al resistance QTL in a Nipponbare (temperate japonica) 9 Kasalath (aus) population (Yokosho et al., 2016). The authors of that study also reported a 1.2-kb insertion in the promoter of OsFRDL4, and demonstrated that a majority of japonica varieties tested carried the insertion, while most indica varieties did not, and suggested that the insertion was associated with higher levels of OsFRDL4 expression under Al stress. Yokosho et al. (2016)

| ART1 cis element
The binding site of ART1 was experimentally defined as the short, degenerate sequence motif GGNVS, found in the promoter region of genes described as ART1-regulated. This motif was further validated for OsSTAR1, a known target of ART1, using EMSA (Tsutsui et al., 2011;Yokosho et al., 2011 TF network analysis using more in-depth transcriptome data, will likely improve our functional understanding of ART1 and help to better refine its cis-regulatory element.
3.8 | ART1-regulated genes and transgressive variation for rice Al resistance 3.9 | ART1 function beyond Al resistance?
While a majority of genes differentially expressed between each NIL and its parent under control conditions was localized to the chromosomal introgression, a small number of them were distributed on other chromosomes. Therefore, the possibility that ART1 also regulates gene expression in the absence of Al cannot be discarded. It is important to note that expression of ART1 itself is constitutive and is not up-regulated by Al (Figure 2c); the signaling mechanism that activates ART1 binding to the promoter of Alresponsive genes in response to Al stress has not been elucidated.
In addition, it is also possible that ART1 could have additional func-

| Phenotyping for Al resistance
Al resistance was phenotyped as described by Famoso et al. (2010).
Individual root seedling digital images were obtained to quantify total root growth (TRG) values according to Clark et al. (2013) using the software RootReader2D (www.plantmineralnutrition.net/rr2d. php). Total root growth (TRG) values from each genotype grown under control and Al stress conditions were used to estimate relative root growth (RRG) indices as described (Famoso et al., 2010). RRG values were used for further statistical analysis.

| Genotyping
For PCR analysis total DNA was extracted from fresh leaf tissue  Imai et al. (2013). Samples were genotyped using the 6K Infinium array (Illumina, www.illumina.com/). DNA was extracted using a modified CTAB protocol described by Imai et al. (2013). The 6K Infinium arrays were performed as described by Arbelaez et al. (2015).

| Quantitative real-time PCR (RT-qPCR)
Total RNA was isolated using the QIAGEN RNeasy Mini Kit according to the manufacturer's instructions (https://www.qiagen.com). Incolumn digestion of genomic DNA was performed using the QIA- replicates were averaged per sample per each assay. RNA from six independent biological replicates was used to measure ART1 expression and confirm that there are no significant differences in expression between genotypes and/or treatments. One biological replicate was used to confirm LOC_Os01g53090 and LOC_Os04g41750 in NILs and parents, and to determine expression levels in the RILs.

| Subcellular localization of ART1:YFP protein fusions
In order to determine the cellular localization of the aberrant protein generated by the art1 mutant, we re-created the mutation in vitro. This clone was also used to generate the yeast one-hybrid construct. To generate the mutant protein, we used the reference genome sequence to identify the new stop codon in the frameshifted protein, which extended into the 3 0 UTR of the ART1 locus.
We then used a series of long oligos to add the required nucleotides to the previously cloned Nipponbare allele, followed by USER-based cloning to introduce the 1-bp deletion (http://www.cb s.dtu.dk/services/PHUSER/). ART1-YFP constructs were made in pGREEN under the control of a 35S promoter. Constructs were infiltrated into 3-to 5-week-old N. benthamiana. After allowing 2 to 4 days for the protein to express, leaf disks were cut from tissue adjacent to the infiltration site and used for confocal microscopy. Confocal images were collected on a Leica TCS-SP5 confocal microscope (Leica Microsystems, Exton, PA) using a 639 water immersion objective. The YFP fusion protein was visualized by excitation with an argon ion laser (512 nm), and emitted light was collected between 525 and 575 nm. Chloroplasts were excited with the blue argon laser (488 nm), and emitted light was collected from 680 to 700 nm. The YFP and chloroplast fluorescence were collected on separate channels and superimposed with a bright field image collected simultaneously with the fluorescence images.
Images were processed using Leica LAS-AF software (version 2.6.0).

| Yeast one-hybrid experiments
A fragment of the STAR1 promoter along with the five ART1 alleles was cloned into pENTR D-TOPO (ThermoFisher, www.thermofishe r.com). The STAR1 promoter fragment (164 bp) was selected according to Yamaji et al. (2009)  Bait strains were generated by homologous recombination of pLacZi-STAR1p into the yeast strain YM4271 according to the manufacturer's protocol (Clontech, www.clontech.com). Integration was verified by PCR. The prey vectors were generated by LR recombination of the pENTR-ART1 constructs into the yeast onehybrid compatible plasmid pDEST22, which fuses the GAL4 activation domain to the N-terminus of the transcription factor. The bait strain was then transformed using the lithium acetate method (Gietz & Schiestl, 2007) with each of the five different pDEST22-ART1 constructs, and assayed for beta-galactosidase reporter gene activity using the liquid ONPG method described in the Clontech yeast protocol handbook (Clontech). Ten independent yeast colonies were screened for each allele. Student t tests assuming unequal variance (and a = 0.05) were performed to determine significant differences in DNA-binding affinity.

| RNA-seq experiment and analysis
Plants were germinated in the dark in moist paper rolls for 3-4 days, and then, 40 uniform seedlings per genotype were grown in control nutrient solution for five days. Subsequently, 20 seedlings per genotype were transferred into a hydroponic solution with 80 lM of Al 3+ activity for 4 hr before total root tissue was harvested for RNA extraction. Four independent biological replicates were performed, for a total of 32 samples (4 genotypes 9 2 treatments 9 4 biological replicates). Total RNA was isolated from root tissue using the QIA- and sequenced in two lanes of an Illumina Hiseq2000 using the single-end reads mode (Pombo et al., 2014). Raw RNA-seq reads were processed to remove adaptor and low-quality sequences using Trimmomatic (Bolger, Lohse, & Usadel, 2014). RNA-seq reads longer than 40 bp were kept. Reads were aligned to a ribosomal RNA database (Quast et al., 2013) using Bowtie (Langmead, Trapnell, Pop, & Salzberg, 2009) and matches were discarded. The resulting, high-quality reads were aligned to the MSUv7.0 rice genome annotation (Kawahara et al., 2013) using HISAT (Kim, Langmead, & Salzberg, 2015).
Following alignments, raw counts for each rice gene were derived and normalized to reads per kilobase of exon model per million mapped reads (RPKM).

| Differential expression analysis
Differentially expressed (DE) genes were identified using the DESeq2 package (Love, Huber, & Anders, 2014), https://biocond uctor.org/packages/release/bioc/html/DESeq2.html). Independent filtering to remove low expression genes previous to DESeq2 analysis was performed by eliminating genes that across all samples their total sum counts is below the threshold estimated as: overall normalized mean counts/number of samples to be analyzed (Bourgon, Gentleman, & Huber, 2010).
Differentially expressed genes (DEGs) were selected according to the following post hoc criteria: (i) twofold or greater change in transcript abundance, up or down, measured in log2 fold-change ratio (log 2 FC ≥1, or ≤À1); (ii) p value of ≤.05, corrected for multiple testing using a false discovery rate (FDR) criteria (Benjamini & Hochberg, 1995); and (iii) minimum of eight normalized, log-transformed counts in at least one sample (Worley et al., 2016). Regularized log-transformed counts (Love et al., 2014) estimated from DESeq2 using the true experimental design option (blind = TRUE; which accounts for different library sizes and stabilizes the variance among counts) were used to generate heat maps and for principal component analysis (PCA). For PCA, the variance across all samples for normalized, log-transformed counts was estimated using the "rowVars" function in R (https://cran.r-project.org/). Genes were ranked from the most to the least variance and the 500 genes with the highest variance across all samples were used to calculate prin-