Zinc (Zn) and iron (Fe) are key elements for plant growth and development. (Zn) is an enzymatic cofactor, participating in electron transport and antioxidant metabolism. On the other hand, it is an integral element of several transcription factors and plays important roles in enzyme activation, gene expression, protein synthesis, photosynthesis and carbohydrate metabolism (Chang et al., 2005; Slaton et al., 2005). (Fe) is required for electron transport in chloroplasts and mitochondria. It is involved in the hormones biosynthesis and integrity of cellular organelles; as well as in the tolerance to oxidative stress and in the N fixation (Hänsch and Mendel, 2009). Zn and Fe deficiency causes many nutritional, biochemical and molecular disorders. A complex control system allows the uptake, distribution and accumulation of these ions, ensuring metal homeostasis. The uptake of metal ions, under deficiency conditions, activates several transport proteins, such as the zinc/iron-regulated transporter-like protein (ZIP) family. Nowadays, members of the ZIP family have been identified in bacteria, archaea, and eukaryotes as well as protists, animal fungi, and plants (Mäser et al., 2001). ZIP gene family in the model species Arabidopsis thaliana has been extensively studied. For most ZIP proteins, the topology with the N-terminal and C-terminal ends is predicted to be located on the outside surface of the plasma membrane, with eight potential transmembrane domains (TMDs). The variable region contains a potential metal-binding domain rich in histidine residues (Guerinot, 2000).

Studies on the function of ZIP genes in A. thaliana and the uptake of metal ions start with the IRTI (iron-regulated transporter) gene, which was the first identified gene encoding a transporter involved in the Fe²⁺ uptake with localization in the plasma membrane (Eide et al., 1996). IRT3 and ZIP4 genes are induced by Zn deficiency and are involved in the Zn increase in xylem (Grotz et al., 1998). IRT1 is also responsible for the uptake of elements such as cadmium (Cd) and Zn. Plants with overexpression of IRT1 hold higher levels of these elements (Connolly et al., 2002). Genes expressed in roots and leaves such as ZIP2, ZIP4, ZIP5 and ZIP9 are up-regulated in Zn-deficient conditions, while IRT2 is involved in Fe uptake. Two, ZIP2 and ZIP4 genes are also transcriptionally regulated by copper (Cu) (Wintz et al., 2003). Other members of this family, such as ZIP1 and ZIP3, are mostly expressed in roots and are implicated in Zn uptake in A. thaliana (Ishimaru et al., 2005). IRT1, a gene which plays a role in the Fe, manganese (Mn), Zn and Cd uptake (Cohen et al., 1998; Vert et al., 2002) also mediates the accumulation of nickel (Ni) in A. thaliana (Nishida et al., 2011). After that, a lot of ZIP family members have been found in the species Oriza sativa (Ishimaru et al., 2011). The genes OsZIP1, OsZIP4, OsZIP5 and OsZIP8 are involved in the Zn transport and distribution. The localization of these proteins is in the plasma membrane of roots and shoots (Ishimaru et al., 2007; Lee et al., 2010; Lee et al., 2010; Yang et al., 2009). Thus the ZIP plays key roles in the uptake and translocation of essential and non-essential metal ions in different plant organs (Eide et al., 1996; Grotz et al., 1998; Wintz et al., 2003; Li et al., 2013; Milner et al., 2013) and has been studied in several plant species of economic importance, such as Zea mays (Li et al., 2013), Glycine max (Moreau et al., 2002), Hordeum vulgare (Tiong et al., 2014), Vitis vinífera (Gainza-Cortés et al., 2012), Solanum tuberosum (Li et al., 2020) and Poncirus trifoliata (Fu et al., 2017).

Information on ZIPs was not systematically analyzed in the T. cacao genome, hampering the understanding of its role in the processes of uptake and distribution of essential and non-essential elements such as Cd, which is toxic and commonly affects this crop (Araújo et al., 2017). T. cacao genome, at the chromosomal level, in its recent high-quality version (Motamayor et al., 2013) provides quality resources to study, at the genomic level, the role of the ZIP gene family.

