Decrypting Molecular Mechanisms Involved in Counteracting Copper and Nickel Toxicity in Jack Pine (Pinus banksiana) Based on Transcriptomic Analysis

The remediation of copper and nickel-afflicted sites is challenged by the different physiological effects imposed by each metal on a given plant system. Pinus banksiana is resilient against copper and nickel, providing an opportunity to build a valuable resource to investigate the responding gene expression toward each metal. The objectives of this study were to (1) extend the analysis of the Pinus banksiana transcriptome exposed to nickel and copper, (2) assess the differential gene expression in nickel-resistant compared to copper-resistant genotypes, and (3) identify mechanisms specific to each metal. The Illumina platform was used to sequence RNA that was extracted from seedlings treated with each of the metals. There were 449 differentially expressed genes (DEGs) between copper-resistant genotypes (RGs) and nickel-resistant genotypes (RGs) at a high stringency cut-off, indicating a distinct pattern of gene expression toward each metal. For biological processes, 19.8% of DEGs were associated with the DNA metabolic process, followed by the response to stress (13.15%) and the response to chemicals (8.59%). For metabolic function, 27.9% of DEGs were associated with nuclease activity, followed by nucleotide binding (27.64%) and kinase activity (10.16%). Overall, 21.49% of DEGs were localized to the plasma membrane, followed by the cytosol (16.26%) and chloroplast (12.43%). Annotation of the top upregulated genes in copper RG compared to nickel RG identified genes and mechanisms that were specific to copper and not to nickel. NtPDR, AtHIPP10, and YSL1 were identified as genes associated with copper resistance. Various genes related to cell wall metabolism were identified, and they included genes encoding for HCT, CslE6, MPG, and polygalacturonase. Annotation of the top downregulated genes in copper RG compared to nickel RG revealed genes and mechanisms that were specific to nickel and not copper. Various regulatory and signaling-related genes associated with the stress response were identified. They included UGT, TIFY, ACC, dirigent protein, peroxidase, and glyoxyalase I. Additional research is needed to determine the specific functions of signaling and stress response mechanisms in nickel-resistant plants.


