The research was conducted in three commercial vineyards located in three regions of central-southern Italy: Campania (CAM), Molise (MOL) and Sicily (SIC), each characterized by different altitudes and latitudes. The vineyard geo-references and elevations are detailed in Online Resource 1: Supplementary Table S1 . For each region, two vineyards were chosen. Two black (Vitis vinifera L.) wine grape cultivars, Aglianico (A) and Cabernet Sauvignon (C), were grown at each vineyard. All the vines, which were grafted onto 140 Ruggieri rootstocks, were approximately fifteen years old and were planted between 2008 and 2010. In Campania and Molise, vines were planted in north−south rows, whereas in Sicily, they were planted on gentle slopes. In Molise, both cultivars were planted at a spacing of 1.20 m (within the row) × 2.90 m (between rows), and both were trained with the simple Guyot method, with a formation height of 60 cm. In Campania, Aglianico vines were spaced at 1.50 m (within the row) × 2.90 m (between rows), and Cabernet Sauvignon vines were spaced at 1.0 m (within the row) × 2.70 m (between rows); these vines were also trained using the simple Guyot method, with a formation height of 60 cm. In Sicily, Aglianico vines were spaced at 1.10 m (within the row) × 1.10 m (between rows), and Cabernet Sauvignon vines were spaced at 1.10 m (within the row) × 1.30 m (between rows). Both cultivars were bush trained at a height of 0.5 m, with two to six main branches, each of which was spur-pruned to one spur and two buds per spur. The experimental design was based on sampling three independent randomized plots of five rows, each containing 30 vines, to ensure representativeness. All the measurements were carried out on seven ‘index’ vines per block (Ferlito et al., 2023).

Berries were collected at the commercial ripening stage during the 2021 growing season (8 September 2021 for the Campania site; 24 September 2021 for the Sicily site; and 20 October 2021 for the Molise site). At least ten berries were randomly selected from each index vine, avoiding those with visible damage and/or signs of pathogen infection, and pooled with berries from the other index plants from different blocks. For subsequent analyses, three independent pools (biological replicates) of 30 whole berries each were selected, immediately frozen in liquid nitrogen, and stored at -80°C. Each analysis included six experimental samples: Aglianico from Campania (A_CAM), Aglianico from Molise (A_MOL), Aglianico from Sicily (A_SIC), Cabernet Sauvignon from Campania (C_CAM), Cabernet Sauvignon from Molise (C_MOL) and Cabernet Sauvignon from Sicily (C_SIC).

Weather course data were collected during the entire growing season by a recording climatic station located inside each vineyard. The daily maximum air temperature (DMAT), relative humidity (RH), air pressure (AP), dew point (DP), total precipitation (TP), evapotranspiration (ET) and hours of sun radiation (HSR) were measured. Weather course indices for WGCNA were calculated as follows: for DMAT, RH, AP and DP, the average values of the thirty days preceding the sampling date were used, whereas for TP, ET and HSR, the sum values of the thirty days preceding the sampling date were used.

The total soluble solids (TSS), pH, titratable acidity (TA), polyphenol (PP) content, total anthocyanin (ANTH) content and maturity index (MI) were measured. Free-run juice was used to determine TSS measured by a digital refractometer with temperature correction (RX-5000 Atago Co., Ltd., Bellevue, WA, USA). The pH and titratable acidity were determined with an automatic titrator (Titrino model 798, Metrohm, Riverview, FL, USA) on 5.0 mL juice samples titrated against 0.1 M NaOH up to pH 8.2. TA was expressed as g/L of tartaric acid equivalents. The MI was estimated at harvest (soluble solids -°Brix and pH) using the following formula (Lana et al., 2021):

The PP content was determined using the Folin–Ciocalteu reagent assay (Singleton et al., 1999) and expressed as mg/kg of grapes. The total ANTH content was measured according to the methods of Lo Cicero et al. (2016) and expressed as mg/kg of fresh weight as reported by Ferlito et al. (2014), Ferlito et al. (2020). Principal component analysis (PCA) of these quality traits was carried out with the R function ‘prcomp’.

