Hidden Markov model (HMM) logos of Aux/IAA (PF02309) and ARF (PF06507) proteins were downloaded from the Pfam database (Finn et al., 2014) and used to scan the E. konishii predicted proteome (Sun et al., 2021) with the HMMER software package (Finn et al., 2011). The resulting Aux/IAA and ARF candidates were further used to generate HMM logos for EkIAAs and EkARFs using hmm-build from the HMMER suite (Finn et al., 2011), before scanning the E. konishii proteome again. Proteins with an E-value lower than 0.01 were retained, and the presence of conserved ARF or IAA domains was confirmed using the Conserved Domains Database (Marchler-Bauer et al., 2011), Pfam (Finn et al., 2014), and the Simple Modular Architecture Research Tool (Letunic and Bork, 2018). The proteins meeting all of the above criteria were used for further study. The number of amino acids, the predicted molecular weight, and the theoretical isoelectric point (pI) were determined using the ExPASy server (http://web.expasy.org/protparam/) (Gasteiger et al., 2003).

TBtools (Chen et al., 2020) was employed to illustrate the exon/intron structures of all EkIAA and EkARF genes. Conserved protein motifs in their encoded proteins were predicted by the MEME program (parameters: number of maximum patterns, 10; maximum width, 50) (http://memesuite.org/tools/meme) (Bailey et al., 2006).

The protein sequences for the 29 ARFs and 34 IAAs from Arabidopsis were obtained from published references (https://www.arabidopsis.org), and the protein sequences of 14 ARFs and 11 IAAs from Amborella trichopoda were identified following the same method described for E. konishii. Full-length protein sequences for all IAAs and ARFs identified in E. konishii, Arabidopsis, and A. trichopoda were used for phylogenetic analysis. The phylogenetic tree was built with the maximum likelihood method on the IQ-TREE web server (http://iqtree.cibiv.univie.ac.at/) (Nguyen et al., 2015).

The upstream sequences (2 kb) of EkIAA and EkARF genes were extracted via TBtools (Chen et al., 2020) and then submitted to the PlantCARE database (Lescot et al., 2002) (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to identify cis-elements.

The chromosomal distribution and location of all EkIAA and EkARF genes were acquired from the E. konishii genome annotation file. Colinear circles for EkIAA and EkARF genes were drawn with TBtools (Chen et al., 2020). Duplication events were confirmed on the basis of coverage (>70% of the entire gene body) and similarity (70%) of the two aligned sequences (Gu et al., 2002) and were considered tandem duplication pairs if they were located within 100 kb (Mehan et al., 2004). Genes located in duplicated regions with 70% similarity were identified as segmental duplications (Mehan et al., 2004). K _a/K _s values were calculated with TBtools (Chen et al., 2020).

The RNA-seq data for three development stages (green, turning, and red fruit) and three tissues (red-winged pericarp, branch, and leaf) were downloaded from the National Center for Biotechnological Information (NCBI) Sequence Read Archive (SRA) under the accession numbers PRJNA548305 and PRJNA548305, respectively. The RNA-seq reads were mapped to the E. konishii reference genome via Salmon algorithm, and the transcripts per million reads (TPM) for AUX/IAA and ARF genes were extracted for further analysis. The heatmaps were drawn using TBtools (Chen et al., 2020).

To study the effect of auxin on the regulation of anthocyanin and terpenoid biosynthesis in E. konishii, a comprehensive correlation analysis was first performed using the correlation test in R between anthocyanin contents and the expression levels of anthocyanin biosynthetic genes. Biosynthetic genes whose expression was positively correlated with anthocyanin contents were selected for further correlation analysis between their expression levels and those of EkARF genes. Pearson’s correlation coefficients (r, p-value < 0.05) were used to define five correlation levels: no correlation (|r| ≤ 0.2), weak correlation (0.21 ≤ |r| ≤ 0.35), moderate correlation (0.36 ≤ |r| ≤ 0.67), strong correlation (0.68 ≤ |r| ≤ 0.90), and very strong correlation (0.91 ≤ |r| ≤ 1), with r > 0 indicating positive correlations and r < 0 negative correlations (Prion and Haerling, 2014).

Fruits at the green stage, turning stage, and red stage from E. konishii were harvested as materials for RT-qPCR. Total RNA was extracted using the RNAprep Pure kit (Tiangen, China), and then 1 µg of total RNA per sample was subjected to reverse transcription using the PrimeScript RT Reagent Kit (Takara, China) with gDNA Eraser (Takara, China). The specific primers for EkARF5.1 (F: 5′-GCA​ACC​TCC​AAC​TCA​AGA​GC-3′, R: 5′-GAC​GCC​TCA​CAC​CCA​CTA​AT-3′) were designed by Primer 5 software and synthetized by Sangon Biotech (Shanghai, China). UBC23 (F: 5′-AGC​CAC​ATA​ATC​TCC​GTG​TAA​G-3′, R: 5′-GCT​GAC​CAT​GTT​CGA​GTA​GTT-3′) was used as an internal reference (Yuan et al., 2018b). The reaction mixture consisted of 10 µl of 2× GoTaq qPCR Master Mix (Promega, United States), 0.4 µl of each gene-specific primer, 1 µl of cDNA (10× dilution), 0.4 µl of dye, and 7.8 µl of nuclease-free water. The reaction conditions were as follows: 95°C for 2 min followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Relative gene expression levels were calculated by the comparative ΔCt method. Three biological replications were assessed per sample.

The full-length EkARF5.1 coding sequence was cloned into the pGBKT7 vector to generate BD-EkARF5.1, which was then introduced into yeast strain Y2HGold. The resulting colonies were grown on synthetic defined (SD) medium lacking Trp and His for 2 days to observe the transcriptional activation activity of EkARF5.1, using empty pGBKT7 as a negative control.

To test the binding of EkARF5.1 to AuxREs, seven repeats of the AuxRE element (TGTCTC) were inserted into the multiple cloning site of the pAbAi vector to generate the 7×TGTCTC-pAbAi vector, which was integrated into the Y1HGold genome to construct the bait reporter strain. The AD-EkARF5.1 clone was generated by subcloning the full-length EkARF5.1 coding sequence into the pGADT7 vector. AD-EkARF5.1 was then transformed into the bait reporter strain. The transformants were spotted onto SD medium lacking Leu or the same medium containing 100 ng/ml of the antibiotic aureobasidin A (AbA) and allowed to grow for 48 h to analyze binding activity.

Constructs consisting of the β-GLUCURONIDASE (GUS) reporter 7×TGTCTC:GUS (cloned in pMDC164) and 35S:EkARF5.1 (cloned in pDMC32) were introduced into Agrobacterium (Agrobacterium tumefaciens) strain GV3101 for infiltration in Nicotiana benthamiana leaf epidermal cells. Agrobacterium harboring 7×TGTCTC:GUS and 35S:EkARF5.1 were infiltrated into the abaxial side of N. benthamiana leaves with a syringe as described previously (Liu et al., 2014). N. benthamiana leaves were stained for GUS activity 3 days after infiltration.

Publicly available datasets were analyzed in this study. This data can be found here: The E. konishii chromosome-level genome assembly and annotation data (Accession No. PRJCA005268/GWHBCHS00000000) were available from National Genomics Data Center at BioProject/GWH (https://bigd.big.ac.cn/gwh c).

To identify EkIAA and EkARF genes, we searched the E. konishii genome using HMMER v3 in two rounds (see Materials and Methods for details). We analyzed the resulting protein sequences using the Conserved Domains Database, Simple Modular Architecture Research Tool, and Pfam database, which resulted in 34 Aux/IAA and 29 ARF candidate proteins (Figure 1; Supplementary Table S1). We numbered the E. konishii IAA and ARF genes based on their homologs in Arabidopsis (Remington et al., 2004; Okushima et al., 2005); the full list is provided in Supplementary Table S1, along with their gene IDs, their coding sequences, genomic DNA and protein sequences, the lengths of the coding and protein sequences, and the predicted molecular weights and isoelectric points (pI) of the proteins.

We observed a broad variation in the lengths and biochemical properties of EkIAA and EkARF proteins. EkIAA proteins ranged from 91 (EkIAA1.4) to 849 (EkIAA29.3) amino acids, with predicted molecular weights from 10.5 to 97.2 kDa (Supplementary Table S1). The predicted pI values of EkIAA proteins varied from 4.7 (EkIAA32/34) to 9.3 (EkIAA33.2). Similarly, EkARF proteins ranged in length from 432 (EkARF17.2) to 1,110 (EkARF24.1) amino acids with predicted molecular weights from 48.2 to 122.7 kDa (Supplementary Table S1). The predicted pI values of EkARF proteins varied from 5.41 (EkARF5.1) to 8.42 (EkARF10.2) (Supplementary Table S1).

To better understand their evolutionary history, we subjected all IAAs and ARFs identified in the model plant Arabidopsis, the early angiosperm A. trichopoda, and E. konishii to phylogenetic analysis with the MEGA-X software package (Kumar et al., 2018) and the IQ-TREE web server (Nguyen et al., 2015). In both protein families, individual members clustered into five branches, indicating that IAAs and ARFs are highly differentiated (Figure 2A), as previously reported in Arabidopsis (Remington et al., 2004) and poplar (Kalluri et al., 2007). EkIAA proteins were equally divided among groups Ⅰ, Ⅲ, and Ⅴ (Figure 2A). EkIAA5 and EkIAA15 from group I appeared to have undergone gene duplication, while group I had no clear E. konishii orthologs for Arabidopsis IAA6 or IAA19 (Figure 2A). As illustrated by the size of groups Ⅱ and Ⅳ, EkIAA genes have undergone gene duplication, especially IAA27 and IAA29 (Figure 2A). As with IAAs, the phylogenetic analysis of ARFs also divided the proteins into five groups (Figure 2B), as previously reported in Arabidopsis (Okushima et al., 2005). Group I consisted of six EkARF and five Arabidopsis ARF members, the latter having been reported to exhibit transcriptional activation activity (Guilfoyle and Hagen, 2007). Of note, ARF5, ARF6, and ARF8 all showed gene duplication in the E. konishii genome, while orthologs for ARF7 and ARF19 appeared to be lacking. Some, but not all, group II members showed signs of gene duplication (Figure 2B), for example, ARF1 and ARF9. Notably, IAAs and ARFs displayed the same distribution across groups in Arabidopsis and E. konishii, with the exception of group Ⅴ ARFs, which indicates that the IAA and ARF gene families in E. konishii are likely conserved. Within each group, several EkIAA and EkARF members had experienced duplication, with the exception of group Ⅴ IAA members (Figure 2A). In addition, we noted the absence of clear orthologs for several IAAs and ARFs in several groups (Figure 2), indicative of their independent evolution in E. konishii since the divergence from Arabidopsis.

We determined the genomic positions of all EkIAA and EkARF genes along the linkage groups (LGs) of the E. konishii genome. Both groups of genes were unevenly distributed in the E. konishii LGs (Figure 3). For example, LG02 alone harbored eight EkIAA genes, whereas no EkIAA gene mapped to LG07 or LG11 (Figure 3). Several EkIAA genes clustered in close proximity on LG01, LG2, LG03, and LG08 (Figure 3). The pattern for the EkARF genes was similar: six EkARF genes mapped to LG03, five to LG05, three each to LG06 and LG07, one each to LG09 and LG12, and none to LG02 (Figure 3).

During genome evolution, gene duplication and neofunctionalization are driven by tandem and segmental duplication (Cannon et al., 2004; Freeling, 2009). To elucidate the expansion of the EkIAA and EkARF gene families in E. konishii, we studied their segmental and tandem duplications. We identified 20 instances of segmental duplication events, involving 11 EkIAA and 20 EkARF genes, but no tandem duplication events in either gene family (Figure 3). For the 10 pairs of segmental duplicated EkIAA genes and the 10 pairs of segmental duplicated EkARF genes identified above, we calculated the ratios of nonsynonymous to synonymous substitutions (K _a/K _s) to evaluate their molecular evolutionary rates (Supplementary Table S2). The ratios for all duplicated pairs were less than 1 (Supplementary Table S2), suggesting that duplicated pairs of genes underwent purifying selection during evolution, thus raising the possibility that the biochemical characteristics of these EkIAAs and EkIAAs may not have changed very much since the initial duplication event.

Protein motifs are critical for protein function and structure maintenance (Smith-Gill, 1991). Accordingly, we looked for conserved functional motifs in the predicted EkIAA and EkARF proteins with the MEME web server tool. We identified four domains conserved in EkIAA proteins (motifs 1–4), corresponding to IAA domains Ⅳ, Ⅲ, Ⅱ, and Ⅰ, respectively (Figure 4). Of the 34 EkIAA proteins, 23 (61.8%) contained all four conserved domains (domains I–IV) (Figure 4). Some EkIAA proteins lost one or more domains: For example, EkIAA1.3, EkIAA29.1, EkIAA29.2, EkIAA29.3, and EkIAA29.4 lack domain I; EkIAA20/30, EkIAA32/34, EkIAA33.1, and EkIAA33.2 lack domain II; EkIAA1.5 lacks domains I and IV; EkIAA27.2 lacks domains III and IV; and EkIAA1.4 and EkIAA27.3 have only domain I (Figure 4).

To understand the evolution of EkIAA genes, we examined their exon-intron structures (Figure 4A). Most EkIAA genes consisted of five exons and four introns, as in Arabidopsis, although introns in EkIAA genes were larger than those in their Arabidopsis counterparts (Figure 4). Three EkIAA genes (EkIAA27.2, EkIAA33.1, and EkIAA33.2) comprised only two exons and one intron, with another six EkIAA genes (EkIAA1.1, EkIAA1.2, EkIAA1.4, EkIAA1.5, EkIAA5.1, and EkIAA5.2) having three exons and two introns. Six EkIAA genes (EkIAA20/30, EkIAA27.3, EkIAA29.1, EkIAA29.2, EkIAA29.4, and EkIAA32/34) had four exons and three introns. EkIAA9.2 had six exons and five introns, while EkIAA29.3 had by far the most exons (20) and 19 introns (Figure 4). As the presence of conserved domains in the EkIAA proteins and the EkIAA gene structure are similar to those in their Arabidopsis orthologs (Remington et al., 2004), these results suggested that the EkIAA family in E. konishii is conserved.

We analyzed the conserved motifs and gene structure of EkARF proteins and EkARF genes, respectively (Supplementary Figure S1). EkARF proteins belonging to the same clade in the phylogenetic tree had the same functional motifs (Supplementary Figure S1). The DNA binding domain was represented by motifs 1, 9, and 10, while motifs 3, 6, and 8 matched the variable middle transcriptional regulatory region (MR). Motifs 7 and 5 formed part of the C-terminal dimerization domain (CTD) (Supplementary Figure S1). Of the 29 EkARFs, 21 (72.4%) contained all three functional domains, with only eight EkARFs (EkARF3, EkARF10.1, EkARF10.2, EkARF16.1, EkARF16.2, EkARF16.3, EkARF17.1, and EkARF17.2) lacking the CTD (Supplementary Figure S1). This variation in functional protein motifs may reflect mutations or deletions in the gene structure. Most EkARF genes contained 14 exons and 13 introns (Supplementary Figure S1). However, EkARF17.1 consisted of only two exons and one intron, three EkARF genes (EkARF16.1, EkARF16.2, and EkARF16.3) were composed of three exons and two introns, EkARF17.2 had five exons and four introns, and EkARF3 contained 11 exons and 10 introns, which is consistent with the observed variation in protein domains (Supplementary Figure S1). Overall, EkIAA and EkARF genes appeared to be relatively conserved during evolution, but those derived from segmental duplication have experienced some structural divergence.

To explore the transcriptional regulation of EkARF and EkIAA genes and predict their functions, we analyzed the cis-regulatory elements in their promoters. We extracted 2,000 bp of upstream sequence, which we submitted to the PlantCARE online tool (Lescot et al., 2002). We then counted the number of phytohormone-, environment-, and flavonoid-responsive elements and noted their locations (Figure 5). cis-elements in EkIAA and EkARF promoters exhibited a similar pattern (Figure 5). Indeed, phytohormone-responsive and environmental stress–related cis-elements were present in all promoters of the EkIAA and EkARF gene family in E. konishii.

We then used publicly available RNA-seq datasets to analyze the expression patterns of EkIAA and EkARF genes at four fruit developmental stages, including green, turning red fruit, and red-winged pericarp stages, as well as in branches and leaves (Yang et al., 2020); the results are summarized as heatmaps in Figure 6. EkIAA gene family members showed varying expression patterns. Most EkIAA genes were highly expressed in green fruits, with the exception of EkIAA15.2, EkIAA33.1, EkIAA33.2, and EkIAA29.3, of which the first two were expressed specifically in branches, whereas the latter two genes were specifically expressed in red-winged pericarp (Figure 6A). EkIAA27.5 and EkIAA32/34 were highly expressed during the red fruit stage. EkARF genes were highly expressed during the fruit maturation stage (Figure 6B), of which EkARF1.2, EkARF2.2, EkARF4.2, EkARF5.1, EkARF5.2, EkARF16.3, and EkARF17.1 showed high expression levels in the red fruit stage. Most EkIAA and EkARF genes were highly expressed in fruits, hinting at their potential involvement in fruit maturation, including the accumulation of associated secondary metabolites.

To assess the extent of functional diversification within these two families, we focused on duplicated gene pairs (10 EkIAA pairs and 10 EkARF pairs) and calculated the Pearson’s correlation coefficients of their expression profiles. Several duplicated gene pairs did in fact exhibit differential expression across the samples tested. Six EkIAA gene pairs (EkIAA12/13.1/-12/13.2, EkIAA14.1/-14.2, EkIAA14.1/-16.1, EkIAA14.1/-16.2, EkIAA14.2/-16.1, and EkIAA16.1/-16.2) and three EkARF gene pairs (EkARF5.1/-5.2, EkARF8.1/-8.2, and EkARF24.1/-24.2) showed similar expression patterns within the pairs, as evidenced by their high correlation coefficients (Figure 6B). We also identified one EkIAA pair and one EkARF pair with distinct expression levels between duplicated copies: EkIAA27.5 (expressed at high levels in all samples) and the duplicated copy EkIAA27.4 (expressed at relatively low levels), with a correlation coefficient of 0.15 (Figure 6C); and EkARF11/18.1 (expressed at low levels in all tissues) and EkARF18.2 (highly expressed in branches), with a correlation coefficient of 0.45. (Figure 6D). These results suggested that the functions of duplicated genes may have diverged following the initial duplication event.

Medicinal compounds such as triterpenes, phenolic acids, and flavonoids have been isolated from Euscaphis fruits, leaves, and roots (Liang et al., 2018). The accumulation of anthocyanin and terpenoid secondary metabolites coincides with E. konishii fruit maturation (Yuan et al., 2018b; Liang et al., 2019). Generally, genes with similar expression patterns may have related roles as they belong to the same regulatory pathway or are regulated by the same upstream factors. Thus, co-expression detected from our correlation analyses may provide cues as to gene regulation or function. Because anthocyanin contents, and the key anthocyanin biosynthetic genes, were previously well characterized during Euscaphis fruit development (Yuan et al., 2018b), we performed a correlation analysis between anthocyanin levels and the expression estimates for anthocyanin biosynthesis genes, EkIAA, and EkARF genes. We determined that anthocyanin contents were positively and strongly correlated with the expression of five genes encoding key enzymes (CHALCONE SYNTHASE8 [CHS8], CHALCONE ISOMERASE2 [CHI2], FLAVANONE 3-HYDROXYLASE1 [F3H1], F3H2, and F3H3) and very strongly correlated with that of another eight genes (CHS2, CHI3, F3H4, FLAVONOID 3′-HYDROXYLASE [F3′H], DIHYDROFLAVONOL 4-REDUCTASE1 [DFR1], LEUCOANTHOCYANIDIN REDUCTASE [LAR], FLAVONOL 3-O-GLUCOSYLTRANSFERASE3 [UFGT3], and UFGT7) (Figure 7A). Because ARF proteins regulate the expression of their target genes by binding to their cognate cis-elements in promoters, we looked for EkARF genes co-expressed with the anthocyanin biosynthetic genes listed above (Figure 7A), leading to the identification of seven such genes (EkARF1.2, EkARF2.2, EkARF4.2, EkARF5.1, EkARF5.2, EkARF16.3, and EkARF17.1).

Triterpenoids accumulate to high levels in Euscaphis fruits and are important raw materials for natural products, food additives, and chemical products. Genes encoding the enzymes involved in the biosynthesis of triterpenoids have been described in the Euscaphis genome (Huang et al., 2019; Liang et al., 2019), prompting us to test for correlations between their expression patterns and those of ARFs genes (Figure 7A). This analysis highlighted seven EkARF genes (EkARF1.2, EkARF2.2, EkARF4.2, EkARF5.1, EkARF5.2, EkARF16.3, and EkARF17.1) whose expression was positively and strongly correlated with genes involved in triterpenoid accumulation (Figure 7). These results suggested that ARF genes contribute to secondary metabolite biosynthesis in Euscaphis.

ARF proteins are transcription factors that can bind to AuxREs (TGTCTC) to regulate the expression of their target genes. We noticed at least one AuxRE either upstream, downstream, or within intronic regions of anthocyanin and triterpenoid biosynthetic genes (Supplementary Table S3). To assess the role of ARFs in anthocyanin and triterpenoid biosynthesis, we selected EkARF5.1 for further characterization. RT-qPCR showed that EkARF5.1 is highly expressed during the red fruit stage, which is consistent with the RNA-seq data (Figure 7B). We fused EkARF5.1 to the GAL4 activation domain (AD) (AD-EkARF5.1) and introduced the resulting construct into yeast strain Y1HGold carrying a reporter consisting of seven copies of the AuxRE sequence driving the expression of Aureobasidin Resistance 1, conferring resistance to the antibiotic aureobasidin A (AbA). Whereas yeast colonies carrying AD or AD-EkARF5.1 grew on synthetic defined medium (Figure 7C), only yeast cells harboring the AD-EkARF5.1 construct survived growth on AbA-containing medium (Figure 7C), supporting the notion that EkARF5.1 binds to the AuxRE. We then tested the transactivation activity of EkARF5.1 by fusing full-length EkARF5.1 to the GAL4 DNA binding domain (BD) to generate the BD-EkARF5.1 fusion protein; the resulting encoding construct was introduced into yeast strain Y2HGold, which harbors the His3 gene driven by a GAL4-responsive promoter (Figure 7C). Only yeast cells carrying the BD-EkARF5.1 construct survived on synthetic medium lacking histidine, unlike yeast cells carrying the empty GAL4 DB vector (Figure 7C). These results indicated that EkARF5.1 has transactivation activity in yeast cells. Finally, we tested EkARF5.1 in N. benthamiana leaf epidermal cells by co-infiltrating a construct overexpressing EkARF5.1 and a β-GLUCURONIDASE (GUS) reporter construct whose expression is driven by seven copies of the AuxRE. We detected GUS activity in plant cells only when EkARF5 was co-expressed (Figure 7D). These results demonstrate that EkARF5.1 may play a role during anthocyanin and triterpenoid biosynthesis.