Introduction
Nickel and copper are micronutrients that become toxic to plants at high concentrations.The mining and processing of these heavy metals can contaminate soils and the atmosphere, which negatively impacts human health and the environment [1][2][3].The Greater Sudbury region has an extensive history of copper and nickel pollution and will continue to export these metals [4][5][6][7].Comparing the responses of plants to copper and nickel is a crucial step in streamlining the remediation strategy and ensuring efficiency and long-term sustainability.As redox heavy metals, copper and nickel both play important roles in plant physiology [8,9].Copper is a cofactor for electron carriers such as the cytochrome B 6 F complex, plastocyanin, and cytochrome oxidase [10][11][12].It is, therefore, Plants 2024, 13, 1042 2 of 40 integral to photosynthesis and cellular respiration.Nickel plays an important role in nitrogen cycling as a cofactor for urease, which is essential for the proper functioning of the ornithine-urea cycle [13,14].
Copper toxicity has many downstream effects and shares symptoms with nickel toxicity.An overabundance of either metal can decrease the concentration of iron, manganese, and zinc, which poses risks to processes that correspond to these metals [15][16][17][18][19][20].Copper or nickel toxicity can inhibit photosynthesis and cause chlorosis by different means.Excess copper reduces electron transfer in photosystem II by displacing the iron cofactor in plastoquinone QA, inhibiting plastoquinone QB, and preventing the transfer of electrons from tyrosine to P680+ [21,22].Excess nickel hinders electron transfer by altering the protein conformation of plastoquinone QB and replacing calcium in the oxygen-evolving complex [23,24].Both excess copper and nickel have been shown to displace magnesium in chlorophyll, leading to decreased chlorophyll content, reduced function, and a decreased capacity to transfer energy to the reaction center [25][26][27][28].Excess copper also facilitates the Haber-Weiss reaction and Fenton-like reactions, generating superoxide radicals, hydroxyl radicals, and hydrogen peroxide [29][30][31].Comparatively, excess nickel increases ROS generation by decreasing the activity of antioxidative enzymes such as SOD, CAT, APX, POD, and GSH-Px [32,33].An abundance of both metals can also cause severe water loss by altering stomata morphology and altering the rate of transpiration on the surface of leaves [34][35][36][37].Studies on various plants have shown that an overabundance of either metal can lead to necrosis and overall decreased plant growth [38].
Genes associated with copper or nickel resistance have been reported in various species.Membrane-bound transporters such as ZIP and COPT can export copper from the cytosol to the outside of the cell, reducing the initial uptake of copper into the root area [39][40][41][42].Upregulation of the IREG2 transporter can increase the transport of nickel into the vacuoles of root cells, contributing to vacuolar sequestration [43].NRAMP transporters may be involved in unloading excess copper or nickel into the xylem prior to root-toshoot translocation [44,45].Reports of elevated HMA transporter expression were found to enhance the translocation of excess copper to aerial tissue in some plants, increasing copper resistance [41,46,47].HMA9, in particular, may be involved in the loading of excess copper in both the xylem and phloem [48].HMA1 and HMA6/PAA1 can transport excess copper to the stroma of the chloroplasts, while HMA8/PAA2 could transport copper to the thylakoid lumen and plastocyanin [49][50][51][52][53][54].The delivery of copper to components in the chloroplasts suggests an increase in photosynthesis activity to counteract stress and ensure the homeostasis of metabolites [54,55].HMA7 may also counteract stress by transporting excess copper to ethylene receptors in the ER, which are involved in the modulation of growth and development [41,56].Excess copper may also elicit increased CCH expression to reduce the transport of copper to younger leaves, thereby safeguarding newly developed tissue [41,46,57].Chelators such as MT2a, MT2b, nicotianamine, and histidine have been reported to be present in xylem sap and were correlated with the improved translocation of copper and nickel [58][59][60][61].In particular, nicotianamine-metal complexes are directly translocated to the aerial component of the plant and may coordinate with YSL transporters to facilitate the loading of copper and nickel into the xylem [45,58,[62][63][64].Glutathione-Stransferases play several roles in copper and nickel resistance, which include increasing antioxidative activity, hormone signal transduction, and the production of the chelator glutathione [65,66].
The discovery of genes in metal-treated plants is crucial to understanding the mechanisms of metal resistance and planning an effective remediation strategy.Many pine species have been reported as having different responses and strategies toward various heavy metals [67][68][69].The successful utilization of Pinus banksiana in the regreening program in the Greater Sudbury Region prompts the investigation of the genetic response to copper toxicity compared to nickel toxicity [70].Transcriptome analysis will evaluate the genetic similarities or differences in gene expression between nickel-resistant and copper-resistant plants.The survival of this species in metal-contaminated soil also provides a unique op-Plants 2024, 13, 1042 3 of 40 portunity to discover genes that are exclusive to each metal, which is especially important when considering the limitations of the candidate gene identification process.The objectives of this study were to (1) extend the analysis of nickel-resistant and copper-resistant P. banksiana, (2) assess the differential gene expression in nickel-resistant genotypes (RGs) compared to copper-resistant genotypes (RGs), and (3) identify mechanisms specific to each metal.

Nickel and Copper Toxicity
Significant differences were observed in the damage ratings among plants for the groups treated with 1600 mg/kg Ni.In fact, 25% were considered resistant, 40% moderately susceptible, and 35% susceptible.Less damage was observed in the seedlings treated with potassium sulfate.In fact, no significant differences were observed in the seedlings treated with the salt used as a control compared to the water control.For copper, segregation in response to Cu was observed in a group treated with 1300 mg/kg of Cu.Overall, 20% were classified as resistant, 35% as moderately susceptible, and 45% as susceptible.Only resistant and susceptible genotypes, along with samples treated with potassium sulfate and water, were selected in triplicate for the transcriptome analysis.

Transcriptome Analysis
The transcriptome shotgun assembly project has been deposited in the NCBI BioProject database with the accession number PRJNA962116.Differentially expressed genes between copper and nickel-resistant genotypes are described in Table 1.The DEG between copper SG and nickel SG is presented in Table 2. Heatmaps of each pairwise comparison included every DEG between each group, providing an encompassing representation of dynamic gene expression patterns (Figure 1).Both heatmaps showed a high degree of uniformity in gene expression for individuals within a genotype (Figure 1).A heatmap of the DEGs between copper RG and nickel RG revealed a contrast in gene expression for each genotype (Figure 1).Copper RG had a slightly higher number of upregulated genes in comparison to nickel RG, indicating an upregulation of protein production and corresponding processes (Figure 1).A slightly higher number of downregulated genes in nickel RG may suggest the modulation or negative regulation of proteins and associated processes (Figure 1).For copper SG in comparison to nickel SG, the opposite pattern of gene expression was observed, as the majority of genes were upregulated in nickel SG and downregulated in copper SG (Figure 1).The volcano plots showed the degree of spread between the upregulated and downregulated genes in each pairwise comparison (Figure 2).For copper RG in comparison Plants 2024, 13, 1042 4 of 40 to nickel RG, the majority of DEGs had a fold change of less than 5 and a lower false discovery rate (FDR) in comparison to genes with a higher gene expression fold change or FDR (Figure 2).DEGs with the highest fold change had a lower FDR in comparison to other DEGs (Figure 2).For copper SG, in comparison to nickel SG, there was a relatively small number of DEGs that surpassed the FDR threshold (Figure 2).In addition, there was a larger number of downregulated genes with a relatively even distribution (Figure 2).The volcano plots showed the degree of spread between the upregulated and downregulated genes in each pairwise comparison (Figure 2).For copper RG in comparison to nickel RG, the majority of DEGs had a fold change of less than 5 and a lower false discovery rate (FDR) in comparison to genes with a higher gene expression fold change or FDR (Figure 2).DEGs with the highest fold change had a lower FDR in comparison to other DEGs (Figure 2).For copper SG, in comparison to nickel SG, there was a relatively small number of DEGs that surpassed the FDR threshold (Figure 2).In addition, there was a larger number of downregulated genes with a relatively even distribution (Figure 2).DEGs from the transcriptomes of copper RG compared to nickel RG were annotated and allocated to terms within the biological processes category using Omicsbox/Blast2GO (Figure 3).Table 3 and Table S2 illustrate top upregulated genes when copper RG was compared to nickel RG.Table 4 and Table S3 illustrate the top downregulated genes for copper RG compared to nickel RG.SG was compared to nickel SG, and Table 6 and Table S4 show the downregulated genes for copper SG vs. nickel SG.
DEGs from the transcriptomes of copper RG compared to nickel RG were annotated and allocated to terms within the biological processes category using Omicsbox/Blast2GO (Figure 3).Tables 3 and S2 illustrate top upregulated genes when copper RG was compared to nickel RG.Tables 4 and S3 illustrate the top downregulated genes for copper RG compared to nickel RG.Table 5 shows the top upregulated genes when copper SG was compared to nickel SG, and Tables 6 and S4 show the downregulated genes for copper SG vs. nickel SG.Overall, 449 DEGs were annotated and organized into the following subcategories of biological processes: DNA metabolic process (19.80%), response to stress (13.15%), response to chemicals (8.59%), signal transduction (7.68%), response to biotic stimulus (6.01%), catabolic process (5.85%), protein modification process (5.43%), transport (5.19%), response to endogenous stimulus (4.79%), carbohydrate metabolic process (3.10%), lipid metabolic process (2.76%), cell death (2.55%), and cell differentiation (2.19%) (Figure 3a) Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other" (12.88%).
(Figure 5a).The upregulated genes were annotated and organized into the following subcategories of molecular function: hydrolase activity (29.41%), nucleotide binding (23.53%),RNA binding (11.76%), kinase activity (11.76%), transporter activity (11.76%), signal receptor activity (5.88%), and signaling receptor binding (5.88%).processes for copper RG compared to nickel RG.The top upregulated and downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the biological process category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other."from copper RG compared to nickel RG were annotated and allocated to terms within the biological process category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other."The top downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the molecular function category using Omicsbox/Blast2GO (Figure 5b).The downregulated genes were annotated and organized according to the following subcategories of molecular function: transferase activity (30%), hydrolase activity (25%), nucleotide binding (20%), protein binding (10%), DNA binding (5%), DNA-binding transcription factor activity (5%), and transporter activity (5%).
The top upregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the cellular compartment category using Omicsbox/Blast2GO (Figure 6a).The upregulated genes were annotated and organized into the following cellular compartment locations: plasma membrane (21.49%),cytosol (16.26%), chloroplast (12.43%), extracellular region (11.31%),mitochondrion (9.48%), Golgi apparatus (5.90%), endoplasmic reticulum (5.47%), and nucleoplasm (3.37%).Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other" (14.28%).The top upregulated and downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the cellular compartment category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other." The top downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the cellular compartment category using Omicsbox/Blast2GO (Figure 6b).The downregulated genes were annotated and organized into the following cellular compartment locations: nucleus (40%), extracellular region (20%), plasma membrane (13.33%),mitochondrion (6.67%), cell wall (6.67%), chloroplast (6.67%), and vacuole (6.67%).

Physiological Mechanisms in P. banksiana Dealing with Soil Nickel and Copper Contamination
Plants growing on metal-contaminated soils are classified as resistant and have adapted to this stressed environment.The physiological copying mechanism of P. banksiana has been described by Moarefi and Nkongolo [71].In general, plants cope with metal contamination by using either avoidance and/or tolerance strategies.Metal-avoider plants use a strategy that prevents the entry of metal ions into their cells.Tolerant plants detoxify metal ions that have entered their cells by crossing the plasma membrane or the biomembranes of internal organelles.Three categories have been identified among tolerant plants.
Excluders maintain a low level of metal in their aerial tissues when growing on metalcontaminated soils.Metals enter root cells, but they prevent their movement from the root to aerial tissues.Indicators accumulate metals in their above-ground tissues, but the levels found in these areal tissues are reflective of the metal concentration in the soil (the metal concentration in their above-ground tissues is similar to the levels found in the soil).In accumulators/hyperaccumulators, metals that have entered the roots are translocated and accumulated in their areal tissues [72,73].Based on this classification, P. bnkasiana was classified as a metal indicator [71].

Nucleus (40%)
Extracellular region (20%) Plasma membrane (13.33%)Mitochondrion (6.67%) Cell wall (6.67%) Chloroplast (6.67%) Vacuole (6.67%) Figure 6.Percent distribution of the (a) top 100 upregulated genes in cellular compartment location for copper resistant (RG) in comparison to nickel resistant (RG) and the (b) top 100 downregulated genes in cellular compartment location for copper RG in comparison to nickel RG.The top upregulated and downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the cellular compartment category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other."

Physiological Mechanisms in P. banksiana Dealing with Soil Nickel and Copper Contamination
Plants growing on metal-contaminated soils are classified as resistant and have adapted to this stressed environment.The physiological copying mechanism of P. banksiana has been described by Moarefi and Nkongolo [71].In general, plants cope with metal contamination by using either avoidance and/or tolerance strategies.Metal-avoider plants use a strategy that prevents the entry of metal ions into their cells.Tolerant plants detoxify metal ions that have entered their cells by crossing the plasma membrane or the biomembranes of internal organelles.Three categories have been identified among tolerant plants.
Excluders maintain a low level of metal in their aerial tissues when growing on metalcontaminated soils.Metals enter root cells, but they prevent their movement from the root to aerial tissues.Indicators accumulate metals in their above-ground tissues, but the levels found in these areal tissues are reflective of the metal concentration in the soil (the metal concentration in their above-ground tissues is similar to the levels found in the soil).In accumulators/hyperaccumulators, metals that have entered the roots are translocated and accumulated in their areal tissues [72,73].Based on this classification, P. bnkasiana was classified as a metal indicator [71].

DEG Analysis Reveals Different Patterns of Gene Expression between Copper RG and Nickel RG
The number of DEGs between genotypes can provide information on the overall response of the copper-resistant genotype in comparison to the nickel-resistant genotype.The hundreds of DEGs between copper RG and nickel RG suggest that some mechanisms are differentially regulated in response to metal stress.Although the majority of the transcriptome was similar, a small number of genes can often orchestrate significant processes and changes pertaining to heavy metal resistance [74,75].The presence of 449 DEGs at high stringency most likely contributed to differences between the response to copper and the response to nickel.The different number of upregulated and downregulated genes between genotypes infers the use of different physiological mechanisms to achieve a similar endpoint phenotype.The presence of more upregulated genes in copper RG suggests more direct mechanisms that confront copper toxicity with the increased production of proteins and associated processes.In contrast, the presence of more downregulated genes in nickel RG indicates the negative regulation of certain processes or signaling pathways.For example, some phytohormones and pathways may require a controlled response with limitations on protein production as a means to optimize the heavy metal stress response and prevent cellular toxicity [76].For copper SG compared to nickel SG, there were only 41 DEGs at high stringency, indicating a high similarity in gene expression between susceptible plants.This finding suggests that susceptible plants have less targeted specificity when responding to copper or nickel stress.Susceptible plants often lack the expressed genes needed to ameliorate heavy metal toxicity [76][77][78].The similarity between the genotypes suggests that gene expression was not directly involved with the detoxification of heavy metals but was instead managing downstream stress mechanisms and tissue damage.After a certain point, the cascading effect of heavy metals can manifest itself as a blanket stress effect, causing nonspecific tissue damage and cytotoxicity [37,79].The higher number of upregulated genes in nickel SG may indicate attempts to increase processes that mitigate nickel stress, counteract cascading downstream effects, or restrict the accumulation of toxic secondary metabolites induced by excess nickel.In contrast, the higher number of downregulated genes in copper SG may suggest the negative regulation of proteins and metabolites that would be otherwise toxic at higher concentrations [66].Susceptible plants may utilize shared stress mechanisms associated with reducing tissue damage, managing the turnover of metabolites, and responding to necrosis [80][81][82][83].No biomarkers specific to copper or nickel susceptibility were identified in susceptible genotypes, although general stress response mechanisms were identified.
The DNA metabolic process had the highest proportion of gene expression, which suggests different levels of DNA regulation between copper and nickel.The DNA metabolic process encompasses processes such as DNA methylation, modification, repair, processing, maintenance of repetitive sequences, and telomere maintenance [84].Excess copper causes double-stranded breaks, DNA oxidation, and an overall reduction in DNA content [85][86][87].Excess nickel has also been shown to elicit DNA damage and reduce DNA content, although the specific mechanisms are not well documented in plants [88,89].The high proportion of genes associated with nuclease activity indicates different levels of DNA repair mechanisms, which may be used to mitigate the oxidative stress on DNA induced by heavy metals [90].It implies a possible function related to DNA methylation and the epigenetic modulation of gene expression.Multiple plants differ in global methylation patterns elicited by excess copper and nickel [88,[91][92][93].Different methylation patterns may also persist in later developmental stages to facilitate long-term adaptation to metal-contaminated sites [91].DNA methylation can orchestrate significant changes to plant physiology, especially at earlier stages of plant development [94].Genes were ranked based on Log2 FC from each pairwise comparison.
Considering that the samples were replicated three times and only the DEG detected at high stringency cut-offs, based on the false discovery rate (FDR) analysis, were selected, it is highly unlikely that the reported gene expression variations were due to chance.
The response to stress and other related terms had a high proportion of gene expression that indicated physiological differences in stress mitigation toward copper and nickel.Both excess copper and nickel can alter the homeostasis of other heavy metal ions, decrease photosystem II functionality, reduce chlorophyll function, and alter phytohormone distribution during early development [15,16,18,21,23,[95][96][97][98].The response to stress term is also highly expressed in the top upregulated and downregulated genes for copper RG and nickel RG, with nickel RG having a comparatively higher allocation of gene expression.Differences in the physiological impact of copper and nickel are likely responsible for this significant difference.

Analysis of the Top Upregulated Genes in Copper-Resistant Plants Compared to Nickel-Resistant Genotypes Reveals Mechanisms That Are Associated with Copper Resistance in Pinus banksiana
Describing the top upregulated genes of copper-resistant plants compared to nickelresistant plants allowed for the identification of mechanisms associated with the response to copper.Analysis of the top upregulated genes also corroborated previous findings on copper resistance, establishing mechanisms that are specific to copper and are potentially involved in copper detoxification and tolerance [99].Additionally, another transporter associated with copper resistance was identified, demonstrating the efficacy of transcriptome analyses on multiple pairwise comparisons for describing gene expression.
Two heavy metal transporters were identified as exclusive to copper in comparison to nickel.One of the most expressed genes with several isoforms encodes the pleiotropic drug resistance protein 1 (NtPDR1), which is an ATP-binding cassette (ABC) transporter belonging to the ABC G subfamily [100,101].ABC transporters are membrane-bound and use ATP hydrolysis to facilitate the active transport of entities [102,103].The conserved protein regions of ABC transporters include two nucleotide-binding folds associated with ATP hydrolysis and two hydrophobic transmembrane domains involved in the determination of substrate specificity [102,104].Upregulation of AtPDR family transporters in Arabidopsis thaliana increased cadmium resistance and lead resistance by exporting these metals from the cytosol to the apoplast [105,106].It has been proposed that PDR may serve a similar function in Pinus banksiana by exporting excess copper from the cytosol and reducing the overall concentration of copper in the cell.PDR may also be involved in the general stress response, as demonstrated by the upregulation of Ospdr9 in response to hormones such as jasmonic acid and ABA [107,108].
One of the top upregulated genes encoded for the Yellow Stripe 1 (YS1) transporter is a transmembrane transporter that facilitates the movement of heavy metals between intercellular and intracellular compartments (Table 3a) [64,109,110].YSLs have a conserved oligopeptide transporter domain that binds to metal complexes and oligopeptides [111].YS1 has been shown to import various heavy metals complexed with nicotianamine (NA) from the rhizosphere into the root cells [109,112].YSL3 can initiate the loading of metal-NA complexes from the apoplast into the xylem and mediate the root-to-shoot translocation of iron and cadmium [64,113,114].In Brachypodium distachyon, BdYSL3 can load copper-NA into the phloem and redistribute copper-NA complexes to the flowers and reproductive organs [110].Similarly, OsYSL16 in Orzya sativa redistributes copper-NA complexes from the phloem to the seeds and younger leaves [115].The dynamic response of YSL family proteins to alterations in heavy metal content may differ between species, highlighting the different strategies that may be utilized in response to excess heavy metals.In Oryza sativa, YSL6 conferred manganese resistance by increasing the import of manganese into the symplast, which reduced manganese accumulation and subsequent damage in the apoplast [116].The upregulation of YSL7 in transgenic tobacco conferred resistance to excess iron, cadmium, and nickel by increasing the translocation of these metals to the shoots and seeds, which may have allowed for effective sequestration [117].The role of YSL in mediating xylem loading and translocation has yet to be described in coppertreated plants.In response to excess copper, AtYSL2 transports copper-NA complexes laterally within the xylem and phloem [118].Upregulation of YSL3.1 in transgenic tobacco and Oryza sativa led to copper accumulation in younger leaves and conferred copper resistance [119].Nicotianamine is the most commonly reported chelator for copper due to having a higher binding affinity in comparison to other metals, although it is also possible for copper to form a complex with phytosiderophore or be transported as an ion in rare circumstances [109,111].It is, therefore, possible for YSL upregulation to be associated with lateral transport in the vasculature or the redistribution of copper from the phloem to younger leaves in Pinus banksiana [110,118].Redistribution of copper to the younger needles may serve as a sync for copper sequestration, especially if older needles are weaker or more prone to potential toxicity.It is unlikely for YSL to redistribute excess copper to the seeds, as this area is commonly targeted by heavy metals to hinder development and germination [120,121].
Another highly expressed gene was the heavy metal-associated isoprenylated plant (HIPP), which is a metallochaperone that captures and transports heavy metal ions to proteins located in other cellular compartments (Table 3a) [122].HIPPs comprise a heavy metal-associated (HMA) domain associated with directly binding to free metal ions and a cysteine-rich isoprenylation site that interacts with other proteins, enzymes, and membranes [122,123].Variation in the isoprenylation site may modify cellular compartment localization, the site of protein targeting, and overall signal transduction [122,[124][125][126].As a divalent ion, copper can bind to other metallochaperones, such as ATFP3 and ATFP6 in Arabidopsis thaliana, supporting the potential role of HIPP in copper chelation and homeostasis [122,126,127].HIPP may be involved in copper resistance, as it has been shown to confer heavy metal resistance in yeast and Arabidopsis thaliana [128,129].HIPP conferred resistance to excess copper and cadmium in Camellia sinensis exposed to excess copper or cadmium [130].The in vitro binding of HIPP to copper, cadmium, and zinc also corroborates this proposed function [126,127].Knockout of OsHIPP17 in Oryza sativa reduced the expression of copper transporters and cytokinin-related genes, demonstrating crosstalk with other transporters and signaling proteins [131].In Pinus banksiana, HIPP may potentially bind to and trap copper, reducing the effects of copper toxicity on plant physiology [132].
Among the top DEGs were two genes associated with photosystem II, encompassing the GO terms response to light stimulus and response to stress.One of the top upregulated genes encodes the blue copper protein (BCP), which catalyzes the reduction and oxidation of substrates, subsequently facilitating the transport of electrons from an electron donor to an electron acceptor [133,134].Excess copper can reduce chloroplast function and inhibit photosynthesis by replacing ions in enzymes involved in the electron transport chain [21,23,95].Increased plastocyanin expression can recover the frequency of electron transfer reactions and improve photosynthesis, thereby counteracting copper-induced photodamage and maintaining a usable pool of carbohydrates for other metabolic needs.In Glycine max, excess copper upregulated plastocyanin expression as part of a concerted effort to increase the rate of photosynthesis [135].Alternatively, the BCP gene may also encode proteins in the stellacyanin subfamily, which could be involved in facilitating redox reactions with substrates in the cell wall [136,137].In addition to having a low redox potential, stellacyanin has a redox-active copper-binding domain, which can aid in lignin synthesis or mediate the stress response [136][137][138].The upregulated expression of BCP in some plants resulted in higher lignin production, which could act as a counteractive response to heavy metal stress by strengthening the cell wall and enhancing its defense properties [139,140].The lignin production function of the BCP gene is supported by the expression of other top regulated genes in Pinus banksiana, which are involved in enhancing the structural integrity of the cell wall.Further research is needed to reveal the molecular structure, function, and subfamily classification of the blue-collar protein.
Another top upregulated gene that was associated with photosynthesis function is a gene encoding the 22 kDa Photosystem II protein, which facilitates the non-photochemical quenching of energy released by pigment-mediated fluorescence (Table 3a) [141,142].Light absorption causes chlorophyll a, chlorophyll b, and xanthophyll to fluoresce, releasing energy that can be utilized in photosynthesis reactions and O 2 radical production [143][144][145][146]. Non-photochemical quenching captures and dissipates the released energy as heat, preventing the production of O 2 radicals, which damage components of the photosynthetic process [141][142][143].Excess copper has been shown to decrease chlorophyll concentration, alter chlorophyll morphology, and reduce enzyme activity associated with photosystem II [95].The copper-mediated generation of ROS can further exasperate these symptoms by targeting thylakoid membranes and chloroplast membranes [29][30][31].The upregulation of the 22 kDa photosystem II protein could increase non-photochemical quenching, thereby reducing ROS-mediated damage and subsequent photodamage [146,147].The 22 kDa Photosystem II protein and the blue copper protein may be upregulated to prevent photodamage caused by copper and protect photosynthesis machinery.Both proteins have a similar function to the light-inducible protein, which was identified in a transcriptome analysis between copper-resistant and susceptible plants [99].
The high proportion of genes associated with the biosynthetic process and the carbohydrate metabolic process play a crucial role in cell wall metabolism (Table 3a).Many of the top upregulated genes have annotated functions that correspond to the protection and reinforcement of the cell wall.One of the top upregulated genes encoded for Shikimate O-hydroxycinnamoyltransferase (HCT), which is involved in lignin synthesis [148,149].HCT catalyzes the conversion of p-coumaroyl-CoA, shikimate, or quinic acid to 4-coumaroylshikimate [148][149][150].HCT also catalyzes the conversion of 5-ocaffeoylshikimic to caffeoyl-CoA and the subsequent conversion of caffeoyl-CoA to shikimate [151].The initial set of reactions in the phenylpropanoid pathway directs substrates and intermediates through a reaction cycle that eventually yields the lignin monomers S-lignin, H-lignin, and G-lignin [152].HCT can thus act as a crucial control point for lignin synthesis [153,154].In addition to altering the proportions of monolignin subunits, silencing HCT can lead to an overall reduction in lignin content, an increase in cellulose, a decrease in cell wall thickness, and an increase in lumen diameter [154][155][156][157][158]. In Zoysia japonica, upregulation of HCT increased the ratio of the S-lignin monomer and elevated the activity of lignin-associated enzymes via phytohormone signaling [151].In transgenic tobacco, overexpression of HCT increased lignification and the number of secondary cell layers, which provided increased stem strength and rigidity [159].Lignin is an important contributor to the formation of secondary cell walls in fibers and vasculature [160].The hydrophobic domains in lignin provide waterproofing properties that are essential for separating anatomical structures [161].Lignin also crosslinks with xylan and cellulose, providing additional mechanical strength, rigidity, and thickness to the cell wall [162,163].Some researchers have proposed that lignin may be able to block and reduce the transport of heavy metals to other areas of the plant [164].Polysaccharides, or functional groups, in the cell wall can bind to copper and act as a sync for copper sequestration, which may prevent the intracellular transport of copper to other organelles, such as the chloroplast and vacuole [165].The identification of top upregulated genes involved in protecting the photosynthesis machinery and the chloroplast corroborates this proposed function.
A gene encoding the cellulose synthase-like E6 (CslE6) protein was also among the top upregulated genes related to cell wall metabolism (Table 3a).Csl family proteins catalyze the addition of UDP-glucose to the beta-1,4 glycan polymer within cellulose [166].Csl family proteins have also been shown to catalyze the synthesis of hemicelluloses such as mannans, glucomannans, xylans, xyloglucans, and homogalacturonan [167][168][169][170][171]. Cellulose comprises the most abundant portion of the primary and secondary cell walls, whereas hemicellulose comprises a considerable portion of the secondary cell wall [172].Cellulose functions as a load-bearing microfibril that imparts tensile strength and turgor pressure to the cell wall, thereby maintaining the structure and morphology of the cell wall [173][174][175].Hemicellulose crosslinks with cellulose and lignin, providing structural support and reinforcement for the cell wall [176].In particular, the presence of hemicellulose lignin matrices provides resilience against shearing forces [177].In cadmium-tolerant Arabidopsis thaliana, cadmium has been found to bind to cellulose and hemicellulose [178].In some plants, Csl can maintain cell division for normal development and provide thickness and density to the cell wall [171,179].Overexpression of CslD3 in defective Arabidopsis has been shown to increase cell elongation and rescue the integrity of the cell wall [180].The upregulation of Csl genes can lead to higher cellulose production, mechanical strength, and overall mass in the primary and secondary cell walls [181].
Another top upregulated gene that was associated with cell wall metabolism encoded mannose-1-phosphate guanylyltransferase (MPG).MPG catalyzes the production Plants 2024, 13, 1042 20 of 40 of GDP-mannose, which is involved in the N-glycosylation of proteins embedded in the cell wall [182].N-glycosylation is a post-translational modification that attaches a sugar residue to the free amine group of asparagine within the sequence asparagine-X-Serine/Threonine [183].N-glycosylation can alter various properties such as polarity, solubility, reactivity, and cell compartment localization, which can alter protein conformation and overall function [183][184][185][186]. GDP mannose is also a precursor for GDP-fucose, which is a component of the cell wall that is involved in adhesion, polysaccharide crosslinking, and lignification [187][188][189][190]. Additionally, GDP mannose is a precursor for ascorbic acid, which is a multifunctional compound involved in stress signaling, defense processes, and crosstalk with phytohormones [188,191,192].In Arabidopsis thaliana, knockout of the MPG gene leads to callose deposition in the primary cell wall, incomplete construction of the cell wall, and a considerable decrease in cellulose content [182].Upregulation of MPG in Pinus banksiana could potentially mitigate copper toxicity by producing GDP-mannose, which contributes to the reinforcement of the cell wall.
A top upregulated gene associated with cell wall metabolism is encoded for polygalacturonase (pectinase).Polylgalacturonase catalyzes the hydrolysis of alpha-1,4 glycosidic bonds in polygalacturonan (pectin) into monosaccharide subunits [193].Pectin comprises a large portion of the primary cell wall in plants, although some studies suggest a minimal presence in wood tissue [194,195].Pectin plays various roles within the cell wall, which include mediating cell-to-cell adhesion, facilitating ion transport, and contributing to porosity [196][197][198].The variable sidechains of pectin can link to cellulose and hemicellulose, contributing to the overall mechanical strength of the cell wall [199].The modification of pectin and the reorganization of cell wall components is a mechanism that can be utilized against heavy metal-induced stress [200,201].Changes in the proportion of cell wall constituents may impart certain strategic advantages, depending on the structure of the cell wall and metabolic requirements [202].Downregulation of pectinase in various plants also improved tolerance to multiple stressors by decreasing cell expansion and cell separation and increasing cell density [193,203,204].Moy et al. found that polygalacturonate was among the most downregulated genes in copper-resistant trees in comparison to coppersusceptible trees, indicating an association with copper tolerance.In a previous study on Pinus banksiana, other cell wall synthesis-related genes were found to be among the top upregulated genes in copper RG compared to copper SG, indicating an association with copper resistance [99].
Information from the top upregulated genes and a previous study comparing copper RG with copper SG can be used to propose a model of copper resistance [99].In comparison to nickel-treated plants, a lower proportion of genes were associated with the response to stress and signal transduction, indicating that Pinus banksiana may prioritize the direct detoxification of copper ions as opposed to counteracting symptoms caused by copper ions.NtPDR could export copper from the cytosol to the apoplast, at which point copper ions may interact with the cell wall [105,106].HCT, CslE6, and MPG are involved with the reinforcement of the cell wall.In addition, the downregulation of beta-glucosidase and polygalacturonase in copper RG further supports this function [99].Copper may then be sequestered in the cell wall and bind to ligands to decrease mobility and accumulation.Alternatively, cell wall synthesis could be driven by a higher amount of ROS mitigation in copper-resistant plants, as ROS may induce oxidative damage.Copper also seems to target photosynthesis machinery and inhibit photosynthesis.In copper-resistant plants, the early light-induced protein (ELIP) was upregulated, supporting the function of photosynthesis protection for the blue copper protein and the 22 KDA protein.Metabolic and protein-based assays should be performed to validate this proposed model.

Analysis of the Top Downregulated Genes for Copper-Resistant Plants Compared to Nickel-Resistant Genotypes Reveal Mechanisms That Are Associated with the Response to Excess Nickel
In nickel-resistant trees, the highest proportion of the top regulated gene expression was associated with the response to stress and signal transduction, indicating the imple-mentation of distinct strategies in response to elevated levels of nickel.Stress response mechanisms may be used to alleviate symptoms of nickel toxicity, such as tissue damage, decreased photosynthesis, reduced plant biomass, and water loss [18,23,97].Mechanisms that were associated with photosystem II and mitigating photodamage were more exclusive to copper than nickel.The differences in stress mechanisms could be explained by nickel and copper targeting different components of photosystem II.Excess copper inhibits plastoquinone QA and QB function, whereas excess nickel inhibits plastoquinone QB and the functionality of the oxygen-evolving complex [21,23,24].Differences in ROS production may explain a higher proportion of genes allocated to cell wall metabolism in copper RG relative to nickel RG.Excess copper can increase the production of hydrogen peroxide and hydroxyl radicals from the Haber-Weiss reaction and Fenton-like reactions [29][30][31]205].Comparatively, excess nickel may indirectly cause ROS stress by regulating the activity of antioxidative enzymes such as SOD, catalase, and peroxidase [32,33].Due to nickel being less involved in plant physiology and being a cofactor for crucial enzymes, it is possible that excess nickel negatively impacts fewer areas of plant physiology and is thus relegated to fewer categories of biological responses.There is also a lack of identified transporters specific to nickel, which severely limits the physiological strategies that Pinus banksiana can utilize to detoxify nickel ions.
Among the top downregulated genes was a gene that encoded a UDP-glycosyltransferase (UGT), which catalyzes the transfer of glucose to abscisic acid [76].The conversion of ABA to ABA-GE is a negative feedback mechanism that negatively regulates ABA synthesis and signaling [206,207].As a phytohormone, ABA performs crosstalk with other stressrelated hormones, modulating various aspects of stress mitigation such as transpiration rate, cell death, stomatal closure, photosynthesis, germination, and growth [208][209][210][211][212][213][214][215].However, multiple instances of ABA application were found to decrease growth and cell replication in lieu of mediating physiological changes [215][216][217].Application of ABA to cadmium-treated plants reduced root-to-shoot translocation and decreased IRT1, which is a nonspecific transporter that also mediates nickel uptake [218,219].Reduction in ABA can thus be seen as a means to recover growth from nickel-mediated stress or indirectly influence nickel transport.In aluminum-stressed Glycine max, high UGT expression leads to the alteration of cell wall components, reducing callose, hemicellulose, glucose, and xylose [220].In transgenic Arabidopsis thaliana, UGT interacted with metallothionein, facilitated hormone crosstalk, and alleviated oxidative stress in response to heat stress [221].Overexpression of UGT in various species subjected to different abiotic stressors leads to decreases in auxin signaling, increases in seed germination, higher growth rates, and elevated flavonoid content [222][223][224][225].The impact of UGT upregulation in nickel-treated Pinus banksiana may depend on various physiological parameters, necessitating future measurement and correlation.
A gene encoding TIFY was one of the most downregulated genes in copper RG in comparison to nickel RG.TIFY are transcription factors that constitutively suppress jasmonatemediated signaling and are modulated by the ubiquitin-proteasome pathway [226,227].Among the other jasmonate-related proteins identified in nickel RG in comparison to water, TIFY was the only protein to be highly upregulated in the nickel-resistant genotype compared to the copper-resistant genotype.TIFY responds to a variety of stressors and can be induced by the jasmonate and ABA signaling pathways [228,229].In Glycine max, upregulated TIFY led to increases in peroxidase and catalase activity to counteract saltinduced stress [230].In several plants subjected to various abiotic stressors, upregulated TIFY led to an increase in proline content and a reduction in MDA content, which is a toxic byproduct of lipid peroxidation [231][232][233][234]. Upregulation of TIFY in nickel could be a counteractive response to jasmonate production, as jasmonate is a phytohormone that reduces cell replication, photosynthesis, and overall growth in exchange for an elevated stress response [235][236][237].Jasmonate inhibitors have been found to increase heavy metal tolerance by decreasing the impact of heavy metals on photosystem II and preventing chlorophyll loss [234,238,239].Depending on the species, increased methyl jasmonate or jasmonic acid production may forego growth and development in exchange for increased heavy metal tolerance [240][241][242][243].
Among the top downregulated genes was a gene encoding for 1-aminocyclopropane-1carboxylate synthase 7 (ACC synthase 7), which catalyzes the synthesis of 1-aminocyclopropane-1-carboxylate (ACC) [244].ACC undergoes a rate-limiting oxidation reaction to produce ethylene, which is a phytohormone involved in stress processes, growth processes, and overall plant development [245,246].Ethylene is a multifaceted phytohormone that may induce large alterations in stress tolerance when overexpressed or underexpressed.In Arabidopsis thaliana, ethylene can transcriptionally regulate the reduction in root growth and is also involved in photosynthesis recovery [247][248][249].Ethylene mediates crosstalk with several hormones and signaling molecules, which include jasmonates, ABA, and nitrogen oxide [249][250][251].In some instances, ethylene can indirectly alter antioxidative enzyme activity or influence glutathione-S-transferase synthesis [249,[252][253][254].In Nelumbo nucifera, the upregulation of ACC synthase increased both ACC and ethylene content, which led to higher sensitivity to cadmium-induced stress [254].The increase in ethylene also caused the inhibition of catalase, ascorbate peroxidase, and glutathione reductase.In cadmium-treated Nicotiana, ethylene was decreased in the roots and increased in the shoots, resulting in higher shoot growth and photosynthesis recovery [255].The endpoint phenotype for ethylene regulation seems to depend on the specific phytohormone crosstalk and the growth or photosynthetic requirements of the plant responding to stress [249,252].Ethylene may also induce different outcomes based on the type of stressor involved [256].Being the precursor to ethylene, the upregulation of ACC in nickel-resistant plants may indicate prioritization of photosynthesis recovery, ROS detoxification, and the stress response in exchange for tissue growth and development.
A top downregulated gene encodes for the dirigent protein, which is involved in the stress response, defense, and reorganization of the cell wall [257].Dirigent proteins are involved in the production of lignins and lignans, which alter cell wall composition to counteract nickel-induced stress [258].Dirigent proteins specifically mediate the phenoxy radical coupling of coniferyl alcohol to stereoselective pinoresinol [259].Other genes involved in the synthesis of other cell wall components are not specific to nickel-resistant plants in comparison to copper-resistant plants.The upregulation of dirigent proteins and lignin may increase the mechanical strength and rigidity of the secondary cell wall [159,162,163].The upregulation of lignin has also been correlated with increased cadmium deposition within the cell wall [260].In Arabidopsis thaliana, DIR is involved in Casparian strip formation, which reinforces the impermeability of the cell wall [261,262].Casparian strip formation may increase the blocking or reduction in heavy metal transport.In addition to being a reserve precursor for lignin formation, lignan is also a secondary metabolite that may be involved in the regulation of oxidative stress and flavonoid production [263][264][265].Reduced lignan has been shown to strongly mediate ROS scavenging in the xylem, serving a strong protective function for the vasculature [264,266].
Another top downregulated gene is encoded for peroxidase, which is an enzyme involved in ROS scavenging [267].Although not specific to any given stressor, peroxidase expression is a strong biomarker for oxidative stress, most notably in response to peroxide [268].Peroxidase facilitates the electron transfer from hydrogen peroxide to a given substrate [267].Peroxidase can quench radical peroxide and oxygen species, which leads to decreased levels of ROS [269].In some trees, the upregulation of peroxidase is associated with increased lignin production, as peroxidase can oxidize monolignol radicals while simultaneously reducing peroxide [160].The dual functionality of targeting oxidative stress and cell wall reinforcement could significantly contribute to heavy metal tolerance [270][271][272][273]. Peroxidase and peroxide levels may also negatively regulate auxin levels, which can impact plant growth and development [274].Decreased auxin levels can impede apical meristem cell division, apical meristem differentiation, shoot branching, root branching, and overall plant growth [274][275][276][277].
A gene encoding for glyoxalase I was among the top downregulated genes.Glyoxalase I catalyzes the conversion of methylglyoxal (MG) to the inert compound S-Dlactoylgluathion, which is subsequently converted to D-lactate by glyoxalase II [278,279].The conversion of MG to an inert compound serves dual roles by diminishing MG toxicity in plant tissue and generating D-lactate, which is a substrate involved in glucose metabolism [280].As a byproduct of glycolysis, MG is involved in signaling pathways that modulate various areas of plant physiology.At low concentrations, MG can regulate ROS content, control stomatal closing, and activate potassium channels [281,282].Excess nickel induces the overproduction of MG [283].At higher concentrations, MG can reduce photosynthesis activity, impede root growth, decrease seed germination, cause chlorosis, and prevent ethylene synthesis [284][285][286][287][288][289].Excess MG can also be detrimental to the function of critical enzymes and proteins, as MG facilitates the glycation of arginine and lysine [290].The high expression of glyoxalase I in nickel RG could suggest a counteractive measure to reduce MG and the corresponding side effects causing cytotoxicity.

Analysis of the Top Upregulated Genes for Copper-Susceptible Plants Compared to Nickel-Susceptible Plants
Comparing the top regulated genes between copper-susceptible plants and nickelsusceptible plants may show tolerance mechanisms that are specific to either metal or reveal biomarkers indicative of copper or nickel vulnerability.Alternatively, DEGs in the susceptible genotype may not be associated with the response to heavy metals but could instead be associated with the management of cell death and tissue necrosis.
Among the top upregulated genes was a gene encoding inositol oxygenase, which catalyzes the oxidation of myto-inositol to D-glucuronic acid [291].D-glucuronic acid is a critical component for hemicelluloses such as xylan and xyloses and can increase the level of binding and stability between components of the cell wall [292][293][294].Upregulating the synthesis of UDP glucuronic acid can establish a pool of precursors to be utilized for building cell wall carbohydrates [295].In Arabidopsis thaliana, mutants that did not produce D-glucuronic acid had altered hemicellulose, which resulted in a lower binding capacity and the accumulation of aluminum in the cell wall [296].Many studies have shown the upregulation of inositol oxygenase to various stresses, with some studies demonstrating an improved stress response correlated to increased ROS mitigation activity [297,298].
A top upregulated gene encodes for a monosaccharide transporter, which facilitates the transport of glucose, xylose, and 3-O-methylglucose across different compartments [299,300].In Oryza sativa, the expression of monosaccharide transporters in leaf cells, xylem cells, and sclerenchyma cells suggests a function that may pertain to cell wall reinforcement [299].In response to stress, the upregulation of monosaccharide transporters may imply the accumulation of constituents in a pool to be utilized for either cell wall development or fulfilling metabolic requirements [299][300][301].In many plants, increased expression of monosaccharide transporters led to increased production of protein, lipids, and overall development [302].Fulfilling metabolic requirements could induce further growth and reproduction to counteract stress-induced damage [301,302].
A top downregulated gene encodes for a calmodulin-related protein, which binds calcium and is a critical component of the calmodulin signaling pathway [303].As a calcium-binding sensor, the activity of calmodulin is dependent on cytosolic calcium levels, which are subject to alteration by a variety of different stressors [304].In response to various stressors, calmodulin can facilitate crosstalk with signaling proteins, interact with secondary metabolites, interact with ion channels, and modulate enzymes to alter various aspects of physiology [305][306][307][308][309].Many ion channels, phosphatases, and kinases are calmodulin-dependent, although the corresponding functions are too broad to attribute a specific function toward heavy metals [310][311][312][313][314]. Similar to other general stress responses, increased calmodulin is correlated with increased ROS scavenging [308].In transgenic Nicotiana tabacum, overexpression of calmodulin conferred nickel tolerance, possibly by providing binding sites for nickel, which obstructs the channel or alters functionality [315].It should be pointed out that the validation of gene expression was performed by analyzing each treatment in triplicate and by applying a stringent pairwise comparison test.Additional validation by RT-qPCR can be considered for further study.Considering the scope of this study, such analysis will be time-consuming and not cost-effective.

Plant Growth and Treatment
The protocols for Pinus banksiana seedling growth, incubation, storage conditions, and metal treatment were similar to previous studies [316].The College Boreal Plant Center in Sudbury, Ontario, provided six-month-old Pinus banksiana seedlings.These seedlings were transplanted into square plantar pots containing a 1:1 composition of sand with a mixture of 79% Sphagnum moss, 17% perlite, and 5% composted peat moss.The seedlings were incubated in a growth chamber for a month, receiving consistent water and fertilizer to ensure adequate growth.Plants were watered as needed and fertilized twice a week with equal amounts of nitrogen, phosphorus, and potassium (20-20-20).All treatments were applied in a randomized block design to account for variability among individual seedlings.Fifteen seedlings were given 1300 mg of copper sulfate per 1 kg of soil, administered through an aqueous solution in the soil.Another fifteen seedlings were treated in the same manner with 1600 mg of nickel sulfate per 1 kg of soil.The concentration of each treatment was representative of metal-contaminated sites in the Greater Sudbury Region [6].Fifteen seedlings were given deionized water, which represented the negative control.Fifteen seedlings were treated with 1300 mg/kg of potassium sulfate as a salt control for the copper sulfate treatment.Fifteen seedlings were given 1600 mg/kg of potassium sulfate to represent the salt control for the nickel sulfate treatment.
After two weeks of incubation with continued monitoring, damages were recorded to identify the most resistant and susceptible seedlings using the 1 to 9 scale, where 1 to 3 represent resistant plants, 4 to 6 represent moderately resistant/susceptible plants, and 7 to 9 represent susceptible plants.Details of the rating scale are described in Table S1.Needles from resistant genotypes (RGs) with no symptoms and green needles from susceptible genotypes (SG) were harvested and wrapped in individual aluminum foils for RNA extraction.For longer-term storage, the needles were flash-frozen using liquid nitrogen and stored in a freezer at −80 • C.

Transcriptome Analysis of Pinus banksiana
Total RNA extraction, sequencing and transcriptome assembly, annotation of Pinus banksiana using BLAT matching, quantification of gene expression, and quality control (QC) analysis were described in Moy et al., 2023 [316].Likewise, differential gene expression (DGE) analysis of pairwise comparisons and the identification of top upregulated and downregulated genes followed the procedure detailed in Moy et al., 2023 [316].Three biological samples of each group, which include copper-resistant genotypes (RGs), Ni RGs, copper susceptible (SG), and Ni SG, along with water-treated samples and potassium sulfate-treated plants, were used for transcriptome analysis.

RNA Extraction
RNA extraction was conducted in accordance with the protocol outlined in the NORGEN BIOTEK Plant/Fungi Total RNA Purification Kit (https://norgenbiotek.com/product/plantfungi-total-rna-purification-kit, accessed on 1 September 2022).RNA quality was examined using agarose gel electrophoresis, and RNA quantification was performed using the Qubit™ RNA BR assay kit (Qubit RNA BR Assay Kit by Invitrogen Life technologies Corp., Eugene, OR, USA).The extracted RNA was stored at −80 • C in a freezer until use.

RNA Sequencing and De Novo Transcriptome Assembly
For every sample, messenger RNA was purified from the total RNA by using poly(dt) oligomers to target and bind to poly(A) tails.RNA fragmentation was performed to account for the size limitation of the sequencing platform.Reverse transcriptase was added to catalyze the reverse transcription of mRNA into cDNA.RNAse was added to prevent the unnecessary ligation or bonding of certain nucleotide pairs and strands.Second-strand synthesis was performed, followed by 3 ′ end ligation with adaptors and adenosine caps.Polymerase chain reaction (PCR) was used to synthesize and amplify the cDNA libraries of every sample.The cDNA libraries were sequenced using the Illumina sequencing platform at Seqmatic in Fremont, CA, USA.FASTQC files were generated from the cDNA libraries for each sample.The FASTQC program was used to assess the quality of the raw data from the files and generate characteristics for each sequence, which included average sequence length, percentage of guanine/cytosine content, the total percentage of deduplicated sequences, and sequences flagged as having poor quality.The Cutadapt program was used to delete adaptor sequences and low-quality bases from the raw sequence read data.The Bowtie2 algorithm in Trinity was used to map the RNA sequence raw reads to the Trinity transcript assembly, generating sequence alignment map (SAM) files, which were then converted to BAM (binary form of SAM) files.Transcript assembly was performed by inputting RNA sequence data from every sample into the TRINITY program, which quantified the number of genes based on the number of detected genes and corresponding isoforms.

BLAT Matching and Annotation of Pinus banksiana Genes
A two-way BLAST-like alignment tool (BLAT) was used to match and align the transcripts from the assembly to the reference genome of Pinus taeda.The characteristics of each transcript were provided, and they included the transcript ID, gene ID, and the corresponding log (E-value) for sequence similarity in relation to the reference genome.Other attributes that were acquired from the BLAT alignment included transcript sequence size, query sequence size, and the percentage of net match for each attribute.Every transcript was mapped to protein sequences in the UniProt database, generating corresponding UniProt IDs.Protein matches with the highest level of similarity were used to annotate genes and assign gene ontology information, such as the functions and cellular localization of the gene product.

Quantification of Gene Expression and Quality Control (QC) Analysis
The RNA-Seq by expectation-maximization (RSEM) abundance estimation method was used to quantify the expression level of each gene/transcript, as well as related isoforms.A quality control assessment was conducted for each read count to validate and assess the number of counts expressed for each gene.Raw reads were filtered and selected for counts of at least 1, 2, 10, 50, or 100.Genes with one read were considered noise.Genes with two or more counts were used as an estimate for the number of genes expressed.Genes with 10 or more counts were considered an adequate indication of the number of genes that had enough reads for downstream statistical analysis.For each treatment group, genes with counts per million (CPM) values of one or higher in at least two samples were included in the downstream analysis.Genes with a CPM value of less than one in at least two samples were considered unexpressed and removed.Normalization factors for raw counts were generated using a trimmed mean of M-values (TMM) from the edge R package (version 4.2) to normalize sample sizes and remove variations between samples.
The normalized read counts underwent log scale transformation using the voom method (log2 scale) from the R limma package.Boxplots of the transformed expression values were generated to show the mean distribution of every sample.Deviation from the mean distribution in a particular sample may indicate variations among experimental conditions, sample contamination, or the batch effect.Samples that deviated significantly from the mean distribution within the same objective group were excluded.
Multidimensional scaling plots were generated to show the amount of clustering for sample groups based on the leading log fold change (logFC) of normalized data.Groups of samples that deviated significantly from other groups of samples were considered differentially regulated.Samples that deviated significantly from the other samples within the same group were considered outliers and were excluded from the downstream analysis.
A heatmap was generated from the logFC of 5000 genes to show the visual relationship of differential gene expression between the samples.Samples that did not have a similar logFC pattern of gene expression as other samples within the same group were considered outliers and were excluded from the downstream analysis.The proportion of raw reads expressed by the top 100 upregulated and downregulated genes was also assessed in every sample to identify potential bottlenecking issues.

Differential Gene Expression (DGE) Analysis of Pairwise Comparisons
The cut-off for pairwise comparisons was calculated to be equivalent to 10 raw counts.A CPM of 0.413 was determined as the minimum threshold required for each pairwise comparison as calculated from the average of the total counts in all of the included samples.Genes that had a CPM higher than the cut-off in at least three samples were included in downstream analysis, whereas genes that did not fulfill these parameters were excluded.Out of 435,293 total genes, 150,739 genes were included in the differential gene expression analysis.The pairwise comparisons of transcript expression were performed at highstringency and low-stringency cut-offs based on the false discovery rate and on p values (0.01) analyses, respectively.The pairwise comparisons were comprised of copper-resistant samples in comparison to nickel-resistant samples and copper-susceptible samples in comparison to nickel-susceptible samples.Differential gene expression expressed as logFC values was evaluated using the R limma package.To assess the interference of sulfate ions on the treatment regimen, pairwise comparisons of expressed genes were also conducted between each sample and the potassium sulfate control.Every gene for each pairwise comparison was annotated using Trinotate and Trinity.Gene ontology was performed by assigning GO terms and gene IDs from available databases to the set of genes for a given pairwise comparison.Genes that could not be annotated were filtered out of the set of annotated genes.Each gene set was run through a plant slim function using the Omicsbox (formally known as BLAST2GO) program.For each pairwise comparison, gene ontology charts were generated to functionally categorize biological processes, metabolic functions, and cellular components.For each functional category, sequences were distributed and organized by the NodeScore of each assigned GO term.

Analysis of Top Differentially Regulated Genes
The top 100 upregulated genes and downregulated genes were ranked for copper RG in comparison to nickel RG.Genes were ranked based on LogFC and the fulfillment of highstringency parameters.This process was repeated for copper SG in comparison to nickel SG.UniProt annotation and a review of the current literature were conducted to characterize genes associated with copper detoxification or tolerance mechanisms.Genes associated with copper or nickel resistance were considered candidate genes for metal resistance.Gene ontology charts were generated to functionally categorize biological processes, metabolic functions, and cellular component localization for the top 100 regulated genes using the aforementioned process in DGE analysis.Charts comprised of the top 50 genes were provided for each pairwise comparison group.

Conclusions
The transcriptome analysis of Pinus banksiana treated with copper and nickel provided a deeper insight into the distinct genetic response of Pinus banksiana to excess Cu and Ni ions.At high stringency, there were 449 DEGs in copper RG in comparison to nickel RG, indicating a similar level of gene expression for the majority of the transcriptome.Out of the 449 genes, various terms in all three main categories were identified as differentially expressed.The decreased proportion of terms associated with DNA processes in the top DEGs suggests the presence of DNA methylation and regulation in the modulation of gene expression.There were not enough DEGs between copper SG and nickel SG to facilitate a GO categorical comparison.Annotation of the top upregulated genes in copper RG compared to nickel RG identified genes and mechanisms that were specific to copper and not to nickel.For biological processes, the biosynthetic process, response to chemicals, and carbohydrate metabolic process comprised the highest proportion of genes.NtPDR, AtHIPP10, and YSL1 were identified as genes associated with copper resistance.Genes encoding the blue copper protein and the 22 kDa Photosystem II protein were associated with the mitigation of photodamage and photosynthesis protection.Various genes related to cell wall metabolism were identified and included genes encoding for HCT, CslE6, MPG, and polygalacturonase.Annotation of the top downregulated genes in copper RG compared to nickel RG identified genes and mechanisms that were specific to nickel and not copper.For biological processes, the majority of gene expression was associated with the response to stress and signal transduction.Various regulatory and signaling-related genes associated with the stress response were identified, and they included UGT, TIFY, ACC, dirigent protein, peroxidase, and glyoxyalase I. Transcriptome analysis revealed the different strategies that are used by Pinus banksiana toward copper and nickel.The identified genes associated with copper resistance can be further researched and implemented in various industries.Additional research is needed to determine the specific functions of various signaling and stress response mechanisms in nickel-resistant plants.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants13071042/s1, Figure S1: Percent distribution of differentially expressed genes (DEGs) for cell compartment localization when copper resistant (RG) and nickel resistant (RG) were compared; Table S1: Damage rating scale and plant classification based on reaction to nickel and copper treatments; Table S2: Top 108 upregulated genes in copper RG in comparison to Nickel RG in Pinus banksiana; Table S3: Top 112 downregulated genes in copper RG in comparison to Nickel RG in Pinus banksiana; Table S4

Figure 1 .
Figure 1.Heatmap of differentially expressed genes in Pinus banksiana seedlings for (a) the copper resistant compared to nickel resistant genotypes and (b) the copper susceptible compared to nickel susceptible genotypes.Differentially expressed gene values for each pairwise comparison were based on the Log2 normalized fold change.Red data represent different levels of upregulation, whereas blue data represent different levels of downregulation.The labels Nir57, Nir30, and Nir5 represent seedlings from the nickel-resistant genotypes.The labels Cur872, Cur33, and Cur67 represent seedlings from the copper-resistant genotypes.The labels Nis15, Nis31, and Nis58 represent seedlings from the nickel-susceptible genotypes.The labels Cus42, Cus16, and Cus542 represent seedlings from the copper-susceptible genotype.

Figure 1 .
Figure 1.Heatmap of differentially expressed genes in Pinus banksiana seedlings for (a) the copper resistant compared to nickel resistant genotypes and (b) the copper susceptible compared to nickel susceptible genotypes.Differentially expressed gene values for each pairwise comparison were based on the Log2 normalized fold change.Red data represent different levels of upregulation, whereas blue data represent different levels of downregulation.The labels Nir57, Nir30, and Nir5 represent seedlings from the nickel-resistant genotypes.The labels Cur872, Cur33, and Cur67 represent seedlings from the copper-resistant genotypes.The labels Nis15, Nis31, and Nis58 represent seedlings from the nickel-susceptible genotypes.The labels Cus42, Cus16, and Cus542 represent seedlings from the copper-susceptible genotype.

Figure 2 .
Figure 2. Volcano plot of differentially expressed genes in Pinus banksiana seedlings for (a) the copper resistance genotype compared to nickel resistance genotypes and (b) the copper susceptible genotypes compared to nickel susceptible genotypes.Brown data points represent upregulated gene expression, whereas blue data points represent downregulated gene expression relative to the susceptible genotypes.Grey points indicate no significant difference from the nickel-susceptible genotypes.LogFC is the log-fold change of the copper-susceptible genotypes relative to the nickelsusceptible genotypes.Log10(FDR) is the log10 of the false discovery rate.The barrier between the nonsignificant data points (grey) and the differentially regulated genes (orange or blue) signifies a false discovery rate of 0.05 (two-fold).

Figure 2 .
Figure 2. Volcano plot of differentially expressed genes in Pinus banksiana seedlings for (a) the copper resistance genotype compared to nickel resistance genotypes and (b) the copper susceptible genotypes compared to nickel susceptible genotypes.Brown data points represent upregulated gene expression, whereas blue data points represent downregulated gene expression relative to the susceptible genotypes.Grey points indicate no significant difference from the nickel-susceptible genotypes.LogFC is the log-fold change of the copper-susceptible genotypes relative to the nickelsusceptible genotypes.Log10(FDR) is the log10 of the false discovery rate.The barrier between the nonsignificant data points (grey) and the differentially regulated genes (orange or blue) signifies a false discovery rate of 0.05 (two-fold).

Figure 3 .
Figure 3. Percent distribution of differentially expressed genes (DEGs) for (a) biological processes for copper resistant (RG) compared to nickel resistant (RG) and (b) molecular functions for copper RG compared to nickel RG.DEGs from copper RG compared to nickel RG were annotated and allocated to terms within the biological processes category and molecular function category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other."

Figure 4 .
Figure 4. Percent distribution of the (a) top 100 upregulated genes in biological processes for copper resistant (RG) compared to nickel resistant (RG) and the (b) top 100 downregulated genes in biological processes for copper RG compared to nickel RG.The top upregulated and downregulated genes

Figure 4 .
Figure 4. Percent distribution of the (a) top 100 upregulated genes in biological processes for copper resistant (RG) compared to nickel resistant (RG) and the (b) top 100 downregulated genes in biological processes for copper RG compared to nickel RG.The top upregulated and downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the biological process category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other."

Figure 5 .
Figure 5. Percent distribution of the (a) top 100 upregulated genes in molecular function for copper resistant RG compared to nickel resistant (RG) and the (b) top 100 downregulated genes in molecular function for copper RG compared to nickel RG.The top upregulated and downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the molecular function category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other."

Figure 5 .
Figure 5. Percent distribution of the (a) top 100 upregulated genes in molecular function for copper resistant RG compared to nickel resistant (RG) and the (b) top 100 downregulated genes in molecular function for copper RG compared to nickel RG.The top upregulated and downregulated genes from copper RG compared to nickel RG were annotated and allocated to terms within the molecular function category using Omicsbox/Blast2GO.Terms that had a total percentage of expressed genes lower than 2% were collectively categorized under the term "other."
been shown to nonspecifically facilitate the transport of other heavy metals, such as lead.

Table 1 .
Differentially expressed genes between copper and nickel-resistant genotypes.
FDR = False discovery rate.

Table 2 .
Differentially expressed genes between copper and nickel-susceptible genotypes.
FDR = False discovery rate.
Table 5 shows the top upregulated genes when copper

Table 3 .
Top upregulated genes in copper resistant (RG) in comparison to nickel resistant (RG) in Pinus banksiana.

Table 5 .
Top upregulated genes in copper susceptible (SG) in comparison to nickel susceptible (SG) in Pinus banksiana.

Table 6 .
Top downregulated genes in copper susceptible (SG) in comparison to nickel susceptible (SG) in Pinus banksiana.