Genome-wide Identi cation and Expression Analysis of Autophagy-related Genes (ATGs) in Medicago Truncatula

Background: Autophagy is a highly conserved degradation process of cytoplasmic constituents in eukaryotes. Autophagy is known to be involved in the regulation of plant growth and development, as well as biotic and abiotic stress response. Although autophagy-related genes (ATGs) have been identied and characterized in many plant species, little is known about the autophagy process in Medicago truncatula. Results: In this study, 39 ATGs were identi ﬁ ed in M. truncatula (MtATGs), and the gene structures and conserved domains of MtATGs were systematically characterized. In addition, many cis-elements which are related to hormone and stress responsiveness were identied in the promoters of MtATGs. Furthermore, phylogenetic analysis and interaction network analysis suggested that the function of MtATGs is evolutionarily conserved in Arabidopsis and M. truncatula. Gene expression analysis showed that most MtATGs were largely induced during seed development, but repressed by nodulation. Moreover, MtATGs were up-regulated in response to salt and drought stresses. Conclusion: These results provide a comprehensive overview of the MtATGs, which provided important clues for further functional analysis of autophagy in M. truncatula.


Genome-wide identi cation of ATGs in M. truncatula
To identify ATGs in M. truncatula, the BLASTP algorithm was employed in searching against M. truncatula proteome sequences (MedtrA17_4.0) using the amino acid sequences of Arabidopsis thaliana ATGs (AtATGs) as queries. A total of 39 MtATGs were identi ed in M. truncatula (Table S1) (Table S1).
The chromosomal distribution of MtATGs was displayed using TBtools software (Fig. 1). In total, 38 MtATGs were distributed across all eight chromosomes except for MtATG7, which did not 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) contained the largest number of MtATGs with 11 genes, whereas the fewest was found on chromosome 6, which had only one MtATG gene. Gene duplication is important for plant adaptation to adverse and complex environments. In M. truncatula, 7 pairs of MtATG genes were predicted to be segmentally duplicated genes. As shown in Figure 1

Phylogenetic analysis of MtATGs
To evaluate the evolutionary relationships of MtATGs, we conducted phylogenetic analyses using the amino acid sequences of multi-member subfamilies (MtATG1s, MtATG8s, MtATG9s, MtATG13s, MtATG16s, and MtATG18s) and their orthologs from Arabidopsis. As shown in Figure 2, MtATG1 and MtATG13 family members were both clustered in two branches (Fig. 2a, b). M. truncatula consists of multiple members in ATG9 and ATG16 families, whereas Arabidopsis has only one member. In addition, members of MtATG9 and MtATG16 families were similar and clustered together, while AtATG9 and AtATG16 were clustered in different branches (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 were absent in the N-terminus (Fig. S2). The C-terminal glycine residue of ATG8, which is exposed by ATG4 protease cleavage, is essential for the conjugation of ATG8 to phosphatidylethanolamine (PE) [27]. However, MtATG8b did not contain the C-terminal glycine residue. This result indicated that MtATG8b might function in other biological processes independent of the autophagy pathway. In addition, one MtATG8 member of clade II, MtATG8f, had a C-terminal extension after the Gly residue, while the AtATG8 members of clade II lack the C-terminal extension (Fig.  S2). Eight MtATG18 members were also clustered in two branches like MtATG8 family (Fig. 2f). Clade I of MtATG18 family consisted of MtATG18a, MtATG18b, MtATG18c MtATG18d, and MtATG18e, whereas clade II was made up of MtATG18f, MtATG18g, and MtATG18h (Fig. 2f).
Gene structure and conserved domain distribution analysis Gene structure is closely related to the expression pattern and function divergence of members of the multigene family. Gene structure analysis revealed that all the MtATG genes contained introns, with the exon numbers 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 that functional redundancy exists among these genes. However, differences in exon-intron patterns and exon numbers also existed within some subfamilies, such as MtATG1t, MtATG8b/h, and MtATG18b/f (Fig. 3a).
The conserved domains of MtATGs were detected using the Pfam database [28]. In general, the composition of the conserved domain in MtATGs was comparable to that in Arabidopsis. Furthermore, the MtATG members of the same family have similar domains. For example, all three MtATG1 proteins contain a protein kinase domain (Pkinase) at the N-terminal (Fig. 3b). A similar phenomenon was also observed in the MtATG8, MtATG9, and MtATG13 subfamilies. However, exceptions were also found in the MtATG16 and MtATG18 subfamilies. All the MtATG16 family proteins had a C-terminal WD40 domain, but lacked an N-terminal ATG16 domain in MtATG16c (Fig. 3b). MtATG18 proteins contain WD40 domain except for MtATG18b and MtATG18h, but members of clade II (MtATG18f/g/h) had a C-terminal BCAS3 domain that was absent in members of clade I (Fig. 3b). The differences of gene structure and conserved domains may be related to functional divergence among the different gene products within some MtATG gene families.

Cis-element analysis 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]. A total of 92 putative cis-elements were found among MtATG promoters (Table S2). Among these, the TATA-box and CAAT-box are the most common cis-elements. Many of the identi ed ciselements were involved in hormone responsiveness, such as ABRE (abscisic acid-related), TCA-element (salicylic acid-related), TCCACCT-motif and TGACG-motif (MeJA-related), TGA-element (auxin-related), TATC-box, P-box and GARE-motif (gibberellin-related) (Fig. 4). Among these, cis-elements that respond to MeJA and ABA were the most abundant. In addition, some stress-related elements were mainly related to anaerobic (ARE), defense (STRE and TC-rich repeats), drought (MBS), low temperature (LTR), and wound (WUN-motif) stresses (Fig. 4).

