ABSTRACT

   Cucurbits are cultivated worldwide in regions with high concentrations
   of arsenic (As), a hazardous metalloid, affecting produce quality and
   increasing the consumer exposure. Cucurbitacins are herbivore‐deterrent
   secondary metabolites that contribute to the plant defense response.
   The impact of As exposure on phenotypic and metabolic traits has not
   been studied in members of the Cucurbitaceae family, such as squash (
   Cucurbita pepo L.). To comprehend the effects of As on the root system
   of C. pepo , we assessed phenotype, cucurbitacin content, and
   transcriptome under low and high As concentrations. We report that at
   low dosages, cucurbitacins are decreased, while growth is not
   significantly affected. Conversely, high dosages impact growth and
   development altering root phenotype but cucurbitacin content is not
   significantly different from untreated plants. Furthermore, gene
   ontology enrichment on results of the RNA‐seq analysis indicate that
   high dosages of As affect cellular regulatory processes, with genes
   related to glutathione metabolism being of the most upregulated.
   Additionally, an in‐depth analysis of orthologs members of the heavy
   metal–associated (HMA)‐domain superfamily and As‐related transporters
   suggest a dosage‐dependent participation of key members. WGCNA analysis
   reveals As‐specific gene co‐expression modules, indicating that low As
   levels induce adaptive responses in energy and allantoin metabolism,
   while higher levels trigger intensified oxidative stress responses,
   including upregulation of MYB transcription factors and heat shock
   proteins, which may support tolerance to the metalloid. Overall, As
   influences the root system physiology and metabolism in a
   concentration‐specific manner, highlighting key defense systems and
   genes involved in C. pepo response to As exposure.

   Keywords: arsenic, Cucurbita pepo , Cucurbitacins, plant defense,
   stress mechanisms

1. Introduction

   Cucurbitaceae species are major crops. Among them, squash ( Cucurbita
   pepo L.) was valued at approximately $215 million US dollars in 2023
   (USDA, 2023) due to its fruit, flowers, seeds, and the production of
   its bioactive oil, which is used in traditional medicine systems
   worldwide (Almeida et al. [36]2023; Sharma et al. [37]2024). In
   addition, C. pepo is a representative species to study gene duplication
   events on this family of crops (Montero‐Pau et al. [38]2018).

   Arsenic (As) is a nonessential hazardous metalloid for life, with
   toxicity levels determined by its oxidation state (El‐Ghiaty and
   El‐Kadi [39]2023; Preetha et al. [40]2023). As (V) is the most abundant
   form (Biswas and Ganesan [41]2024), particularly in aerobic soils where
   cucurbits are cultivated, followed by As (III), which is 60 times more
   toxic (Chen et al. [42]2021). In agriculture, the anthropogenic sources
   of As include pesticides and irrigation water (Mukherjee
   et al. [43]2024), whereas the geogenic incidence resides in important
   agricultural places (Zhang et al. [44]2021). For example, Mexico is an
   important producer of cucurbits, but As concentrations can reach
   naturally up to 35 times the World Health Organization's (WHO) maximum
   permissible recommendation of 10 μg L^−1 in water (Mahlknecht
   et al. [45]2023; Osuna‐Martínez et al. [46]2021).

   Cucurbitacins are a diverse group of highly oxygenated triterpenoid
   characteristic of species in the Cucurbitaceae family. Cucurbitacins
   are categorized into 17 groups (A‐T) where cucurbitacins B, D, E, and I
   are the most widespread in the Cucurbitaceae family (Almeida
   et al. [47]2022; Liu et al. [48]2024). Besides their biological
   function as defense compounds against herbivory and infectious
   microorganisms (Dong et al. [49]2021), their production is also highly
   valued for medicinal purposes (Maja et al. [50]2022). Interestingly,
   abiotic stress, for example, cold and drought, induce cucurbitacin
   biosynthesis in Lagenaria siceraria , Cucumis sativus, and Luffa
   acutangula , upregulating genes in the cucurbitacin‐biosynthesis
   pathway (Liu et al. [51]2019; Mkhize et al. [52]2023; Shang
   et al. [53]2014). Abiotic stresses including nutrient deficiency and
   heavy metals have been shown to influence the primary metabolome of
   cucurbits (Lee et al. [54]2024; Sun et al. [55]2020; Zhao
   et al. [56]2016). However, the targeted effects of metalloids (e.g.,
   As) on the metabolome, cucurbitacin content, and the cucurbitacin
   biosynthesis pathway have not been investigated in this plant family
   (Martínez‐Castillo et al. [57]2022).

   As exposure in plants results in evident phenotypic effects such as
   reduced size, decreased biomass, and lower growth ratio (Nabi
   et al. [58]2021), which are linked to physiological and molecular
   changes including oxidative stress and protein folding disruption
   caused by the depletion of antioxidants, for example, glutathione (GSH)
   and phytochelatins (Rajendran et al. [59]2024). Plants have evolved
   specialized proteins to regulate the uptake, translocation, and
   homeostasis of heavy metals (Li et al. [60]2020), and it would follow
   the same for metalloids as some are essential for plant growth (Kesawat
   et al. [61]2020). HMA‐domain‐containing proteins share a conserved
   30‐amino acid domain containing a crucial CxxC, where two cysteine
   residues facilitate the binding of essential and nonessential metals
   such as copper (Cu) and cadmium (Cd) (Takahashi et al. [62]2012; Wong
   et al. [63]2009). While As has a high affinity for cysteine‐rich
   peptides (Vergara‐Gerónimo et al. [64]2021), the role of the HMA
   superfamily in response to As remains poorly understood. In a previous
   study, we identified that members of the P[1B]‐ATPase clade, a
   subfamily of HMA‐domain‐containing proteins, differentially expressed
   when C. pepo was grown in the presence of As (Flores‐Iga
   et al. [65]2023), but other members of the HMA superfamily, such as
   heavy metal isoprenylated plant proteins (HIPP), were not studied.
   Moreover, heavy metalloid uptake, translocation, and extrusion are key
   mechanism for tolerance that warrant further study in Cucurbits,
   especially regarding As. For instance, Nodulin 26‐like transporters are
   recognized to participate in As extrusion and translocation, and ABCC
   transporters participate in the subcellular sequestration of As in
   vacuoles (Kamiya et al. [66]2009; Song et al. [67]2010).

   In this paper, we mechanistically explored the effects of As exposure
   on the phenotype and cucurbitacin content in the root system of C. pepo
   , and employed RNA‐seq to gain further insights into the molecular
   responses of C. pepo roots to As. While at low As concentrations the
   plant growth is not affected, cucurbitacin concentration decreases,
   which might negatively impact herbivory defense. Conversely, at high As
   concentrations, growth is diminished but cucurbitacin production is
   induced along with activation of oxidative and abiotic‐stress‐coping
   mechanisms in a concentration‐dependent manner. Additionally, we
   curated orthologs of the HMA‐domain‐containing superfamily and related
   As transporters to further explore characterized mechanisms employed to
   cope with As stress in cucurbits. Our examination of the underlying
   effects of As in the primary root system of C. pepo suggested
   preliminary insights that could eventually contribute to the
   development of healthier Cucurbitaceae crops in As‐contaminated
   environments.

