Two different accessions of the genus Lactuca were included in this study: a commercial variety, the red-leaf lettuce ‘Romired’, and a wild relative species, Lactuca homblei De Wild. They were selected from a previous drought stress experiment in which two irrigation regimes, control (C, week 1: 1350 mL, weeks 2-3: 2100 mL/each) and deficit irrigation (DI, weeks 1-3: 0 mL), were tested in two consecutive years (Medina-Lozano et al., 2024). Three biological replicates for the two accessions in each of the two conditions (C and DI) from the experiment carried out in winter 2020/2021 were used to proceed with the transcriptomic studies.

Total RNA extraction from lyophilized samples coming from 12 samples (2 accessions x 2 irrigation regimes x 3 biological replicates) was performed using the NZY Total RNA Isolation kit (NZYtech Lda.-Genes and Enzymes, Lisbon, Portugal) as described before (Medina-Lozano et al., 2023). RNA was treated with DNase using the TURBO DNA-free™ kit (Invitrogen, Waltham, MA, USA), following the manufacturers’ instructions. RNA quantity and purity were assessed with the Eukaryotic Total RNA Nanobioanalyzer Assay in a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

The obtained RNA samples from ‘Romired’ and L. homblei were used to perform the RNA-seq. They were processed to build a total of 12 strand-specific cDNA libraries. Sequencing of the libraries was performed in both directions with a NovaSeq 6000 S1 instrument (Illumina, San Diego, CA, USA) using the TruSeq Stranded mRNA protocol (Illumina) to obtain between 36 and 111 strand-specific pair-end reads of 100 base pair (bp) lengths per sample. Sequencing was carried out at the National Centre for Genomic Regulation (CNAG-CRG, Barcelona, Spain).

