Evolutionary origins and functional diversification of Auxin Response Factors

Abstract

The Auxin Response Factors (ARFs) family of transcription factors are the central mediators of auxin-triggered transcriptional regulation. Functionally different classes of extant ARFs operate as antagonistic auxin-dependent and -independent regulators. While part of the evolutionary trajectory to the present auxin response functions has been reconstructed, it is unclear how ARFs emerged, and how early diversification led to functionally different proteins. Here, we use in silico and in vivo analyses to revisit the molecular events that led to the origin and subsequent evolution of the ARFs. We reveal the shared origin of ARFs from preexisting domains, uncovering a protein fold homologous to the ARF DNA-binding fold in a conserved eukaryotic chromatin regulator. Building on this, we reconstruct the complete evolutionary history of ARFs, including the divergence events leading to the appearance of the ARF classes and defining the main molecular targets for their functional diversification. We derive a complete evolutionary trajectory that led to the emergence of the nuclear auxin signalling pathway.

Fig. 1. The ARF DNA-binding scaffold evolved from an ancestral Tudor-like domain.

a PHMMER hits using DD-AD region on the UniProt Reference Proteomes database. Summary of architectures found merged into either ARF or PHIP are shown. b PHIP and ARF architecture comparison and multiple protein sequence alignment of the double Tudor domain regions based on Marchantia polymorpha protein primary sequences of MpPHIP and MpARF1. c Unrooted maximum likelihood gene tree of the Tudor-like regions of ARFs and PHIPs. Bootstrap values are indicated as color-coded bubbles in branch nodes. Scale bar indicates substitutions per residue. d Structural alignment of the DD-AD of PHIP and ARF, represented by Marchantia polymorpha PHIP and ARF1 predicted structure. Blue, MpARF1; Grey, MpPHIP. Key alpha helixes are pointed and DD, B3 and AD domains marked. Red circles represent the position of the hydrophobic cage residues in AD. e Magnification of the hydrophobic cage structure characteristic of Tudor domains in MpPHIP and MpARF1 with key amino acids highlighted. f Arabidopsis thaliana wild-type (WT) and mp/arf5 mutant seedlings used to score complementation in i. Scale bar, 2 mm. g Raincloud plot of half-thallus area measurements of 10-day-old Mparf1 mutant complemented lines with AD-deleted versions grown in mock (DMSO) or auxin (3 μM NAA). n = 19,19,17,20,17,20 (left to right, mock/treatment). Boxplots in Raincloud indicate the following parameters: centrum, median; upper bound, first quartile; lower bound, third quartile; whiskers maximum and minimum refer to highest and lowest values, respectively, within 1.5*inter-quartile range (IQR). Statistical groups are determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA. h Representative pictures of M. polymorpha plants in (g), Scale bar, 5 mm. i Summary of AtARF5 (MONOPTEROS, MP) and MpARF1 complementation of each corresponding mutant per transgenic line with AD deleted versions. j Multiple protein sequence alignment of the AD region on selected PHIP and ARF proteins showing hydrophobic cage residues, highlighting mutations for analyses in k referring to AtARF5 and MpARF1 residue positions. k AtARF5 and MpARF1 complementation of each corresponding mutant per transgenic line with hydrophobic cage point mutants. Source data are provided as a Source Data file.

Fig. 2. RAVs and ARFs are ancestral streptophyte-specific ortholog genes.

a Schematic representation of RAV and ARF proteins primary sequence based on Marchantia polymorpha protein primary sequences of RAV and ARF1. BRD stands for B3 repression domain. b Maximum likelihood square phylogenetic tree of B3 domains in Viridiplantae rooted using REM as outgroup. Bootstrap values are indicated as color-coded bubbles in branch nodes. Scale bar represents distance in substitutions per residue. c Raincloud plot of half-thallus area measurements of 10-day-old Mparf1 mutant complemented a line with MpARF1 B3 domain swapped for MpRAV B3 grown in mock (DMSO) or auxin (3 μM NAA). n = 16,15,24,28,14,19,16,15 (left to right, mock/treatment). Statistical groups are determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA. d Representative pictures of plants in c), Scale bar, 5 mm. e Unrooted maximum likelihood gene tree of the eukaryotic PB1s re-analysed from Mutte & Weijers 2020. Bootstrap values are indicated in key nodes as color-coded bubbles in branch nodes. Scale bar represents distance in substitutions per residue. f Maximum likelihood square phylogenetic tree of ARF/RAV-related PB1 domains rooted using Chlorophyta Pho1-like PB1 sequences as outgroup. Bootstrap values are indicated as color-coded bubbles in branch nodes. Scale bar represents distance in substitutions per residue. g Summary of Marchantia polymorpha PB1 pairwise interaction assays in yeast-two-hybrid combining drop and galactosidase activity assays (see Supplementary Fig. 6, and Supplementary Data 2). h Raincloud plot of half-thallus area measurements of 10-day-old Mparf1 mutant complemented lines with MpARF1 PB1 domain swapped for MpRAV and MpIAA PB1s grown in mock (DMSO) or auxin (3 μM NAA). n = 18,17,27,24,23,28,11,15,24,20 (left to right, mock/treatment). Boxplots in Raincloud indicate the following parameters: centrum, median; upper bound, first quartile; lower bound, third quartile; whiskers maximum and minimum refer to highest and lowest values, respectively, within 1.5*inter-quartile range (IQR). Statistical groups are determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA. Source data are provided as a Source Data file.