2. Methods

2.1. C. pepo Growth and Root Phenotype

   Cucurbita pepo var. golden glory seeds were used, and four subsequent
   steps were employed for seed sterilization. First was a 20‐min wash in
   1% Mild Cream Soap (ABENA, Denmark), then a 1‐min rinse with 70%
   ethanol with a further 15‐min soak in 1.5% NaClO. Finally, the seeds
   were rinsed with sterile water and placed on a plate for surface
   drying. Because As(V) is the most prevalent form of As present in
   aerobic soils (Panda et al. [68]2010; Khalid et al. [69]2017) in which
   C. pepo would be cultivated, Murashige and Skoog media, with 0, 50,
   100, and 200 μM of As(V) ion from Na[2]HAsO[4] (Sigma‐Aldrich, St
   Louis, MO, USA), was used for C. pepo growth in square culture vessels.
   Growth conditions were light at 79 μE, 25°C, 80% humidity, and 16‐h
   photoperiod in Fitotron type SGC120‐H (Weiss Technik UK Ltd., Epinal
   Way, Loughborough, UK).

   Root phenotypes including length, area, diameter, volume and number of
   tips, forks, and crossing were analyzed in WinRhizo
   ([70]http://regentinstruments.com).

2.2. Metabolic Profile Analysis by LC–MS

   Frozen powdered materials (100 mg) of 2‐week‐old C. pepo var. golden
   glory roots grown under ½ MS agar supplemented with 0, 50, 100, and
   200 μM As were extracted (in triplicates) with 600 μL of methanol
   containing 50 μM of protopanaxatriol (Extrasynthese, Lyon, France) as
   internal standard. Samples were ultrasonicated for 20 min at room
   temperature, and further centrifuged for 10 min at 13,000 rpm. The
   supernatant was filtered employing 0.22 μm centrifugal filters
   (UFC30GV, Merck Millipore). Cucurbitacin quantification was performed
   as previously reported (Almeida et al. [71]2022). We focused
   exclusively on roots because they are the site of cucurbitacin
   production at the seedling stage, and it is reported that cucurbitacin
   biosynthetic genes are not expressed in cotyledons (Almeida
   et al. [72]2022; Brzozowski et al. [73]2020).

2.3. RNA‐Sequencing and Analysis

   Total RNA from 12 samples: three biological replicates of 2‐week‐old C.
   pepo var. golden glory roots grown under ½ MS agar supplemented with 0,
   50, 100, and 200 μM As, were isolated with the Spectrum Plant Total RNA
   Kit (Sigma‐Aldrich, St Louis, MO, USA) and treated with RQ1 RNase‐free
   DNase I treatment (Promega, Madison, WI, USA). RNA‐seq library
   preparation of RNA samples involved NEBNext Ultra II RNA Library Prep
   Kit for Illumina with NEBNext Poly (A) mRNA Magnetic Isolation Module
   (NEB, Ipswich, MA, USA). Quality and quantity of total RNA samples and
   their generated libraries were subjected to Bioanalyzer 2100 (Agilent,
   Santa Clara, CA, USA) and Qubit 3.0 fluorimeter (Invitrogen, Waltham,
   MA, USA). Paired‐end sequencing (2 × 75 bp) of libraries were carried
   out with the Illumina NextSeq 500 platform (Illumina, San Diego, CA,
   USA).

   The bcl2fastq tool (Illumina, San Diego, CA, USA) was used to convert
   the BCL files generated from NextSeq 500 to FASTQ files with 2 × 75 bp
   reads. Quality of reads was assessed using the fastqc tool
   ([74]https://www.bioinformatics.babraham.ac.uk/projects/fastqc).
   Adapters and low‐quality reads (Phred score QV < 30) were removed with
   Trimmomatic (Bolger et al. [75]2014). The quality‐filtered reads were
   mapped to the C. pepo genome
   ([76]http://cucurbitgenomics.org/v2/ftp/genome/Cucurbita_pepo/),
   excluding chromosome 0, with STAR RNA‐Seq aligner
   ([77]https://github.com/alexdobin/STAR). BAM and GFF were used to
   generate the count table with HTSeq
   ([78]https://github.com/htseq/htseq). Fragments per kilobase per
   million mapped reads (FPKMs) were calculated with the count table and
   length of each gene (Corchete et al. [79]2020). To identify
   differentially expressed genes (DEGs) between control and 50, 100, and
   200 μM As treatments, DESeq2 was employed (Love et al. [80]2014). DEGs
   with a minimum log2fold change ≥ 1 and an adjusted p value ≤ 0.05 were
   further considered. PyWGCNA (Rezaie et al. [81]2023) was used to
   identify co‐expressed gene clusters through weighted gene co‐expression
   network analysis. Gene annotation, classification, gene ontology (GO),
   and pathway enrichment were determined using the CuGenDB tools
   ([82]http://cucurbitgenomics.org/).

2.4. Identification, Phylogenetics, and Structural Analysis of HMA
Superfamily Members and As Transporters in C. pepo

   To identify HMA‐domain‐containing superfamily, aquaporin and ABCC
   transporters in C. pepo , characterized Arabidopsis protein sequences
   (Kamiya et al. [83]2009; Li et al. [84]2020; Song et al. [85]2010),
   were retrieved from TAIR ([86]https://www.arabidopsis.org/), and
   subjected to BLAST in the Cucurbit Genomics Database v2 (CuGenDBv2;
   [87]http://cucurbitgenomics.org/v2) to collect C. pepo protein
   sequences as well as Cucumis melo , Citrullus lanatus , and Cucurbita
   maxima when mentioned. To ensure the presence of conserved domains
   (HMA, ATPase, Hydrolase, MIP, and ABC), the conserved domain database
   (CDD) tool from NCBI was used
   ([88]https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) with the
   Pfam database as reference. In the case of HMA‐domain‐containing
   proteins, we named identified sequences according to chromosomal
   number. For identified Aquaporins and ABCC transporters, we used the
   orthologous name. MEGA11 (Tamura et al. [89]2021) was used for
   phylogenetic analysis. Multiple sequence alignments were performed
   using the Muscle (Edgar [90]2004) algorithm. Phylogenetic tree was
   built with the maximum likelihood method with the Jones–Taylor–Thronton
   (JTT) matrix–based model, using partial deletions and gamma distributed
   rates, further tested with bootstrap analysis (1000 replicates).
   Phylogenetic trees were exported as newick format for display in iTOL
   version 6.9 ([91]https://itol.embl.de/). TBtools (Chen et al. [92]2023)
   was used to display gene and protein structure as well as to perform
   Ka/Ks analysis. PlantCARE database
   ([93]https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was
   used to survey cis‐acting elements in 2000 bp upstream gene promoter
   sequences.

2.5. RT‐qPCR

   Equivalent concentrations of total RNA used for sequencing were
   subjected to first‐strand cDNA synthesis with SuperScript™ IV
   First‐Strand Synthesis System (Invitrogen, Waltham, MA, USA). Primers
   from selected genes (Table [94]S1) were designed with the CLC Genomics
   workbench software (QIAGEN, Aarhus, Denmark). The actin gene was used
   as endogenous control to normalize the relative gene expression. The
   reaction mixture contained 1 μL of cDNA template, 2 μL of forward and
   reverse primer at 10 μM, 10 μL of PowerTrack SYBR Green Master Mix
   (Applied Biosystems, Waltham, MA, USA), and 7 μL of nuclease‐free water
   for a final volume of 20 μL. Three technical replicates were assessed
   per each biological replicate. The qPCR was run in a QuantStudio 6 Pro
   (Applied Biosystems, Waltham, MA, USA). Melting curves were assessed to
   confirm the presence of a single amplicon and the absence of primer
   dimers. Ct values were further analyzed with the 2−ΔΔCt method (Livak
   and Schmittgen [95]2001).

2.6. Statistical Analysis

   Statistical analysis for phenotyping, cucurbitacins content and RT‐qPCR
   was carried out using R version 4.2.0 ([96]http://www.r‐project.org/).
   Normality was assessed using the Shapiro–Wilk test, and means were
   compared using a one‐way ANOVA with Tukey post hoc test for group
   comparison. The principal component analyses (PCA) were carried out
   with metabolomics profile (Table [97]S2), and RNA‐seq count table
   (Table [98]S3), to build a DESeq object in the DESeq package
   ([99]https://bioconductor.org/packages/DESeq2/).

3. Results

3.1. As Stress Deteriorates Root Growth and Modulates Cucurbitacin Production
in C. pepo

   In order to understand the detrimental effects of As on C. pepo , we
   analyzed the root phenotype and cucurbitacin production in the roots of
   seedlings. High As concentrations primarily affect root development
   (Figure [100]1). The root length, projected and surface area, root
   volume, tips, and forks were significantly decreased with increasing As
   concentrations, especially in the 100‐ and 200‐μM treatments (p
   value < 0.001) (Table [101]S4). Whereas in plants grown at 50 μM, the
   overall root phenotype did not differ significantly from those of
   control plants, but there was an increase in the number of crossings
   and forks 103% (p value < 0.05). Finally, the average root diameter
   increased by 17% in the 200 μM As treatment (p value < 0.001); this was
   the only root phenotype showing a significant increase due to high
   levels of As.

FIGURE 1.

   FIGURE 1
   [102]Open in a new tab

   High arsenic concentrations are detrimental for Cucurbita pepo root
   growth. (a–d) Representative pictures of Cucurbita pepo roots grown
   under different Arsenic concentrations for 2 weeks. Differences between
   control and As treatments on (e) primary root length, (f) root
   projected area, (g) root surface area, (h) root diameter, (i) root
   volume, (j) number of root tips, (k) number of root forks, and (l)
   number of root crossings. Data are means of the biological replicates
   with SEM represented as error bars. “ns,” “*,” “**,” and “***” indicate
   statistical significance between control and treatments by one‐way
   ANOVA with post hoc Tukey HSD at the levels of p values ≥ 0.05, ≤ 0.05,
   ≤ 0.01, and ≤ 0.001, respectively.

   In terms of C. pepo root metabolism, it can be observed that their
   metabolic profile is affected by varying As concentrations, as
   summarized by the principal component analysis (PCA) generated from the
   LC–MS analysis (Figure [103]2a). The 200‐μM treatment deviates most
   from the control, primarily along the first principal component (PC1),
   which accounts for 64.8% of the variance. In the particular case of
   cucurbitacins, there is a sharp decrease of cucurbitacin B (p
   value < 0.001), E (p value < 0.05), and I (p value < 0.01) content in
   the 50‐μM treatment (Figure [104]2b). In the case of cucurbitacin B,
   the 100‐μM As treatment also showed lower content than controls (p
   value < 0.01), although the decrease was less pronounced than at the
   50‐μM treatment. Conversely, in the 200‐μM As treatment, the levels of
   all three evaluated cucurbitacins were not significantly different when
   compared with those in controls. These observations suggest As exposure
   influences C. pepo root metabolism, showing the greatest shift at
   200 μM, while cucurbitacin content is affected only at low As levels.

FIGURE 2.

   FIGURE 2
   [105]Open in a new tab

   Arsenic affects the metabolic profile of Cucurbita pepo roots
   modulating their cucurbitacin content. (a) Distribution of individual
   roots exposed to different concentrations of Arsenic in the space of
   the first two principal components of the principal component analysis
   on LC–MS analysis. (b) Cucurbitacin B, E, and I concentrations in 50,
   100, and 200 μM As treatments. Data are means of the biological
   replicates with their SEM represented as error bars. “ns,” “*,”“**,”
   and “***” indicate statistical significance between control and
   treatments by one‐way ANOVA with post hoc Tukey HSD at the levels of p
   values ≥ 0.05, ≤ 0.05, ≤ 0.01, and ≤ 0.001, respectively.

3.2. Expression Levels of Cucurbitacin Pathway Genes Do Not Explain the
Observed Cucurbitacin Content Regulation in As Treatments

   The observed differences in the C. pepo root metabolism prompted us to
   also investigate the transcriptomic profile, the molecular mechanisms
   underlying the coping responses to As at different concentrations, as
   well as the associated phenotypic changes and cucurbitacin modulation
   observed in C. pepo under As stress, by conducting RNA‐sequencing and
   confirming the expression of specific genes using RT‐qPCR. Similar to
   the metabolic profile, a PCA reveals As influences the C. pepo
   transcriptome in a concentration‐dependent manner (Figure [106]3a). The
   transcriptomic profile shows distinct responses depending on the
   concentration, with the 200‐μM As treatment diverging the most from the
   control, primarily along PC1 (77% variance). Additionally, the
   proximity of the 50 and 100‐μM treatment clusters suggests similar
   transcriptomic profiles at these concentrations. Moreover, we
   identified common top upregulated and downregulated genes across As
   treatments, along with a clear clustering pattern of highly variable
   genes based on the presence or absence of As (Figure [107]S1). Notably,
   upregulated genes primarily encode stress‐related proteins, such as
   members of the heat‐shock family and MYB transcription factors, while
   downregulated genes are associated with growth and pathogenesis‐related
   responses.

FIGURE 3.

   FIGURE 3
   [108]Open in a new tab

   Arsenic impact on the transcriptomic responses of Cucurbita pepo . (a)
   Distribution of individual roots exposed to different concentrations of
   arsenic in the space of the first two principal components of the
   principal component analysis on gene expression profiles derived from
   RNA‐seq analysis; (b) expression of cucurbitacin pathway‐related genes;
   (c) gene ontology (GO) enrichment analysis of biological process (BP),
   cellular component (CC), and molecular function (MF); and (d) pathway
   enrichment analysis of the common 592 DEGs (cut‐off p value ≤ 0.05).
   Data are means of the biological replicates with their SEM represented
   as error bars. “ns,” “*,”“**,”and “***” indicate statistical
   significance between control and treatments by one‐way ANOVA with post
   hoc Tukey HSD at the levels of p values ≥ 0.05, ≤ 0.05, ≤ 0.01, and
   ≤ 0.001, respectively.

   Remarkably, we observed an upregulation in the expression of
   cucurbitacin‐biosynthetic pathway genes at 200 μM As, including
   cucurbitadienol synthase, CpCPQ (Cp4.1LG17g04420), characterized
   cytochromes P450, CpCYP81BF1 (Cp4.1LG12g10100), and CpCYP81BB1
   (Cp4.1LG17g09460), but not for the acyl transferase CpACT1
   (Cp4.1LG17g09410) (Figure [109]3b). Nevertheless, expression of these
   genes was not significantly different at the 50‐ and 100‐μM As
   treatments when compared with controls, explaining why these genes did
   not pass the thresholds of the downstream RNA‐seq analysis to be
   classified as differentially expressed genes (DEGs) and therefore are
   not listed in the DEGs (Table [110]S5–[111]S9). Although the expression
   profiles of cucurbitacin pathway genes do not fully explain the changes
   in cucurbitacin content observed across the treatments, the metabolic
   profiles suggest that As exposure modulates the root metabolism, which
   could indirectly cause the changes in cucurbitacin content.

   Furthermore, gene ontology (Figure [112]3c; Table [113]S10) and pathway
   (Figure [114]3d; Table [115]S11) analyses were conducted to gain
   insights into the mechanisms affected by As. Responses to chemicals,
   toxic substances, and other stresses were the most enriched biological
   processes (BP). Accordingly, the glutathione‐mediated detoxification
   pathway in response to oxidative stress was commonly enriched among the
   upregulated genes in the treatments with As exposure. The extracellular
   region was the only cellular component (CC) enriched, while for
   molecular function (MF), DNA‐binding transcription factor activity and
   glutathione transferase activity were the most enriched. On the other
   hand, we found downregulated common peroxidase genes involved in lignan
   biosynthesis such as justidicin and matairesinol. Additionally, the
   correlation between gene expression and As concentration was assessed
   using WGCNA (Figure [116]S2). After assigning genes into modules, we
   found that fundamental biosynthesis and gene regulation processes, as
   well as the glutathione‐mediated detoxification II pathway, were
   coregulated with the 100 and 200‐μM As treatments, whereas signaling
   and energy process and the urate degradation to allantoin pathway were
   coregulated with the 50‐μM As treatment (Table [117]S12). These results
   suggest that exposure to different concentrations of As leads to the
   upregulation of detoxification and stress tolerance mechanisms, while
   also affecting fundamental processes in a concentration‐specific
   pattern, which correlates with its phenotypic traits.

3.3. Genes Encoding HMA‐Domain‐Containing Proteins and As‐Transporters Are
Differentially Expressed Upon As Stress

   We studied in more detail genes encoding members in the HMA‐domain
   superfamily, aquaporin and ABCC transporter families due to their
   association to metalloid tolerance, As transport, and their GO terms
   matching with those enriched in this study (response to oxidative
   stress, inorganic substance, and detoxification). Moreover, we aimed to
   curate the orthologs of these known metal (loid)‐tolerance related
   genes in this non‐model crop to enable further in‐depth study within
   the Cucurbitaceae family. Importantly, we did not observe the
   modulation of As(V) transporters. For instance, PHT1:1 (AT5G43350) is
   known to participate in As(V) entry in Arabidopsis thaliana (Catarecha
   et al. [118]2007). In our dataset, the putative homologs
   Cp4.1LG03g10070, Cp4.1LG12g01500, and Cp4.1LG15g00730 were not found to
   be differentially regulated across conditions (Table [119]S5–[120]S9).

   Members of the HMA‐domain superfamily fall into different subfamilies,
   including HIPP (Heavy metal‐associated Isoprenylated Plant Protein),
   HPP (Heavy metal‐associated Plant Protein), P[1B]‐ATPases, and ATX1
   (antioxidant protein 1) differing in protein structure and
   functionality (De Abreu‐Neto et al. [121]2013; Li et al. [122]2020; Ye
   et al. [123]2022). Using a comparative genomics approach, we identified
   77 genes encoding HMA‐domain‐containing proteins, 53 HIPP, 14
   P[1B]‐type ATPases, nine HPP, and one ATX1 across the genome of C. pepo
   (Table [124]S13; Figure [125]S3). We further studied the relationships
   of C. pepo orthologs to characterized HMA‐domain‐containing proteins in
   A. thaliana through phylogenetic analysis of amino acid sequences.

   The HIPP family in A. thaliana is known to be involved in the coping
   mechanisms against different metals and metalloids (Gao
   et al. [126]2009; Tehseen et al. [127]2010; Wei et al. [128]2023).
   After curation, the resulting phylogenetic tree of the HIPP protein
   family revealed CpHIPP42 (Cp4.1LG15g05950) and CpHIPP17
   (Cp4.1LG04g04920) as the C. pepo orthologs of AtHIPP26 and CpHIPP36
   (Cp4.1LG13g03670) as an ortholog of AtHIPP22 (Figure [129]4a).
   Construction of a phylogenetic tree including amino acid sequences of
   HIPP homologs from other Cucurbitaceae species (Figure [130]S4)
   suggests CpHIPP17 is an ortholog that arose from a duplication event
   unique to the Cucurbiteae tribe (squash, pumpkin) as no sequences from
   members of the Benincaseae tribe (cucumber, melon and watermelon) or A.
   thaliana group with it. Although being a segmental duplication from
   Chromosome 4 to 15, the Ka/Ks ratio of this event is less than 1
   (0.048), which implies purifying selection, also shown in the protein
   structure. Nevertheless, the promoter of both genes differ, as a survey
   of their cis‐acting elements revealed that salicylic acid
   responsiveness and specific MYB binding sites were found only for
   CpHIPP42 (Table [131]S14). Accordingly, gene expression of these HIPP
   paralogs is different in response to As. Whereas the gene expression of
   CpHIPP17 significantly decreased in the 100 μM As treatment compared
   with controls, for CpHIPP42, the expression increased six‐fold in the
   200‐μM As treatment (Figure [132]4a). As for CpHIPP36, its expression
   was upregulated as As concentrations increased, reaching up to a
   40‐fold increase in the 200‐μM As treatment. Despite these differences
   in regulation, the structure of these three HIPP proteins share similar
   features (Figure [133]4a).

FIGURE 4.

   FIGURE 4
   [134]Open in a new tab

   Expression of HMA‐containing‐protein genes is upregulated under As
   stress in Cucurbita pepo . (a) (top) Maximum‐likelihood tree
   constructed on deduced amino acid sequences of 13 plant HIPPs aligned
   by ClustalW spanning 157 positions using the MEGAX program. The
   statistical significance of each node was tested by the bootstrap
   method using 1000 iterations. Accession numbers and representative
   names are given in the figure. (bottom‐left) Relative expression
   profiles of C. pepo HIPP‐protein genes under different arsenic
   concentrations. (bottom‐right) Protein and gene structure of C. pepo
   HIPP proteins. (b) (top) Maximum‐likelihood tree constructed on deduced
   amino acid sequences of 13 plant HMA‐containing proteins aligned by
   ClustalW spanning 717 positions using the MEGAX program. The
   statistical significance of each node was tested by the bootstrap
   method using 1000 iterations. Accession numbers and representative
   names are given in the figure. (bottom‐left) Relative expression
   profiles of C. pepo HMA‐containing‐protein genes under different
   arsenic concentrations. (bottom‐right) Protein and gene structure of C.
   pepo HMA‐containing proteins. Gene and protein structural features are
   shown in the right side of the expression bar plot for each gene. The
   different colored rectangles represent CDS, UTRs, and structural
   domains mentioned in the color legend. The below scale indicates the
   length in kb and aa, accordingly. Data are means of the biological
   replicates with their SEM represented as error bars. “ns,” “*,” “**,”
   and “***” indicate statistical significance between control and
   treatments by one‐way ANOVA with post hoc Tukey HSD at the levels of p
   values ≥ 0.05, ≤ 0.05, ≤ 0.01, and ≤ 0.001, respectively.

   Next, we studied members of the P[1B]‐ATPase family as they also
   contain HMA domain(s) in their protein structure (Axelsen and
   Palmgren [135]2001). Phylogenetic analysis of a curated set of amino
   acid sequences including functionally characterized P[1B]‐ATPases of A.
   thaliana shows that CpHMA2 (Cp4.1LG05g03530) and CpHMA8
   (Cp4.1LG10g10980) group in the clade of AtHMA5, while CpHMA11 groups
   with that of AtHMA6. Whereas CpHMA8 appears to be the true ortholog of
   AtHMA5, CpHMA2 seems to have arisen from a duplication event in the
   common ancestor of the Benincaseae and Cucurbiteae tribes
   (Figure [136]4b), which was consistent with our previous study
   (Flores‐Iga et al. [137]2023). AtHMA5, AtHMA7, CpHMA8, and CpHMA2 all
   have a similar structure containing three HMA domains in the first 300
   aa of the N‐terminal, while AtHMA6 and CpHMA11 contain a single HMA
   domain in their N‐terminal. In regard to gene expression, in general,
   the three C. pepo P[1B]‐ATPases are upregulated in the presence of As,
   but their profiles vary significantly (Figure [138]4b). Specifically,
   CpHMA8 increases more than five‐fold in all As treatments but being
   most highly expressed at the 50‐μM As treatment. Interestingly, despite
   belonging to different phylogenetic clades, CpHMA2 and CpHMA11 follow a
   similar trend of increased gene expression with ascending As
   concentration in the 50‐ and 100‐μM treatments but having a decrease at
   200 μM As.

   Furthermore, we surveyed aquaporins in the NIP group, which are key As
   transporters involved in its extrusion and translocation (Kamiya
   et al. [139]2009). There are two proteins in C. pepo , named here as
   CpNIP1:1 (Cp4.1LG19g06570) and CpNIP1:2 (Cp4.1LG10g05550), that group
   in a clade with AtNIP1:1 and AtNIP1:2 (Figure [140]5a). CpNIP1:1
   appears to be the result of a duplication event exclusive to the
   Cucurbiteae tribe, whereas CpNIP1:2 groups with other aquaporins of the
   Benincaseae tribe (Figure [141]S5). Interestingly, expression of
   CpNIP1:1 remains unaffected at low As concentrations but increases
   10‐fold in the 200 μM As treatment, while expression of CpNIP1:2
   decreases nearly half under 50 and 100 μM of As (Figure [142]5a).

FIGURE 5.

   FIGURE 5
   [143]Open in a new tab

   Orthologous transporters are regulated on Cucurbita pepo in high As
   conditions. (a) (left) Maximum‐likelihood tree constructed on deduced
   amino acid sequences of eight NIPs aligned by ClustalW spanning 218
   positions using the MEGAX program. The statistical significance of each
   node was tested by the bootstrap method using 1000 iterations.
   Accession numbers and representative names are given in the figure.
   (middle) Relative expression profiles of C. pepo NIP‐protein‐encoding
   genes under different arsenic concentrations. (right) Protein and gene
   structure of C. pepo NIP proteins. (b) (left) Maximum‐likelihood tree
   constructed on deduced amino acid sequences of 10 ABCC transporters
   aligned by ClustalW spanning 936 positions using the MEGAX program. The
   statistical significance of each node was tested by the bootstrap
   method using 1000 iterations. Accession numbers and representative
   names are given in the figure. (middle) Relative expression profiles of
   C. pepo ABCC transporter genes under different Arsenic concentrations.
   (right) Protein and gene structure of C. pepo ABCC transporters. The
   different colored rectangles represent CDS, UTRs, and structural
   domains according to the color legend. The below scale indicates the
   length in kb and aa, accordingly. Data are means of the biological
   replicates with their SEM represented as error bars. “ns,” “*,” “**,”
   and “***” indicate statistical significance between control and
   treatments by one‐way ANOVA with post hoc Tukey HSD at the levels of p
   values ≥ 0.05, ≤ 0.05, ≤ 0.01, and ≤ 0.001, respectively.

   Finally, we studied ABCC transporters because they have a functionally
   characterized role in transporting As‐complexes in A. thaliana (Song
   et al. [144]2010) and are recognized to participate in the subcellular
   sequestration of As (Zhang et al. [145]2021). Our focus was on two
   identified C. pepo ABCC transporters, namely, CpABCC1 (Cp4.1LG15g02900)
   and CpABCC2 (Cp4.1LG19g05820). These C. pepo ABCC transporters
   exhibited distinct expression patterns across As treatments
   (Figure [146]5b), with a significant up‐regulation at 200 μM: with
   CpABCC1 increased 12‐fold, whereas CpABCC2 showed a three‐fold
   (Figure [147]5b).

   In summary, the phylogenetic analysis showed that evolution of the
   HIPP, P1[B]‐ATPase, and aquaporin families in the Cucurbiteae tribe,
   specifically C. pepo , appears to have diverged from that of A.
   thaliana and members of the Benincaseae tribe, possibly as a result of
   duplication events, while the ABCC transporter family remains
   relatively conserved.

4. Discussion

4.1. As Exposure Is Detrimental to Plant Growth and Plant Defense of C. pepo

   Previous studies on model and non‐model plants have established plant
   responses to As stress, a hazardous nonessential metalloid
   (Gracia‐Rodriguez et al. [148]2024). Before this study, little was
   known about the molecular responses of cucurbits to As and how exposure
   to this metalloid influenced C. pepo metabolism. In the Cucurbitaceae
   family, exposure to As significantly diminishes root and shoot length
   in Cucurbita moschata and C. sativus (Hong et al. [149]2011; Mushtaq
   et al. [150]2020). In our study, we noticed As exposure negatively
   impacts the growth of the root system in C. pepo , particularly at 100
   and 200 μM, indicating its vulnerability at high As levels that have
   been previously studied in A. thaliana (Navarro et al. [151]2024;
   Peña‐Garcia et al. [152]2021; Zou et al. [153]2024). Damage to the root
   system likely triggers a cascade of negative effects on plant health.
   For instance, reduced root function hinders water and nutrient uptake,
   ultimately impacting plant metabolism (Kosakivska et al. [154]2021).
   Although members of the Benincaseae tribe, C. melo and C. lanatus
   exhibited a species‐specific pattern of As bioaccumulation in response
   to As stress (Pineda‐Chacón and Alarcón‐Herrera [155]2016), both
   species presented significant biomass decrease and plant weakening. As
   the root diameter is the only trait that increased in size in response
   to As exposure, it warrants further investigations.

   Cucurbitacins are bitter secondary metabolites recognized for medicinal
   purposes (Abdelmonsef et al. [156]2024); within the plant system, they
   serve as a defense mechanism against biotic stress such as herbivores
   and microbial pathogens (Brzozowski et al. [157]2020). Our LC–MS
   analysis indicates that As exposure modulates the metabolome of C. pepo
   , particularly the cucurbitacin profile, dramatically decreasing the
   concentrations of cucurbitacin B, E, and I in the 50 μM As treatment.
   Interestingly, the expression profile of the cucurbitacin biosynthetic
   pathway genes was only modulated in the 200‐μM As, with the exception
   of CpACT1, which catalyzes the final step of cucurbitacin E
   biosynthesis (Ma et al. [158]2022). Our observations suggest that,
   under low As concentrations, cucurbitacin levels are indirectly
   influenced by changes in metabolic pathways or regulatory networks
   rather than by direct transcriptional regulation. Although no
   significant changes in cucurbitacin content were detected at 200 μM,
   the increased gene expression at this concentration suggests a
   metabolomic shift that redirects isoprenoid precursors and/or
   intermediates towards cucurbitacin biosynthesis. According to our PCA,
   this shift was most pronounced at this high As concentration,
   warranting further investigation. A comprehensive untargeted
   metabolomic analysis using MS/MS data would be crucial to gain deeper
   insights into the metabolic changes affecting cucurbitacin
   accumulation. To our knowledge, this is the first report linking the
   levels of defense compounds, such as cucurbitacins, in Cucurbitacea
   crops to As exposure.

   Analysis of our RNA‐seq data showed downregulation of genes associated
   with immune and stress signaling pathway genes, providing further
   evidence that metalloid exposure, specifically As, not only affects
   oxidative stress pathways but also disrupts C. pepo biotic defense
   mechanisms, weakening the plant system (Labidi et al. [159]2023; Labidi
   et al. [160]2021). Interestingly, we found degradation of Allantoin
   correlated with the 50‐μM As treatment in our WGCNA. Allantoine is a
   nitrogen‐rich compound that participates in the resistance mechanism of
   C. melo to powdery mildew, a major fungal disease (Dun
   et al. [161]2024). These results provide insights into how As stress
   influences plant growth and plant defense, potentially contributing to
   an increased susceptibility of crops to infections.

   In addition, upregulated genes showed GO enrichment of toxic substance
   and oxidative stress responses in the biological process category
   across As treatments. Notably, genes related to the
   gluthathione‐mediated detoxification II pathway were enriched, which is
   well‐known for its role in metal and metalloid stress (Kumar and
   Trivedi [162]2018). When we assessed WGCNA to understand the
   co‐expression networking, we found regulatory processes governing
   transcription and biosynthetic pathways being disturbed, particularly
   in the 100‐ and 200‐μM treatments, which concur with previous results
   in A. thaliana (Abercrombie et al. [163]2008; Peña‐Garcia
   et al. [164]2021). Other genes with expression patterns worth
   mentioning are the top up‐regulated DEGS, specifically, three heat
   shock‐like proteins (Cp4.1LG08g07390, Cp4.1LG03g14960, Cp4.1LG01g11230)
   and a MYB transcription factor (Cp4.1LG02g16350). In Arabidopsis, genes
   in the MYB family have been linked to metalloid tolerance, regulating
   key players such as members of the HMA family (Cao et al. [165]2024;
   Chen et al. [166]2021), and we found MYB regulatory sites in the
   promoter of CpHIPP42, which was upregulated as concentrations
   increased; further studies could reveal a link between this MYB factor
   and HMA regulation or other factors in the As response. On the other
   hand, biotic stress signaling pathway genes such as jasmonic acid
   (Cp4.1LG18g02900, Cp4.1LG03g14930, Cp4.1LG11g01710) and pathogenesis
   related signaling (Cp4.1LG15g04120) were found as the most
   down‐regulated genes. Consistently, jasmonic acid biosynthesis genes
   have been found to be down‐regulated when C. pepo is exposed to trace
   heavy metals (Labidi et al. [167]2023).

4.2. HMA‐Domain‐Containing Proteins and As‐Related Genes Undergo Different
Genetic Expressions Under As Stress

   Our study applies comparative genomics to explore metalloid tolerance
   in C. pepo by examining mechanisms from A. thaliana and Oryza sativa ,
   specifically focusing on the HMA‐domain superfamily, whose members play
   key roles in the transport, binding, homeostasis, and detoxification of
   metals such as Cu, Ag, Pb, and Cd (Jamla et al. [168]2021; Li
   et al. [169]2020). In contrast to 55 and 46 members in A. thaliana and
   O. sativa , respectively (De Abreu‐Neto et al. [170]2013; Li
   et al. [171]2020), we identified 77 HMA‐domain‐containing superfamily
   members across the C. pepo genome. The known high affinity to thiol
   groups (Sun et al. [172]2023) made us suspect that As could interact
   with cysteine residues in the HMA‐domain of members in this superfamily
   and regulation of these genes could be influenced. Our prior analysis
   of the P[1B]‐ATPase family in Cucurbitaceae pointed this was a
   hypothesis worthwhile testing as we had observed gene expression
   differences in C. pepo under As stress, particularly in the Cu/Ag clade
   (Flores‐Iga et al. [173]2023). Indeed, initial observations on the gene
   expression of members of this superfamily in the RNA‐seq analysis
   suggested their regulation was different in response to As
   (Figure [174]S6), which was confirmed with RT‐qPCR (Figures [175]4,
   [176]5).

   In the HIPP family, we found that CpHIPP17 and CpHIPP42 are orthologs
   of AtHIPP26, and CpHIPP36 of AtHIPP22. Both AtHIPP26 and AtHIPP22 are
   known and required for heavy metal tolerance such as Cd (Manara
   et al. [177]2020; Tehseen et al. [178]2010). In this study, we observed
   a substantial expression of CpHIPP36 only, especially in the highest As
   concentration. Moreover, intron length and regulation of CpHIPP17 and
   CpHIPP42 is clearly different but the protein structure is not, such
   differences in transcription activity of duplicated genes are
   recognized (Panchy et al. [179]2016). In a survey of their cis‐acting
   elements, salicylic acid responsiveness elements were unique for the
   promoter region of CpHIPP42. Salicilyc acid is a ROS‐stimulated plant
   hormone that has been reported to activate antioxidant defense systems
   (Emamverdian et al. [180]2020).

   In the P1B‐ATPase family, CpHMA8 and CpHMA2 were found to be orthologs
   of AtHMA5 and AtHMA7, respectively. AtHMA5 participates in Cu
   detoxification in roots, while AtHMA7 is essential for Cu transport to
   ethylene receptor ETR1 (Andrés‐Colás et al. [181]2006; Puig
   et al. [182]2007; Migocka [183]2015). CpHMA11 is orthologue of AtHMA6,
   which contains a lone HMA domain and is involved in delivery of Cu to
   the chloroplastic compartment (Mayerhofer et al. [184]2016). Elevated
   transcription of CpHMA8 and CpHMA11 across As treatments suggests that
   gene structure, specifically HMA domain quantity, does not determine
   transcriptional response, but instead transcription may be influenced
   by functional factors (Abercrombie et al. [185]2008), which is also
   observable in different tissues at later developmental stages of C.
   pepo (Flores‐Iga et al. [186]2023). AtATX1 (AtCuSOD) is a Cu chaperone
   for superoxide dismutase (SOD) with a role of homeostasis under heavy
   metal stress, for example, Cu (Shin et al. [187]2012). Besides
   delivering Cu to Cu/Zn SOD, it participates in Cu trafficking for
   homeostasis networking with Cu/Ag members of the P[1B]‐ATPase subfamily
   (Zhang et al. [188]2018). Interestingly, expression levels of CpCuSOD
   were higher in all As treatments. Whether the higher expression of
   CpCuSOD is caused directly by As, P[1B]‐ATPase induction or indirectly
   due to ROS production by As remains to be determined. The genome‐wide
   identification of HMA members and gene expression analysis in this
   study narrow down the members that are likely activated in the primary
   stage due to As exposure. The HMA investigation needs further
   functional characterization to understand what are the triggering
   requirements of different HMA‐domain‐containing proteins to act
   synergistically towards heavy metal and metalloid tolerance and its
   related oxidative stress, which is proposed as an auxiliary mechanism
   for plants to cope with a wide range of metals and metalloid stressors
   (Rono et al. [189]2022).

   A common mechanism of As tolerance is the reduction of As (V) into a
   higher affinity and likely extrusive form, As (III) (Chao
   et al. [190]2014; Wang et al. [191]2018). Non‐enzymatic and enzymatic
   antioxidants are responsible for this mechanism, which is observed in
   our RNA‐seq analysis with the upregulation of glutathione
   detoxification pathway and glutathione metabolic processes. In this
   sense, aquaporins are recognized to participate in As exclusion and
   translocation, particularly of the As (III) form (Sun
   et al. [192]2024). Here, high expression levels of CpNIP1:2 were
   observed in the 200 As μM treatment, but not in CpNIP1:1. AtNIP1:1 and
   AtNIP1:2 are determinants of As transport and permeability across root
   cells, and higher expression of their encoding genes has been reported
   in earlier reports (Ji et al. [193]2017; Kamiya et al. [194]2009).
   Similarly, OsNIP1:1 participates in root‐shoot transport, and
   accumulation (Sun et al. [195]2018). Therefore, the studied aquaporins
   may play a supportive role in the plant's response to high cellular As
   (III) toxicity after the first As(V) to As (III) reduction step
   (Dhankher et al. [196]2006; Zou et al. [197]2024), potentially
   facilitating its release from the root system after prolonged
   intracellular accumulation.

   We ultimately investigated C. pepo orthologues of AtABCC1 and AtABCC2,
   which are well‐characterized vacuolar transporters that facilitate the
   cellular sequestration of As by transporting As (III)‐phytochelatin and
   As (III)‐GS complexes (Song et al. [198]2010). We observed that CpABCC1
   and CpABCC2 have sequence divergence in the leading ABC membrane
   superfamily domain, and their expression levels differ. Notably, only
   CpABCC2 exhibits high transcription rates, especially at 200 μM
   (Figure [199]5b), consistent with observations that the expression of
   these transporters is solely induced when chelate action and vacuolar
   As transport are coordinated (Briat [200]2010). Divergence in the
   function of CpABCC1 reflected by difference in its structure could be a
   factor contributing to the lack of transcription induction for its gene
   under As stress. To study the effect on gene activation, further
   investigation is needed to explore the regulatory mechanisms behind the
   transcriptional control of CpABCC2, such as upstream signaling pathways
   and cis‐binding regulators providing insights into how this plant
   orchestrates its stress response at the molecular level.

5. Conclusion

   Cucurbita pepo plants exposed to As stress exhibit root damage and
   altered cucurbitacin levels, possibly weakening the plant to biotic
   stress. Specifically, our GO enrichment analysis of DEGs suggests that
   to cope with As toxicity, the plant activates oxidative stress
   responses with regulatory and metabolomic disruption at differing As
   levels. This could include a possible increase of glutathione
   production and other thiol‐rich peptides. Furthermore, WGCNA analysis
   of As‐concentration‐specific gene co‐expression modules indicate that
   low As levels stimulate energy and signaling pathways, along with
   allantoin metabolism, while higher As concentrations intensify
   oxidative stress responses. This includes upregulation of MYB
   transcription factors and heat shock proteins, which may play a role in
   regulating As tolerance. In addition, by curating C. pepo orthologs of
   A. thaliana , we were able to gain insights into possible downstream
   mechanisms of As transport, extrusion, and sequestration, highlighting
   potential pathways for As detoxification.

   Collectively, these preliminary findings suggest a shift from pathogen
   defense towards metalloid stress management in As ‐polluted
   environments, providing a foundation for future research aimed at
   understanding the underlying mechanisms of the responses observed in in
   this study and developing As‐tolerant cucurbits in the future.

Author Contributions

   G.F.‐I.: writing – original draft preparation, writing – review and
   editing, visualization, software, investigation, formal analysis, data
   curation, conceptualization. C.L.‐O.: writing – review and editing,
   validation, formal Analysis, conceptualization. P.N.: software. P.N.:
   project administration, funding acquisition. U.K.R.: writing – review
   and editing, supervision, project administration, funding acquisition.
   N.B.: writing – review and editing, supervision, project
   administration, funding acquisition. A.A.: writing – original draft
   preparation, writing – review and editing, investigation, validation,
   supervision, software, project administration, methodology, funding
   acquisition, formal analysis, conceptualization.

Conflicts of Interest

   The authors declare no conflicts of interest.

Peer Review

   The peer review history for this article is available in the
   [201]Supporting Information for this article.

Supporting information

   Data S1. Peer review.
   [202]PLD3-9-e70074-s009.pdf^ (108.7KB, pdf)

   Table S1. Primer sequences of selected Cucurbita pepo genes evaluated
   by RT‐qPCR under arsenic stress.
   [203]PLD3-9-e70074-s010.xlsx^ (10.6KB, xlsx)

   Table S2. Metabolomic profile of Cucurbita pepo roots under different
   As treatments (raw data).
   [204]PLD3-9-e70074-s005.xlsx^ (32.1KB, xlsx)

   Table S3. RNA‐seq count table of Cucurbita pepo roots under different
   As treatments.
   [205]PLD3-9-e70074-s014.xlsx^ (1.9MB, xlsx)

   Table S4. Data of Cucurbita pepo phenotyping under different As
   concentrations.
   [206]PLD3-9-e70074-s015.xlsx^ (11.3KB, xlsx)

   Table S5. FPKM values of control and As treatment RNA‐seq samples from
   Cucurbita pepo .
   [207]PLD3-9-e70074-s003.xlsx^ (4.4MB, xlsx)

   Table S6. DEGs between control and 50 μM As treatment in Cucurbita pepo
   .
   [208]PLD3-9-e70074-s008.xlsx^ (125.6KB, xlsx)

   Table S7. DEGs between control and 100 μM As treatment in Cucurbita
   pepo .
   [209]PLD3-9-e70074-s012.xlsx^ (117.4KB, xlsx)

   Table S8. DEGs between control and 200 μM As treatment in Cucurbita
   pepo .
   [210]PLD3-9-e70074-s011.xlsx^ (130.7KB, xlsx)

   Table S9. Common DEGs between control and As treatments in Cucurbita
   pepo .
   [211]PLD3-9-e70074-s007.xlsx^ (48.4KB, xlsx)

   Table S10. enriched GO terms of common and specific DEGs in different
   As treatments.
   [212]PLD3-9-e70074-s001.xlsx^ (15.8KB, xlsx)

   Table S11. Enriched pathways of common and specific DEGs in different
   As treatments.
   [213]PLD3-9-e70074-s002.xlsx^ (10.7KB, xlsx)

   Table S12. Enrichment of genes found in modules from WGCNA in Cucurbita
   pepo to As stress.
   [214]PLD3-9-e70074-s004.xlsx^ (68.1KB, xlsx)

   Table S13. Identified HMA domain containing proteins in Cucurbita pepo
   to As stress and their log2(FPKM) expression.
   [215]PLD3-9-e70074-s006.xlsx^ (15.9KB, xlsx)

   Table S14. Cis‐acting elements present in CpHIPP17 and CpHIPP42.
   [216]PLD3-9-e70074-s016.xlsx^ (30.1KB, xlsx)

   Figure S1 Overview of RNA‐seq analysis in Cucurbita pepo under
   different arsenic treatments. Volcano plots highlighting top
   differentially expressed genes from 50 μM As (a), 100 μM As (b), and
   200 μM As (c) with thresholds: minimum log2 fold change = 1, p value
   ≤ 0.05 and lcshrink() function for low counts. Heatmap of the 25 common
   top and bottom genes with the highest expression variation (d).
   Dendograms explain the clustering pattern of samples and gene
   expression values. (e) Venn diagram showing the number of DEGs shared
   and unique for each of the three As treatments.

   Figure S2. WGCNA analysis of Cucurbita pepo exposed to As. WGCNA
   preprocessing results (a). In the left panel, scale‐free fit index. In
   the right panel, mean connectivity. The panels display the relationship
   and function with respect to the soft‐thresholding power (X‐axis),
   respectively. Clustering dendrogram of genes (b), with dissimilarity
   based on topological overlap, together with assigned merged module
   colors and the original module colors. Module–trait relationship
   heatmap generated from WGCNA for C. pepo ‐As concentration (c). Rows
   represent different As concentrations, and columns represent modules.
   Values in the squares indicate correlation module‐treatment according
   to Pearson correlations, p values are in parenthesis (null hypothesis;
   probability to find the current results if the correlation coefficients
   were zero). Color fade indicate positive (red) and negative (blue)
   correlation.

   Figure S3. Phylogenetic analysis and gene expression heatmap of the
   HMA‐containing proteins superfamily in Cucurbita pepo . HMA‐containing
   proteins were clustered with a Maximum‐likelihood dendogram. Heatmap
   was generated with normalized values of expression (Table S14). From
   inside to outside, expression of HMA‐containing proteins in 50, 100,
   and 200 μM As treatments are displayed.

   Figure S4. Maximum likelihood phylogenetic tree of orthologous HIPP
   members with CuSOD (ATX1) sequences used as outgroup. Arabidopsis
   thaliana, Cucurbita pepo, Cucurbita maxima, Cucumis sativus, Citrullus
   lanatus, and Cucumis melo sequences were used to understand the
   orthologs and duplication events between the Cucurbitae and Benincaseae
   tribes.

   Figure S5. Maximum likelihood phylogenetic tree of orthologous
   Aquaporin members with AtNIP5:1 and AtNIP6:1 used as outgroup.
   Arabidopsis thaliana , Solanum lycopersicum , Cucurbita pepo ,
   Cucurbita maxima , Citrullus lanatus , and Cucumis melo sequences were
   used to understand the orthologs, conservation and duplication events
   between the Cucurbitae and Benincaseae tribes.

   Figure S6. Gene expression RNA‐seq values for genes of interest.
   P1B‐ATPase (a), HIPP (b), aquaporin (c), and ABCC (d) families. Data
   are means of the biological replicates with their SEM represented as
   error bars. “ns,” “*,”“**,”and “***” indicate statistical significance
   between control and treatments by one‐way ANOVA with post hoc Tukey HSD
   at the level of p values ≥ 0.05, ≤ 0.05, ≤ 0.01, and ≤ 0.001,
   respectively.
   [217]PLD3-9-e70074-s013.pdf^ (2MB, pdf)

   Funding: The authors received no specific funding for this work.

   Gerardo Flores‐Iga and Aldo Almeida contributed equally to this work.

Data Availability Statement

   The RNA sequencing project is available in the NCBI database under the
   BioProject accession number PRJNA1177526. Scripts to analyze and
   display the data can be found at
   [218]https://github.com/gerapskx/Cucurbita‐pepo_Arsenic.

References