Genome-Wide Identiﬁcation of the Auxin Response Factor ( ARF ) Gene Family and Their Expression Analysis during Flower Development of Osmanthus fragrans

: Auxins have long been implicated in many aspects of plant growth and development. Auxin response factors (ARFs) are important proteins in auxin-mediated pathways and they play key roles in plant physiological and biochemical processes, including ﬂower development. Endogenous indoleacetic acid (IAA) levels were measured and ARFs were studied in the ﬂowers during the developmental stages in order to further elucidate the role of auxin in ﬂower development of Osmanthus fragrans . A systematic analysis of OfARFs was conducted by carrying out a genome-wide search of ARFs . A total of 50 ARF genes ( OfARF s) were detected and validated from the Osmanthus fragrans genome. Furthermore, a comprehensive overview of the OfARF s was undertaken, including phylogenetic relationship, gene structures, conserved domains, motifs, promoters, chromosome locations, gene duplications, and subcellular locations of the gene product. Finally, expression proﬁling, while using transcriptome sequencing from a previous study and quantitative real-time PCR (qRT-PCR), revealed that many OfARF genes have di ﬀ erent expression levels in various tissues and ﬂower developmental stages. By comparing the expression proﬁles among the ﬂower developmental stages, and the relationship between ARFs and endogenous IAA levels, it can be supposed that OfARF s function in ﬂower development of O. fragrans in an auxin-mediated pathway.


Introduction
Flowering is a crucial developmental stage in the life cycle of plants, especially for ornamental flowering plants [1]. Previous studies have indicated that flowers react differently to the various endogenous phytohormone types and levels. For example, the levels of indole-3-acetic acid (IAA), gibberellic acid (GA 3 ), and zeatin riboside (ZR) are downregulated, while the levels of abscisic acid (ABA) are upregulated during florescence and senescence in the petals of Paeonia suffruticosa [2]. Auxins, a class of phytohormone previously reported to be an essential factor in plant growth and development, play a pivotal role in the development process of flowers, especially in floral primordium initiation and were elucidated. Finally, the endogenous concentrations of IAA at different flowering periods were determined, and then correlated with the expression profile of OfARF genes in various O. fragrans tissues, including the flower, during the developmental process, with the goal of providing a sound basis for analysis of the auxin-mediated pathway in flowers and to further functional studies of OfARF genes.

