Genome-wide analysis of the Dof gene family in durian reveals fruit ripening-associated and cultivar-dependent Dof transcription factors

DNA binding with one finger (Dof) proteins constitute a ubiquitous plant-specific transcription factor (TF) family associated with diverse biological processes, including ripening. We conducted a genome-wide analysis of durian (Durio zibethinus Murr.) and identified 24 durian Dofs (DzDofs), 15 of which were expressed in fruit pulp. Gene expression analysis revealed differential expression of DzDofs during ripening in two commercial durian cultivars from Thailand, Monthong and Chanee. Comparing the expression levels of fruit pulp-expressed DzDofs between cultivars revealed ten potential cultivar-dependent Dofs, among which DzDof2.2 showed a significantly greater fold increase at every ripening stage in Chanee than in Monthong. The prediction of DzDof2.2’s function based on its orthologue in Arabidopsis revealed its possible role in regulating auxin biosynthesis. We observed significantly higher auxin levels during ripening of Chanee than Monthong which concurred with the greater expression of auxin biosynthetic genes. Transient expression of DzDof2.2 in Nicotiana benthamiana significantly upregulated the expression levels of auxin biosynthetic genes. Higher expression levels of DzDof2.2 in Chanee would enhance auxin levels through transcriptional regulation of auxin biosynthetic genes. Higher auxin levels in Chanee could activate auxin-mediated transcription, contributing to its faster ripening compared to Monthong through earlier initiation of the ethylene response (auxin-ethylene crosstalk).

TFs in climacteric ripening. In a study by Zhang et al. 12 , the role of a kiwifruit Dof (AdDof3) in starch degradation during ripening was investigated. AdDof3 was shown to physically interact with the AdBAM3L (β-amylase) promoter, a key gene involved in starch degradation.
Durian (Durio zibethinus Murr.) is an economically important tropical fruit crop native to Southeast Asia, which has recently begun to be distributed globally. Durian is distinctive for its strong and unique odour, formidable spiny husk, and overwhelming flavour. It plays a vital role as an important export commodity in the Southeast Asian region. Thailand is the top exporter of durian across this region. There are more than 200 durian cultivars, but only a few cultivars are commercially grown and are in high demand in the international market. In Thailand, two cultivars, Monthong (D. zibethinus Murr. cv. 'Monthong') and Chanee (D. zibethinus Murr. cv. 'Chanee') are cultivated commercially in several provinces, particularly in the eastern and southern parts of Thailand 15 . Notably, these two cultivars exhibit distinctive sensory characteristics from one another. Monthong is a less sweet cultivar with a mild odour and light yellowish pulp, whereas Chanee has an overpowering aroma, sweeter taste, and dark yellowish pulp. Generally, Thai farmers harvest durian at the commercially mature stage, which occurs around 90 and 105 days after anthesis for Chanee and Monthong, respectively. Harvesting is done at the mature stage to facilitate fruit handling and transportation, protect the fruits from pathogen infection on the tree, and enhance their sensory characteristics. Chanee ripens four days after harvesting at the mature stage, whereas Monthong has a longer post-harvest ripening period that requires five days after harvesting at the mature stage. Given that durian is a climacteric fruit, its shelf life is restricted once the ripening process starts. This phenomenon could negatively affect the handling and transportation of durian and cause major economic losses. Therefore, gaining a deeper understanding of the ripening process in durian is of the utmost importance to control its post-harvest ripening. Despite the rising importance of durian as an economic fruit crop, genomic research on this species is almost non-existent, especially with regard to its Dof TFs. Our study aimed to conduct a genome-wide analysis to identify the Dof TFs in the durian genome, which has recently been published 16 . Further, we performed some routine in silico analyses, including gene structure analysis, identification of conserved motifs, and phylogenetic analysis. Finally, we analysed the expression patterns of durian Dofs (DzDofs) during post-harvest ripening in both the Monthong and Chanee cultivars and identified ripening-associated and cultivar-dependent DzDofs. The findings of our study deepen our knowledge of the role of Dof TFs during post-harvest ripening of climacteric durian fruits, and how differential expression patterns of some Dofs could contribute to differences in the ripening behaviours of the Monthong and Chanee cultivars.