In the present study, we were the first to identify and characterize members of the ZIP gene family based on phylogenetic relationships, physico-chemical properties, transmembrane domains, conserved motifs, exon/intron organization, substitution rates, chromosomal locations and cis-regulatory elements in cacao. We also explored post-translational modifications, 3D modeling, prediction of ZIP transport protein pocket sites and protein-protein interaction (PPI) network. The results reported in the present study will increase our knowledge about the evolution, function, structure and interaction of ZIPs, being the first step towards future studies on the processes of uptake and translocation of essential and non-essential elements in cacao.

The ZIP family is a group of metal ion transporter genes found in various kingdoms including plants (Guerinot, 2000). This family is involved in the transport of various essential and non-essential metal ions. In the present paper, we performed a general analysis of the ZIP gene family in cacao, including analysis of its phylogeny, chromosomal location, gene structure, conserved motifs, prediction of three-dimensional structure and protein-protein interactions (PPIs). A total of 11 ZIP genes were identified in the cacao genome ( Table 1 ), which is fewer than the ZIP members identified in the model species A. thaliana (15) and a similar number as in Solanum tuberosum (12) (Li et al., 2020). Eventually the ZIP family is relatively small, so far the smallest number of ZIP members has been found in Zea mays (Li et al., 2013). Sequence analysis showed that the TcZIP proteins bear similarity to the previously studied ZIP families (Guerinot, 2000). Specifically, TcZIP members encode proteins with lengths from 309 to 435 aa, in agreement with the predicted range in known plants (Guerinot, 2000). Studies in perennial plants such as Poncirus trifoliata reveal that ZIP members encode proteins with 334 to 419 aa (Fu et al., 2017). Most of the TcZIPs were predicted to be located in the plasma membrane, similarly to AtIRT1, OsZIP4, OsZIP5, HvIRT1 and VvZIP3 (Vert et al., 2002; Ishimaru et al., 2005; Pedas et al., 2008; Lee et al., 2010; Gainza-Cortés et al., 2012; Krausko et al., 2021), and also predicted in the ZIPs members of the Poncirus trifoliata and Zea mays genomes (Li et al., 2013; Fu et al., 2017). However, two TcZIP genes are predicted to be located in the chloroplast. Certainly, ZIP genes can be expressed in various parts of plants, and they can be located in the membrane of internal organelles such as chloroplasts (Ajeesh Krishna et al., 2020). Furthermore, between 6-9 MTs were identified in TcZIP proteins, similarly to ZIPs Poncirus trifoliata and Zea mays (Li et al., 2013; Fu et al., 2017).

On the phylogeny analysis, the 11 TcZIP genes identified are grouped into four subgroups, similar to ZIP genes of A. thaliana (Mäser et al., 2001). In each group, most genes have characteristics of the exon/intron structure and relatively conserved motif compositions in recent paralogs, except for the exon/intron structures of two paralog pairs, in which some differences still exist ( Figure 1 ). For instance, structural analysis showed that TcZIP4 and TcZIP9 have different numbers of introns, in contrast to HD-zip gene family in Glycine max that showed relatively conserved exon/intron structure among the paralog pairs. This is due to the few insertions and deletions that accumulated along the evolution (Schmutz et al., 2010). Studying the structure of the exon-intron gene may provide important clues to the evolution of the gene (Freeling, 2009), because genes in the same group can have similar functions. The interspecies phylogeny results show that they are grouped into four distinct groups, in which the first group has the highest number of proteins. Although the groups identified are different from those found by (Li et al., 2020), our results showed consistency in the groups identified in the phylogenetic relationships made in the Vitis vinifera perennial species (Gainza-Cortés et al., 2012). Thus, the discrepancy with that presented by Li et al. (2020) may be due to the numbers of genomes used; in this one ZIP members of five plant genomes were used.

