Genome-wide analysis of autophagy-related genes in Medicago truncatula highlights their roles in seed development and response to drought stress

Autophagy is a highly conserved process of degradation of cytoplasmic constituents in eukaryotes. It is involved in the growth and development of plants, as well as in biotic and abiotic stress response. Although autophagy-related (ATG) genes have been identified and characterized in many plant species, little is known about this process in Medicago truncatula. In this study, 39 ATGs were identified, and their gene structures and conserved domains were systematically characterized in M. truncatula. Many cis-elements, related to hormone and stress responsiveness, were identified in the promoters of MtATGs. Phylogenetic and interaction network analyses suggested that the function of MtATGs is evolutionarily conserved in Arabidopsis and M. truncatula. The expression of MtATGs, at varied levels, was detected in all examined tissues. In addition, most of the MtATGs were highly induced during seed development and drought stress, which indicates that autophagy plays an important role in seed development and responses to drought stress in M. truncatula. In conclusion, this study gives a comprehensive overview of MtATGs and provides important clues for further functional analysis of autophagy in M. truncatula.

Protein sequence alignment and analysis of the phylogenetic relationship. The phylogenetic analysis of MtATGs was performed using the MEGA7 software 28 . The amino acid sequences of MtATGs and AtATGs in different gene families were aligned independently using the ClustalW algorithm with the default parameters. An unrooted phylogenetic tree was constructed with the neighbor-joining statistical method, and the following parameters were used: uniform rates are used as rates among sites, gaps/missing data are treated as pairwise deletion, and the bootstrap analysis was performed with 1000 replicates to obtain a support value for each branch.
Identification of cis-elements. The 1.5 kb genomic DNA sequence upstream of the initiation codon of each MtATG was retrieved from the M. truncatula genome (MedtrA17_4.0). The assumed cis-elements of MtATGs were predicted using the PlantCARE web servers (http:// bioin forma tics. psb. ugent. be/ webto ols/ plant care/ html/) 29 .
Construction of the protein-protein interaction (PPI) network. The PPI networks were constructed using the STRING database (http:// www. string-db. org). A total of 39 MtATGs were selected as input, and the PPI network of the MtATGs was constructed with a medium confidence (0.4).
Analysis of the expression profiles using microarray data. The M. truncatula microarray data were downloaded from the MtGEA v3 database (https:// mtgea. noble. org/ v3/) 30 . Expression values were normalized using the z-score method, and plotted using GraphPad Prism 8. Plant materials and growth conditions. Medicago truncatula (cv. Jemalong A17) seeds were scarified with sulfuric acid, and vernalized on wetted filter paper at 4 °C for 7 days. Seedlings were grown in a greenhouse at 24 °C, 16-h light/8-h dark cycle, with humidity ranging from 60 to 80%. Different plant tissues (roots, stems, leaves, petioles, buds, flowers, and pods) were harvested from multiple plants. Material for the seed developmental was collected from pods at 5 different stages. For drought stress, 7-day-old seedlings were treated by withholding watering for 2 days. For mannitol treatment, 2-weeks-old seedlings were transferred to liquid 1/2 MS medium supplemented with 300 mM mannitol for additional 2 days. All plant samples were frozen immediately in liquid nitrogen after harvest and stored at − 80 °C until use. Plant material collections in this study complied with relevant institutional, national, and international guidelines and legislation.
RNA isolation and quantitative PCR (qPCR) analysis. Total  www.nature.com/scientificreports/ eight in different groups (three in the MtATG1 family, eight in the MtATG8 family, two in the MtATG9 family, three in the MtATG13 family, three in the MtATG16 family, and eight in the MtATG18 family) ( Table 1). The chromosomal distribution of MtATGs determined using the TBtools software is shown in Fig. 1. In total, 38 MtATGs were found to be distributed across all eight chromosomes except for MtATG7, which could not be mapped to any chromosome according to data from MedtrA17_4.0 (Fig. 1). The number of MtATGs located on each chromosome varies dramatically. Chromosome 4 (Chr4) contains the maximum number (11) of MtATGs, whereas chromosome 6 has only one MtATG gene. Gene duplication is important for adaptation of plants to adverse and complex environments. In M. truncatula, 7 pairs of MtATGs were predicted to be segmentally duplicated. As shown in Fig The subcellular localization of the MtATGs was predicted using the CELLO system (http:// cello. life. nctu. edu. tw). Most of the MtATGs were predicted to localize to the nucleus, plasma membrane, and cytoplasm, followed by extracellular space, chloroplast, and mitochondria (Table 1, Supplementary Figure S1). Furthermore, some MtATG families exhibited different subcellular localization. For example, MtATG8 proteins were predicted to be mainly cytoplasmic or nuclear, but were also found to localize to the mitochondria and extracellular space ( Table 1). The prediction was the same for MtATG18 family members, which were localized to both the plasma membrane and nucleus (Table 1). Taken together, the diverse subcellular localization of MtATGs implies that they have distinct functions.  Fig. 2, members of the MtATG1 and MtATG13 families were clustered in two branches ( Fig. 2A,B). There are two ATG9s and three ATG16s in M. truncatula, whereas only one ATG9 and ATG16 in Arabidopsis (Fig. 2C,D). ATG8 plays a central role in autophagy by promoting autophagosome formation and cargo recruitment. As in Arabidopsis, eight MtATG8 members were clustered into two distinct groups in M. truncatula: MtATG8a, MtATG8b, MtATG8c, MtATG8d, and MtATG8e were grouped into clade I, whereas MtATG8f, MtATG8g, and MtATG8h were clustered in clade II (Fig. 2E). MtATG8 proteins showed high identity with ATG8 proteins from Arabidopsis, except for MtATG8h, in which half of the amino acids from the N-terminus were absent (Supplementary Figure S2). The C-terminal glycine residue in ATG8, which is exposed upon protease cleavage by ATG4, is essential for the conjugation of ATG8 to phosphatidylethanolamine 33 . However, MtATG8b did not contain the C-terminal glycine residue. This result indicates that MtATG8b might function in other biological processes independent of autophagy. In addition, one MtATG8 member of clade II, MtATG8f, had a C-terminal extension after the Gly residue, whereas the AtATG8 members of clade II lack the C-terminal extension (Supplementary Figure S2). Eight MtATG18 members were also clustered in two branches like the MtATG8 family

Analyses of gene structures and distribution of conserved domains.
Gene structure is closely related to the expression pattern and function divergence of members of multigene families. Gene structure analysis revealed that all the MtATGs contain introns, with the number of exons ranging from 2 to 17 (Fig. 3A).
In addition, similar exon-intron patterns and the same number of exons were observed in some ATG subfamilies, such as MtATG1a/b, MtATG8a/c/d/e/f/g, MtATG13a/b/c, MtATG18a/c/d/e, and MtATG18g/h (Fig. 3A). The similar gene structures suggest functional redundancy among these genes. However, differences in exonintron patterns and exon numbers were also seen within some subfamilies, such as MtATG1t, MtATG8b/h, and MtATG18b/f (Fig. 3A). The conserved domains of MtATGs were detected using the Pfam database 25 . In general, the composition of the conserved domains in MtATGs is comparable to that in Arabidopsis. Furthermore, members of the same MtATG families have similar domains. For example, all three MtATG1 proteins contain a protein kinase domain (Pkinase) at their N-terminus (Fig. 3B). In addition, almost all MtATG8 proteins (except MtATG8h) are similar in length and have identical ATG8 domains (Fig. 3). A similar phenomenon was also observed in the MtATG9 and MtATG13 subfamilies. However, exceptions were also found in the MtATG16 and MtATG18 subfamilies. All the MtATG16 family proteins have a C-terminal WD40 domain, but lack an N-terminal ATG16 domain in MtATG16c (Fig. 3B). MtATG18 proteins contain the WD40 domain except for MtATG18b and MtATG18h, but members of clade II (MtATG18f/g/h) have a C-terminal BCAS3 domain that is absent in members of clade I (Fig. 3B). The differences in the gene structure and conserved domains may be related to functional divergence among the different gene products within some MtATG families.

Analysis of cis-elements in the promoter regions of MtATGs.
Cis-elements regulate genes through interactions with their corresponding transcription factors. To further understand the gene regulation network of MtATGs, cis-elements were identified using the online tool PlantCARE 29 . Ninety-two putative cis-elements www.nature.com/scientificreports/ were found among MtATG promoters (Supplementary Table S3). Among these, the TATA-box and CAAT-box are the most common cis-elements. Many of the identified cis-elements, such as ABRE (abscisic acid-related), TCA-element (salicylic acid-related), TCC ACC T-motif and TGACG-motif (MeJA-related), TGA-element (auxin-related), TATC-box, and P-box and GARE-motif (gibberellin-related), are involved in hormone responsiveness (Fig. 4). Among these, cis-elements that respond to MeJA and ABA were found to be the most abundant. In addition, some stress-related elements are mainly related to anaerobic (ARE), defense (STRE and TCrich repeats), drought (MBS), low temperature (LTR), and wound (WUN-motif) stresses (Fig. 4). The diversity of cis-elements in the promoter regions of MtATGs provided evidence for their potential biological functions in response to phytohormone, abiotic and biotic stresses.   www.nature.com/scientificreports/ MtATG18b) were inspected by qPCR. Most of the selected genes were highly expressed in roots, which was very similar to those of microarray analysis (Fig. 6B). Consistent with previous studies, most of the MtATGs were upregulated during seed development (Fig. 7A). In particular, MtATG2, MtATG3, MtATG4, MtATG5, MtATG6, MtATG13a, and MtATG18b, were highly expressed in the late stage of seed development (Fig. 7A). In contrast, a few MtATGs, including MtATG7 and MtATG8b, were downregulated after pollination (Fig. 7A). To validate the results of the microarray data, seeds were collected from pods at 5 different stages of seed developmental (Fig. 7B). As shown in Fig. 7C, the expression levels of five selected genes (MtATG2, MtATG4, MtATG5, MtATG8a, and MtATG18b) were considerably increased, only MtATG4 showed no gene expression change during seed development. These results were very similar to those of microarray analysis, and indicate that autophagy is essential for seed development in M. truncatula.

Analysis of the protein-protein interaction network of
Expression of MtATGs in response to drought stress. To investigate the putative roles of autophagy in the response of M. truncatula to drought stress, the expression profiles of MtATGs were analyzed using microarray data from the MtGEA database 35,36 . Generally, most MtATGs were upregulated after drought treatment (Fig. 8A). Specifically, 26 of 34 MtATGs (e.g., MtATG1t, MtATG8d, MtATG9a, and MtATG18b) were continuously upregulated when plants were subjected to drought stress by withholding watering, and the transcripts of MtATGs rapidly dropped to their basal levels after resuming the watering (Fig. 8A). Interestingly, MtATG8g showed an opposite trend: the expression level of MtATG8g dramatically decreased under drought stress compared with other MtATGs (Fig. 8A). To validate the results of the microarray data, six genes (MtATG1a, MtATG2, MtATG4, MtATG5, MtATG8a, and MtATG18b) were selected for independent validation by qPCR. The expression levels of most of the selected genes were significantly higher after 2 days of drought treatment (Fig. 8B). To examine autophagy activity under drought stress, antibodies against ATG8a were used to detect ATG8 protein www.nature.com/scientificreports/ by western blotting. ATG8 proteins are lipidated with phosphatidylethanolamine (PE) to promote autophagosome formation in response to drought treatment, whereas no changes in the level of ATG8-PE were detected under control condition (Fig. 8C). Furthermore, MDC staining showed that the number of autophagosomes was significantly increased after drought treatment (Fig. 8D). These results suggested that autophagy might play a crucial role in M. truncatula response to drought stress.

Discussion
In this study, 39 ATGs were identified in M. truncatula. These ATGs are similar to orthologous genes in Arabidopsis. For example, phylogenetic analysis revealed that ATG families in M. truncatula are very similar to those in Arabidopsis. In addition, the PPI network analysis shows that the interaction pattern of MtATGs is also similar to that of ATGs in Arabidopsis. These results indicate that the autophagy pathway is highly conserved across different plant species. However, the number of members in some ATG families differs among plant species. For example, the ATG8 family contains eight genes in M. truncatula, but nine in Arabidopsis, seven in rice, and thirteen in wheat 6,9,37 . In addition, the gene structure and conserved domains of some MtATG families, such as MtATG16 and MtATG18 subfamilies, also differ from those of other plants. Furthermore, different types of www.nature.com/scientificreports/ cis-elements were identified in the promoters of MtATGs in the same gene family. These results suggest that M. truncatula may have species-specific autophagy mechanism. Hence, it is necessary to illustrate the conserved and specific functions of MtATGs in future studies. Autophagy has been shown to play crucial roles in the growth and development of plants 4 . In this study, we found that all ATGs were expressed in the tested tissues of M. truncatula, but their expression levels varied among different tissues. The tissue-specific expression of MtATGs suggests that different functions are required in different tissues. Seed development consists of embryo morphogenesis and seed maturation 38 . In rice, autophagy has been shown to be involved in the regulation of starch and sugar metabolism during seed maturation 39 . In Norway spruce (Picea abies), autophagy is also involved in embryogenesis in which it regulates vacuolar cell death of the embryo suspensor 40 . Furthermore, autophagy plays an important role in microspore embryogenesis in Brassica napus 41 . The seed weight in autophagy-defective mutants of Arabidopsis and maize was reported to For drought stress treatment, soil-grown plants were subjected to drought stress by withholding watering (Drought) for 14 days, followed by rewatering. Scale bar represents the fold change (log2 value) relative to the corresponding control. (B) qRT-PCR validation of MtATGs expression under drought stress. For drought stress treatment, 7-day-old seedlings were subjected to drought stress by withholding watering for 2 days. MtACTIN was used as a reference gene. Error bars represent SD of three independent experiments. Significant differences were indicated with an asterisks (*), P < 0.05. (C) Analysis of ATG8 lipidation by western blot. Two week seedlings were transferred to liquid 1/2 MS medium with or without 300 mM mannitol, and whole seedlings were collected at 0, 1, and 2 day after treatment. The anti-ATG8a antibodies were used for immunoblotting. (D) MDC staining of root cells with or without drought treatment. Two-wk-old seedlings were transferred to liquid 1/2 MS medium with or without 300 mM mannitol for 2 days followed by staining with MDC. The labeled autophagosomes (arrows) were visualized by epifluorescence microscopy. Scale bar: 50 μm. www.nature.com/scientificreports/ be lower than in the wild-type plants 7,42 . In the present study, we found that most of the MtATGs were induced during seed development and were highly expressed at the late stage of seed development, which indicates that autophagy is necessary for seed development in M. truncatula. Overall, autophagy plays crucial roles in the growth and development of plants through a pathway that is conserved across different species. Autophagy has been demonstrated to promote plant survival by maintaining cellular homeostasis under drought stress 43,44 . In A. thaliana, the transcriptional level of ATG18a was rapidly upregulated by mannitol treatment 45 . In O. sativa, the expression levels of OsATG6 genes were also induced by drought stress 46 . Moreover, ATG genes were upregulated by drought stress in many other plant species, such as barley 47 , pepper 48 , apple 49 , and banana 50 . Besides changes in gene expression, the Arabidopsis autophagy-defective mutants (atg5, atg7, and RNAi-ATG18a) showed more sensitivity to drought treatment than the wild type 45,51 . Inhibition of autophagy by 3-MA or knockdown of ATG6 sensitized wheat seedlings to drought stress 52 . Furthermore, virus-induced gene silencing of ATG8d or ATG18h significantly reduced drought tolerance in tomato 53 . However, overexpression of MdATG5 or MdATG18a enhanced tolerance to drought stress in apple trees 54,55 . In addition, overexpression of SiATG8a from foxtail millet improved drought tolerance in Arabidopsis 56 . Recently, it was reported that autophagy improves drought tolerance in M. truncatula through degradation of the aquaporin MtPIP2;7, which interacts with the cargo receptor MtCAS31 56 . Consistent with previous studies, our results reveal that the promoter of many MtATGs contain the drought-related MBS cis-element. Furthermore, the transcriptional levels of most of the MtATG genes, especially those of the MtATG8 family, significantly increased after drought treatment. The lipidation of ATG8 protein and accumulation of autophagosome are enhanced in M. truncatula during drought stress. Our findings indicate that autophagy is is largely induced by drought stress in M. truncatula, and can be considered an adaptive response under drought stress.

Conclusion
This study provided comprehensive analysis of ATGs in M. truncatula. In total, 39 ATGs were identified in M. truncatula. Members of the same ATG family showed similar gene structure and conserved domains. Analysis of cis-elements implied that MtATGs have potential biological functions in response to phytohormone, abiotic and biotic stresses. Phylogenetic and interaction network analyses suggested that the function of MtATGs is evolutionarily conserved in Arabidopsis and M. truncatula. The expression pattern of MtATGs indicates that autophagy possibly participates in seed development and plays an important role in plant responses to drought stress. In conclusion, this study gives a detailed overview of MtATGs and their expression patterns. The results obtained in this study provide useful information for further functional characterization of autophagy in M. truncatula. www.nature.com/scientificreports/