Identification of 24 DzDof transcription factors.
Our genome-wide attempts to identify DzDofs revealed a total of 24 genes encoding transcription factors of the DzDof family in Musang King genome, which was equal to the number of previously annotated DzDofs in the durian genome 16 (Table 1). The sizes of the fulllength coding sequences of DzDof genes ranged from 519 to 1560 bp. The deduced amino acid lengths of proteins made from DzDofs varied from 172 to 519 amino acids, with the isoelectric point (pI) of the resultant proteins predicted to range from 4.88 to 9.65 and their molecular masses from 19.19 to 56.59 kD (Table 1). Notably, DzDof2.4 and DzDof5.1 showed the highest levels of sequence similarity (89.85%), whereas DzDof1.2 and durian cyclic Dof1 (DzcDof1) were the least similar (22.82%) (Supplementary Table S1). Multiple sequence alignment of DzDofs revealed a highly conserved DNA binding domain of 50 amino acids at the N-terminus, structured as a Cys2/Cys2 Zn 2+ finger (Fig. 1). This conserved motif, designated as the Dof domain, is the defining characteristic of Dof TFs. Moreover, a nuclear localization signal (NLS) was also detected in all DzDofs by in silico analyses. This sequence, which partly overlapped with the Dof domain, indicated the point of the nuclear localization of DzDofs.
Phylogenetic analysis of DzDofs. To investigate the evolutionary relationships among DzDof TFs and those from eight other representative species (dicots: Arabidopsis thaliana, Solanum lycopersicum, Capsicum annuum, and Vitis vinifera; monocots: Oryza sativa and Sorghum bicolor; non-vascular species: Chlamydomonas reinhardtii and Physcomitrella patens subsp. patens), a neighbour-joining (NJ) phylogenetic tree was constructed by aligning their full-length protein sequences, with 1000 bootstrap replicates. The phylogenetic tree divided the DzDof homologues into nine groups (A to I), which could be further classified into 14 subgroups: A1, A2, B, C, D1, D2, D3, D4, E1, E2, F, G, H, and I (Fig. 2). Group A was the largest clade, containing 56 members, while group B was the smallest clade and harboured only seven members. Interestingly, most of the Dofs from O. sativa and S. bicolor were placed in group F, highlighting the fact that these two species belonged to the grass family. Moreover, the only Dof TF from C. reinhardtii was classified into group F. Many Dofs from S. lycopersicum and C. annuum were clustered as pairs, which might point to some levels of domain consensus before the divergence of tomato and pepper. Fourteen pairs of DzDofs with V. vinifera Dofs (VvDofs) were shown while only one pair of DzDof with Arabidopsis Dof (AtDof) was found. All durian cyclic Dofs (DzcDofs) were grouped together in subclade A1.

Conserved motifs and gene structure analysis. Gene structure analysis revealed no introns in all
DzDofs ( Supplementary Fig. S1). This structural feature (having few (0-2) introns) might be correlated with the important functions of Dof, as similar findings were previously reported for Arabidopsis and rice Dofs 9 . Moreover, to gain better insights into the conserved motifs of DzDofs, the MEME program was used to analyse putative motifs. As shown in Supplementary Fig. S2, ten distinct motifs were identified. Motif 1 (50 amino acids), known as the Dof domain, was detected at the N-terminus of all DzDofs. All members clustered into one group showed similar motif compositions, which suggested functional similarities among the members of the same group ( Supplementary Fig. S3).  7) were not expressed in the fruit pulp at ripe stage. However, their expressions could be detected in the other tissues ( Supplementary Fig. S4).

Differential expression patterns of DzDofs during post-harvest ripening of the Monthong cultivar.
We examined the expression levels of 24 DzDofs at three different stages (unripe, midripe, and ripe) during post-harvest ripening. Using qRT-PCR, we observed that 15 of the 24 DzDofs were expressed in the fruit pulp, which were: DzcDof1, DzcDof2, DzcDof3, DzDof1. 4 6, were expressed at decreasing levels over the course of fruit ripening. In addition, DzDof5.3 was expressed at significantly higher levels at the midripe compared to the unripe stage. However, its expression level declined from the midripe to the ripe stage. Taken together, these 11 DzDofs showed ripening-associated expression patterns, and their expression levels were further examined under ethylene treatment. In contrast, the expression of DzcDof3, DzDof2.2, DzDof2.4, and DzDof5.1 did not vary significantly during post-harvest ripening, and thus these were not considered to be ripening-associated TFs.
Validation of ripening-associated DzDofs under three ripening treatments of the Monthong cultivar. Ethylene plays a major role in climacteric fruit ripening. To further confirm the expression levels of the putative ripening-associated DzDofs under ethylene treatment, their transcript levels were measured under three different ripening conditions: natural, ethephon-induced, and 1-methylcyclopropene (1-MCP)-delayed ripening. As shown in Fig. 4 ethylene treatment, whereas 1-MCP treatment significantly enhanced their expression levels relative to those during natural ripening. These results provide convincing evidence regarding the ripening-associated roles of these seven DzDofs during fruit ripening. However, the expression of DzcDof1, DzDof3.1, DzDof4.6, and DzDof5.3 did not vary significantly under different ripening treatments, and therefore these were not considered to be ripening-associated TFs.  6, showed significant differences between the two cultivars in at least one stage of ripening ( Supplementary Fig. S5). These ten DzDofs were thus considered potential cultivar-dependent TFs. However, the expression levels of the remaining DzDofs did not show any significant differences between the two cultivars during ripening. Among the cultivar-dependent DzDofs, DzDof5.4 had the greatest fold increase at the ripe stage in Chanee compared with that in Monthong. Moreover, DzDof2.2 was the only Dof whose expression level was higher in Chanee, and it showed a significant fold difference between the two cultivars at each ripening stage. This DzDof might be responsible for the different ripening characteristic of Chanee and Monthong, since its expression levels were significantly higher at each ripening stage in the Chanee cultivar than those in the Monthong cultivar. Hence, DzDof2.2 was further selected for functional prediction analysis via orthologue identification.