The phylogenetic study showed putative orthologous ZIP genes, unexpectedly there was also an orthologous pair of monocots and dicots (TcZIP7-OsZIP2), suggesting that the orthologous pair originated from common ancestral genes that existed before the divergence of monocots and dicots. Due to the existence of several paralogous genes in the TcZIP family, it is possible to confirm that the cacao tree underwent two duplication events after the monocotyledonous/dicotyledonous split. This suggests that most of the ZIP genes in cacao expanded in a species-specific manner. Gene duplication plays an important role in the evolution of the organism; the generation of new genes will allow organisms to adapt to different environments. Two main evolutionary mechanisms have been attributed to gene duplication, including segmental and tandem duplication (Horan et al., 2005; Kong et al., 2007). In our analysis, we found five pairs of duplicated ZIP genes, the distribution was preferentially for segmental duplications, suggesting that segmental duplication was the main driver for the expansion of the ZIP gene family in cacao. Certainly segmental duplication is the mechanism frequently discovered in plants, as most plants are diploidized polyploid and retain numerous duplicated chromosomal blocks within their genomes (Cannon et al., 2004; Del Pozo and Ramirez-Parra, 2015). Estimates of the Ka/Ks ratio were utilized to measure the selection rate of genetic pressure. Based on these results, we propose that purifying selection is the main evolutionary force in the ZIP paralog pairs in cacao. Since, based on the theory of natural selection, when the ratio Ka/Ks >1 indicates positive selection, while the Ka/Ks =1 indicates neutral selection, and finally when the ratio Ka/Ks <1 indicates negative or purifying selection (Gaut et al., 1996). Previous studies indicate that the cacao genome has undergone lineage-specific scrambling events from the paleohexaploid ancestor. It is estimated to have occurred 123 million years ago (MYA), and at least 11 chromosomal fusions have occurred to achieve the recent structure of ten chromosomes (n), estimated to have occurred around 59 MYA (Argout et al., 2011). By calculating the duplication dates of the paralog pairs, we showed that all segmental duplication events in the ZIP family in cacao occurred after the evolutionary scenario that gave rise to the ten chromosomes. The results indicated that the gene pair with tandem duplication appeared 5.99 MYA ago, being the most recent event. Almost all cacao chromosomes contain duplicated genes, so we assume that the ZIP gene family is a slowly evolving gene family, in the same way as was found in the TCP gene family in Z. mays (Chai et al., 2017). Evidently segmental duplication played an important role in the expansion of ZIP genes. Furthermore, the cis-acting elements that control many biological processes and responses to different stimuli through a regulatory system allow genes to have important unique combinations (Ibraheem et al., 2010). In this study ABRE, TGA-element, Box4, GATA-motif, GT1-motif and ATCT-motif, Myoblastosis and Tc-rich are essential elements for the stress response. These cis-acting elements are abundant in TcZIP as well as in potato ZIP members (Li et al., 2020).