Sequences were analysed using the Galaxy tool (The Galaxy Community, 2022). Adapter sequences were removed by processing the reads from the 12 individual datasets using Trimmomatic (Galaxy version 0.38.1) (Bolger et al., 2014). RNA-seq data alignment to the lettuce reference genome Lactuca sativa ‘Salinas’ v8 (Reyes-Chin-Wo et al., 2017) was performed using HISAT2 (Galaxy version 2.2.1+galaxy0) (Kim et al., 2015), with a maximum intron length set at 20,000 bp. The Picard tools (http://broadinstitute.github.io/picard) MarkDuplicates (Galaxy version 2.18.2.2) and FixMateInformation (Galaxy version 2.18.2.1) were used to filter out the optical duplicates and to mate-pairs, respectively. featureCounts (Galaxy version 2.0.1+galaxy2) (Liao et al., 2013) was used to generate read counts using the gene annotation available in the literature (Reyes-Chin-Wo et al., 2017).

Analysis of differential gene expression between treatments (C and DI) within each of the two accessions was conducted using edgeR (Galaxy version 3.36.0+galaxy0) (Robinson et al., 2009). Genes were considered to be differentially expressed when values of |log₂(FC, fold change)|>1 and FDR (False Discovery Rate)<0.05 (adjusted p-value via the Benjamini-Hochberg method).

Venn diagrams were performed with DEG datasets using the R stats package VennDiagram (https://CRAN.R-project.org/package=VennDiagram). GO (Gene Ontology) enrichment analyses were carried out using the tool GOEnrichment from Galaxy platform (Galaxy version 2.0.1) (Faria, 2017), with p-value cut-off < 0.05 and using Benjamini-Hochberg multiple test correction. Enriched GO terms involving three categories, biological processes, cellular components, and molecular functions, were evaluated. GO terms of the DEGs were obtained from predicted data using information available in the literature (Reyes-Chin-Wo et al., 2017). Heatmaps were constructed using gplots (https://CRAN.R-project.org/package=gplots) and ggplot2 (Wickham, 2009) R stats packages.

The selection of genes for expression data validation was based on different criteria. First, genes were filtered out for values of |log₂(FC)|>4 and FDR<0.05 (substantial increase or decrease in expression levels). Among them, those exhibiting high and significant correlation (both positive and negative) with, first, anthocyanin content and, second, drought stress treatment, were selected. Finally, gene functions related to both anthocyanin content and/or response to different stresses were also taken into account for the selection of a total of 19 DEGs.

Correlations between gene expression and both anthocyanin content and treatments were established through weighted gene co-expression network analysis (WGCNA), which was conducted using the R stats package WGCNA (Langfelder and Horvath, 2008). Normalised RNA-seq data of all genes were used for the WGCNA, except for those with a very low expression among the no DEGs (i.e., less than 5 reads per sample in the three biological samples of each group (C and DI)). Data from both species, L. sativa (cultivated lettuce ‘Romired’) and L. homblei, were analysed separately.

Total RNA was extracted from each of the 12 samples described above. Subsequently, mRNA was purified using the Dynabeads mRNA DIRECT™ kit (Invitrogen) and cDNA was synthesized using the NZY M-MuLV First-Strand cDNA Synthesis, separate oligos kit (NZYTech) as described before (Medina-Lozano et al., 2023).

Specific pairs of primers for each of the 19 selected DEGs ( Supplementary Table S1 ) were designed using OLIGO software version 6.45 (Cascade, CO, USA) from a consensus sequence of the two species under study, L. sativa and L. homblei, excluding any ambiguity in the sequences. Real-time qPCR reactions were performed on a StepOnePlus™ System (Applied Biosystems, Waltham, MA, USA) with two technical replicates per each of the three biological replicates. Each reaction contained 1 µL of 1:5 diluted cDNA, 0.40 µM of forward and reverse primers (Integrated DNA Technologies, IDT, Coralville, Iowa, USA), and 1x NZYSupreme qPCR Green Master Mix, ROX plus (NZYTech) in a final volume of 11 µL. The amplification conditions were: 2 min at 95°C, 40 cycles of 5 s at 95°C, 15 s at 52-66°C ( Supplementary Table S1 ) and 30 s at 72°C, followed by the melting curve analysis that ranged from 72°C to 95°C with 0.3°C increment per cycle to verify that a single product was amplified. Non-template controls were included to ensure that contamination with genomic DNA had not occurred.

TRXL3-3 was used as reference gene to normalise qPCR data (Medina-Lozano et al., 2023). Relative expression levels were obtained using the Pfaffl method (Pfaffl, 2001) with some modifications: arithmetic instead of geometric mean was calculated due to the presence of zero values in the raw data (either genes completely shut down as a consequence of the DI or the other way round, unexpressed genes in C conditions that were activated with the DI). This explains values different from 1 in C samples and why they have been represented separately from the DI data in qPCR results.

Student t-test was used to assess whether the differences between the means from the qPCR expression data of samples under C and DI conditions were statistically significant. Data transformations (1/(1+x)² or 1/√(x+1)) were applied when needed to achieve a normal distribution. Alternatively, Wilcoxon test was used with non-normally distributed data. Statistical analyses were conducted using the software JMP v5.1.2 for Windows (SAS Institute Inc. Cary, NC).

Detection of polymorphisms was carried out using the sequences of the 12 samples aligned to the lettuce reference genome (Reyes-Chin-Wo et al., 2017) and processed as explained in subsection 2.2. Firstly, variant calling was performed using FreeBayes package (Galaxy Version 1.3.6+galaxy0) (Tange, 2011; Garrison and Marth, 2012) from Galaxy platform. Then, VCFfilter (Galaxy Version 1.0.0_rc3+galaxy3) (Garrison, 2015) was used to remove polymorphisms with a total read depth at the locus < 10, QUAL < 20, and number of alternative alleles in called genotypes > 0. In addition, those polymorphic sites exhibiting different genotypes among the total number of samples within accessions and/or more than two different genotypes in comparison with the reference genome in more than 70% of the cases, were filtered out with Excel. Any possible ambiguous polymorphism was also eliminated.

The effect of each polymorphism was annotated and predicted using the SnpEff eff tool (Galaxy Version 4.3+T.galaxy1) and a snpEff database created using the SnpEff build tool (Galaxy Version 4.3+T.galaxy4) (Cingolani et al., 2012) from the annotation dataset and the FASTA file of L. sativa ‘Salinas’ v8 (Reyes-Chin-Wo et al., 2017).

To investigate the involvement of anthocyanins at molecular level in the response mechanism to drought stress of Lactuca spp., a transcriptomic analysis via RNA-seq was performed using plants belonging to the CWR L. homblei and to the red commercial lettuce variety ‘Romired’ coming from a previous experiment carried out in winter 2020/2021 (Medina-Lozano et al., 2024). Samples of both accessions showed an accumulation of anthocyanins under DI in comparison to C conditions, though the differences only resulted statistically significant in the case of the wild species L. homblei (Medina-Lozano et al., 2024). In particular, three different anthocyanins were identified: cyanidin 3-O-(6’-O-malonylglucoside) was the predominant one and was detected in both accessions and treatments; peonidin 3-O-glucoside appeared under both treatments in the commercial variety, but only under DI in the CWR; and cyanidin 3-(6’’-acetylglucoside), exclusively identified under DI conditions in the commercial variety.

After processing the data from L. homblei, the clean reads ranged from 40.76 to 51.75 Gb and the percentage of uniquely mapped sequences to the reference genome ranged from 32.39% to 37.40%. In the case of ‘Romired’, the clean reads ranged from 35.43 to 110.62 Gb and the uniquely mapped sequences from 81.04% to 84.81% ( Table 1 ). RNA-seq data from both accessions were aligned to the lettuce reference genome L. sativa ‘Salinas’ v8 (Reyes-Chin-Wo et al., 2017). However, L. homblei belongs to the tertiary lettuce gene pool (PGR (Plant Genetic Resources) Lettuce, https://www.pgrportal.nl/en/lettuce-geneticresources-portal.htm), so it is quite distant from L. sativa, what might explain its lower values in terms of uniquely mapped sequences.

A total of 6,179 DEGs were identified when L. homblei plants under C and DI treatments were compared (3,113 upregulated and 3,066 downregulated genes), whereas a total of 5,329 DEGs were obtained in ‘Romired’ plants for the same treatments (1,747 upregulated and 3,582 downregulated genes) ( Figure 1A ). A total of 2,272 DEGs were common to both accessions: 847 genes were upregulated, 1,347 downregulated, and 78 exhibited an opposite behaviour in the two accessions ( Figure 1B ). Attending to the differences, the CWR showed a total number of DEGs higher than the cultivated species (42.3% of exclusive DEGs in L. homblei vs. 33.1% in ‘Romired’). The same happened in the case of the upregulated genes, where the disparity was the largest, a 56.5% of the DEGs was exclusively upregulated in L. homblei, which was more than twice the upregulated DEGs only in ‘Romired’ (22.4%). In the case of the downregulated genes, we observed the opposite, the number was higher in the cultivated species than in the wild relative (42.2% vs. 32.4%, respectively) ( Figure 1B ).

To deeply explore the DEG functions in the drought response of Lactuca spp., analyses of GO enrichment were conducted using the GO annotations found in Reyes-Chin-Wo et al. (2017). The three main GO categories, biological processes, cellular components and molecular functions, were studied within the upregulated and downregulated genes ( Figure 2 ). Within the upregulated genes, the number of enriched GO terms in biological processes was higher in L. homblei than in ‘Romired’ ( Figure 2A ). In particular, an important number of L. homblei DEGs belonged to significantly enriched GOs that were involved in transmembrane transport and different metabolic processes, though the most significantly upregulated DEGs were those in enriched GOs related to gene expression regulation and response to abiotic stimulus, water included. The response to water resulted to be also among the enriched GO terms in ‘Romired’, but the number of genes, and especially the significance level, were lower than in the CWR L. homblei. In the case of cellular components ( Figure 2A ), a similar number of genes were part of the GO terms membranes and lipid storage bodies in both species, with a higher significance in ‘Romired’. In addition, the GO term cellular anatomical entity was enriched exclusively in the cultivated accession. Finally, for the molecular function category ( Figure 2A ), the most significantly enriched term was the transcription regulator activity in the CWR L. homblei, and the endopeptidase inhibitor activity in the commercial lettuce ‘Romired’. Interestingly, anthocyanins (among other flavins) could be involved in the oxidoreductase activity in which a flavin group acts as acceptor, being actually the only enriched GO common to both species.

Within the downregulated DEGs of biological processes ( Figure 2B ), genes involved in the carbohydrate metabolism were the most represented in L. homblei, while in ‘Romired’ were those implied in protein phosphorylation, that in fact, appeared only in this accession. However, several enriched processes, as well as their significance levels, were common or similar in both species, such as those related to cellular division and multiplication (DNA replication, nucleosome assembly, and microtubule-based processes). In the cellular component category ( Figure 2B ), the enriched GO terms found in ‘Romired’ appeared also enriched in L. homblei: nucleosome, microtubule associated complex, and external encapsulating structure, with similar significance and number of genes, except for the external encapsulating structure that resulted more significant in L. homblei. In fact, many more GO terms were enriched in L. homblei, with the nucleus and the protein-containing complex being the most represented ones. On the contrary, we found many more downregulated DEGs with enriched GO, as well as more terms and with a higher significance, in the cultivated (‘Romired’) than in the wild species (L. homblei) in the molecular function category ( Figure 2B ).

Some DEGs were assigned to more than one GO term, either because a gene can participate in different biological processes and molecular functions and be part of different cellular components or because GO is loosely hierarchical, with genes belonging to both ‘parent’ and ‘child’ terms. Thus, counting genes only once within each category, we obtained that the number of DEGs with enriched GO terms was very similar between L. homblei and ‘Romired’ within the upregulated genes (746 vs. 794, respectively), while it was lower in the CWR L. homblei than in the commercial variety ‘Romired’ within the downregulated ones (1,591 vs. 1978, respectively). Even so, the percentages of common DEGs in the two species was considerably lower in the case of upregulated genes than in downregulated: 3.23%, 18.25%, and 1.14% vs. 12.81%, 22.98%, and 23.00% in biological processes, cellular components and molecular functions, respectively.

Two heatmaps constructed using normalised expression data of DEGs in L. homblei and ‘Romired’ confirmed the effect of the drought stress treatment in Lactuca spp. plants in terms of gene regulation ( Figure 3 ). A hierarchical clustering conducted with all the DEGs allowed us to identify two separate groups in the two accessions, as expected, the upregulated and the downregulated ones. This clustering also divided clearly the two conditions (C and DI) in both species, what was even more evident for L. homblei ( Figure 3A ). In addition, both heatmaps showed again that the number of upregulated genes under water deficit was clearly higher in the CWR L. homblei ( Figure 3A ), whereas those downregulated were more numerous in the commercial variety ‘Romired’ ( Figure 3B ). These results show that the wild species was activating more mechanisms in response to drought stress. Figure 3 also shows two heatmaps constructed using the data of the anthocyanin content variation as a consequence of the drought stress for L. homblei and ‘Romired’ (Medina-Lozano et al., 2024), and the data of the two treatments themselves (C and DI). Content of all detected anthocyanins was higher under DI treatment than in C conditions in both Lactuca spp. Similar to what happened with the upregulated genes, the accumulation of anthocyanins in response to water stress was higher (and only significant) in L. homblei ( Figure 3A ). Actually, in L. homblei all the DI replicates showed a higher content than the C replicates. This was especially remarkable in the case of peonidin 3-O-glucoside, which was present under DI conditions and in the 3 biological replicates, but not under C conditions ( Figure 3A ). However, differences were not so clear (and not significant) between DI and C replicates in ‘Romired’, despite mean anthocyanin content being higher in DI than in C conditions, as commented above ( Figure 3B ). Even though, one of the anthocyanins was also identified exclusively under DI conditions in ‘Romired’, as observed in L. homblei, but in this case, it was cyanidin 3-(6’’-acetylglucoside) ( Figure 3B ).

Nevertheless, the molecular mechanisms underlying the anthocyanin accumulation as a consequence of water deficiency have been barely studied in lettuce, unlike in other crops like grapevine (Castellarin et al., 2007) or canola (Chen et al., 2022b). To gain a more comprehensive understanding of the process in Lactuca spp., we mapped the expression profiles of the DEGs identified in the drought experiment which participate in the biosynthesis pathway of the detected anthocyanins (the general phenylpropanoid pathway and the flavonoid pathway, this last one leading specifically to the anthocyanin biosynthesis) ( Figure 4 ). We found a total of 36 DEGs involved in the pathway in either L. homblei, ‘Romired’ or both. Different expression profiles were observed between both species. Our results confirmed that, in these routes, more DEGs were activated in the CWR L. homblei than in the commercial variety ‘Romired’, 17 vs. 7 upregulated genes, respectively, what was concordant with the higher accumulation of anthocyanins in the wild relative (Medina-Lozano et al., 2024). Not all the isoforms of the genes coding for the enzymes catalysing each step were upregulated under DI. The activation happened mainly in the first steps of the pathway, that is, in the early biosynthesis genes (EBGs), especially at the beginning of anthocyanin-specific route (flavonoid pathway). This becomes glaringly obvious in the first step which is catalysed by the chalcone synthase (CHS), whose gene isoforms are all strongly and significantly upregulated only in L. homblei (triggering of the anthocyanin synthesis in the wild species). This pattern is not so obvious in the preceding genes from the general phenylpropanoid pathway as they participate in the biosynthesis of many other compounds apart from anthocyanins. Furthermore, the late biosynthesis genes (LBGs) were mostly upregulated in the CWR L. homblei but not in all the isoforms as observed in CHS (EBG), except in the last step which leads to the synthesis of the specific major anthocyanin (cyanidin 3-O-(6’-O-malonylglucoside)) where most of the genes coding for the isoforms were significantly activated in L. homblei. The final steps to produce the two minor anthocyanins, cyanidin 3-(6’’-acetylglucoside) and peonidin 3-O-glucoside, are not clearly described in the literature. They might be catalysed by some acetyltransferases and O-methyltransferases, respectively, as suggested by Ino and Yamaguchi (1993) and Hugueney et al. (2009), respectively. It is possible that the genes coding for these enzymes were activated under drought as those anthocyanins were detected in ‘Romired’ and L. homblei, respectively, only under stress conditions though they have not been characterised in L. sativa yet.

Selection of anthocyanin-related genes potentially involved in the response to drought stress was based on different criteria. First, we searched for important changes in the expression levels. Second, we selected DEGs with high (positive and negative) and significant values of correlation with both anthocyanin content and drought stress treatment, obtained through a WGCNA. WGCNA allows to identify genes correlated with certain traits (anthocyanins and irrigation treatment in our case) to reveal putative genes with particular interest (Horvath and Dong, 2008). Lastly, we paid attention to gene function, so that DEGs were related to stress and/or anthocyanin content. Finally, 19 genes were selected for validation through real-time qPCR. Remarkably, genes meeting all these criteria resulted to be differentially expressed exclusively in L. homblei.

L. homblei |log₂(FC)| values between C and DI ranged from 4.03 to 6.23 ( Table 2 ). Both up- and downregulated genes were included in the selection. The higher accumulation of anthocyanins and/or the activation of stress response may result either from the upregulation of activators or from the downregulation of repressors. This was also observed in a previous study that characterised four genes related to anthocyanin content in lettuce (Su et al., 2020). From the WGCNA results, we obtained absolute correlation values with anthocyanins ranging from 0.81 to 0.93, and with treatment, from 0.86 to 0.99 ( Table 2 ). Both positive and negative correlations were also considered here. By contrast, in the case of ‘Romired’, we found 5 out of the 19 genes showing a significant correlation with the anthocyanin content, but none of them exhibited a significant change of expression level nor a significant correlation with the stress treatment (data not shown).

The expression data of the 19 selected genes obtained from the RNA-seq analysis were validated by real-time qPCR. In L. homblei, the 13 downregulated and the 6 upregulated genes according to the RNA-seq analysis showed concordant expression profiles with the qPCR results ( Figure 5 ). Significant, and very significant differences were observed between C and DI treatments for seven and two genes, respectively, according to qPCR data ( Figure 5 ). In the case of ‘Romired’, the selected genes did not show any differential expression in the RNA-seq analysis, as mentioned above, nor by qPCR. Even so, the expression of 15 out of the 19 genes followed the same profile using the two different techniques ( Supplementary Figure S1 ). Therefore, we were able to confirm the reliability of the results from the RNA-seq analysis.

The expression profiles of the selected DEGs obtained with both techniques (RNA-seq and qPCR) were also analysed by hierarchical clustering both using the mean values ( Figure 6 ) and all data ( Supplementary Figure S2 ). According to the RNA-seq data, two clearly differentiated expression patterns could be observed for L. homblei: the expression levels were noticeably lower in C than in DI conditions in the case of upregulated genes, and vice versa in the case of downregulated genes ( Figure 6A ; Supplementary Figure S2A ). In contrast, for ‘Romired’ the differences in expression were not clear in most of the selected genes ( Figure 6A ; Supplementary Figure S2A ), as was expected since they did not result to be differentially expressed in the RNA-seq analysis. Comparable patterns, especially for L. homblei samples, were observed when qPCR data were represented ( Figure 6B ; Supplementary Figure S2B ).

Turning the attention to gene function, it could be confirmed that those 19 DEGs with large changes in expression levels and high correlations with treatment and anthocyanins were indeed related to stress responses and/or to anthocyanin content ( Table 3 ). Specifically, most gene products of the selected DEGs have been described to participate in the response to one or more types of stresses. Several of the DEGs are involved in the response to biotic stresses, like resistance to bacteria (Lsat_1_v5_gn_1_21441 (Ramírez et al., 2013)), virus (Lsat_1_v5_gn_3_1101 (Tsitsekian et al., 2023)), or fungi (Lsat_1_v5_gn_3_20640 (Wu et al., 2022)), but most have been described to act in abiotic stress responses. In particular, genes related to water deficit and/or the stress-responsive hormone ABA (abscisic acid) stood out, such as Lsat_1_v5_gn_2_116640 (Clément et al., 2011), Lsat_1_v5_gn_2_43400 (Palusa et al., 2007), and Lsat_1_v5_gn_3_20640 (Tahmasebi et al., 2019). Furthermore, not only genes reported to be generally activated under water stress conditions were included in the selection, but also some described as negative regulators, which were in fact downregulated (inhibition of suppressors) in our samples subject to the drought treatment ( Tables 2 , 3 ), like Lsat_1_v5_gn_1_50480 (Lee et al., 2022), Lsat_1_v5_gn_2_90361 (Qu et al., 2016), and Lsat_1_v5_gn_6_67540 (Na and Metzger, 2020). Some of the selected DEGs have been related to other abiotic stresses such as salt (Lsat_1_v5_gn_1_127541 (Bhattarai et al., 2021) and Lsat_1_v5_gn_3_1101 (He et al., 2020)), heat (Lsat_1_v5_gn_2_47181 (Zhang et al., 2022)), and nutrient deficiency (Lsat_1_v5_gn_9_80621 (Zhou et al., 2021)), which makes sense as especially salt and heat stresses often occur simultaneously with drought.

Genes related to flavonoid (e.g., anthocyanins) accumulation were mostly upregulated in our samples ( Tables 2 and 3 ). Two of those genes encode putative transcription factors (TFs): a zinc finger protein (Lsat_1_v5_gn_5_26000 (Zhang et al., 2021)) and a NAC TF (Lsat_1_v5_gn_5_7401 (Wei et al., 2020)), both having been described to induce anthocyanin-related genes. Among the DEGs identified in this work, other two have been previously described to cause the increase of flavonoid or anthocyanin content under abiotic stress. Specifically, Lsat_1_v5_gn_2_15680 was found to play important roles in flavonoid accumulation under drought stress in tea (Li et al., 2020), and Lsat_1_v5_gn_8_165301 was involved in the increase of anthocyanin content when overexpressed in Arabidopsis thaliana under ABA-mediated abiotic stress responses (Wang et al., 2018). Another selected gene (Lsat_1_v5_gn_8_157561) is not well characterised but could possibly be related to the anthocyanin content as it contains a MYB domain and, in plants, MYB TFs have been described as one of the major transcriptional regulators of anthocyanin pathway, both activators and repressors (Yan et al., 2021; Cao et al., 2024). This specific gene was downregulated in our samples when exposed to water stress, so it might be a transcriptional repressor of anthocyanin biosynthesis. Lsat_1_v5_gn_7_92980 encodes a ribonuclease III-like protein 3 (RTL3) that cleaves doble-stranded RNA (Comella et al., 2008) and might be also related to anthocyanin content. Proteins with ribonuclease III domains have been described to participate in the regulation of seed coat in soybean and fruit colour in peach through the production of siRNAs (small interfering RNAs) and the increase of transcription levels of genes implied in anthocyanin regulation (Zhu et al., 2012; Jia et al., 2020). Finally, two genes that encode proteins involved in DNA damage repair were also selected (Lsat_1_v5_gn_1_109200 (Fujimori et al., 2014), and Lsat_1_v5_gn_5_10141 (Dangel et al., 2014)). This is not surprising as the generation of ROS is a potential cause of DNA damage under drought stress and a fine-tuned regulation of DNA repair is required to tolerate it (Shim et al., 2018).

In silico search and prediction of polymorphisms were performed to get an overview of the variation in the sequences of the total number of DEGs detected and to explore more deeply the structural variation of the set of 19 selected DEGs.

A total of 235,600 polymorphisms were found in the whole set of DEGs (9,236) in L. homblei and ‘Romired’ compared to the reference genome (both shared and species-exclusive). Most polymorphisms were detected in L. homblei, as expected since it is a wild species that is very distant from the cultivated L. sativa used as reference. The predominant types of polymorphisms were SNPs (Single Nucleotide Polymorphisms) (89.22%), followed by MNPs (Multiple Nucleotide Polymorphisms) (9.97%), and, in a much smaller extent, by indels (insertions-deletions) (0.81%). The most abundant polymorphism effects were synonymous (63.20%) and missense (27.96%). We also identified intron (6.14%) and splice region (2.23%) variants, as well as others that were present in less than 0.1%, so they are not detailed here.

In the subset of 19 DEGs selected, a total of 404 polymorphisms with 408 predicted effects were identified in L. homblei ( Table 4 ), in contrast to the 11 polymorphisms with 12 predicted effects found in those same 19 genes non-differentially expressed in ‘Romired’ ( Supplementary Table S2 ). Considering only the 19 DEGs in L. homblei, the proportions of both polymorphism types and effects were almost the same than those found in the whole set of DEGs, 87.87% of polymorphisms were SNPs, 11.14% were MNPs, and 0.99% were indels. Once again, we found that the predominant predicted effect was synonymous (70.10%), followed by missense (27.44%) type ( Table 4 ). We also identified, though in a reduced number of genes, effects in splice regions (0.98%) and introns (0.49%), disruptive and conservative in-frame deletions (0.49 and 0.25%, respectively), and a frameshift variant (0.25%). The impact of the polymorphisms was frequently low, which makes sense considering that most of them were predicted to have a synonymous effect. However, a polymorphism with high impact was detected. It was a 2-bp insertion that theoretically causes a frameshift mutation in the Lsat_1_v5_gn_3_1101 gene of the wild species which is responsible for the appearance of a premature stop codon. A conservative in-frame deletion was also found in this same gene. According to our results, this gene was downregulated in L. homblei whereas in ‘Romired’ was not differentially expressed, in which showed low expression levels in both C and DI conditions. This gene codes for a zinc-finger protein and appears in the literature as a putative TF that intervenes in the response to abiotic stress (He et al., 2020). A possible effect of one or both polymorphisms might be that the truncated protein acts as a repressor in the wild species under C conditions but stops inhibiting its target(s) as a consequence of its own downregulation under water stress. Other polymorphisms with possible important effects were the putative disruptive in-frame deletions found in Lsat_1_v5_gn_2_116640 and in Lsat_1_v5_gn_8_165301 genes, whose predicted impact was moderate. Lsat_1_v5_gn_2_116640 encodes a 70-kDa heat shock cognate protein. Heat shock proteins (HSPs) were initially described in relation to heat tolerance (Ritossa, 1962), although nowadays they are well known to be expressed in response to a great diversity of environmental stressors besides heat (reviewed in Ul Haq et al. (2019)). According to the RNA-seq analysis, this gene (Lsat_1_v5_gn_2_116640) showed a considerable increase in expression in L. homblei in response to drought, whereas in ‘Romired’ there was no significant change, with the values under both C and DI being similar to those in C plants of L. homblei. Therefore, the disruptive in-frame deletion in this gene could be inducing the activation of this HSP when L. homblei plants are subject to drought stress. Lsat_1_v5_gn_8_165301 encodes a Phospholipase A1 phospholipid-inositol phosphatase 2 (PLIP2) that has been described to be involved in the accumulation of anthocyanins under ABA-mediated abiotic stress responses (Wang et al., 2018). This gene also exhibited a highly significant upregulation in L. homblei and no change of expression in ‘Romired’. In this case, its expression levels in ‘Romired’ under C and DI conditions were similar to those of L. homblei under DI. Thus, the disruptive in-frame deletion found in L. homblei sequence might be causing the gene to be activated only under stress in the wild plants.

The number of exclusive up- and downregulated genes was more than twice in L. homblei and 1.3 times higher in ‘Romired’, respectively. In general, genes related to regulation within biological process (e.g., response to water) and molecular function categories were more abundant and more intensively upregulated in the wild species whereas genes responsible for cellular components were more commonly and significantly upregulated in the cultivated species. In the case of the downregulated DEGs, the most represented terms in both species were those related to catalytic activities. In general terms, basal and growth-related processes were deactivated in both species, which probably contributes to redirect resources to guarantee plant survival.

Interestingly, activation of responses seemed to be species specific, and it looks like the CWR was triggering more mechanisms of response to drought stress as the number of upregulated genes was clearly larger in L. homblei under DI. In contrast, the higher number of downregulated DEGs common to both species, many of them implied in basal processes, could be due to the deactivation of basal metabolism processes to designate more resources to water deficit tolerance, previously described in different plant species subject to water stress (Shao et al., 2009). The results from the GO enrichment analysis are in agreement with those found in other studies that assessed different stresses in lettuce, in which response to stimulus, biological regulation, metabolic processes, binding and catalytic activities, as well as membrane components, were the most represented terms (Wang et al., 2017; Zhou et al., 2023). Interestingly, in a transcriptomic analysis carried out to identify genes involved in lettuce anthocyanin accumulation, the most represented GO terms were the same (Zhang et al., 2016).

Both anthocyanin contents and the number of upregulated genes were clearly larger in L. homblei under DI. Therefore, results point to a relationship between gene expression profiles (for some DEGs) and changes in the accumulation of these antioxidant compounds. In addition, not only the most abundant anthocyanins showing the biggest change in quantity in response to water stress but also the minor anthocyanins only identified under DI (peonidin 3-O-glucoside and cyanidin 3-(6’’-acetylglucoside) in L. homblei and ‘Romired’, respectively) could play a role in the response to drought. That is more plausible in the case of L. homblei where the differences in anthocyanin content between C and DI conditions were significant (Medina-Lozano et al., 2024).

Interestingly, to activate the anthocyanin biosynthesis route in L. homblei seems to be enough to upregulate the isoforms of the gene controlling the first step of the specific pathway branch (i.e., CHS), even when the preceding genes from the general phenylpropanoid pathway (e.g., 4-coumarate-CoA ligase (4CL)) could be downregulated or not significantly differentially expressed, as they are involved in the biosynthesis of many other phenylpropanoids apart from anthocyanins. Something similar has been shown in a previous study on the expression of those genes and the anthocyanin content of poplar leaves (Tian et al., 2021).

All these differences between both Lactuca spp. might reflect a lager plasticity of the wild species to adapt to environmental changes. The great genetic diversity of wild species allows them to counteract the effects of different stresses more effectively (Jordanovska et al., 2020), whereas the cultivated species could have lost these mechanisms through domestication. In fact, the common DEGs to both accessions which show an opposite sense in the change of expression could consist of genes that have acquired a different mode of action as L. homblei belongs to the lettuce tertiary gene pool, the most genetically distant from L. sativa. Alternatively, they could be artefacts, either methodological (e.g., library preparation) or statistical or even both.

The fact that the genes with the strongest change of expression and correlation with anthocyanin content and drought were differentially expressed only in L. homblei might reveal, once again, the existence of tolerance mechanisms in the wild species that are not present in the cultivated one. This is in agreement with the wild species showing the highest increase (and the only resulting statistically significant) of anthocyanins in a previous study on drought stress with the same accessions, among others (Medina-Lozano et al., 2024).

The reliability of the results from the RNA-seq analysis was confirmed as the candidate DEGs were validated by real-time qPCR, being all differentially expressed only in the wild species, in which the expression profiles obtained with the two techniques coincided. Besides, most gene products of the selected DEGs have been described to participate in the response to one or more types of stresses which makes sense as some stresses often occur simultaneously.

The fact that the genes with a high change in the level of expression identified in this study resulted to be related to both the stress response and the anthocyanin content could indicate that these compounds are playing an important role for plants to cope with the drought conditions, as was also proposed before in purple-stem Brassica napus L. (Chen et al., 2022b).

Talking about the polymorphisms found in all DEGs and in the candidate genes, our results are in agreement with studies carried out in other crops that also used transcriptomic data, in which the number of SNPs was also much higher than the number of indels, and the most abundant polymorphism effects were synonymous and missense variants too (Iquebal et al., 2017; Muñoz-Espinoza et al., 2020). The impact of the polymorphisms identified in the candidate DEGs frequently resulted low as most of them were predicted to have a synonymous effect. However, a few of them showed a high or moderate predicted impact what could be a reflection of the drastic changes in the gene expression profiles (either activation or inhibition) in the wild species when subject to drought stress. Among the rest of the polymorphisms found, the missense variants could have an impact on the function of the resultant protein, though this one has been predicted to be moderate.

In silico tools are truly useful for obtaining information of functional and structural variants on the transcriptome and their possible correlation with phenotypic changes (Yazar and Özbek, 2021). However, further experimental approaches like functional analyses are essential to verify in the future the polymorphism effect found in putative candidate genes involved in the anthocyanin accumulation and the response to drought stress.