Differential expression levels of fruit pulp-expressed
Orthologue identification and gene function prediction of the candidate cultivar-dependent DzDof2.2. Characterized orthologue identification is widely used as a powerful tool for functional prediction in genes [17][18][19] . It has been shown that orthologous proteins have similar biological functions in different species 17 . The protein sequence of DzDof2.2 was used as a query to find orthologous proteins in Arabidopsis. As a result, AtDof2.4, AtDof5.1, AtDof3.6 and AtDof1.1 (OBP2) were found to be the orthologues of DzDof2.2. This was consistent with the phylogenetic tree ( Fig. 2), which clustered DzDof2.2 with the mentioned AtDofs together in the clade F. According to the functional characterization of AtDof1.1 13 , significantly higher auxin levels were measured in the leaves and seedlings of AtDof1.1 over-expressing lines compared to those in wild-type plants. www.nature.com/scientificreports www.nature.com/scientificreports/ Accordingly, it might be possible that DzDof2.2 could also regulate auxin biosynthesis during ripening in durian. This finding prompted us to further measure endogenous auxin levels in the fruit pulps of both cultivars at different stages of ripening.

Higher indole-3-acetic acid (IAA) and auxin biosynthetic gene expression levels during the ripening of Chanee compared to Monthong.
Since the cultivar-dependent DzDof2.2 could play a role in ripening by regulating auxin biosynthesis, endogenous levels of IAA were measured in the pulps of durian fruits at five different stages before (immature and mature) and during post-harvest ripening (unripe, midripe, and ripe) of both the Monthong and Chanee cultivars. As shown in Fig. 5A, during the immature and mature stages, IAA was detected at low levels, which were not significantly different between the two cultivars. However, after being harvested (from the mature to unripe stages), IAA levels significantly increased in both cultivars, with there being a significantly greater increase in Chanee than in Monthong. Moreover, IAA levels were significantly higher during the post-harvest ripening of Chanee compared to those in Monthong. The higher IAA levels in Chanee could contribute to its faster post-harvest ripening compared to that of Monthong through enhancing climacteric ethylene biosynthesis (auxin-ethylene crosstalk). Assessing transcript levels by qRT-PCR confirmed the higher expression levels of both ethylene biosynthetic genes, e.g. ACC synthase (ACS) and ACC oxidase (ACO), during ripening of Chanee compared to Monthong ( Supplementary Fig. S6), suggesting a higher endogenous ethylene content in Chanee. Moreover, the expression levels of both genes peaked at the ripe stage, indicating the highest levels of climacteric ethylene biosynthesis at this stage during post-harvest ripening of durian fruit ( Supplementary Fig. S6).
To assess the correlation between DzDof2.2 expression levels and IAA content, we analysed the expression levels of DzDof2.2 before ripening at the immature and mature stages, in addition to the previous results assessed during post-harvest ripening ( Supplementary Fig. S5), in both cultivars. At these two stages, it was found that DzDof2.2 was expressed at low levels in both the Monthong and Chanee cultivars, without any significant difference between the two cultivars (Fig. 5B). However, there was a significant increase in the expression levels of www.nature.com/scientificreports www.nature.com/scientificreports/ DzDof2.2 from the mature to the unripe stage in both cultivars, with a significantly greater fold change in Chanee compared to that in Monthong (Fig. 5B). The consistency of the expression levels of DzDof2.2 with auxin levels before and during ripening could support the possible role of DzDof2.2 in controlling auxin biosynthesis.
From the analysis of our in-house transcriptome data of Chanee fruit pulp, we identified two candidate auxin biosynthetic genes, including L-tryptophan aminotransferase 1 (TAA1) and the indole-3-pyruvate monooxygenase (YUCCA4) of the indole-3-pyruvic acid (IPyA) pathway. Interestingly, TAA1 and YUCCA4 had low expression levels during the immature and mature stages (before post-harvest ripening), with no significant differences between the two cultivars ( Fig. 6). However, there was a dramatic increase in the expression levels of both genes after harvesting (from the mature to the unripe stage) in both the Monthong and Chanee cultivars. This phenomenon was consistent with the increased auxin levels observed after harvesting, which could occur as a stress response to mechanical injury during fruit harvesting. Notably, TAA1 and YUCCA4 were expressed at significantly higher levels during the post-harvest ripening of Chanee compared to that of the Monthong cultivar (Fig. 6). This observation was consistent with the higher IAA levels observed during ripening in the Chanee than in the Monthong cultivar. Interestingly, we scanned the promoter regions of TAA1 and YUCCA4 and identified 40 and 34 cis-regulatory motifs (AAAG or its reverse-oriented sequence CTTT, known as the binding site for Dof TFs) 1,2 , respectively ( Supplementary Fig. S7). This observation could pinpoint the possible role of Dof TFs in the transcriptional regulation of auxin biosynthetic genes.