Predicting the 3D structure and identifying the pocket site of proteins provides the structural and microenvironment basis for proteins to perform functions such as ligand binding, enzymatic activity and DNA interaction (Edelsbrunner et al., 1998; Tian et al., 2018). Analysis of the three-dimensional structure showed that α-helices conformations are predominant and β-sheets are practically absent, in the same way as reported in ZIP protein structures in other species (Ajeesh Krishna et al., 2020), showing that α-helices play an important role in ZIP proteins. All predicted TcZIP models have structural similarity to the recent Zrt-/Irt-like protein (ZIP) crystal structure of Bordetella bronchiseptica (Zhang et al., 2020), which has a key function in the transport of bivalent transition metal ions (Zn, Fe, Cd, Mn, Ni, Co e Cu) (Antala and Dempski, 2012; Jeong and Eide, 2013). Pocket sites and predicted cavities were diverse in the structures of TcZIP proteins, which provides valuable insights into protein function based on metal-binding sites. The amino acid residues LEU, ALA, VAL, GLY, SER, LYS and HIS were frequently predicted within the pockets for almost all candidate proteins. Among the residues found, GLU and HIS were identified in previous studies as part of the substrate binding sites (Antala et al., 2015). Another study showed that the residues GLU, HIS, ASP and CYS coordinate iron binding in various wheat proteins (Verma et al., 2017). A later study shows that LYS, PHE and SER residues also contribute to the transport process (Gyimesi et al., 2019). These residues were found in pocket sites of TcZIP proteins, suggesting they are important for transport activity. Post-translational modifications are modifications in the side chain of amino acids in some proteins after their biosynthesis, these modifications regulate several cellular processes (Ramazi and Zahiri, 2021). All ZIP members in cacao have potential phosphorylation sites. Previous studies have revealed that ZIP proteins undergo many phosphorylations. Protein phosphorylation by a kinase (CIPK23) and IDF1-mediated polyubiquitination are important for efficient endosomal sorting and vacuolar degradation of IRT1. Moreover, phosphorylation as a subsequent process of protein-metal binding helps to optimize the process of uptake and protection of plants from harmful metals (Dubeaux et al., 2018). It is also known that phosphorylation at SER residues of ZIP7 in humans dramatically increases its transport activity (Hu, 2021). In addition, potential sites for N-glycosylation were also predicted. Six TcZIP proteins are glycosylated at the N-terminal. In the current study, TcZIP3 and TcZIP7 are the highlighted proteins (three sites both). This indicates that ZIP proteins in cacao have few sites for N-glycosylation. It has been demonstrated in previous studies, proteins from the ZIP gene family are rarely glycosylated, and those that show N-glycosylation exhibit values of less than four sites as observed in ZIP4, ZIP6, ZIP8 and ZIP14 (Taylor et al., 2003; Wang et al., 2004; He et al., 2006; Girijashanker et al., 2008).

Studies of PPIs have revealed valuable information about the possible interactions and processes involved in target proteins. The function of a protein is governed by its interaction with other proteins within a cell, PPIs establish one of the most decisive conditions for sustaining life in living organisms (Yilmaz et al., 2022). A systematic analysis of the system biology was performed to obtain a broader view of the processes that evolve the ZIP transport proteins homologous to those identified in cacao. In this study, we identified bottlenecks and hubs proteins (Verli, 2014) and 109 proteins constituted the interactome. Proteins qualified as bottlenecks proteins (betweenness value above average) are of importance because they are related to their ability to connect several clusters, which, in this case, represent biological processes. The network also contains 53 highly connected proteins, qualified as hubs proteins (above average node degree), since these proteins have the characteristics of having an important regulatory function within the interactome. The ZIP1 protein with both characteristics (bottleneck and hub) allows the junction between cluster 2 related to cation transport and cluster 1 related to DNA metabolic process. Surprisingly, the interaction of ZIP members reveals other processes besides the transport of cations. Previous studies show that the ZIP1 overexpression leads to apoptosis of tumor prostatic epithelial cells, due to the increase in the concentration of Zn (Samavarchi Tehrani et al., 2019), which leads us to believe that this homeostatic Zn dysregulation causes damage to DNA molecule and cell replication. DNA damage interferes with the development of organisms. Molecular mechanisms for metal toxicity in plants have been described previously (Gill, 2014). A complex network of proteins is activated to protect DNA. The DNA repair response includes different pathways in plants (Tuteja et al., 2001) that are similar to those used in other organisms. For example, MSH family proteins found in cluster 1 can form MSH2/MSH6 and MSH2/MSH7 complexes important in recognizing nucleotide mismatches in A. thaliana (Moura et al., 2012), in this way, they mitigate the damage resulting from metal toxicity. Moreover, in cluster 2 with proteins related to the cation transport process, it was shown that the bottleneck proteins ZIP11, ZIP2, ZIP3, and IRT1 interact with other transport proteins.