Identification of ARF Genes in Osmanthus fragrans
A systematic identification of the ARF family members (OfARFs) in O. fragrans was conducted while using the full-length Arabidopsis ARF (AtARF) proteins previously achieved as initial BLAST queries out of the O. fragrans Genomic Database with the e-value of "10 −3 " [33,34]. Moreover, the ARF Hidden Markov Model (HMM), Pfam 06507, ARF (AUX_RESP), was also employed to hmmsearch the ARFs from the O. fragrans genome, with the e-value of "10 −3 " [35]. Furthermore, the HMMSCAN (http//www.ebi.ac.uk/Tools/hmmer/search/hmmscan) and NCBI Conserved Domain Search Service (CDSearch) (http//www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) tools were used to further confirm the position and integrity of the detected conserved domains in individual OfARFs [36,37]. The sequences with broken B3, AUX_RESP, or AUX/IAA domains were deleted. The number of amino acids, molecular weight (MW), and theoretical isoelectric point (pI) of the OfARFs were obtained from the online ProtParam tool (http//www.expasy.org/tools/protparam.html) [38].

Analysis of Gene Structures, Conserved Domains and Subcellular Localizations
The gene structure details of OfARFs were obtained from the O. fragrans genome database [33]. The conserved domains in all of the identified OfARFs were analyzed and confirmed while applying the Pfam tool (http//pfam.sanger.ac.uk/) [39]. The gene structure information and conserved domain distributions of OfARFs were then drawn while using TBtools [40]. The subcellular localization predictions of the OfARFs were performed on the database Subcellular Localization Predictor (CELLO v.2.5; http//cello.life.nctu.edu.tw) [41].

Chromosomal Localization, Gene Duplication and Evolution Analysis of OfARF Genes
The OfARF chromosomal locations were obtained from the O. fragrans Genomic Database [33], and then drafted, using MapInspect software (http//www.plantbreeding.wur.nl/uk/softwaremapinspect. html). The OfARF protein sequences were compared to themselves with the Blast Compare program in Tbtools, with the e-value of 1 × 10 −5 , NumOfHits of 5, NumOfAligns of 5. Synteny analysis of OfARFs was conducted with the blast result and the O. fragrans gff file while using Quick MCScanX Wrapper program in TBtools. The graph of collinearity analysis was drawn by Amaizing Super Circos in TBtools software. Finally, the tool Circos were applied to locate OfARFs, and draft duplication relations of the 23 chromosomes in O. fragrans [42].

Conserved Motifs and Promoter Cis-acting Regulatory DNA Elements Analysis
The 12 conserved protein motifs were obtained while using the MEME program (http://memesuite.org/tools/meme) and the optimum motif width was set to . The PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was employed to analyze the cis-acting regulatory elements of ARF genes in O. fragrans [43].

Phylogenetic Analysis of OfARF Genes
The Arabidopsis ARF (AtARF) protein sequences were downloaded from the online PlnTFDB as described in a previous report [34,44]. The OfARF proteins were sequence aligned by the DNAMAN (v8.0) software with its default settings. The sequences of O. fragrans and A. thaliana were aligned while using the Muscle program that was built in MEGAX software. The phylogenetic relationship tree was constructed and drafted with MEGAX, using the neighbor-joining method [45]. The bootstrap was set with 1000 replicates. Additionally, two separate phylogenetic trees were constructed to explore the relationships among the OfARFs themselves and between OfARFs and AtARFs. The images of phylogenetic trees were drawn while using Fig Tree v1.4.2 (http://tree.bio.ed.ac.uk/). The gene structures of OfARFs were predicted using TBtools software with GFF3 file and ID numbers of OfARFs.

Plant Material Collection, RNA extraction and cDNA Synthesis
In the present study, the flower samples were collected from a mature O. fragrans 'Rixianggui' on the campus of Nanjing Forestry University, Nanjing, Jiangsu Province, China. All of the samples were collected and immediately frozen in liquid nitrogen and then stored at −80 • C until further use. The flower developmental process in cultivar Rixianggui generally goes through the following stages of flower development: linggeng stage (S1), bud eye stage (S2), initial flowering stage (S3), full flowering stage (S4), and late flowering stage (S5). The total RNAs were extracted while using the TRNzol Universal Reagent (TIANGEN Biotech Co. Ltd., Beijing, China), according to the manufacturer's instructions [46]. The RNA quantity and quality were evaluated while using Agilent Bioanalyzer 2100 (Agilent Technologies, USA) and NanoDrop™ One UV-Vis spectrophotometer (Thermo Fisher Scientific, USA), and Qubit ® 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). The cDNA synthesis was completed while using TransScript one-step gDNA Removal and cDNA Synthesis SuperMix (Transgen, Beijing, China), following the manufacturer's guidelines [47].

Expression Analysis of OfARF Genes
Transcriptome sequencing was applied to assess the expression levels of 50 OfARF genes in roots (young and mature), stems, leaves, and flowers at three developmental stages (initial, full, and late flowering stages). The TBtools software was used to generate the heatmap, while using transcriptome data obtained from a previous study [33]. Expression level was determined according to reads per kilobase of transcript per million mapped reads (RPKM). RPKM values of OfARF genes were normalized by log2 transform. The edgeR package (http://www.r-project.org/) was used to identify differentially expressed genes (DEGs) across initial flowering stage (S1), full flowering stage (S2), and full flowering stage (S3) flower samples. Identified DEGs were then used for enrichment analysis of GO functions. All DEGs were searching against GO terms in the Gene Ontology database (http://www.geneontology.org/). Numbers of the genes were calculated for every term. Enriched GO terms in DEGs when comparing to the genome background were defined by hypergeometric test. Furthermore, in order to verify the expression profiles of OfARF genes at different developmental stages of O. fragrans flowers, twelve OfARF genes from Class II, Class III, and Class IV were selected for quantitative real-time PCR (qRT-PCR) analysis. The genes were selected based on the RPKM values and the phylogenetic relationship of 50 Osmanthus fragrans auxin response factors (OfARFs) and 23 Arabidopsis thaliana auxin response factors (AtARFs). OfARF11a, 12, 13, 14a, and 14b were homologous genes of AtARF1 and AtARF2. OfARF10a and OfARF10b were homologs of AtARF4. OfARF21a, 22a, 22b, 24a, and 24b were homologs of AtARF6 and AtARF8. The gene-specific primers that were used in qRT-PCR were designed while using Primer Premier 5 software (Supplementary Table S1). The OfACT gene (forward, 5 -CCCAAGGCAAACAGAGAAAAAAT-3 ; reverse, 5 -ACCCCATCACCAGAATCAAGAA-3 ), which was the optimal reference gene for different floral developmental stages in O. fragrans, was used as an internal control in each qRT-PCR experiment [48]. The qRT-PCR was performed, while using an ABI 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) with SYBR Green Realtime PCR Master Mix. The total reaction volume was 10 µL, which included 5 µL of SYBR Mix, 3 µL of ddH2O, 1 µL of cDNA sample, 0.4 µL of each gene-specific forward and reverse primer, and 0.2 µL of ROX Reference Dye. The PCR conditions were, as follows: polymerase activation for 3 min. at 94 • C, and then 20 s at 94 • C for 40 cycles, followed by 20 s at 62 • C, and finally 40 s at 72 • C. The amplification specificity of the primers was tested by the melting curve method. Three biological and three technical replicates were conducted for each sample in the qRT-PCR system. The comparative cycle threshold (44Ct) and the 2 −∆∆CT method were used to evaluate the relative expression level of all twelve genes. The error bars represent the standard error from the three biological replicates. Analysis of variance (ANOVA) was used to detect significant differences within the experiment, with pairwise multiple comparisons to identify significant differences between the samples being determined with Duncan's multiple range test at p < 0.05; analyses were carried out while using SPSS v.24 (IBM, Armonk, NY, USA) [46]. In the transcriptome analysis, there were three flowering stages (initial, full, and late flowering stages) ( Figure  S3). However, five flower stages, including linggeng stage (S1), bud eye stage (S2), initial flowering stage (S3), full flowering stage (S4), and late flowering stage (S5-1 & S5-2), were used in the qRT-PCR process in order to analyze the expressing profiles of ARFs in the whole flower development process ( Figure S4).

IAA Content Measurement
The flower samples were collected, and the surfaces of them were washed three times in deionized water. The samples were then blotted dry with a clean paper towel and weighed. After the addition of 500 pg of the 13 C 6 -IAA internal standard (MSD, ISOTOPES, Quebec, Canada), five independent biological replicates, each representing a 50-mg sample, were purified while using a ProElu C18 column (http://www.dikma.com.cn). The IAA concentrations were determined by separation and detection with FOCUS GC-DSQII single quadrupole GC-MS (Thermo Fisher Scientific Inc., USA) [49].

Identification of the OfARF Gene Family
Two strategies were employed to comprehensively identify the ARFs in O. fragrans: BLASTP using the 23 full-length Arabidopsis ARF (AtARF) protein sequences previously identified as queries [34], and HMM search with the AUX_RESP HMM model (PF06507) against the O. fragrans genome database. HMMSCAN (http//www.ebi.ac.uk/Tools/hmmer/search/hmmscan) and NCBI Conserved Domain Search Service (CDSearch) (http//www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) were applied to ensure the position and integrity of the detected conserved domains [38]. The sequences with broken domains (B3, AUX_RESP, or AUX/IAA) were deleted. Finally, a total of 50 ARF genes were identified in O. fragrans, including 39 complete genes and 11 truncated genes. The entire collection of OfARFs that were identified in this study was named OfARF1 to OfARF29b, based on the homology relationship that was shown in the phylogenetic tree. The coding sequence (CDS) lengths of the OfARF genes ranged from 999 bp to 3327 bp, with deduced proteins of 332-1108 amino acids. The molecular weights (MW) of the 50 OfARF proteins were predicted to range from 37.46 to 122.67 kDa ( Table 1). The theoretical isoelectric points (pI) were predicted to range from 5.22-9.05. Among them, seven were basic proteins, and the left forty-three were acidic proteins. The instability indexes were predicted to range from 41.46-68. 19, which suggested that all of the proteins were unstable proteins. The grand averages of hydropathicity (GRAVY) values of them were all negative, which suggested that they were all hydrophilic proteins.

Phylogenetic Analysis and Gene Structure of OfARFs
A phylogenetic tree was constructed by aligning all the 50 OfARF and 23 AtARF protein sequences, while using the neighbor-joining method, to detect phylogenetic relationships. It was revealed that the 50 OfARF proteins were broadly distributed into four classes, namely Class I, Class II, Class III, and Class IV, with well-supported bootstrap values ( Figure 1). Among all of the OfARFs, Class IV had the largest number (16) of OfARFs. In contrast, there were only four members in Class II, with Class I and Class III each containing fifteen members (Figure 2A). The exon-intron components of the genes were obtained by comparing the cDNA sequences with the corresponding genomic DNA sequences to compare the structure of the different OfARF genes ( Figure 2B). The number of exons varied from 2 to 15 (Table 1). Generally, the members that shared similar gene structures tended to be distributed in the same subfamilies.

Analysis of Conserved Domains of Osmanthus fragrans Auxin Response Factor (OfARF) Proteins
The domains of OfARF protein sequences were aligned and analyzed while using the online Pfam tool (http//pfam.sanger.ac.uk/) and SMART database (http://smart.embl-heidelberg.de/). Three conserved domains, namely B3, AUX_RESP, and AUX_IAA, were identified and drafted into four classes according to the phylogenetic tree ( Figure 2C). All of the OfARF proteins, except OfARF7a-OfARF9b, OfARF15, OfARF16b, and OfARF17, contained a C-terminal Aux/IAA domain, which was consistent with AtARF3, 13, and 17. However, when compared with A. thaliana, OfARF25b showed an absence of DNA-binding domains (DBDs), which indicated that there were some differences between A. thaliana and O. fragrans ARF genes.

Motif Characterization and Promoter Analysis of Osmanthus fragrans Auxin Response Factor (OfARF) Proteins
The online software MEME was employed to investigate the motif distributions in all 50 OfARF proteins, with 12 individual motifs identified, to further study the characteristic regions of OfARF proteins. Among all the motifs, motifs 2, 3, 6, and 7 were found in all 50 OfARF proteins. As predicted, those with common motif compositions came mostly from the same subfamily from the phylogenetic

Analysis of Conserved Domains of Osmanthus fragrans Auxin Response Factor (OfARF) Proteins
The domains of OfARF protein sequences were aligned and analyzed while using the online Pfam tool (http//pfam.sanger.ac.uk/) and SMART database (http://smart.embl-heidelberg.de/). Three conserved domains, namely B3, AUX_RESP, and AUX_IAA, were identified and drafted into four classes according to the phylogenetic tree ( Figure 2C). All of the OfARF proteins, except OfARF7a-OfARF9b, OfARF15, OfARF16b, and OfARF17, contained a C-terminal Aux/IAA domain, which was consistent with AtARF3, 13, and 17. However, when compared with A. thaliana, OfARF25b showed an absence of DNA-binding domains (DBDs), which indicated that there were some differences between A. thaliana and O. fragrans ARF genes.

Motif Characterization and Promoter Analysis of Osmanthus fragrans Auxin Response Factor (OfARF) Proteins
The online software MEME was employed to investigate the motif distributions in all 50 OfARF proteins, with 12 individual motifs identified, to further study the characteristic regions of OfARF proteins. Among all the motifs, motifs 2, 3, 6, and 7 were found in all 50 OfARF proteins. As predicted, those with common motif compositions came mostly from the same subfamily from the phylogenetic tree, which suggests functional similarities among the close homologs ( Figure 3). By comparing with the conserved domain distribution, it was concluded that motif 4 and 5 belong to B3 domain, motif 7 and motif 8 belong to Auxin_resp domain, and motif 11 and 12 belong to Aux/IAA domain. Previous studies indicated that AtARF1, 2, 3, 4, and 9 act as transcriptional repressors [13,14], which have middle regions rich in threonine, serine, and proline residues, whereas AtARF5, 6, 7, 8, and 19 were transcriptional activators [15][16][17], with middle regions that were rich in glutamine residues. The detailed sequence analysis of all 50 OfARF proteins identified threonine-, serine-, and proline-rich middle regions in OfARF9a, 9b, 10a, 10b, 11a, 11b, 12, 13, 14a, 14b, 20, and 21, indicating that they are likely to be repressors (Figure 4). The glutamine-rich middle regions were detected in OfARF25a, 25b, 21a, 21b, 23, 24a, 24b, 26a, 26b, and 27, indicating that they are probably activators of the O. fragrans flower development process. Putative promoter sequences (2000-bp upstream of the 5′UTR (untranslated region) of OfARFs were isolated from the O. fragrans genome sequence. Several auxin response elements were detected through searching with the PlantCARE database. Three auxin response elements (two TGA-boxes and one AuxRR element) were found in the promoter region of the OfARF4a gene, and three auxin response elements (two TGA-elements and one AuxRR element) were detected in OfARF4a. Two TGA-elements were detected in OfARF2 and OfARF24b, whereas two AuxRRs were discovered in OfARF7b. Only one TGA-element was detected in each of 10 other OfARF genes, and one AuxRR Putative promoter sequences (2000-bp upstream of the 5 UTR (untranslated region) of OfARFs were isolated from the O. fragrans genome sequence. Several auxin response elements were detected through searching with the PlantCARE database. Three auxin response elements (two TGA-boxes and one AuxRR element) were found in the promoter region of the OfARF4a gene, and three auxin response elements (two TGA-elements and one AuxRR element) were detected in OfARF4a. Two TGA-elements were detected in OfARF2 and OfARF24b, whereas two AuxRRs were discovered in OfARF7b. Only one TGA-element was detected in each of 10 other OfARF genes, and one AuxRR element was present in OfARF23 ( Figure 5). However, no auxin response elements were detected in 34 OfARF genes.
Forests 2020, 11, x FOR PEER REVIEW 11 of 21 element was present in OfARF23 ( Figure 5). However, no auxin response elements were detected in 34 OfARF genes.   Figure 6). A collinearity analysis was performed while using MCScanX software to reveal the duplication events that occurred within the OfARF genes. Thirty-three fragment duplication gene pairs were found among the 50 OfARFs. Fragment duplication genes were the most abundant on chromosome 03, followed by chromosomes 07 and 13, whereas there were no pairs of fragment duplication genes on chromosomes 09 and 23 (Figure 7).

Subcellular Localization Analysis of OfARFs
The online tool CELLO V2.5 was utilized to predict the subcellular localization of OfARFs (http//cello.life.nctu.edu.tw/). To detect the subcellular localization, the full-length coding sequences of OfARFs were employed. Only OfARF3a, OfARF7b, and OfARF8a were predicted to be located in the chloroplast, while all of the remaining 47 proteins were predicted to be located in the nucleus (Supplementary Table S2). The prediction results were consistent with the function of being transcription factors as expected, which should regulate the expressions of the downstream genes.

Expression of OfARF Genes in Different Tissues and Various Flower Development Stages of O. fragrans
The spatial-specific expression profiles of the 41 OfARF genes were determined in different tissues (nine genes not expressed), namely roots, stems, leaves (young and mature), and flowers based on data from the sequencing of transcriptomes that were performed in the previous study ( Figure 8). The expression of most OfARF genes occurred in almost all tissues, which indicated that they might have various functions in different aspects of plant physical and biological events.
The expression profile of ARF genes in flower development, as previously reported, as well as the role of auxin in controlling flower development, prompted us to analyze the expression of OfARF genes during the O. fragrans flower development process [2][3][4]. The transcript abundance level of individual OfARF genes was determined at five developmental stages of O. fragrans flowers by qRT-PCR to investigate the expression profiles. The results indicated that the expression of all twelve genes tested was down-regulated from stage S1 to stage S3, and it was markedly up-regulated from stage S3 to stage S5. Similar expression patterns for all twelve OfARF genes indicate that they play important roles during flower developmental processes in O. fragrans (Figure 9). DEGs were found in S1 vs. S3 with 15 (OfARF2, 4a, 4b, 5b, 7b, 12, 13, 14a, 21a, 22a, 25b, 26b, 27, 28a, 28b) upregulated and 20 (OfARF1a, 1b, 3a, 6a, 6b, 7a, 9a, 9b, 10a, 10b, 11a, 15, 16b, 17a, 17b, 24a, 24b, 25a, 29a, 29b) downregulated, respectively ( Figure S1). These identified DEGs were then used for enrichment analysis of GO functions. In S1 vs. S2, the most enriched term in "biological process" was "response to stimulus" (Counts = 27), the most enriched terms in "cellular component" was "cell", "cellular anatomical entity", and "intracellular" (Counts = 26), and the most enriched term in "molecular function" was "binding" (Counts = 26). In S2 vs. S3, the most abundant term in "biological process" was "cellular process" (Counts = 29), the most abundant term in "cellular component" was "cellular anatomical entity", (Counts = 28), and the most abundant term in "molecular function" was "binding" (Counts = 29). In S1 vs. S3, the most enriched term in "biological process" was "cellular process" (Counts = 30), the most enriched term in "cellular component" was "cellular anatomical entity" (Counts = 30), and the most enriched term in "molecular function" was "binding" (Counts = 29). Totally, the most enriched term in "molecular function" of the three groups was "binding", which is in accordance with the functions of the ARFs that functions by binding to the downstream genes ( Figure S2). The expression profile of ARF genes in flower development, as previously reported, as well as the role of auxin in controlling flower development, prompted us to analyze the expression of OfARF genes during the O. fragrans flower development process [2][3][4]. The transcript abundance level of individual OfARF genes was determined at five developmental stages of O. fragrans flowers by qRT-PCR to investigate the expression profiles. The results indicated that the expression of all twelve genes tested was down-regulated from stage S1 to stage S3, and it was markedly up-regulated from stage S3 to stage S5. Similar expression patterns for all twelve OfARF genes indicate that they play important roles during flower developmental processes in O. fragrans (Figure 9). The flower development process in cultivar Rixianggui was marked by stages S1 (linggeng stage), S2 (bud eye stage), S3 (initial flowering stage), S4 (full flowering stage), and S5 (late flowering stage). Changes in the expression levels during different flower developmental stages depicted above are relative to RNA accumulation levels. The error bars represent the standard error from three biological replications. Significant differences detected with Duncan's multiple range test at p < 0.05 are indicated by different lowercase letters above the error bars.
After comparison, some differences would easily been seen between the expression profiles in the heatmap and qRT-PCR results. This is perhaps due to the different flower stage distributions between the transcriptome analysis and the qRT-PCR. In the transcriptome analysis, there were three Figure 9. Expression profiles of 12 Osmanthus fragrans auxin response factor genes (OfARFs) by qPCR. The flower development process in cultivar Rixianggui was marked by stages S1 (linggeng stage), S2 (bud eye stage), S3 (initial flowering stage), S4 (full flowering stage), and S5 (late flowering stage). Changes in the expression levels during different flower developmental stages depicted above are relative to RNA accumulation levels. The error bars represent the standard error from three biological replications. Significant differences detected with Duncan's multiple range test at p < 0.05 are indicated by different lowercase letters above the error bars.
After comparison, some differences would easily been seen between the expression profiles in the heatmap and qRT-PCR results. This is perhaps due to the different flower stage distributions between the transcriptome analysis and the qRT-PCR. In the transcriptome analysis, there were three flowering stages (initial flower, full flower, and late flower) (Figure 3). However, there were five flower stages: Linggeng stage (S1), bud eye stage (S2), initial flowering stage (S3), full flowering stage (S4), and Late lowering stage (S5-1 & S5-2) (Figure 4). The initial stages in the transcriptome analysis should be equivalent to the stages of S1, S2, and S3 in qRT-PCR, and the full flower stages should be equivalent to the stage of S4 in qRT-PCR, and late flowering stages should be equivalent with the stage of S5-1 in the qRT-PCR. For example, the expression profiles of OfARF10a, 10b in the heatmap showed a slow decline from the initial flowering stage to the full flowering stage, but a sharp decline from the full flowering stage to the late flowering stage, which were different with the results in qRT-PCR. This is perhaps because the expression levels of OfARF10a, 10b in the earlier time (linggeng stage, bud eye stage) were higher than the levels in the initial flowering stage, and there should be a sharp increase of expression level during the later time of the late flowering stage.

Endogenous Indoleacetic Acid (IAA) Concentration Measurement
Endogenous IAA levels were measured to examine the role of auxin concentration in different stages of flower development. The data showed that the endogenous IAA levels were higher in the flowers at the early stages of development than in the flowers at later stages. Specifically, the endogenous IAA level declined from S1 to stage S3, and then rose significantly from S3 to S5 ( Figure 10).

Endogenous Indoleacetic Acid (IAA) Concentration Measurement
Endogenous IAA levels were measured to examine the role of auxin concentration in different stages of flower development. The data showed that the endogenous IAA levels were higher in the flowers at the early stages of development than in the flowers at later stages. Specifically, the endogenous IAA level declined from S1 to stage S3, and then rose significantly from S3 to S5 ( Figure  10).

Discussions
Auxins play an important role in the development of plant reproductive organs, especially the flowers, as has been previously reported. These organs react differently to the diverse levels of auxins. It was hypothesized that the endogenous IAA concentration affected flower development, so the endogenous IAA concentration in the flowers was measured at different development stages in this study [4]. The O. fragrans ARFs were identified and characterized in order to investigate the putative functional genes involved during the development stages of flowering.
The auxin-mediated pathway is mainly comprised of the products of two gene families, the ARFs and the Aux/IAAs. As an important element in the auxin-signaling pathway, ARFs bind directly to downstream target genes and regulate their expression during flower development, and they are also involved in the reproduction of various plant species [50]. The most widely accepted model is that Aux/IAA inhibits transcriptional activities of ARF proteins by dimerizing with ARFs

Discussions
Auxins play an important role in the development of plant reproductive organs, especially the flowers, as has been previously reported. These organs react differently to the diverse levels of auxins. It was hypothesized that the endogenous IAA concentration affected flower development, so the endogenous IAA concentration in the flowers was measured at different development stages in this study [4]. The O. fragrans ARFs were identified and characterized in order to investigate the putative functional genes involved during the development stages of flowering.
The auxin-mediated pathway is mainly comprised of the products of two gene families, the ARFs and the Aux/IAAs. As an important element in the auxin-signaling pathway, ARFs bind directly to downstream target genes and regulate their expression during flower development, and they are also involved in the reproduction of various plant species [50]. The most widely accepted model is that Aux/IAA inhibits transcriptional activities of ARF proteins by dimerizing with ARFs and by recruiting additional proteins that can inhibit gene expression. When auxin concentration increases in the cell, the ubiquitin-proteasome pathway degrades Aux/IAA proteins. Thus, the ARF protein is available to homodimerize with itself or to heterodimerize with another ARF protein to regulate transcription [51]. Identification and analysis of OfARFs, as achieved in the current study, should allow us to reveal the mechanisms operating in auxin-mediated flower development of O. fragrans.
In the current study, the recently developed reference genome database of O. fragrans was employed to identify the members of the OfARF gene family [33]. The number of OfARF genes (n = 50) was greater than that of AtARF genes (n = 23), Oryza sativa ARF genes (n = 25), and Populus trichocarpan ARF genes (n = 39). Protein domain analysis provided useful information on the biological function of ARFs. A typical ARF contains three domains, namely DBD, MR, and CTD. Aux/IAAs bind to CTDs and form heterodimers. The presence of numerous OfARFs without CTDs indicated that some auxin-response genes in O. fragrans might be regulated in an auxin-independent way.
A phylogenetic tree was built to analyze the relationship of ARFs between O. fragrans and A. thaliana. The results showed that six sister gene pairs with high bootstrap values (≥99%) were identified between O. fragrans and Arabidopsis, which suggested that OfARFs were highly homologous to AtARFs. Many AtARFs have been already been reported, which reveals useful information on the respective biological functions in O. fragrans. In Arabidopsis, transcription factors AtARF1 and AtARF2 have been reported to function in flower initiation and senescence. OfARF11a, 12, 13, 14a, and 14b, which were homologous genes of AtARF1 and AtARF2, were highly expressed in floral development stages S1 and S5, indicating that they play key roles in early flower formation and floral senescence, whereas AtARF4 was previously reported to be involved in flower patterning [14]. OfARF10a, OfARF10b, which are both homologs of AtARF4, showed high expression levels in the early flower developmental stages. AtARF8 regulates a complex process by promoting petal growth, filament elongation, and flower bud opening [16]. The expression of OfARF24a and 24b, homologs of AtARF8, showed the greatest expression levels at the initial flowering stage, indicating a putative function in flower initiation. In Solanum lycopersicum, SlARF1-10 and 14-17 have been reported to function in regulation of cell division and enlargement during flower development [29]. In Carica papaya L., CpARF1 and CpARF2 have been reported to participate in flower bud formation. Furthermore, CpARF1, CpARF2, and CpARF4-related auxin expression regulation was decreased, while CpARF6-mediated auxin expression regulation was activated in the mature flowers [52].
It is well known that auxins may play important roles in the early flower development stages [4]. However, there are still no reports that endogenous IAA plays roles in O. fragrans flower development. The developmental pathways that are regulated by ARF proteins have been reported in a previous study [51,53]. In the absence of auxin, ARF proteins bind to Aux/IAA proteins that inhibit transcription by chromatin remodeling. Auxins induce Aux/IAA turnover, freeing ARF proteins to activate transcription. In the current study, the endogenous IAA concentrations decreased from S1 to S3, and then increased from S3 to S5, which suggested that a high level of endogenous IAA contributes to both the initiation and senescence of flowers. ARF activity is inhibited at lower auxin levels, which results in the repression of auxin-responsive genes, so that the ARF is released and can regulate the transcription of its target auxin response genes when the auxin levels are higher, as has been previously reported [4]. The transcript abundance level of OfARF genes was measured at five developmental stages of O. fragrans flowers. The results indicated that all twelve selected genes were down-regulated from S1 to S3 and they were markedly upregulated from S3 to S5. Similar expression patterns for all twelve OfARF genes indicated that they might play important roles during flower developmental processes in O. fragrans. The relationship between the IAA levels and ARF gene expression profiles during flower development showed a significant positive correlation, which is in accordance with the earlier hypothesis of a positive link between IAA concentration, flower development, and ARF gene expression (Supplementary Table S3).

Conclusions
This is the first study on the genome-wide identification and analysis of the ARF family members in O. fragrans. The comprehensive analysis of the ARF gene family in O. fragrans has revealed structural characterization, molecular evolution, expression profiling, and functional properties of the OfARFs, especially the functions that are related to flower development. Studies of the phylogenetic relationship and chromosomal distribution of the OfARF genes have provided valuable insights into the evolutionary aspects of the O. fragrans genome. Expression analysis showed that all twelve selected genes might play important roles during flower developmental processes in O. fragrans. The relationship between the IAA levels and ARF gene expression profiles during flower development showed a significant positive correlation, especially OfARF11a, 12, 13, and 14a, which could be taken as the target candidate genes for the future transgenic improvement. The results of an expanded expression analysis of all identified OfARF genes relative to endogenous IAA concentration will help to focus molecular genetic studies, leading to a better understanding of the functions of the OfARF genes in O. fragrans and future practical applications of this knowledge, such as extending the flowering period of O. fragrans. The comprehensive identification and subsequent characterization of OfARF genes that are described above would provide new insight into the potential role of some ARF genes in mediating O. fragrans flower development via the auxin-mediated pathway.