Analysis of variance (ANOVA) was performed with the R function ‘aov’ on the differences among quality traits for each genotype. A post hoc analysis based on the Tukey honestly significant differences (Tukey HSD) test was performed with the ‘TukeyHSD’ R function at significance levels (p values) of 0.05, 0.01 and 0.001.

Ten berries from each biological replicate, which were kept frozen by the continuous addition of liquid nitrogen, were ground using a precooled mortar and pestle. Total RNA was isolated from ground whole berries, excluding seeds, using a slight modification of the method described by Japelaghi and coworkers (Japelaghi et al., 2011). Briefly, because grape berries are rich in water and secondary metabolites that might interfere with RNA yield and quality, a scale-up of the ground tissue was needed: in detail, 1 g of sample was added to 10 ml of extraction buffer. Next, a chloroform:isoamyl alcohol (24:1) extraction was performed three times. RNA degradation and DNA contamination were monitored by 1% agarose gel electrophoresis. The RNA purity and concentration were assayed using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) (Russo et al., 2021; Strano et al., 2014). Before sequencing, RNA integrity was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) (Santoro et al., 2022).

Library preparation, clustering and sequencing were performed by Novogene Co., Ltd., UK (25 Cambridge Park, Milton Road, Cambridge, CB4 OFW, United Kingdom). One µg of RNA was used as input material for library preparation (eighteen libraries: three biological replicates × two varieties × three sites). The sequencing libraries were generated using the NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (New England Biolabs, Ipswich, MA, USA), as reported in Sicilia et al. (2019). After cluster generation, the library preparations were sequenced on the Illumina HiSeq2000 platform to generate paired-end reads with a size of 2×150 bp. Raw reads in fastq format were first processed using in-house Perl scripts (Novogene Co., Ltd., UK). Clean data were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads (Q score below 5 for more than 50% of the bases) (Online Resource 2: Supplementary Table S2 ). At the same time, the Q20, Q30, GC content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on high-quality clean data.