Fig. 3. PHIP and RAV are not related to ARF function in M. polymorpha.

a Pictures of 10-day-old wild-type (Tak-1) and genomic edited phip plants grown in mock (DMSO) or auxin (3 μM NAA). Scale bar, 5 mm. b Raincloud plot of thallus area measurements (in halves) of 10-day-old wild-type and phip mutant plants grown in (a). n = 12,14,11,11 (left to right, mock/treatment). Statistical groups are determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA. c Quantification of gametangiophore formation after sexual organ induction (continuous far-red light, cFR) in 10-day-old plants. Data represents mean and error bars each day after transfer to cFR, and standard deviation of five plants. Asterisks indicate statistical difference compared to the wild-type group (Tak-1) by Kruskal–Wallis test (*p < 0.05; **p < 0.01; ***p < 0.001). d, g Expression analysis of MpWIP by qRT-PCR in 10-day-old plants treated for 1 h with 3 μM NAA or DMSO (Mock), using MpSAND and MpEF1α as reference genes. Dots represent average of two technical replicates of biological replicates (n = 4), but for Mparf1 (n = 3); bars represent the average of the biological replicates. Statistical groups are determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA. e Pictures of 10-day-old wild-type (Tak-1) and genomic edited Mparf1;rav plants grown in mock (DMSO) or auxin (3 μM NAA). Scale bar, 5 mm. f Raincloud plot of thallus area measurements (in halves) of 10-day-old wild-type and rav mutant plants grown in (e). n = 14. Letters (black, mock; blue, NAA) indicate statistical groups as determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA. Source data are provided as a Source Data file. Boxplots in Raincloud indicate the following parameters: centrum, median; upper bound, first quartile; lower bound, third quartile; whiskers maximum and minimum refer to highest and lowest values, respectively, within 1.5*inter-quartile range (IQR).

Fig. 4. ARF originated and diverged from a single ABC-class.

a Phylogenetic tree of ARF proteins using DD-to-AD protein sequences, rooted using PHIP sequences. Bootstrap values are indicated as color-coded bubbles in branch nodes, with the main ARF class ancestral nodes highlighted. Scale bar represents distance in substitutions per residue. b Reconstruction of the evolutionary pathway of ARF and related proteins in streptophytes as a schematic summary of presence/absence extracted from data in phylogenetic analyses. Right tree depicts the known phylogenomic tree of plant lineages with ARF major origin and duplication events marked. Expansion indicates lineage-specific family expansions inferred from phylogenetic data as a qualitative approximation. Source data are provided as a Source Data file.

Fig. 5. A-class ARFs evolved through fast neofunctionalization.

a Pictures of 10-day-old plants grown in mock (DMSO) or auxin (3 μM NAA) in a Mparf1 complementation assay. Full length CDS of referred ARFs are expressed under the endogenous MpARF1 promoter in a Mparf1 background. Scale bar, 5 mm. b Heatmap representation of Mparf1 thallus phenotype complementation with full length ARFs of 10-day-old plants grown in mock (DMSO) or auxin (3 μM NAA) extracted from thallus area measurements. Asterisks indicate statistical differences with Mparf1 as determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA (see Supplementary Fig. 9). Asterisks indicate statistical difference compared to the mutant in mock conditions. Upper tree indicates phylogenetic representations among the ARFs used in the complementation assay. c, Yeast transactivation assays of full length ARFs fused to the yeast Gal4 DNA binding domain showing the quantification of the UAS:LacZ reporter activation as colorimetric-measured β-galactosidase activity. Dots represent the average of three semi-technical replicates in independent biological replicates (independently transformed lines, n = 4). d Dual luciferase transactivation assay in Arabidopsis protoplasts using a LUC gene under the control of 5xGal UAS motif as reporter, 35S:REN as ratiometric control, and the yeast Gal4 DNA binding domain fused to ARFs as effectors, fused additionally to mNeonGreen for protein expression normalization. e Qualitative summary of interspecies PB1 pairwise interaction assays obtained from yeast-two-hybrid quantitative β-galactosidase activity assays (extracted from Supplementary Data 2). In (c, d), statistical groups are determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA. Cm Chlorokybus atmophyticus, Me Mesotaenium endlicherianum, Sp Spirogyra pratensis, Sm Spirogloea muscicola, Pm Penium margaritaceum, Co Coleochaete orbicularis, Kn Klebsormidium nitens, Aa Anthoceros agrestis, Tl Takakia lepidozioides. Source data are provided as a Source Data file.

Fig. 6. C-class ARFs function is specific and deeply conserved.