Protein-protein interaction network analysis of MtATGs
To investigate the protein-protein interaction (PPI) between MtATGs, all 39 MtATGs were submitted to STRING (Search Tool for the Retrieval of Interacting Genes database). A total of 22 MtATGs were found to form a complex interaction network that can be divided into four major modules according to the functional classi cation in Arabidopsis (Fig. 5). In the first module, MtATG1a, MtATG11, MtATG101, and three MtATG13 members (MtATG13a, MtATG13b, MtATG13c) interacted with each other and function as the ATG1 kinase complex. The second module consisted of two members of the PI3K complex, MtATG6 and MtVPS34. MtATG2 and six MtATG18 family members (MtATG18a, MtATG18b, MtATG18c, MtATG18f, MtATG18g, and MtATG18h) made up the third module, which played a role in autophagic membrane recruitment. The last module served as a ubiquitin-like conjugation system and is composed of MtATG4, MtATG5, MtATG12, and four MtATG8 members (MtATG8a, MtATG8d, MtATG8f, and MtATG8g). This interaction pattern of MtATGs was similar to that of Arabidopsis, suggesting that ATG proteins are likely evolutionarily conserved in Arabidopsis and M. truncatula.

Gene expression patterns of MtATGs during growth and development
To investigate the possible roles of the MtATGs in growth and development, the expression patterns of these genes were determined among different tissues and different developmental stages in seeds and root nodules. MtATG microarray data were acquired from the MtGEA database (https://mtgea.noble.org/v3/). All of the MtATG genes were expressed in tested tissues, indicating that autophagy is critical for the growth and development of plants (Fig. 6a). However, MtATG genes showed signi cantly distinct tissue-speci c expression patterns across tissues. Speci cally, the expression levels of many MtATG genes, such as MtATG4, MtATG8b, MtATG8g, MtATG9a, MtATG13a, MtATG13c, MtATG18b, MtATG18c, MtATG18e, MtATG18h, MtATG101, VPS15, and VPS34, were signi cantly higher in roots than in other tissues (Fig. 6a). In addition, some MtATG genes (MtATG1a, MtATG1t, MtATG2, MtATG7, MtATG9b, MtATG10, and MtATG18f) were highly expressed in leaves, whereas other MtATG genes (MtATG3, MtATG8a, MtATG8e, MtATG8f, and MtATG11) were highly expressed in owers (Fig. 6a). The results revealed that different MtATG genes might function in different tissues.
M. truncatula can develop root nodules through its symbiotic association with rhizobia, which is an important strategy to enhance nitrogen acquisition in legumes [31]. To investigate the role of MtATGs in root nodule symbiosis, we analyzed their expression pro les during nodule development at different days post inoculation (dpi). The expression of most MtATGs rapidly decreased at the stages of early lumps (4 dpi), and showed a slight increase in immature nodules (10 dpi), before declining to relatively low expression levels in mature nodules (14 dpi) (Fig. 6c). In the stable nodule organs (28 dpi), the expression of MtATGs was even lower. However, some genes, including MtATG9a, MtATG9b, MtATG13c and MtATG18g, were found to be most highly expressed at 28 dpi (Fig. 6c). These results indicated that autophagy is inhibited during nodulation.

Gene expression of MtATGs in response to abiotic stresses
Previous studies indicated that autophagy is essential for plant response to abiotic stresses, including those of heat, hypoxia, salt and drought [32][33][34][35]. To investigate the putative roles of autophagy in response to various abiotic stresses in M. truncatula, the expression pro les of MtATGs under salt and drought stresses were analyzed using microarray data from the MtGEA database. As shown in Figure 7, MtATGs displayed similar expression patterns in response to salt and drought stresses. Most MtATGs were up-regulated after salt and drought treatments, which indicated that autophagy may play a positive role in response to salt/drought stresses. Speci cally, 26 of 34 MtATG genes (e.g., MtATG1t, MtATG8d, MtATG9a, and MtATG18b, among others) were continuously up-regulated when plants were subjected to drought stress by withholding watering, and the transcripts of MtATGs rapidly dropped to their base levels after watering recovery (Fig. 7a). Interestingly, MtATG8g showed the opposite trend: the expression level of MtATG8g dramatically decreased under drought stress compared to other MtATGs (Fig. 7a). Furthermore, two members of the MtATG13 family, MtATG13b and MtATG13c, were also slightly downregulated when watering was stopped (Fig. 7a). Similar to drought stress, the expression of most MtATGs was induced under high salinity conditions. Twenty-one of 34 MtATG genes were continuously induced by salt stress (Fig. 7b). However, the expression levels of several MtATG genes (MtATG5, MtATG8c, MtATG18c, and MtATG18e) were decreased after NaCl treatment (Fig. 7b). These results indicated that autophagy played an essential role in plant response to salt and drought stresses.

Discussion
In this study, 39 ATGs have been identified in M. truncatula ( Fig. 1; Table S1). The ATG genes of M. truncatula are similar to orthologous genes in Arabidopsis. For example, most MtATG8s are similar in length and have identical ATG8 domains (Fig. 3). This result indicates that the autophagy pathway is highly conserved across different plant species. However, members of some MtATG families differ among plant species, such as ATG1, ATG4, ATG8, ATG9, ATG13, and ATG16 families ( Fig. 2; Table S1). In addition, MtATG8h showed distinctive gene and protein structures compared with the homologue in Arabidopsis, in which half of N-terminal amino acids and UTRs are absent ( Fig. 3a; Fig. S2). 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 further studies.
Autophagy has been shown to play crucial roles in plant growth and development [4]. In this study, we found that all ATG genes were expressed in the tested tissues of M. truncatula, but that their expression levels varied among different tissues, and especially differed over the course of the nodulation process (Fig. 6). Nodulation is a speci c process for nitrogen xation in leguminous plants [36]. A previous study reported that loss of PI3K function severely inhibited nodule formation in M. truncatula [37]. In addition, downregulation of ATG6 impaired the symbiosis between rhizobium and P. vulgaris [38]. These results revealed that autophagy is required for mycorrhization and nodulation in leguminous plants. On the contrary, silencing of TOR, a negative regulator of plant autophagy, caused a drastic reduction in nodule number in the common bean [39]. Consistent with this, we found that the transcript levels of most MtATGs decreased during nodule development (Fig. 6c). This result suggests a new role of autophagy in the nodulation of leguminous plants. Autophagy plays an important role in plant immunity. On one hand, autophagy protects plants from pathogen infection by limiting the spread of hypersensitive responses to programmed cell death (PCD) [40,41]. On the other hand, autophagy mediates selective degradation of viral components to promote plant survival [41]. In turn, increasing evidence has shown that pathogens have developed diverse strategies to modulate host autophagy and evade autophagic clearance for their own survival [41,42]. Therefore, we hypothesized that rhizobium represses gene expression of MtATGs to alleviate plant immune response and promote nodule formation, and then improve nitrogen xation during nodulation.
Seed development consists of embryo morphogenesis and seed maturation [43]. In rice, autophagy has been shown to be involved in the regulation of starch and sugar metabolism during seed maturation [44]. In Norway spruce (Picea abies), autophagy is also involved in embryogenesis in which it regulates vacuolar cell death of the embryo suspensor [45]. Furthermore, autophagy plays an important role in microspore embryogenesis in Brassica napus [46]. Previous studies have reported that autophagydefective mutants of Arabidopsis and maize exhibited lower seed weights [7,47]. 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 (Fig. 6b), indicating that autophagy is necessary for seed development in M. truncatula. Overall, autophagy plays crucial roles in plant growth and development through a pathway that is conserved across different plant species.
Autophagy has been demonstrated to promote plant survival through maintaining cellular homeostasis under stress conditions, including under nutrient starvation, oxidative stress, and drought, salt, and other abiotic stresses [48]. An early study reported that autophagy could be induced by salt and drought stresses. In A. thaliana, the transcriptional level of ATG18a was rapidly unregulated by NaCl and mannitol treatments [49]. Furthermore, ATG8a and NBR1 proteins rapidly accumulated after NaCl treatment [50]. In O. sativa, the expression levels of OsATG6s were also induced by drought stress [51]. In addition to gene expression, the autophagy-defective plant, RNAi-AtATG18a, showed more sensitivity to drought and salt treatment than the wild type [52]. However, overexpression of MdATG18a enhances tolerance to drought stress in apple trees [53]. A recent study reported that autophagy improves drought tolerance in M. truncatula through degradation of the aquaporin MtPIP2;7, which interacts with the cargo receptor MtCAS31 [54]. Consistent with previous studies, our results revealed that the promoter of many MtATG genes contained the drought-related MBS cis-element (Fig. 4). Furthermore, the transcriptional levels of most of MtATG genes, especially those of the MtATG8 gene family, signi cantly increased after NaCl treatments (Fig. 7a). Our results suggest that autophagy is involved in the positive regulation of plant resistance to drought and salt stress in M. truncatula, and that autophagy plays a conserved role in drought and salt stress response across different plant species.

Conclusions
This study provides a comprehensive overview of the ATGs in M. truncatula. A total of 39 MtATGs were identified, and their chromosomal locations, evolutionary relationships, cis-elements, gene and protein structures, and protein-protein interaction were systematically characterized. Gene expression analysis showed that MtATGs played important roles in seed development, nodulation, and response to salt and drought stresses. Our work provides an important reference for further research on the functional analysis of autophagy in M. truncatula.

Identi cation of MtATGs
The identi cation of putative MtATGs was conducted using a bi-directional blast analytical strategy, and was performed using the BLASTP program that is integrated into BioEdit software. First, the protein sequences of published autophagy-related genes in Arabidopsis were used to search against M. truncatula proteome sequences (MedtrA17_4.0) with the e-value cutoff at 1×E-5. Then, all output M. truncatula protein sequences were aligned back to Arabidopsis proteome sequences. Only the M. truncatula genes that shared the highest similarities to the AtATGs from the second BLAST analysis were considered putative MtATGs. To further verify that the candidate MtATGs are indeed MtATGs, the protein domain architectures were analyzed in the Pfam database (http://pfam.xfam.org) [28]. The chemical features of the MtATG proteins, including their molecular weights and theoretical isoelectric points, were obtained using the online tool ExPASy (http://web.expasy.org/compute_pi/). Subcellular localization of MtATGs was predicted using the CELLO v2.5 system (http://cello.life.nctu.edu.tw). The gene and protein structures of MtATGs were extracted from the annotation le of the M. truncatula genome (MedtrA17_4.0) and visualized with the integrating bioinformation analysis toolkit Tbtools [55].

Chromosomal location and gene duplication analysis
MtATGs were mapped to the M. truncatula chromosomes based on their physical positions in the the M. truncatula genome (MedtrA17_4.0). To investigate the synteny relationships of related genome regions in M. truncatula, putative orthologous genes were identi ed using the BLASTP program, and the results were used to generate a synteny map with the MCScanX program [56]. The genome locations of MtATGs and the duplicated gene pairs were visualized using Tbtools.
Protein sequence alignments and phylogenetic relationship analysis Phylogenetic analysis of MtATGs was performed using MEGA7 software [57]. 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 neighborjoining statistical method, and the following parameters were used: uniform rates were used as rates among sites, gaps/missing data were treated as pairwise deletion, and the bootstrap analysis was performed with 1000 replicates to obtain a support value for each branch.
Protein-protein interaction (PPI) network Construction PPI network were constructed using the STRING v11 database (http://www.string-db.org). A total of 39 MtATGs were selected as input, and the PPI network of the MtATGs were constructed with medium con dence (0.4).
Expression pro le analysis using microarray data The M. truncatula microarray data were downloaded from the MtGEA v3 database (https://mtgea.noble.org/v3/) [58] for datasets generated for different organs [59], drought stresses condition [60], and salt stresses condition [61]. Expression values were normalized using the z-score method or transformed to log2 fold change, and plotted using GraphPad Prism 8.

Funding
This work was supported by the National Natural Science Foundation of China (31700236). The funding organization provided the nancial support to the research project, but had no role in the design of the study, data collection and analysis, and writing of the manuscript.

Availability of data and materials
All data analyzed during this study are contained in the materials and methods section of the manuscript and its additional les.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.        The number of cis-elements in promoter of MtATGs. The assumed cis-elements of MtATGs were predicted using the PlantCARE web servers, and the number of cis-elements in each promoter of MtATGs are visualized for heatmap using GraphPad Prism 8.   and nodules were dissected from roots at indicated days after being inoculated (DPI) with S. meliloti.
Scale bar represents the relative expression value after z-score normalization.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.