Abstract

   Cadmium (Cd) is among the world’s major health concerns, as it renders
   soils unsuitable and unsafe for food and feed production.
   Phytoremediation has the potential to remediate Cd-polluted soils, but
   efforts are still needed to develop a deep understanding of the
   processes underlying it. In this study, we performed a comprehensive
   analysis of the root response to Cd stress in A. thaliana, which can
   phytostabilize Cd, and in A. halleri, which is a Cd hyperaccumulator.
   Suitable RNA-seq data were analyzed by WGCNA to identify modules of
   co-expressed genes specifically associated with Cd presence. The
   results evidenced that the genes of the hyperaccumulator A. halleri
   mostly associated with the Cd presence are finely regulated (up- and
   downregulated) and related to a general response to chemical and other
   stimuli. Additionally, in the case of A. thaliana, which can
   phytostabilize metals, the genes upregulated during Cd stress are
   related to a general response to chemical and other stimuli, while
   downregulated genes are associated with functions which, affecting root
   growth and development, determine a deep modification of the organ both
   at the cellular and physiological levels. Furthermore, key genes of the
   Cd-associated modules were identified and confirmed by differentially
   expressed gene (DEG) detection and external knowledge. Together, key
   functions and genes shed light on differences and similarities among
   the strategies that the plants use to cope with Cd and may be
   considered as possible targets for future research.

   Keywords: heavy metal stress, cadmium stress, phytoremediation, root
   biology, Arabidopsis halleri, Arabidopsis thaliana