a Pictures of Mparf3 mutant complemented with different ARF classes. MpARF1, MpARF2, and MpARF3 coding sequences are expressed under the control of the endogenous MpARF3 promoter in the Mparf3 mutant. Upper row represents 20-day-old apical notches right after excision from adult plants. Lower row are the same plants after seven days of re-growth. White-shaded area in lower row is the same area occupied at day 0. Red dotted line indicates excision, while filled arrowheads point to apical notches present at excision. Scale bar, 5 mm. b Dot plot of notch-driven thallus growth measured as projected area fold-change 14 days after excision (percentage of times day 0 area). ARF coding sequences are expressed under the control of the endogenous MpARF3 promoter in the Mparf3 mutant. Mparf3 plants expressing CmARF do not produce visible notches, preventing comparable excision. Dots represent thallus area growth percentage average of three notches derived from a single plant; n = 5,5,6,8,7,7,6,8,8 (left to right). Statistical groups are determined by non-pooled Welch’s t-test and Benjamini–Hochberg adjustment (adjusted p < 0.01 for non-overlapping letters). c Pictures of regenerated plants showing the phenotype of a MeARFc-complemented Mparf3. White dotted squares in first and third images indicate the zoomed are shown in the second and fourth pictures, respectively, highlighting a mature gemma cup. Cm Chlorokybus melkonianii, Me Mesotaenium endlicherianum, Sp Spirogyra pratensis, Sm Spirogloea muscicola. Source data are provided as a Source Data file.

Fig. 7. ARF classes diverged through sequential domain specialization.

a Pictures of 10-day-old plants grown in mock (DMSO) or auxin (3 μM NAA) in a Mparf1 complementation assay. Chimeric MpARF1 with heterologous ARF^DBD sequences as indicated are expressed under the endogenous MpARF1 promoter in a Mparf1 background. Scale bar, 5 mm. b Heatmap of Mparf1 thallus phenotype complementation with MpARF1-ARF^DBD chimeras of 10-day-old plants grown in mock (DMSO) or auxin (3 μM NAA) extracted from thallus area measurements. Asterisks indicate statistical differences with Mparf1 as determined by Tukey’s Post-Hoc test (p < 0.05) following one-way ANOVA (see Supplementary Fig. 12). c Electrophoretic mobility assay showing MBP-SmARF3 DNA-binding domain interaction with a bipartite ARF-specific binding site (IR7, inverted repeat, 7 bp spacing between two Auxin Response Elements, AuxRE), with one or two mutated AuxRE (IR7^m and IR7^2m, respectively). Marchantia B-class ARF (MpARF2) DBD is shown as positive control. Right-most lanes show a MBP-SmARF3 protein titration experiment indicating concentration-dependent complex-formation ([SmARF3^DBD]=100, 50, 25, 5 μM, left to right). This assay has been performed twice with equivalent results. d Summary MpARF3 domain swaps with MpARF1 and MpARF2 in Mparf3-related phenotypes complementation extracted from experiments shown in Supplementary Fig. 13. Full complementation (Yes) is assigned when phenotype is comparable to that of the wild-type (Tak-1) instead of Mparf3. Partial complementation indicates a quasi-wild-type phenotype. e, f Dot plots of notch-driven thallus growth measured as projected area doubling 14 days after excision (percentage of times day 0 area). Domain-swapped MpARF3 chimeras are expressed under the control of the endogenous MpARF3 promoter in the Mparf3 mutant. Dots represent thallus area growth percentage average of three notches derived from a single plant. e MpARF3-ARF^DBD chimeras; n = 5. f MpARF3-MeARFc^MR-PB1 chimera; n = 4. Statistical groups are determined by non-pooled Welch’s t-test and Benjamini–Hochberg adjustment (adjusted p < 0.01 for non-overlapping letters). Sm Spirogloea muscicola, Pm Penium margaritaceum. Source data are provided as a Source Data file.

Fig. 8. Origin and evolutionary history of ARF proteins.

a Homologous domain origin leading to ARF proteins. Cartoon adaptation from known structures of domains, connected by an artistic representation of regions without a known structure. ARF homologous domains are highlighted in the related protein cartoons (Dark rosy-brown shade, cTudor from PHIP; Tan-brown, PB1, and light rosy-brown, B3 domains from RAV). b Hypothetical step-wise domain acquisition originating ARF and RAV domain architectures from an ancestral ARF/RAV gene. (I) and (II) assume a RAV-like ancestral RAV/ARF (AP2 + B3 containing). (I) Duplication of ancestral copy giving rise to ARF and RAV lineages; ARF-independent acquisition of PHIP-cTudor scaffold. (II) Acquisition of cTudor in ancARF/RAV and subsequent duplication and loss of this domain in the RAV lineage. (III) Roadmaps from a non-AP2 containing ancRAV/ARF; subsequently AP2 domain acquired or lost at different points. BRD stands for B3 repression domain. c ARF evolution from the ancestral ARF. Major divergence points associated to molecular sub- and neo-functionalization events are indicated.