P-type ATPases (HMAs) are another family of metal transporters found in various kingdom including plants (Baxter et al., 2003). HMA proteins play an important role in Zn, Cd, Pb and Cu homeostasis and have been expressed in various tissues (Hussain et al., 2004; Mills et al., 2003), which demonstrates that they are involved in compartmentalization and detox (Andrés-Colás et al., 2006). Thus, AtHMA3 plays a role in the detoxification of biological (Zn) and non-biological (Cd, Co, and Pb) heavy metals, participating in their vacuolar sequestration, a unique function for a P1B-2 ATPase in a multicellular eukaryote (Morel et al., 2009). The Yellow Stripe-Like (YSL) protein can mediate the transport of metals in the form of a complex. Heterologous expression of BjYSL7 revealed the transport of Fe, Cd and Ni (Wang et al., 2013) using a metal-nicotianamine complex. Similarly SnYSL3 revealed the transport of Fe, Cu, Zn and Cd (Feng et al., 2017). IAR (IAA-Ala–resistant) gene encode a protein involved in auxin metabolism or response. The protein encoded by the IAR gene has several His-rich regions, with structural features shared with the ZIP family of transporters (Lasswell et al., 2000).

Natural resistance-associated macrophage proteins (NRAMP) are a family of metal ion transport proteins identified in yeast, insects, mammals and plants (Cellier et al., 1995) and have substrate specificity for various metal ions, including Cd (Nevo and Nelson, 2006). Recently, several members of NRAMP have been identified in cacao, with NRAMP5 being a potential Cd transporter in cacao (Ullah et al., 2018). Metal transporter CDF proteins encode proton antiporters that effuse heavy metals out of the cytoplasm (Gaither and Eide, 2001). The expression of CDF genes such as MTP1 confers Zn tolerance, as well as knockdown in mtp1 mutants produces Zn-sensitive phenotypes, demonstrating the role of MTP1 in Zn homeostasis (Kobae et al., 2004; Desbrosses-Fonrouge et al., 2005). Studies show that the MTP5 and IAR1 genes transport metals in an antagonistic way, modifying the intracellular ionic composition to influence the activity of the metallohydrolase responsible for hydrolyzing the conjugates with of indole-3-acetic acid (IAA), influencing the homeostasis of IAA (Rampey et al., 2013).

Knowing about the proteins of our interactome, we can infer that there is a complex network of interaction that mediates metals homeostasis in cacao, whose ZIP family is a part of this complex. The activities of import to cell, distribution to organelles and export from the cell, need a controlled and balanced activity between several transporters. Transcriptional control is an important factor in the regulation of cellular homeostasis. Although little is known about the transcriptional profiles in cacao and the main transcription factors that mediate metal homeostasis, this analysis of ZIP transporters PPI contributes to the understanding of a complex network of transport proteins that interact and possibly coordinates from the entry of metallic ions to their use in different tissues and deposition in cacao beans.

In summary, 11 ZIP genes were identified using the cacao genome for the first time. All sequences analyzed suggest that they have basic characteristics of previously studied members of the ZIP family. Among genes identified, three including TcZIP1, TcZIP2 and TcZP4 are the most distant from their counterparts in the model species. Segmental and tandem duplication have been identified as the main patterns contributing to the expansion of the ZIP gene in cacao. All TcZIP play important roles in the uptake and distribution of transition metals, including the toxic element cadmium, according to pocket site prediction analyses. This study provides essential and comprehensive information for ZIP genes in cacao, but we propose that further work is needed to delve into exact subcellular and tissue localization, gene expression between different genotypes, and KO analysis (gene knockout).

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material .

DP and A-AA performed the original design preparation, methodology, validated the data, wrote and edited the manuscript and visualized the data. CP and BS performed the formal analyses, curated data and revised the manuscript. A-AA and CP performed resource administration, project administration and financing acquisition and oversaw the data. All authors contributed to the article and approved the submitted version.