Functional characterization of DzDof2.2 through transient expression in Nicotiana benthamiana.
To assess the function of DzDof2.2 in planta, it was transiently expressed in the leaves of N. benthamiana and the expression levels of major auxin biosynthetic genes from N. benthamiana, NbTAA1 and NbYUCCA3, were analysed. We observed a significantly higher expression levels of both NbTAA1 and NbYUCCA3 in the leaves transiently expressing DzDof2.2 when compared to the control (Fig. 7). This observation confirmed the transcriptional regulation of auxin biosynthetic genes by DzDof2.2 TF. Consistently, our in silico analysis identified cis-regulatory binding sites (AAAG/CTTT) for Dof TFs in the promoter regions of both NbTAA1 and NbYUCCA3 (Supplementary Fig. S8).   Comparisons were made between the two cultivars at each stage. An asterisk ( * ) above the bars indicates a significant difference between cultivars (Student's t-test, P < 0.05). (B) Expression levels of DzDof2.2 were measured at five different stages (immature, mature, unripe, midripe, and ripe) in both the Monthong and Chanee cultivars. Comparisons were made between the two cultivars at each stage. Three independent biological replicates were used. Three technical replicates were also performed for each biological replicate. An asterisk ( * ) above the bars indicates a significant difference between cultivars (Student's t-test, P < 0.05).
(76 BraDof) 21 , cucumber (36 CsDof) 9 , and potato (35 StDof) 4 . Consistent with previous studies, the alignment of DzDof protein sequences revealed a highly conserved Dof domain, which was observed at the N-terminal regions of all 24 DzDofs (Fig. 1), suggesting that Dof TFs have been evolutionarily conserved among plants. Apart from the Dof domain, nine other distinct motifs were identified that were differentially distributed among DzDofs. In addition, as predicted by in silico analyses, all of the 24 deduced DzDofs harboured NLSs to help localize them to the nucleus. Digital expression profiles of DzDofs from Musang King cultivar revealed a total of 15 ripe fruit pulp-expressed Dofs (Supplementary Fig. S4). These 15 DzDofs were identical to those expressed during post-harvest ripening of both Monthong and Chanee cultivars, suggesting their functions in the pulp tissue. Dof TFs have been shown to play crucial roles in various physiological phenomena, such as the regulatory networks of plant responses to diverse biotic and abiotic stressors, seed germination, morphology, light responses, and fruit ripening 3,11,22 .
In our study, among the 15 fruit pulp-expressed Dofs in Monthong cultivar, seven ripening-associated DzDofs, including DzcDof2, DzDof1.4, DzDof1.5, DzDof2.1, DzDof2.5, DzDof5.4, and DzDof5.6, were identified based on their differential expression during post-harvest ripening in the Monthong cultivar (Fig. 3). The expression levels of DzcDof2, DzDof2.1, DzDof5.4, and DzDof5.6 were confirmed to be induced under ethephon and suppressed under 1-MCP treatment (Fig. 4). These ethylene-inducible DzDofs were expressed at increasing levels during ripening, and thus could act as either putative transcriptional activators or repressors of ripening. Zhang et al. 12 identified an ethylene-induced Dof TF in kiwifruit (AdDof3) that acted as a ripening-activator TF and enhanced starch degradation during fruit ripening. On the other hand, an ethylene-induced banana Dof (MaDof23), which was expressed at increasing levels during ripening, was a transcriptional repressor. This Dof acted antagonistically to a ripening-activator ethylene response factor (MaERF9) to regulate banana fruit Figure 6. Expression levels of auxin biosynthetic genes. Expression levels of two candidate auxin biosynthetic genes, L-tryptophan aminotransferase 1 (TAA1) (left panel) and the indole-3-pyruvate monooxygenase (YUCCA4) (right panel) of the indole-3-pyruvic acid (IPyA) pathway, were measured at five different stages (immature, mature, unripe, midripe, and ripe) in both the Monthong and Chanee cultivars. Comparisons were made between the two cultivars at each stage. Three independent biological replicates were used. An asterisk ( * ) above the bars indicates a significant difference between cultivars (Student's t-test, P < 0.05).

Figure 7. Expression levels of auxin biosynthetic genes in Nicotiana benthamiana transiently expressing
DzDof2.2. Expression levels of two candidate auxin biosynthetic genes from N. benthamiana, L-tryptophan aminotransferase 1 (NbTAA1) and the indole-3-pyruvate monooxygenase (NbYUCCA3) of the indole-3pyruvic acid (IPyA) pathway, were measured in N. benthamiana leaves infiltrated with pGWB2-DzDof2.2 (treatment) and pGWB2-GFP (control). Comparisons were made between the control and treatment. Three independent biological replicates were used. An asterisk ( * ) above the bars indicates a significant difference (Student's t-test, P < 0.05). www.nature.com/scientificreports www.nature.com/scientificreports/ ripening 11 . Notably, the expression levels of the remaining ripening-associated DzDofs (DzDof1.4, DzDof1.5, and DzDof2.5) were dramatically suppressed under ethephon and induced under 1-MCP treatment. It is interesting that these ripening-associated DzDofs were expressed at decreasing levels during ripening. To the best of our knowledge, this is the first report describing this particular effect of 1-MCP treatment on Dof TFs. The 1-MCP treatment delays climacteric ripening by inhibiting ethylene-dependent signalling by competitively binding to ethylene receptors 23,24 . Besides this inhibitory effect, 1-MCP also downregulates the expression of ethylene biosynthetic genes [25][26][27] (e.g., 1-aminocyclopropane-carboxylase (ACC) synthase (ACS) and ACC oxidase (ACO)), ethylene perception genes 27-30 (e.g., ethylene receptors ETR and ERS), and ethylene signalling genes 31,32 (e.g., ethylene response factor (ERF)). It is possible that 1-MCP upregulates the expression of these three DzDofs by inhibiting ethylene biosynthesis and signalling. This observation suggests that these DzDofs could be negatively regulated by ethylene, and act as putative transcriptional repressors of ripening. Consistent with our observations, the expression of a kiwifruit Dof (AdDof4) was suppressed by ethylene treatment and was thus considered to be a ripening-repressor TF 12 .
In addition, we compared the expression patterns of DzDofs between the two cultivars and found that DzDofs in the Chanee cultivar showed similar expression patterns, including identical ripening-associated DzDofs, to those in the Monthong cultivar. The similar expression patterns of DzDofs between Monthong and Chanee cultivars suggest that these DzDofs might have similar regulatory roles during fruit ripening in these two cultivars. However, although the expression patterns were similar between the two cultivars, comparisons of the expression levels of fruit pulp-expressed DzDofs between cultivars revealed a total of 10 potential cultivar-dependent Dofs. Among these, the expression levels of DzDof2.2 underwent a significantly greater fold change at each ripening stage in the Chanee than in the Monthong cultivar. Although the expression level of DzDof2.2 did not vary significantly during the unripe, midripe, and ripe stages, its expression level harboured a significant increase from mature (before ripening) to unripe stage (post-harvest ripening) in both cultivars with a greater fold increase in Chanee than in Monthong (Fig. 5B). Therefore, we hypothesized that DzDof2.2 would play a role during ripening of durian fruit and the difference in its expression levels between the two cultivars might contribute to the faster ripening and different sensory characteristics of the Chanee cultivar from those of the Monthong one.
To examine the possible function of DzDof2.2, we found its orthologues in Arabidopsis. Among the four identified orthologues, AtDof1.1 (OBP2) is involved in regulating glucosinolate biosynthesis. In addition, transgenic Arabidopsis lines over-expressing OBP2 had higher auxin concentrations in their leaves and seedlings compared to those in the wild type 13 . This phenomenon could suggest a role of DzDof2.2 in regulating auxin biosynthesis. Higher expression levels of DzDof2.2 could enhance auxin biosynthesis, leading to the higher auxin levels during ripening in the Chanee than in the Monthong cultivar (Fig. 5A). Notably, the increased expression levels of DzDof2.2 were consistent with the increased auxin levels from the mature to the unripe stage (Fig. 5). We further confirmed that the higher level of auxin in Chanee co-occurred with the increased expression of auxin biosynthetic genes. Four tryptophan-dependent auxin biosynthetic pathways have been well-characterized in Arabidopsis, including the indole-3-pyruvic acid (IPyA), tryptamine (TAM), indole-3-acetaldoxime (IAOx), and indole-3-acetamide (IAM) pathways 33 . Among these pathways, IPyA is considered the major auxin biosynthetic pathway in many plant species 34 . The main biosynthetic enzymes of the IPyA pathway in durian, DzTAA1 and DzYUCCA4, were expressed at significantly higher levels during post-harvest ripening of the Chanee cultivar compared to those in the Monthong cultivar. Moreover, there was a significant increase in the expression levels of both genes from the mature to the unripe stage, with a higher fold increase in Chanee compared to that in Monthong (Fig. 6). This expression profile was consistent with the changes in expression levels of DzDof2.2 and IAA levels. It has been well-documented that auxin plays a major role during the ripening of climacteric fruits, such as peach, tomato, and papaya 35,36 . Concomitant with the increase in climacteric ethylene biosynthesis during ripening of tomato and peach, elevated levels of auxin have also been measured 35 . Auxin can stimulate climacteric ethylene biosynthesis through upregulating the ACC synthase gene family (auxin-ethylene crosstalk) [37][38][39] . Consistently, we observed significantly higher expression levels of both ethylene biosynthetic genes (ACS and ACO) during post-harvest ripening of Chanee compared to Monthong, suggesting higher ethylene content in Chanee than in Monthong. It has been well documented that the expression levels of ethylene biosynthetic genes are positively correlated with the ethylene content during climacteric ripening 40,41 . Besides this well-known role of auxin in upregulating ethylene biosynthesis, auxin can also play a role of its own during climacteric ripening by transcriptionally modulating other ripening-related auxin-responsive genes 35 . Taken together, our results demonstrate that auxin is actively involved in durian fruit ripening.
Interestingly, the promoter regions of DzTAA1 and DzYUCCA4 were enriched with cis-regulatory binding sites (AAAG/CTTT) for Dof TFs. This observation could suggest the possible transcriptional regulation of these two genes by Dof TFs. Our results obtained from transient expression of DzDof2.2 in N. benthamiana leaves provided convincing evidence of the possible role of DzDof2.2 in transcriptional regulation of auxin biosynthetic genes during ripening in durian (Fig. 7). Consistent with the promoters of auxin biosynthetic genes from durian, we identified the occurrences of Dof TF binding sites in the promoter regions of NbTAA1 and NbYUCCA3 (Supplementary Fig. S8).
It is notable that, among the 10 identified cultivar-dependent DzDofs, DzDof5.4 had the greatest fold increase at the ripe stage in Chanee relative to that in Monthong (Supplementary Fig. S5). OBP4 (AtDof5.4) was the closest orthologue of DzDof5.4. Interestingly, exogenous IAA induced the expression of AtDof5.4. This might explain why the expression levels of DzDof5.4 peaked at the ripe stage in both cultivars, with a dramatically greater fold increase in Chanee compared to that in Monthong. The higher IAA levels during ripening, the greater expression levels of DzDof5.4 were thus observed in our study.
To summarize, here, we reported a comprehensive analysis of durian Dof TFs. Out of 24 identified DzDofs, 15 were fruit pulp-expressed, among which seven DzDofs were identified as being ripening-associated for both the Chanee and Monthong cultivars. Notably, comparing the expression patterns of fruit pulp-expressed DzDofs (2019) 9:12109 | https://doi.org/10.1038/s41598-019-48601-7 www.nature.com/scientificreports www.nature.com/scientificreports/ between the Monthong and Chanee cultivars revealed 10 potential cultivar-dependent TFs, among which DzDof2.2 showed a significantly greater fold increase at every ripening stage in Chanee than in Monthong. We suggest that DzDof2.2 could transcriptionally control the expression of auxin biosynthetic genes. Higher expression levels of DzDof2.2 in Chanee could enhance its auxin levels during ripening through transcriptional regulation of auxin biosynthetic genes. Higher auxin levels in Chanee would upregulate ethylene biosynthesis by transcriptionally activating the ACC synthase gene family, and lead to an earlier ethylene response (auxin-ethylene crosstalk), and thus faster post-harvest ripening compared to that in the Monthong cultivar. Our results provide compelling evidence of the role of Dof TFs during the post-harvest ripening of durian. The results obtained herein should provide an important foundation for future studies involving gene cloning and the functional characterization of ripening-associated and cultivar-dependent DzDofs.

Materials and Methods
Plant materials and treatments. Durian (Durio zibethinus Murr.) cv. Monthong and Chanee fruits were collected from a commercial orchard located in Trat province in the eastern part of Thailand in early April 2017. Fruit samples of similar size and weight (~3-4 kg each) were harvested at two different stages: immature (at 70 days (for Chanee) and 85 days (for Monthong) after anthesis) and mature (at 90 days (for Chanee) and 105 days (for Monthong) after anthesis). Immature samples were peeled immediately after harvesting, whereas some mature samples were peeled immediately, and some were kept at room temperature (30 °C) until peeling. Five types of samples were used in this study: immature, mature, unripe (harvested at the mature stage and kept at room temperature for one day (both cultivars) and then peeled), midripe (harvested at the mature stage and kept at room temperature for two days (Chanee) or three days (Monthong) and then peeled), and ripe (harvested at the mature stage and kept at room temperature for four days (Chanee) or five days (Monthong) and then peeled). After peeling, two of the central pulps were collected from each fruit sample, following the method described by Pinsorn et al. 15 . Briefly, the first pulp was collected along with a seed, and was used to measure fruit firmness with a texture analyser (TA-XT2i; Stable Micro Systems, Godalming, UK) to assure that samples of the two cultivars were compared at the same ripening stages. A puncture test with a 6-mm probe, at a test speed of 2 mm/s and testing distance of 5 mm, was carried out on five random points in each pulp. Midripe and ripe pulps had a mean ± SD firmness of 3.4 ± 0.81 and 1.55 ± 0.45 N, respectively, in both cultivars. After this test, the second pulp was collected without a seed, immediately frozen in liquid nitrogen, and stored at −80 °C for further use.
To further investigate the role of ethylene during post-harvest ripening, three different ripening treatments were applied as follows: natural, ethephon-induced, and 1-methylcyclopropene (1-MCP)-delayed ripening. Mature durian samples were collected and treated with either ethephon (48% 2-chloroethylphosphonic acid; Alpha Agro Tech Co., Ltd., Thailand) or 1-MCP (0.19% 1-MCP tablet; BioLene Co., Ltd., China) for ethephon-induced and 1-MCP-delayed ripening, respectively. For ethephon treatment, the ethephon solution was applied to the upper area of each fruit stalk. For 1-MCP treatment, each fruit was placed inside a closed 20-L chamber. Then, one tablet of 1-MCP was placed into a beaker inside the chamber. Water (5 mL) was added to the beaker to generate gaseous 1-MCP and the chamber was immediately closed for 12 h at room temperature (30 °C) while the control samples were kept under similar conditions without 1-MCP. After treatment, control and treated samples were kept at room temperature (30 °C) until the ethephon-treated samples ripened. All samples from the three ripening conditions were then peeled. After that, two central pulps were collected from each sample and processed as mentioned previously. In this study, for each type of sample (immature, mature, unripe, midripe, and ripe, as well as natural ripening, ethephon-induced ripening, and 1-MCP-delayed ripening), three biological replicates were used. Each biological replicate was defined as one durian fruit harvested from a separate tree.
For the agroinfiltration experiment, N. benthamiana seeds were sown in pots containing peat moss and were grown under controlled conditions (temperature 25 °C and 16/8 h light/dark photoperiod; artificial light of 4,500 Lux). Two-week-old plants were transplanted individually into pots and were grown under similar conditions.
In silico identification and analysis of durian Dofs. The whole-genome sequence of durian cv. 'Musang King' was used to identify putative durian Dofs. First, the conserved Dof domain (PF02701) determined based on a Hidden Markov Model (HMM) was obtained from the Pfam protein database (http://pfam.xfam.org). The HMM profile of the Dof domain was then used as a query to search against in the Musang King genome database 16 . The amino acid sequences of DzDofs were further obtained and checked for the presence of the highly conserved motif (Dof domain) at the N-terminus. The isoelectric point (pI) and molecular weight of each of the identified DzDof proteins were then calculated using the web-based tool ExPASy (http://web.expasy.org/ compute_pi/) 42 . Nuclear localization signal (NLS) analysis was also done using a web-based tool, found at: http:// cello.life.nctu.edu.tw 9 .
Phylogenetic analysis of DzDofs. To investigate phylogenetic relationships among DzDofs and Dofs from eight other representative species (dicots: Arabidopsis thaliana, Solanum lycopersicum, Capsicum annuum, and Vitis vinifera; monocots: Oryza sativa and Sorghum bicolor; non-vascular species: Chlamydomonas reinhardtii and Physcomitrella patens subsp. patens), a neighbour-joining (NJ) phylogenetic tree was constructed by aligning their full-length protein sequences with 1000 bootstrap replicates using MEGA X software (with a JTT model and pairwise gap deletion) 43 .

Gene structure analysis and identification of conserved motifs. To obtain a deeper understanding of
DzDof structures, Gene Structure Display Server 2.0 (with default parameters) (http://gsds.cbi.pku.edu.cn) was used to characterize the exon-intron structures of DzDof-family genes. Motifs were identified using the MEME program (http://meme-suite.org) 44 , with the following parameter settings: motif length = 6-100; motif sites = 2-120; maximum number of motifs = 10; and the distribution of one single motif was "any number of repetitions". (2019) 9:12109 | https://doi.org/10.1038/s41598-019-48601-7 www.nature.com/scientificreports www.nature.com/scientificreports/ Tissue-specific expression analysis. In order to analyse the expression levels of DzDofs in different tissues, Illumina reads were obtained from a public repository database (SRA, Sequence Read Archive) as follows: SRX3188225 (root tissue), SRX3188222 (stem tissue), SRX3188226 (leaf tissue), and SRX3188223 (aril/pulp tissue). FASTQC software was used to analyse the total reads quality parameters (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). To filter low quality reads, the trimming process was performed using Trimmomatic V0.32 (http://www.usadellab.org/cms/?page=trimmomatic) (https://github.com/seqan/flexbar). Thereafter, transcriptome assembly was performed using Trinityrnaseq-Trinity-v2.4.0 (https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity) 45 . For expression profiling analysis, abundance estimation tool in Trinity package was used to align the input reads to the de novo assembled transcriptome to obtain raw counts of each contig. The raw read counts were then merged in a single read count matrix. This matrix was used as input for normalization to generate a Trimmed Mean of M-values (TMM) normalized matrix. Thereafter, the normalized total read counts (RC) were used to generate heatmap using MetaboAnalyst 4.0, an open source R-based program 46 . For the heatmap construction, the values were sum normalized, log2 transformed and autoscaled.
RNA isolation and quantitative real time polymerase chain reaction (qRT-PCR). Total RNA was isolated from durian fruit pulp and purified using PureLink Plant RNA Reagent (Thermo Fisher Scientific ™ ) according to the manufacturer's instructions. Samples were treated with DNase I (Thermo Fisher Scientific ™ ) to remove genomic DNA. The quality and quantity of RNA samples were assessed using agarose gel electrophoresis and an Eppendorf BioPhotometer D30 with A260/280 and A260/230 ratios between 1.8 to 2.0 and 2.0 to 2.2, respectively. One microgram of total RNA was reverse-transcribed to cDNA using a RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific ™ ), in accordance with the manufacturer's recommended protocol.
All primers used in this study were designed using Primer 3 online (http://primer3.ut.ee/). The primers used are presented in Supplementary Table S2. qRT-PCR was performed to investigate patterns in the expression of DzDofs at the immature and mature stages (before ripening) and unripe, midripe and ripe stages (post-harvest ripening) of fruits of both the Monthong and Chanee cultivars. The PCR reaction was performed in a total volume of 20 μL containing 1 μL of diluted cDNA (1 ng of cDNA), 10 μL of 2x QPCR Green Master Mix LRox (biotechrabbit, Berlin, Germany), and 200 nM of each gene-specific primer. Amplification was carried out with a CFX95 Real-time System (Bio-Rad Laboratories Inc., California, USA) under the following conditions: initial activation at 95 °C for 3 min, followed by 40 cycles of denaturation at 95 °C for 15 s, then annealing at 56-60 °C for 30 s, and finally extension at 72 °C for 20 s. A melting curve analysis for testing specific products was carried out by heating the products from 55 °C to 95 °C at increments of 0.5 °C. A single product was obtained for each DzDof. All qRT-PCR experiments were performed with three independent biological replicates. The elongation factor 1 alpha (EF-1α) gene was used as a reference gene (See Supplementary Table S3 for primers). The relative expression level of each gene compared with that of the reference gene was calculated using the 2 -ΔCt method 47 .

Identification of orthologues.
Orthologues of candidate DzDofs in Arabidopsis were identified using the web-based tool EnsemblPlants, release 40 (https://plants.ensembl.org/index.html), using the default parameters 48 .
Indole-3-acetic acid (IAA) measurement. Endogenous IAA levels in fruit pulps from both cultivars were measured. IAA was extracted from fruit pulp samples following the method described by Sinha and Basu 49 . Briefly, each pulp sample was freeze-ground into a fine powder. After this, 500 mg of the powder was mixed with 10 mL of chilled 80% ethanol, vortexed vigorously for 3 min, and kept at 4 °C for 2 h with occasional vortexing (at least five times). Afterwards, the mixture was centrifuged at 10,000 rpm (4 °C) for 15 min, and the supernatant was then used for IAA measurements. The IAA concentration (A 535 ) was quantitatively measured via the colorimetric assay of Gordon and Weber 50 using a standard curve prepared from serial dilutions of commercial IAA stock solutions (Duchefa Biochemie, The Netherlands), ranging from 7.5 to 500 μM and with an R 2 of 0.99.
In silico identification and promoter analysis of auxin biosynthetic genes. Based on the auxin biosynthetic pathways identified in Arabidopsis, the indole-3-pyruvic pathway was selected as the candidate pathway for auxin biosynthesis 34 in durian. The conserved domain of each biosynthetic enzyme (based on a HMM) was obtained from the Pfam protein database (http://pfam.xfam.org/) as follows: using the PF00155 for aminotransferase (TAA1) and the PF00743 for flavin-containing monooxygenases of the YUCCA family (YUCCA4). The HMM profile of each enzyme was then used as a query to search against our in-house transcriptome database and the Musang King genome (i.e., TAA1 (XM_022878297.1) and YUCCA4 (XM_022900772.1)). For qRT-PCR, primers were designed as mentioned previously (Supplementary Table S3). In order to analyse cis-regulatory elements in the promoter regions of auxin biosynthetic genes, the 2000-bp upstream promoter regions of these genes were scanned for the Dof domain (AAAG/CTTT) using the online tool PLACE (http://www.dna.affrc.go.jp/) 36 .

expression analysis of ACC synthase (ACS) and ACC oxidase (ACO) in Monthong and Chanee cultivars.
In order to compare the expression levels of ethylene biosynthetic genes during post-harvest ripening of Monthong and Chanee, the HMM profile of ACS and ACO (PF00155 from the Pfam protein database (http:// pfam.xfam.org/)) was used as a query to search against our in-house transcriptome database and the Musang King genome (i.e., ACS (XM_022901720.1) and ACO (XM_022903266.1). To perform qRT-PCR, primers were designed as mentioned previously (see primers in Supplementary Table S3).  Supplementary Table S4). Then, the PCR amplicon was introduced into pDONR207 vector using Gateway ® BP Clonase ® II (Thermo Fisher Scientific, USA). The putative DzDof2.2 from pDONR207-DzDof2.2 was inserted in to the pGWB2 expression vector using Gateway ® LR Clonase ® II (Thermo