De novo transcriptome assembly was performed with Trinity software (version 2.6.6) with the parameters min_Kmer_Cov=3 and min_glue=4 (Grabherr et al., 2011). Hierarchical clustering was carried out with Corset (version 4.6) to remove redundancy (parameter -m 10) so that the longest transcript of each cluster was selected as the Unigene (Davidson and Oshlack, 2014). The assembly assessment and gene prediction were evaluated by comparing the Unigenes to the set of Embryophyta genes using the BUSCO quality assessment tool coupled with the OrthoDB (version 9.0) database of orthologues (Simão et al., 2015). The gene functional annotation was obtained using seven different databases: National Centre for Biotechnology Information (NCBI) non-redundant protein sequences (Nr, Diamond software, version 0.8.22, e-value threshold 1e-5) (Buchfnk et al., 2015), NCBI non-redundant nucleotide sequences (Nt, NCBI blast software, version 2.9.0, e-value threshold 1e-5), Protein Family (Pfam, hmmscan software, HMMER version 3.1, e-value threshold 0.01) (Finn et al., 2011), Cluster of Orthologous Groups of Proteins (KOG/COG, Diamond software, version 0.8.22, e-value threshold 1e-5) (Buchfnk et al., 2015), SWISS-PROT (Diamond software, version 0.8.22, e-value threshold 1e-5) (Buchfnk et al., 2015), Kyoto Encyclopedia of Genes and Genome (KEGG, Diamond and KAAS software, version 0.8.22, e-value threshold 1e-5) (Buchfnk et al., 2015; Moriya et al., 2007) and Gene Ontology (GO, blast2GO software, version b2g4pipe_v2.5, e-value threshold 1e-6) (Götz et al., 2008). The Unigene coding sequences were aligned to the Vitis vinifera genome PN40024 v4 using the Grapedia BLAST tool (https://grapedia.org) PN40024.v4.1_REF_prot database with an e-value threshold of 1e-5, followed by annotation as described by Grimplet et al. (2012). The annotated sequences not corresponding to the plant species were filtered out.

The gene expression level was estimated with RSEM software (version 1.2.28) by mapping each clean read back onto the assembled transcriptome, and the read counts for each gene were obtained from the mapping results (Li and Dewey, 2011). Furthermore, the read counts of each gene were used as input data for DESeq2 (version 1.26, padj ≤ 0.05) to obtain differentially expressed genes (DEGs) (Santoro et al., 2022). Six comparisons were made to identify the set of DEGs of each variety at the three sites (A_SIC vs. A_CAM; A_MOL vs. A_CAM; A_SIC vs. A_MOL; C_SIC vs. C_CAM; C_MOL vs. C_CAM; and C_SIC vs. C_MOL). An adjusted p value cut-off of 0.05 and a |log₂fold change| (log₂FC) threshold of 1 were the criteria used to filter significantly up- and downregulated genes. A correlation analysis was performed to demonstrate experimental repeatability and to reveal differences in gene expression among samples. PCA was performed using R language, considering the read counts of each sample as input data, including the biological replicates.

Total RNA (2.5 µg) was reverse transcribed using a SuperScript™ Vilo™ cDNA synthesis kit from Thermo Fisher Scientific according to the manufacturer’s instructions. To validate the reliability of RNAseq in determining gene expression levels, real-time qRT−PCR was performed for 10 random DEGs using PowerUp SYBR Green Master mix from Thermo Fisher Scientific and the Rotor-Gene Q 2plex detection system (Qiagen, Venlo, Netherlands) (Online Resource 3: Supplementary Table S3 ). All the genes were normalized to Vitis vinifera ubiquitin-60S ribosomal protein L40-2 (LOC100253716), which is a suitable housekeeping gene (Cheng et al., 2021b). All reactions were performed in duplicate, and the fold change values were calculated using the 2^−ΔΔCT method.

Both GO and KEGG enrichment analysis of DEGs separated for each variety were performed using the ShinyGO V0.80 online tool (Ge et al., 2020). To define the genotype-specific transcriptomic response at each site, fragments per kilobase million (FPKM) values were used to identify genes specifically expressed at each site by each genotype. Six different gene lists were constructed (A_CAM, A_MOL, A_SIC, C_CAM, C_MOL, and C_SIC) by filtering genes that were expressed by a specific genotype at the site (FPKM ≠ 0) and had FPKM values at the other sites equal to zero. Genes of interest (expressed in only one site/genotype) were selected from the ten highest FPKM values in each cluster list.

A co-expression analysis of the gene expression data, weather course indices and quality traits was conducted. The co-expression analysis was performed using the Weighted Gene Co-expression Network Analysis (WGCNA) package in R (Langfelder and Horvath, 2008) to obtain hierarchical clustering and identify co-expressed genes (hub genes). Specifically, an adjacency matrix was created using the FPKM values of all the DEGs. The pickSoftThreshold() function was used to determine the optimal soft-thresholding power (Langfelder and Horvath, 2008). For each analysis, the lowest power for which the scale-free topology fit index was 0.90 was used (Online Resource 4: Supplementary Figure S1 ). The specific WGCNA parameters were set as follows: Soft powers β = 20 were selected using the function pickSoft Threshold; WGCNA ‘mergeCutHeight’ was set at 0.25. The adjacency matrix was transformed into a topological overlap matrix (TOM) and the corresponding dissimilarity matrix (1-TOM). Afterwards, a hierarchical clustering dendrogram of the 1-TOM matrix was constructed to classify genes with similar expression levels into different gene co-expression modules. Modules were merged by using a criterion of MEDissThres = 0.25. Finally, the relationships between each module and either the weather course indices or quality traits that were significantly different among sites for each genotype were estimated by Pearson’s correlation using the module eigengene values. Modules with high correlation coefficients and a correlation padj ≤ 0.05 were selected for subsequent analysis.

The weather course indices obtained as described in the Materials and Methods section are reported in Table 1 . The highest DMAT, AP, DP, ET and HSR values were recorded in Campania, as well as the lowest values of RH and TP. As expected, the highest TP was recorded in Sicily because of the seasonal rainfall pattern that is typical of the eastern slopes of Mt. Etna, where the vineyards were located. In addition, the lowest ET value was also recorded in Sicily. Finally, the highest RH and lowest DMAT, AP, DP and HSR values were recorded at the Molise site, thus indicating opposite weather conditions compared with those at the Campania site.

Table 2 shows the mean values of the quality parameters of grape berries at the commercial ripening stage at the three sites. On average, Cabernet Sauvignon had a higher PP content than that of Aglianico at all three sites and the highest ANTH content, except at the CAM site, where the lowest value was recorded. Moreover, both Cabernet Sauvignon and Aglianico presented higher contents of TSS at the CAM and MOL sites than at the SIC site. Aglianico presented a significantly greater TA at the MOL site ( Table 2 ). Interestingly, at all three sites, the Cabernet Sauvignon genotype exhibited significantly higher MI values than those of the Aglianico genotype ( Table 2 ). Figure 1 shows a PCA in which all the measured qualitative traits (PP, pH, TSS, TA and ANTH) were used as the variables. A total of 77.2% of the variability was explained by the first two components (dim1 and dim2), and the three variables that best differentiated the samples were TA, pH and PP. Interestingly, both the genotypes cultivated in Sicily clustered separately from the same genotypes cultivated at the other two sites. PP, TSS and MI were the traits that best separated the Sicilian site from the other sites. In CAM and MOL, the two genotypes were similar in terms of quality traits (TSS and MI), except for A_MOL, which had higher TA values than those of A_CAM, C_MOL and C_CAM.

RNAseq was used to comprehensively identify the transcriptional profiles of Aglianico and Cabernet Sauvignon berries at the three sites in central-southern Italy. After library sequencing, the raw reads were filtered to remove adapter-based or poor-quality reads. A total of 598 million clean reads, representing 98.3% of the total reads, were obtained, with a Q30 value and a GC content of 91.54% and 46.62%, respectively ( Table 3 ). The clean read de novo assembly yielded 226,498 transcripts and 97,442 Unigenes with N50 lengths of 2753 bp and 1881 bp, respectively ( Table 3 ), which is consistent with previously reported N50 values (Santoro et al., 2022, 2023; Sicilia et al., 2022) and indicates good contiguity of the transcriptome. To assess assembly consistency, filtered unique reads were mapped back to the reconstructed transcriptome, and the average read mapping rate determined using Bowtie2 alignment software was 77.3% (mean value between 75.25 and 78.43) ( Table 3 ). Among the 1440 groups searched with BUSCO, 63.1% (909 groups) were complete single-copy orthologues ( Table 3 ). Functional annotation of the Unigenes was conducted by performing BLAST searches against public databases, such as the NCBI, Pfam, KOG/COG, SWISS-PROT, KO, and GO databases ( Table 4 ). A total of 94,058 Unigenes were annotated in at least one database, corresponding to 96.52% of the total Unigenes. Among them, 64,056 (65.73%) and 87,354 (89.64%) Unigenes showed identity with the sequences in the Nr and Nt databases, respectively. The distributions of Unigenes homologous to the sequences in the KO, SWISS-PROT, Pfam, GO, and KOG databases were 21.61%, 45.03%, 41.87%, 41.87% and 18.50%, respectively ( Table 4 ). With respect to the GO annotation, most of the Unigenes were annotated as “cellular process” (GO:0009987), “metabolic process” (GO:0008152), “biological regulation” (GO:0065007), “localization” or “regulation of biological process” (GO:0050789) within the Biological Process category. “Cellular anatomical entity” (GO:0110165), “Intracellular” and “Protein-containing complex” (GO:0032991) were annotated from the Cellular Component. Finally, “Binding”, “Catalytic activity” and “Transporter activity” (GO:0005215) were annotated from the Molecular Function category ( Figure 2 ). In the KOG functional annotation, the three most annotated categories were “Posttranslational modification, protein turnover, chaperones”, “General function prediction only” and “Translation, ribosomal structure and biogenesis” (Online Resource 5: Supplementary Figure S2 ). Finally, in the KEGG pathway functional annotation, the pathways with the highest percentage of annotated genes were “Translation”, “Signal transduction” and “Carbohydrate metabolism” ( Figure 3 ). The functional annotation analysis revealed that a total of 39,947 out of 94,058 annotated Unigenes belonged to the grape fungal community, as characterized by Iorizzo et al. (2023). These Unigenes were removed, and the subsequent analyses were conducted on 57,495 Unigenes.

With the aim of deciphering the transcriptomic profile of each sample, the Unigene read counts of the two genotypes cultivated at the three sites were used as inputs for PCA and explained 57% of the total variability ( Figure 4 ). PCA revealed a sharp separation of Aglianico (circles) from Cabernet Sauvignon (triangles), suggesting that the differences in the transcriptomic profiles are due mainly to the genotype ( Figure 4 ). Moreover, the green and red triangles grouped at the bottom of the figure indicate that Cabernet Sauvignon is less influenced by the growing site (CAM and MOL) than the Aglianico samples are, which shows a greater variance in PC2. However, the weather variable at the SIC site greatly influences the plant transcriptomic profiles of both genotypes (light blue circles and triangles) ( Figure 4 ).

The characterization of the Aglianico and Cabernet Sauvignon transcriptomes was accomplished by identifying those Unigenes with expression levels that changed as a function of the growing site. According to the experimental design, a total of 11,887 DEGs among the two genotypes and the three sites were identified, 1937 of which were differentially expressed in A_CAM vs. A_MOL (847 upregulated and 1090 downregulated), 3397 in A_SIC vs. A_MOL (1622 upregulated and 1775 downregulated), 4711 in A_SIC vs. A_CAM (2221 upregulated and 2490 downregulated), 1914 in C_CAM vs. C_MOL (888 upregulated and 1026 downregulated), 2328 in C_SIC vs. C_MOL (1655 upregulated and 673 downregulated) and 2504 in C_SIC vs. C_CAM (1860 upregulated and 644 downregulated) ( Figure 5 ). The greatest number of DEGs was found in the Aglianico genotype ( Figure 5 ). Furthermore, the comparisons that included the SIC site presented the greatest number of DEGs for both genotypes ( Figure 5 ). Validation of the RNAseq experiment was performed by measuring the expression levels of ten selected DEGs by quantitative real-time PCR (qRT−PCR) (Online Resource 3: Supplementary Table S3 and Online Resource 6: Supplementary Figure S3 ). The results revealed high congruence between the RNAseq and qRT−PCR results (coefficient of determination R² = 0.93), which confirms the high reliability of the RNAseq analysis in the quantification of gene expression.

The enrichment of AGL and CAB DEGs in the GO and KEGG databases was performed to identify the principal pathways and biological processes involved in the transcriptome reprogramming of each genotype. Among the AGL enriched GO terms of the Biological Process, “Oligopeptide transmembrane transport”, “Amide transport” and “Organic hydroxy compound metabolic process” were the three most represented categories ( Figure 6A ). Differently, “Ethylene-activated signaling pathway”, “Organic acid catabolic process” and “Hormone-mediated signaling pathway” were the three most represented categories in CAB ( Figure 6B ). The main KEGG pathways in AGL were associated with the “Tropane piperidine and pyridine alkaloid biosynthesis” and “Isoquinolone alkaloid biosynthesis” categories, followed by “Sphingolipid metabolism” and “Tyrosine metabolism” ( Figure 6C ). As regards the KEGG pathways in CAB, “DNA replication”, “Homologous recombination” and “Mismatch repair” are the most represented, followed by “Glycine serine and threonine metabolism” and “Phenylalanine tyrosine and tryptophan biosynthesis” ( Figure 6D ).

The transcriptomic response of each genotype at the three distinct sites was evaluated by recruiting clusters that were expressed in a single genotype per site (FPKM ≠ 0) and were completely unexpressed at the other sites (FPKM = 0). In particular, a total of 185 genes were specifically expressed in A_CAM, 90 genes in A_MOL, 171 genes in A_SIC, 226 genes in C_CAM, 165 genes in C_MOL, and 206 genes in C_SIC, all of which represent a specific transcriptomic signature of the genotype at that particular site. Most of these genes were not annotated and encoded hypothetical and/or uncharacterized proteins (data not shown). In Table 5 , some genes of interest are listed along with their cluster ID, annotation, accession number, e-values and specific site. In A_CAM, probable genes encoding serine/threonine-protein kinase WNK5 (cluster-39.0) (Su et al., 2024) and DExH-box ATP-dependent RNA helicase (cluster-27627.0) involved in ribosome biogenesis (Liu and Imai, 2018) were specifically expressed. The gene encoding vegetative cell wall protein gp1 (cluster 27903.60085), which functions as a major component of the outer cell wall layer (Niu et al., 2020), was among the genes specifically expressed in A_MOL, along with the gene encoding UDP-glycosyltransferase 85A2-like (cluster-1161.0), which is involved in the key step of secologanin biosynthesis (Jin et al., 2022). The genes encoding phenylalanine N-monooxygenase (cluster-27903.31192), which is part of the pathway that converts phenylalanine to glucosinolate; glucotropaeolin (Wittstock and Halkier, 2000); and myb-related protein 315 (cluster-27903.28077), which confers cold and drought tolerance (Fang et al., 2024), were exclusively expressed in A_SIC. Clusters specific to C_CAM encode proteins homologous to the homeobox-leucine zipper protein ATHB-15 (cluster-27903.3260) (Li et al., 2020) and an endo-1,4-beta-mannosidase 5-like (cluster-27903.633) (Rahmani et al., 2017). The transcription factor JUNGBRUNNEN 1 (cluster-17927.0), a regulator of drought tolerance in tomatoes (Thirumalaikumar et al., 2018), and MLP-like protein 43 (cluster-27903.61696), which confers drought tolerance (Fujita et al., 2021), are among the genes specifically expressed in C_MOL. Finally, the specific expression of an APETALA2/ethylene responsive factor (cluster-20558.0), which is involved in berry firmness (Wong et al., 2016), and an SRG1 protein (cluster-21223.0) (Cui et al., 2018) were detected in C_SIC.

Co-expression networks of weighted genes associated with the weather course parameters and quality traits (DMAT, RH, AP, DP, TP, ET, and HSR values as weather course parameters; TA, TSS, ANTH, MI and PP values as quality traits) of each site were constructed by WGCNA based on 11,887 DEGs. These genes were grouped into 20 co-expressed modules ( Figure 7A, B ). Each set of highly correlated genes corresponded to a branch of the tree ( Figure 7A ). With respect to the weather course parameters, DMAT (positive correlation, r² = 0.9), RH (negative correlation, r² = -0.9), AP (positive correlation, r² = 0.86), DP (positive correlation, r² = 0.89) and HSR (positive correlation, r² = 0.92) were significantly correlated with the “lightcyan1” module (80 genes), and both TP (negative correlation, r² = -0.85) and ET (positive correlation, r² = 0.91) were significantly correlated with the “indianred3” module (1078 genes) ( Figure 7B ). With respect to the quality traits, ANTH, TA and MI were not significantly correlated with any module (data not shown). Conversely, co-expression modules that were highly correlated with TSS and PP were identified; in particular, the “deeppink” module (1400 genes) was significantly correlated with TSS (positive correlation, r² = 0.9), and two modules, namely, “chocolate3” (113 genes) (positive correlation, r² = 0.91) and “plum1” (557 genes) (negative correlation, r² = -0.81), were significantly correlated with PP ( Figure 7B ). The genes in these latter modules were filtered according to the highest intramodular connectivity (hub genes, module membership [MM] > 0.65 and gene significance [GS] > 0.65), as these genes might represent points of biological interest in defining plant adaptability to the environment (Online Resource 7: Supplementary Figure S4 ).

To describe the biological functions of the hub genes, the functional categories annotated in the PN40024 v4 Vitis vinifera genome were considered. Although the modules included a substantial number of genes, most of them encode unknown or uncharacterized proteins, resulting in very few enriched genes in each functional category. As shown in Table 6 , TSS was positively correlated with the deeppink module, with a total of 326 genes involved in several functional categories, such as cellular component organization and biogenesis (10 genes), nucleic acid metabolism (18 genes), protein metabolism and modification (14 genes), regulation of transcription (31 genes), hormone signaling (14 genes), signaling pathways (12 genes), and transport (19 genes). PP was positively correlated with the chocolate3 module (32 genes), containing genes involved in cellular component organization and biogenesis (3 genes), amino acid metabolism (4 genes), regulation of transcription (6 genes), and signaling pathways (2 genes). Interestingly, PP was also negatively correlated with the plum1 module (121 genes), including genes involved in amino acid metabolism (6 genes), protein metabolism and modification (9 genes), phenylpropanoid metabolism (3 genes), terpenoid metabolism (3 genes), regulation of transcription (9 genes), hormone signaling (5 genes), and transport (10 genes). A substantial number of genes are included in the indianred3 module (255 genes), which is negatively correlated with TP; most of the genes are involved in cytoskeleton organization and biogenesis (2 genes), lipid metabolism (3 genes), and the regulation of transcription (4 genes) ( Table 6 and Online Resource 8: Supplementary Table S4 ). ET was positively correlated with the indianred3 model (347 genes), with genes involved in lipid metabolism (3 genes), protein metabolism and modification (3 genes), and signaling (6 genes) ( Table 6 and Online Resource 8: Supplementary Table S4 ). As previously described, most of the weather course traits were significantly correlated with the lightcyan1 module (DMAT, RH, AP, DP and HSR with 18, 16, 21, 15 and 31 genes, respectively) ( Table 6 and Online Resource 8: Supplementary Table S4 ). These genes are related mainly to the regulation of transcription (4 genes), signaling (4 genes) and transport (5 genes). Eight genes are shared among the weather course traits: a cation transporter (autoinhibited Ca²⁺-ATPase 10, ACA10), a plasma membrane protein [Vitvi11g01176]; a cation exchanger (cation exchanger 7, CAX7) [Vitvi06g00999_t001]; SHATTERPROOF 2 [Vitvi12g00019_t002]; a transcription factor (DExH-box ATP-dependent RNA helicase, DExH8 [Vitvi07g00315_t001]); a disease resistance protein (PREDICTED: TMV resistance protein N isoform X1 [Vitvi18g01746_t001], putative disease resistance protein [Vitvi13g04668_t001]); and a serine/threonine kinase (Ribosomal-protein S6 kinase p70 [Vitvi17g00727_t001]). As shown in Figure 8 , the expression levels of these genes in both genotypes and at all three sites increased with increasing HSR and DMAT and decreasing RH, indicating a decreasing expression pattern from the Campania site to the Molise site.

A total of 58 genes encoding transcription factors were retrieved since all the modules were significantly correlated with a trait (Online Resource 9: Supplementary Table S5 ). Most of these genes belong to the MYB, bHLH, AUX/IAA, AP2/ERF, C2C2 and WRKY families ( Figure 9 ). Table 7 lists the genes whose expression greatly changed among sites. Specifically, MYB61 regulates a specific set of target genes by binding DNA to the AC within the cis-element 5’-ACCTAC-3’. Higher expression of this gene was observed in both genotypes grown at the CAM site, which had higher ET and HSR values and lower TP values. Cluster-27903.36947 encodes multiprotein-bridging factor 1C (MBF1C), a transcriptional coactivator that stimulates transcriptional activity by bridging regulatory proteins and TBP, recruiting TBP to promoters occupied by DNA-binding regulators. This gene was more highly expressed in Aglianico than in Cabernet Sauvignon at all sites. Cluster-27903.38074 encodes Arabidopsis pseudo-response regulator 7 (APRR7) involved in the positive and negative feedback loops of the circadian clock. This gene was highly expressed at the MOL site, where both the DMAT and HSR were lower than those at the other sites. Cluster-27903.17293 encodes a CPC (CAPRICE) MYB transcription factor that determines the fate of epidermal cell differentiation. This gene was more highly expressed in Aglianico than in Cabernet Sauvignon at sites with low TP (Campania and Molise).