Identification of CsARGs in tea plants
Based on three different identification methods, a total of 35 CsARGs were identified from the two published tea plant genomes (‘ShuChaZao’ and ‘YunKang10’). Among these genes, four genes (CsATG1s, CsATG8s, CsATG18s and CsVTI13s) were determined to have isoforms (Table 1). In this study, we identified a UV radiation resistance protein/autophagy-related protein 14, known as CsATG14, which has not been fully studied in plants to date and may be deficient in Arabidopsis. In addition, as there is no single strict criterion to identify all the paralogs of these four genes, and as the genome assembly and gene annotation of the two reported tea plant genomes have not been fully completed to date compared to the Arabidopsis, rice and tobacco genomes, the products of the partial paralogs of these four genes, as mentioned in Arabidopsis, rice, and tobacco, were not determined in the tea plant genomes. Bioinformatic analysis results showed that as a type of biological macromolecule, the CsARG ORF lengths varied from 285 to 7410 bp, the corresponding numbers of deduced amino acids ranged from 94 to 2469 aa, and the molecular weights ranged from 10.51 to 276.87 kD. The theoretical isoelectric points (pIs) were predicted to range from 4.52 to 9.41. The prediction of subcellular location results suggested that most CsARGs were located in the nucleus, and some of them were also predicted to be located in the cytoplasm, chloroplasts and mitochondria. The results of signal peptide prediction showed that none of these CsARGs contained signal peptides. In addition, CsATG9 was predicted to contain 5 TMHs, CsATG18b, and three transport vesicle-soluble NSF attachment receptor (v-SNARE) proteins, namely, CsVTI12, CsVTI13a and CsVTI13b, were predicted to contain 1 TMH.
Phylogenetic analysis of CsARGs in tea plant
To explore the evolutionary relationships and classification of CsARGs in tea plant, a total of 177 ARG proteins from tea plants, Arabidopsis, Setaria italic, Oryza sativa, and Nicotiana tabacum were aligned to construct a phylogenetic tree. As shown in Fig. 1, except for CsATG14, which had only been identified in tea plants, all of the CsARG proteins were highly clustered together with the homologous proteins derived from the other four species, and almost all of the CsARGs showed the closest relationship with NtARGs. Meanwhile, we observed that the bootstrap values among the different ARG proteins in each subtree were nearly 100%, except for the ATG8 subfamily, which suggests that ARG protein sequences are highly conserved and may exhibit similar functions among different species.
Gene structure, protein domain distribution and cis-acting element analysis
Understanding the exon-intron structure is beneficial for exploring the evolution of multiple gene families [57]. To investigate how the differences in exon-intron structure were generated, both the genomic and ORF sequences of CsARGs were uploaded into GSDS v2.0 to predict the exon-intron structure. As shown in Fig. 2a, the numbers of exons in the CsARG family varied, with members within the CsATG8s or CsVTI13s subfamilies exhibiting similar exon-intron structures.
To further dissect the functions of CsARG proteins, the protein domains of each CsARG were analyzed by the SMART program. As shown in Fig. 2b, CsATG1s encode serine/threonine protein kinases, which contain catalytic domains involved in protein phosphorylation occurring during the progression of autophagy [58]. CsATG9 contains 5 transmembrane helix regions (88-110, 155-177, 320-342, 403-425 and 438-457), as detected by the TMHMM v2.0 program, which play a unique role in autophagosome formation derived from the endoplasmic reticulum (ER) in plants [59]. CsATG16 contains a coiled coil region and 7 WD40 domains, which form a conserved Atg12-Atg5-Atg16 complex during the autophagy process. Each member of the CsATG18s subfamily is largely a β-propeller and is formed by 2 or 3 WD40 domains. CsATG20, also known as Snx42, contains a PX domain, which plays a central role in efficiently inducing nonselective autophagy. CsVPS15 encodes a serine/threonine protein kinase, which is formed by 5 WD40 domains and is regulated by a phosphoinositide 3-kinase (PI3K), CsVPS34. At present, CsVPS34 is characterized as a central regulator in mediating vesicular trafficking and cellular homeostasis [60]. As an ATG8-interacting protein, CsATI1 contains a transmembrane region that may help the protein complex to move to the ER network and reach the lytic vacuole. CsVTI1s, including CsVTI12, CsVTI13a and CsVTI13b, all contain a coiled coil region, a t-SNARE domain and a transmembrane region, and these domains may be involved in the trafficking of cargoes to the vacuole. CsTOR, as a conserved phosphatidylinositol kinase-related protein kinase, contains a specific rapamycin-binding domain, a P13Kc-catalyzing domain and a FATC domain, which suggests that it is involved in mediating redox-dependent structural and cellular stability.
To elucidate the regulatory mechanisms governing CsARGs in response to growth and development, stress defenses, and hormone signaling, a 2000-bp 5’-upstream noncoding region sequence of each CsARG was isolated to predict cis-elements. As Fig. 2c shows, the distributions, numbers and types of cis-elements vary among the promoter sequences. Nevertheless, most of the promoters contain a number of MYB- and MYC-binding sites, except for the promoter of CsATG12, which lacks MYC-binding sites. In addition, most of the promoters of CsARGs contain ABA- and MeJA-responsive elements, and some of them contain GA, SA, auxin, cold, drought, defense and stress (TC-rich repeats)-responsive elements. In addition, all of the promoters of CsARGs contain numerous light-responsive elements, including G-boxes, MREs, Box-4, and AE-boxes (not shown in Fig. 2c). Overall, each CsARG may play a vital role in responding to circadian variation, hormones, and biotic and abiotic stresses.
Protein-protein interaction networks of CsARGs
To investigate the interactions among CsARGs in tea plants, Arabidopsis was used as the reference species to construct PPINs. As Fig. 2d shows, 34 CsARGs were matched to 33 AtARGs, and those 34 CsARG proteins formed 333 protein-protein association patterns. The ARGs were determined to be closely related to one another, except for ATG14, which has not been fully studied in Arabidopsis. Among these proteins, 23 ARGs, including ATG9, ATG7, and ATG1c, were reported or predicted to interact with more than 20 ARGs, suggesting that the occurrence of autophagy requires interactions among numerous ARG proteins.
Conserved domain and motif distribution analysis of CsATG8s
As members of the UBQ superfamily, ATG8s coupled with their conjugation system are key components for autophagy. In our study, to clearly understand the regulatory mechanisms governing these proteins, the bioinformatic characteristics of CsATG8 subfamily proteins were further explored. As Fig. 3 shows, a total of 5 CsATG8s were identified in tea plants based on homologous alignment analysis. Phylogenetic analysis results showed that CsATG8s were subdivided into 3 clades. Among these proteins, CsATG8a/c and CsATG8g/f were clustered into one group, and CsATG8i was aligned closely with MdATG8i (Fig. 3). Motif distribution analysis results showed that CsATG8s contain motifs 1-4, and CsATG8i contains an additional motif 7 at the C-terminus (Additional file 2). All of these CsATG8s proteins contain conserved GABARAP domains, four putative tubulin binding sites, three ATG7 binding sites, and a conserved glycine (G) residue. In addition, we found that the conserved G residues in both CsATG8a and CsATG8g were directly exposed at the C-terminus (Fig. 3), suggesting that the functions of CsATG8a and CsATG8g involved in autophagy may be distinct from those of the other three proteins. These proteins may not require the cysteine protease Atg4 to cleave the C-terminus but may directly bind to the E1-like enzyme Atg7.
Expression profiles of CsARGs in various tea plant tissues
To confirm the tissue specificity of CsARGs, the roots, stems, mature leaves, tender leaves, and seeds of tea plants were obtained for qRT-PCR analysis. The results showed that the transcription of all CsARGs was detected among the above mentioned tissues, although the mRNA level of each CsARG varied across the various tissues (Fig. 4). In addition, we found that most CsARGs exhibited higher transcription abundances in stems and seeds, suggesting that autophagy plays important roles in the development of stems and seeds in tea plants. Moreover, we found that the CsATG3/7/101, CsVPS15/34, CsATI, CsVTI12/13b, CsATG8s and CsATG18s subfamily genes were highly expressed in different tissues. Notably, CsATG8c was significantly expressed in mature leaves and seeds, and CsATG8i was dramatically expressed in stems and seeds. In brief, our results found that the expression patterns of the CsARGs varied across different tissues, but some of them showed a degree of tissue specificity, suggesting that CsARGs mediated the growth and development of tea plants.
Differential expression of CsARGs in response to hormone treatments
To elucidate the comprehensive roles of CsARGs under ABA and GA treatment conditions, we analyzed the expression pattern of each CsARG. Under ABA treatment conditions, we observed that multiple CsARGs were highly induced after 12 h and/or 2 d of ABA treatment. Among these genes, the expression levels of CsATG6/11/12/14/16/18b/18c/18f/20/101 and CsVTI12/VTI13b were induced more than 2-fold for at least at one time point. Meanwhile, some genes, such as CsATG3/11/18c/101/VPS15/TOR, were directly upregulated during the entire ABA treatment period (Fig. 5). In contrast, CsARGs exhibited contrary expression profiles under GA stress compared to ABA stress. Most CsARGs initially decreased but significantly increased after 2 d of GA treatment. Among these genes, the expression levels of 13 genes, including CsATG12/14/16/18a/18b/18c/18g/18h/20/ATI/VTI12/VTI13a/VTI13b, were increased more than 2-fold after 2 d of GA treatment. Furthermore, 6 genes, CsATG7/8c/8f/101/VPS15/VPS34, were negatively downregulated throughout GA stress periods. These results indicated that CsARGs are required for the response to hormone treatments in tea plants.
Expression patterns of CsARGs in response to different abiotic stresses
Similarly, to explore the temporal expression patterns of CsARGs under abiotic stress conditions, the related expression level of each CsARG was determined by qRT-PCR. As shown in Fig. 6, the expression levels of all CsARGs were regulated to different degrees under various abiotic stress conditions.
During the CT period, the expression levels of almost all CsARGs were upregulated at different processing time points, except for 2 genes, CsATG18c/h, which were downregulated throughout the CT period. In addition, 20 induced CsARGs showed the highest expression levels after 12 h of CT. Specifically, the expression levels of CsATG3/4/6/8i/10/11/18a/18g/101/VTI12/VTI13a/VTI13b/VPS34 were more than 5-fold higher than those at 0 h of CT. In addition, CsATG5/ATG18g/ATI expression was gradually induced within 2 d of CT. Under DT conditions, 15 CsARGs were induced at different DT time points, and most of these genes showed the highest expression levels after 12 h of DT. The expression levels of CsATG2/3/6/8i/11/16/VPS34 were more than 2-fold higher than those at 0 h of DT. Moreover, CsATG18h and CsATG101 were gradually upregulated as the DT time extended. Conversely, the expression of the remaining CsARGs, such as CsATG4/7/8a/10/14, was slightly deduced or not affected by DT. Within NT periods, many CsARGs were also induced to different degrees. In contrast to CT, most of the induced CsARGs showed the highest expression levels after 1 d of NT but were deduced following the extended treatment time. In particular, the expression levels of CsATG9/10/18a/18g/101/VTI13 were increased more than 2-fold after 1 d of NT compared to the control. However, CsATG8g/12/16/ATI was significantly upregulated after 12 h of NT. Taken together, these results demonstrated that CsARGs play central roles in responding to abiotic stress in tea plants.
Differential expression of CsARGs in different tea plant cultivars during CA periods
The differential expression patterns of CsARGs were compared in different tea plant cultivars (the cold-resistant cultivar ‘Longjing43’ and the cold-susceptible cultivar ‘DaMianBai’) within CA periods. Eleven CsARGs, which were confirmed to be notably induced under CT conditions, as shown in Fig. 7, were selected for analysis by qRT-PCR. The results showed that these 11 genes presented different expression patterns in ‘Longjing43’ and ‘DaMianBai’ during the CA period in 2018-2019. Specifically, these 11 genes showed contrary expression patterns from November 14th to December 13th in ‘Longjing43’ and ‘DaMianBai’. With the exception of CsVPS34, the other CsARGs were continuously upregulated in ‘Longjing43’, but all of them were gradually reduced in ‘DaMianBai’ from November 14th to December 13th. In contrast, the transcription levels of these 11 genes were all increased in ‘DaMianBai’ from Dec 13th to Jan 17th, but many CsARGs, such as CsATG16/18g/101/VTI12, were decreased in ‘Longjing43’. Notably, we observed that the transcription level of CsATG101 was lower in ‘LongJing43’ than in ‘DaMianBai’ throughout the CA period. These results indicated that these CsARGs play important roles in the response of tea plants to cold resistance, but their regulatory mechanisms may vary among different cultivars.