1. Introduction

   Cadmium (Cd) is a non-essential heavy metal and one of the most
   widespread pollutants in both terrestrial and marine environments with
   largely demonstrations about its high cytotoxicity toward all organisms
   [[30]1]. Plants can remediate Cd-metal polluted soils by diverse
   strategies generally termed phytoremediation, which are promising,
   clean, and easy to apply [[31]2]. To address phytoremediation, plants
   may act on the heavy metals through different schemes such as uptake
   (phytoextraction) or stabilization in the root system
   (phytostabilization) [[32]3]. Among all, phytostabilization is of
   particular interest, not only for its feasibility and effectiveness but
   also for its being peculiar to a wide range of plants which include
   commercial species and the model Arabidopsis thaliana [[33]4]. In A.
   thaliana, the obvious effect of Cd toxicity is a reduction of plant
   growth related to the inhibition of photosynthesis, respiration, and
   nitrogen metabolism, as well as to a reduction in water and nutrient
   uptake [[34]5]. The root system, as the first organ-sensing soil heavy
   metal, is strongly affected by Cd, which in A. thaliana inhibits
   primary root growth via affecting the root apical meristem (RAM) stem
   cell niche, while lateral root formation is somehow stimulated due to
   changes in root radial patterning [[35]1]. The most recent findings
   evidenced that many signaling molecules and plant growth regulators are
   involved in Cd sensing and downstream responses, which encompass
   modifications at both transcriptional and post-translation levels
   [[36]6,[37]7,[38]8,[39]9]. However, the relationship between all these
   factors remains unclear, differing with species, plant organs, heavy
   metal concentrations, and treatment durations (reviewed in [[40]8]).
   Recently, Arabidopsis halleri gained attention for its constitutivee
   tolerance to zinc (Zn) and Cd and for its close phylogenetic
   relationship to the model plant A. thaliana; thus, it is now emerging
   as a promising model to study heavy metal hyperaccumulation traits
   [[41]10,[42]11,[43]12]. Despite this evolutionary closeness between A.
   halleri and A. thaliana, the main mechanisms for tolerating metal
   stress are specific to one or the other species and therefore may be
   based on some different and particular biological processes, or the
   activation of them according to species-specific patterns.

   A. halleri, as a metal(loid) hyperaccumulator [[44]13], is attracting
   the research toward the use of this trait in the remediation of
   metal-polluted soils [[45]10,[46]14]. Different mechanisms underlying
   Cd and Zn tolerance and accumulation were found in the natural A.
   halleri population, highlighting contrasting Cd/Zn accumulation
   capacities and differences in the ability to adjust micronutrient
   homeostasis and flavonol content upon high Zn or Cd exposure
   [[47]15,[48]16]. Comparative transcriptomic studies of A. halleri and
   A. thaliana also identified ten genes candidates in metal homeostasis
   with putative roles in metal hyperaccumulation or metal hypertolerance,
   which encode various transmembrane transporters of metal cations and
   isoforms of a metal chelator biosynthetic enzyme
   [[49]17,[50]18,[51]19]. Specifically, hypertolerance in A. halleri
   relies on the activity of genes involved in metal uptake, xylem
   loading, root-to-shoot translocation, and metal sequestration in the
   leaf vacuole or metal accumulation in the veins
   [[52]15,[53]18,[54]20,[55]21,[56]22,[57]23]. Transcript levels of Heavy
   Metal ATPase 4 (HMA4), which encodes a plasma membrane P1B-type metal
   ATPase mediating cellular export for Zn and Cd xylem loading in the
   root, was substantially higher in A. halleri than in A. thaliana, as
   was also observed for several other candidate genes such as
   zinc-regulated transporter and iron-regulated transported (IRT)-like
   proteins (ZIPs). Additionally, the cell wall has been described in A.
   halleri to play a role in the adaptation to metalliferous soils
   [[58]15,[59]16] with several genes, being involved in its response to
   Cd stress, interacting among each other and with polysaccharides to
   structure the cell wall [[60]24,[61]25].

   However, to our knowledge, key connectors/biomarkers specifically
   related to A. thaliana phytostabilization or A. halleri
   hyperaccumulation abilities in response to Cd stress remain unclear.
   Thus, further insights are required mainly concerning the genetic
   network and functional relationships among genes/pathways involved in
   root growth and development under Cd stress. In particular, revealing
   of how proteins cooperate and participate in complexes or pathways to
   accomplish vital tasks is pivotal, especially for distinguishing
   functions and peculiarities which underlie a specific trait such as
   phytostabilization or hyperaccumulation. The in wet approaches are
   common practices to detect interactions among biological entities, but
   they are frequently insensitive to weak and transient interactions
   [[62]26].

   Nowadays, taking advantage of large-scale omics data and annotations in
   databases, bioinformatics-based high-throughput analyses are widely
   used to infer biological knowledge by applying holistic approaches from
   human health to agricultural sciences [[63]27]. In this sense,
   network-based strategies have great potential since networks can be
   combined with enrichment, clustering, data-guided prioritization, and
   other approaches for biomarker discovery and confirmation [[64]28].
   These holistic approaches resolve some problems associated with the
   traditional methods based on the comparison of group pairs, such as
   differentially expressed gene (DEG) analysis, and thus have more
   advantages than traditional methods, being able to comprehensively
   explore the associations between genes and target traits [[65]29].
   Among a variety of methods applied to infer networks, weighted gene
   co-expression network analysis (WGCNA) has shown its potential in
   developing an understanding of the relationships among genes and in
   focusing on specific modules of co-expressed genes and key genes
   correlated with specific phenotypes and, thus, have been used to select
   targets or to characterize traits [[66]12,[67]30,[68]31]. Additionally,
   expression data can be analyzed by traditional methods based on
   pairwise comparisons to identify the statistically significant
   variations of the expression patterns to profile traits, uncover
   targets, or focus on specific biomarkers in plants [[69]32,[70]33]. In
   this perspective, the identification of DEGs is a strong and relatively
   simple analysis whose application for validating key genes derived via
   WGCNA is a robust practice that has become very common [[71]34,[72]35].

   Thus, gaining knowledge from available RNA-seq data and based on
   co-expression analysis, the present study aimed to deepen the
   understanding of A. thaliana and A. halleri root responses to Cd stress
   by performing an in silico identification and comparison among key
   pathways and genes underlying the phytostabilization and
   hyperaccumulation traits. Specifically, RNA-seq data were used to
   derive highly correlated modules of co-expressed genes associated with
   Cd responses. By analyzing the interaction network between the
   co-expressed genes of the two Arabidopsis species (halleri and
   thaliana) and the Cd-related genes, the ultimate goal was to identify
   the hub genes and preliminarily explain the reasons for the different
   abilities of the species in terms of Cd responses. Moreover, the
   detection of DEGs under Cd stress was used to confirm the key genes as
   reliable markers, and data retrieved from databases, including
   literature, were used for additional confirmation. Key functions and
   genes can be considered as further in wet strategies or to be furtherly
   explored by in silico techniques to strengthen phytoremediation
   knowledge and related technologies.

2. Results

2.1. Identification of Functionally Active Genes of A. thaliana Root

   Various studies related to RNA-seq analysis were selected by a
   literature search to identify genes expressed in the root of A.
   thaliana. Other genes were retrieved by exploiting the functional
   annotation by GO terms as reported in the TAIR database. This led to a
   comprehensive list of 11,897 functionally active genes in A. thaliana
   root ([73]Figure 1), of which only 3502 (29.44% of the total) were
   present in two or more datasets.

Figure 1.

   [74]Figure 1
   [75]Open in a new tab

   The functionally active genes of A. thaliana root: the Venn diagram
   shows the overlap among the dataset of genes expressed in the root of
   A. thaliana (set1, set2, set3, set4, and set5) or functionally
   annotated by GO terms including “root” as reported in the highly
   curated database TAIR (TAIR). The table shows the references to the