A systems view on DNA damage response kinetics in Tetrahymena

Known DNA damage repair factors are differentially regulated in response to genotoxic stressors

To study DNA damage response kinetics, we exposed Tetrahymena thermophila to six well-established DNA damaging agents: 254-nm ultraviolet light (UV), cis-diamineplatinum (II) dichloride (CPT), hydrogen peroxide (HP), methyl methanesulfonate (MMS), ionizing radiation (IR), and hydroxyurea (HU). Treatment conditions were determined either by establishment of EC50 or based on previous DNA damage studies of Tetrahymena (Supplemental Table S1; Campbell and Romero 1998; Loidl and Mochizuki 2009; Sandoval et al. 2015). To obtain transcriptome and proteome expression profiles, we harvested samples at 0, 1, 2, 3, 4, 6, and 8 h after treatment (0 h–8 h) in quadruplicate and performed mRNA sequencing (RNA-seq) and high-resolution quantitative mass spectrometry (MS) measurement (Fig. 1A). Transcriptomes and proteomes were measured in sets of three, with two treatments paired with a nontreated control. With this, we could calculate changes in transcript expression and protein intensity in each drug treatment condition over a matched control, thereby correcting for any potential batch effect (log₂ fold change values, Supplemental Fig. S1). After stringent filtering across treatments and time points, 20,443 transcripts and 6551 protein groups were quantified, with 99.6% of proteins associated with transcriptome data (Fig. 1B).

Figure 1.

Screen to explore kinetics of DNA damage response in Tetrahymena. (A) Schematic of the screen workflow. Cells were treated with a mutagenic agent, and samples were harvested incrementally over 8 h. At each time point, samples were collected for processing for both RNA sequencing and quantitative mass spectrometry. (B) Venn diagram depicting the overlap between the identified transcripts and proteins. (C) Bar graph of enriched (full color) and detected (transparent color) DNA damage repair factors, where enriched proteins have L2FC > 1 or transcripts have L2FC > 2. Detected proteins or transcripts are present and found in the respective data sets, but not enriched. (D) Heat map of hierarchical clustering of DNA repair genes of interest. Light blue box highlighting the group of core damage response transcripts. NER: nucleotide excision repair, ICL: inter crosslink repair, BER: base excision repair, DSBR: double-strand break repair, MMR: mismatch repair, DDR: DNA damage response, UV: ultraviolet light, CPT: cisplatin, HP: hydrogen peroxide, MMS: methyl methanesulfonate, IR: ionizing radiation, HU: hydroxyurea.

Open in new tabDownload PowerPointLink to figure

To confirm that each treatment induced the anticipated DNA damage, we examined the differential upregulation of primary DNA repair factors of NER (Tatum and Li 2011; Schärer 2013), BER (Kelley et al. 2003; Beard et al. 2019), MMR (Bowen et al. 2013; Kunkel and Erie 2015; Chakraborty and Alani 2016), DSBR (Li and Heyer 2008; Mathiasen and Lisby 2014; Pannunzio et al. 2018; Scully et al. 2019), ICL (Lehoczký et al. 2007), and general DDR (Ciccia and Elledge 2010; Pizzul et al. 2022). Our list of these factors, containing 173 DNA repair genes, is not intended to be comprehensive but highlights established key genes in the different DNA repair pathways as well as non-pathway-specific overall DNA damage response (Supplemental Table S2). Across all treatments, between 81.4% and 100% of detected pathway-associated transcripts (fold change (FC) > 2) and between 82.4% and 96% of detected pathway-associated proteins were upregulated (FC > 1). In the transcriptome and proteome data set, 81.4% and 88.9% of the detected global, non-pathway-specific, DNA damage response transcripts and proteins were upregulated, respectively. For each treatment, we examined the amount of upregulated transcripts and proteins with its commonly described DNA repair pathway. Across all treatments, between 68.4% and 96.6% of the transcripts and between 33.3% and 92% of the proteins associated with the respective pathway were upregulated (Fig. 1C).

Selected members of the primary DNA repair pathways were further examined for each treatment over time (Fig. 1D). Hierarchical clustering revealed a group of 10 transcripts (MSH6L3, RAD51, RAD4, SNML1, RAD5, RLP1, RAD5L4, RFA1, TTHERM_00316410 [Rad3 homolog], and TTHERM_00391570 [Rev3 homolog]) with a minimum 2.3 log₂ fold change across all treatments (blue box, Fig. 1D). This indicated an unexpected core DNA damage response that involves parts of MMR, DSBR, and NER, as well as expected general responders, regardless of damage origin.

Genotoxic stressors induce core and specific dynamic gene expression responses

We examined expression dynamics of both the transcriptome and proteome over time by calculating the Gini coefficient, a measurement of inequality, for every quantified transcript and protein (Fig. 2).

Figure 2.

Mutagenic treatments cause global dynamic response. (A) Volcano plots plotting dynamic transcripts for each treatment. The x-axis contains the maximal positive or negative fold change during the duration of the time course, and the y-axis contains a Gini score evaluating the dynamicity throughout the entire time course. (B) Example line plots of the expression profiles dynamic (MSH6L3 and RAD51) and stable transcripts (TTLL6B). (C) Upset plots of overlapping dynamic transcripts between treatments. (D) Volcano plots plotting dynamic proteins for each treatment. The x-axis contains the maximal positive or negative fold change during the duration of the time course, and the y-axis contains a Gini score evaluating the dynamicity throughout the entire time course. (E) Line plots of the expression profiles dynamic (MSH6L3 and RAD51) and stable proteins (BTU1). (F) Upset plots of overlapping dynamic proteins between treatments.

Open in new tabDownload PowerPointLink to figure

To categorize dynamic and stable transcripts, we applied a filter at the 60^th quantile of the calculated Gini coefficients (Gini score > 0.042) (Fig. 2A). In addition to this dynamicity filter, to select transcripts generally differentially regulated in response to the DNA damaging agents, we filtered based on log₂ fold change (L2FC), requiring the expression change over the time course to reach a magnitude of at least 1 L2FC, and a significance of adjusted P-value < 0.05 (FDR). Surpassing these thresholds, of the overall 20,443 detected transcripts, we classified 8815 as dynamic. Among them were the previously characterized DSBR gene RAD51 and the MSH6 homolog, MSH6L3 (Campbell and Romero 1998; Marsh et al. 2000, 2001), whereas, for example, TTLL6B, a tubulin-related ligase, was found to be a stable transcript (Fig. 2B).

There were 57 overlaps of shared dynamic transcripts among treatments, of which 42 were significantly more than expected (P-value ≤ 0.04, Fisher's exact test). The seven groups which contained overlaps of dynamic transcripts in five or more treatments were significant (P-value < 0.001) (Fig. 2C; Supplemental Tables S3, S4). Of the 78 transcripts in these overlaps, 39 had no S. cerevisiae homolog. This indicates both a strongly conserved and unique DNA damage response in Tetrahymena. Among these transcripts with no yeast homolog were three PARP and PARP-correlated genes, PCP3 (PARP12, TTHERM_00467770), PARP3 (TTHERM_00030430), and PCP5 (PZN1, TTHERM_00773650). Of the 15 PARP and PARP-correlated genes, 12 members were dynamic in one or more treatments. Whereas PCP3, PARP3, and PCP5 act as core responders, there are also damage-specific PARP responses. For example, PARP1, PARP7, and PCP1 have a dynamic expression profile only in response to UV. Previously, biologically specific PARP responses have also been observed in the interface between RNA interference pathways and DNA repair (Lee et al. 2021), indicating there is a conserved and unique DNA damage response both globally and within protein families.

We performed the same analysis with the 6551 proteins quantified across these six treatments. Here, the Gini coefficient threshold for dynamicity was also set at the 60^th quantile (Gini score > 0.021 with a minimum L2FC > |1|, and P-value < 0.05 [Welch t-test]) (Fig. 2D). We found a total of 2582 proteins to be dynamically regulated. RAD51 and MSH6L3 were dynamic in each treatment (Fig. 2E). Of the 57 overlapping groups of dynamically upregulated proteins, there was significantly more overlap than expected in 38 groups (P ≤ 0.05, Fisher's exact test) (Fig. 2F; Supplemental Tables S5, S6). There were 33 dynamic proteins that were present in overlaps containing five or more treatments. Of these 33 core responders, 15 of these proteins have no homolog in S. cerevisiae. The dynamic proteins across all treatments had overrepresentation of the GO terms “DNA repair” and “cellular response to damage.” Thus, as with the dynamic transcripts, both a global and damage-specific protein response was observed.

Early transcriptional activity is critical to the DNA damage response

To cluster the 8815 dynamic and differentially regulated transcripts, we used self-organizing maps (SOMs), an unsupervised machine learning approach. The transcripts detected across all six treatments fell into 15 distinct clusters (T1–T15); 689 transcripts could not be assigned to any cluster (Fig. 3; Supplemental Fig. S2A). Each cluster exhibits unique temporal and intensity of response. Clusters T1–T6 have the highest fold change at either 0 h or 1 h, followed predominantly by downregulation throughout 2 h to 8 h (Fig. 3A). These clusters represent early responders in all treatments that are progressively downregulated as repair progresses. Conversely, clusters T8, T10, T12, T13, and T14 show the most pronounced differential regulation CPT-treated cells, underscoring a treatment-specific response. Next, each transcript was mapped to its respective homolog(s) in S. cerevisiae, and functional enrichment analysis was performed using Gene Ontology (GO) of the S. cerevisiae homologs in each cluster (Fig. 3B,C).

Figure 3.

Transcript clusters reveal dynamic DNA damage response. (A) Average expression profiles for clusters of dynamic transcripts. Using an unsupervised machine learning technique, we clustered dynamic transcripts based on their shared expression profiles into 15 different clusters. These line graphs are the average expression profiles for each treatment. (B) Heat map of functional enrichment analysis using KEGG. Each row contains an over-represented KEGG term with a gradient representing the adjusted P-value. (C) Heat map of functional enrichment analysis using GO analysis. Each gene was mapped to its respective homolog(s) in S. cerevisiae. Each row contains an overrepresented GO term with a gradient representing the adjusted P-value.

Open in new tabDownload PowerPointLink to figure

Cluster T10 uniquely showed overrepresentation of genes related to “DNA repair,” “cellular response to DNA damage stimulus,” “DNA replication,” and “DNA metabolic process.” Within the average expression profile for CPT-treated cells in this cluster, there is an immediate strong and continual upregulation until 8 h. This response could be due to the known long half-life of CPT (Evans et al. 1982). At 0 h, the average expression profile of IR- and UV-treated cells shows strong upregulation, followed by gradual downregulation throughout the remaining time points. This also reflects these specific treatments, as both UV and IR treatments had only one initial application and were not sustained during the time course, whereas the average fold change of HP-, MMS-, and HU-treated cells decreases over time. All clustered DDR genes are included in Supplemental Figure S3. Within the 15 clusters, nine clusters are enriched for “protein phosphorylation” (T1, T2, T4, T5, T8, T10, T12–T14). The phosphorylation of DNA repair proteins has been well established to regulate all major repair pathways (Lanz et al. 2019; Huang and Zhou 2020). Overall, our data suggest that protein phosphorylation is critical to immediate and sustained DNA damage response as a whole in Tetrahymena.

Further general responses to DNA damage were found also in the other clusters. All histones (T7), 20S ribosomal proteins (T2), and dense core granules (T3) clustered together and showed primarily similar regulation trends for each treatment (Supplemental Fig. S4A–C). The degree and time point of decline depended on the actual treatment, but each gene family or complex subunits behaved similarly within each treatment. In contrast, the three families of previously studied chromatin remodelers in Tetrahymena—the Poly-(ADP-ribose) polymerases (PARPs)/PARP-associated proteins, histone acetyl transferases (HATs), and histone deacetylases (HDACs)—are differentially regulated and span across multiple clusters (Supplemental Fig. S4D,E; Slade et al. 2011; Chalker et al. 2013; Ashraf et al. 2019; Saettone et al. 2019; Wahab et al. 2020). The PARP and PARP-correlated proteins mediate DNA repair by chromatin modifications via ADP-ribosylation, as well as direct binding, modification, and recruitment of DNA repair proteins (Sousa et al. 2012; Morales et al. 2014). PARP7, PARP8, and PCP3 (T2), PCP1 (T7), PARP6 (T9), and PARP2 and PARP5 (T10) all showed unique responses to DNA damage (Supplemental Fig. S4D). The histone acetylases and deacetylases are critical to changing chromatin architecture to facilitate DNA repair (Wahab et al. 2020). The histone acetylases (HATs) HAT1 (T10) and MYST2 (T2) and histone deacetylases THD4, THD17, and THD18a (T2, T3, T6, respectively) also showed greatly differential regulation (Supplemental Fig. S4E). This indicates that each of the PARPs, HATs, and HDACs in Tetrahymena has a particular role. Together, the transcriptome data show specific transcriptional regulation kinetics of the DNA damage response dependent on the genotoxic stressor.

Protein expression over time reveals specific trends involved in DDR

As for the transcriptome, we used self-organizing maps clustering 2582 dynamically expressed proteins into seven distinct clusters (P1–P7); 202 proteins could not be assigned (Fig. 4A; Supplemental Fig. S2B). Each protein included in these seven clusters was mapped to its respective homolog(s) in S. cerevisiae, and functional enrichment analysis was performed using GO (Fig. 4B,C). Cluster P6 uniquely showed overrepresentation of genes related to “DNA repair,” “cellular response to DNA damage stimulus,” “DNA replication,” and “DNA metabolic process.” In P6, there were 15 known DNA damage factors, which included ATR1, RAD53/Chk1, TKU80, RAD3, DNA2, and three members of the RFC complex. ATR1 and TKU80 are critical in both DNA damage response and conjugation in Tetrahymena (Loidl and Mochizuki 2009; Lin et al. 2012). Additionally, in this cluster, two MMR proteins were identified, MSH3L6 and TMLH1. Together with MSH2, MSH3 is part of the MMR MutSβ complex, which repairs larger insertions and deletions. Alternatively, if MSH6 and MSH2 form an MMR recognition heterodimer MutSɑ, one to two base pair mismatches and indels are repaired. MSH6 is found in P3, indicating a differential regulation of these portions of the MMR recognition complex. TMHL1 and PMS2 interact with the recognition complex to initiate cleavage events. Our clustered data suggest that TMLH1 expression profiles are more comparable to MSH3 rather than MSH6. All clustered DDR proteins are included in Supplemental Figure S5.

Figure 4.

Protein clusters reveal a specific and dynamic DNA damage response. (A) Average expression profiles for clusters of dynamic proteins. Using an unsupervised machine learning technique, we clustered dynamic transcripts based on their shared expression profiles into seven different clusters. These line graphs are the average expression profiles for each treatment. (B) Heat map of functional enrichment analysis using KEGG. Each row contains an overrepresented KEGG term with a gradient representing the adjusted P-value. (C) Heat map of functional enrichment analysis using GO analysis. Each gene was mapped to its respective homolog(s) in S. cerevisiae. Each row contains an overrepresented GO term with a gradient representing the adjusted P-value.

Open in new tabDownload PowerPointLink to figure

Similar to transcriptional regulation, the chromatin remodelers within the PARP family were differentially regulated across clusters. However, this family of genes is not just being regulated at the transcriptional level. For example, PARP7, PARP8, and PCP3 all fell within T2, whereas now they are found within P6, P7, and P4, respectively (Supplemental Fig. S6A).

Other gene families not directly involved in DNA repair exhibited similar time-dependent transcriptional regulation but grouped into distinct clusters, displaying divergent expression patterns at the protein level. For example, all dynamic dense core granules clustered together in T3 but were part of three different protein clusters (P1 and P3, P7, Supplemental Fig. S6B). All histones were transcriptionally regulated in a highly similar manner (T7), whereas the histone proteins had different protein expression profiles (P6 and P7, Supplemental Fig. S6C). Also, the DNA repair-related PARP family was differentially regulated. Whereas the dynamic transcripts also showed a differential regulation across all clusters, there was clustering of some of the PARP transcripts expression profiles in the same SOM. However, there were some instances of similar protein and transcriptional regulation, such as for the 20S proteasome (T2 and P7, Supplemental Fig. S6D) and transcription related factors (T6 and P7, Supplemental Fig. S6E). Generally, the differences between transcriptome and protein expression profiles unveil additional regulation at the protein level.

Cell culture

The Tetrahymena thermophila wild-type strain SB210 (Tetrahymena Stock Center) was used throughout the study. Cultures were grown in a medium of 2% proteose peptone (BD Biosciences), 0.2% yeast extract (BD Biosciences), 12 µM FeCl₃, and 1× penicillin/streptomycin/amphotericin B (Hyclone) at 30°C at 100–150 rotations per minute.

Collection of Tetrahymena for mass spectrometry and RNA sequencing

Tetrahymena were grown to a concentration between 1.5 × 10⁵ and 3 × 10⁵ cells/ml in 500-mL cultures. Samples were treated with six different conditions, grouping two treatments with one nontreated group (grouped as MMS, HP; CP, UV; and HU, IR). Details of treatments are described in Supplemental Table S1. Cells were harvested samples at 0, 1, 2, 3, 4, 6, and 8 h after the initial treatment. To collect samples for later quantitative mass spectrometry, 5 × 10⁴ cells were centrifuged at 9400g for 5 min. The supernatant was removed and cells were washed with 1 mL 10 mM Tris-HCl [pH 7.5], and centrifuged at 9400g for 5 min. The supernatant was discarded leaving a total ∼15 µL of cells and Tris, and 5 µL of 4× LDS (Thermo Fisher Scientific) and 2 µL of 1 m DTT (Sigma-Aldrich) were added. Samples were heated to 90°C for 10 min. Samples were stored at −20°C until mass spectrometry sample preparation. To collect samples for later RNA sequencing (RNA-seq), 5 mL of cells were collected and centrifuged at 1400g for 3 min. The supernatant was decanted and cells were washed with 5 mL 10 mM Tris-HCl [pH 7.5]. Cells were centrifuged at 1400g for 3 min and the supernatant was removed. The cell pellet was resuspended in 600 µL Buffer RLT (Qiagen, RNeasy mini kit), flash frozen in liquid nitrogen, and stored at −80°C until RNA sequencing sample preparation.

Mass spectrometry sample preparation

An LDS sample was loaded on a 4%–12% NuPage NOVEX Bis-Tris gel (Thermo Fisher Scientific) and run for 10 min at 180 V, in 1× MES buffer (Thermo Fisher Scientific). Samples were processed as previously described (Scherer et al. 2020). In short, the gel was stained and fixed with Coomassie Brilliant Blue G250 (Sigma-Aldrich), and initial destaining of the gels was done overnight with water. Gel pieces were cut, further destained with 50% EtOH/50 mM ammonium bicarbonate (ABC), dehydrated with acetonitrile (VWR), reduced with 10 mM DTT (Sigma-Aldrich), alkylated using iodoacetamide (Sigma-Aldrich), subsequently again dehydrated with acetonitrile (VWR), and digested with 1 µg of MS-grade trypsin (Sigma-Aldrich) at 37°C overnight. The peptides were eluted from the gel pieces and loaded onto activated C18 material (Empore) StageTips (Rappsilber et al. 2007) and stored at 4°C until elution and measurement.

Mass spectrometry measurement

Peptides were eluted from the StageTips using 80% acetonitrile/0.1% formic acid and concentrated prior to loading on an Easy-nLC-1200 system coupled to an Orbitrap Exploris 480 mass spectrometer (Thermo Fisher Scientific). The peptides were loaded on a 50-cm column (75 µm inner diameter, New Objective) in-house packed with ReproSil-Pur 120 C18-AQ (Dr. Maisch GmbH). We used a 103-min gradient from 3% to 40% acetonitrile with 0.1% formic acid at a flow of 250 nL/min. The mass spectrometer was operated in positive ion mode with a top 20 MS/MS data-dependent acquisition strategy of one MS full scan (scan range 300–1650 m/z; 60,000 resolution; normalized AGC target 300%; max IT 28 ms) and up to 20 MS/MS scans (15,000 resolution; AGC target 100%, max IT 40 ms; isolation window 1.4 m/z) with peptide match preferred using HCD fragmentation.

Mass spectrometry data analysis

Raw files were analyzed using MaxQuant (version 1.6.10.43). As a search space, the T. thermophila protein database was used (June 2014, TGD). Oxidation and acetylation were set as variable modifications, carbamidomethylation as a fixed modification. Fast LFQ was used to calculate and normalize intensities. The minimum ratio count used was 2. Match between runs was used to match within each time point per treatment and to the time points right before and after, with a match time window of 0.7 min, match ion mobility window of 0.05, an alignment time window of 20 min, and alignment ion mobility of 1. Matching of unidentified features was deactivated. For protein quantification, label minimum ratio count 2 and unique + razor peptides were used.

RNA sample preparation and sequencing

Previously obtained samples were thawed on ice. RNA isolation was performed with an RNeasy mini kit (Qiagen) in accordance with the manufacturer's instructions with the addition of the optional DNase I on column digestion. This digestion was carried out with 3 units of DNase I (Qiagen) per sample, and samples were digested for 15 min at room temperature on the column. NGS library prep was performed with Lexogen`s QuantSeq 3´mRNA-Seq Library Prep kit FWD following Lexogen`s standard protocol (015UG009V0252). Libraries were prepared with a starting amount of 300 ng and amplified in 14 PCR cycles. Libraries were profiled in a High Sensitivity DNA on a 2100 Bioanalyzer (Agilent Technologies) and quantified using the Qubit dsDNA HS Assay kit, in a Qubit 2.0 Fluorometer (Life Technologies). All libraries from the two treatments and coordinating nontreatment were pooled together in equimolar ratio and sequenced on 1 NextSeq 500 high output flow cell, SR for 1 × 84 cycles plus seven cycles for the index read.

RNA-seq analysis

All demultiplexed, raw sequencing files of each treatment set were analyzed together. Initial analysis was done through a modified version of the NGS pipeline by the bioinformatics core facility of the IMB (available at GitLab, https://gitlab.rlp.net/imbforge/NGSpipe2go). For reference, “subread2rnatypes,” “genebodyCov2,” “rMATS,” and the GO enrichment analysis were removed from the pipeline. In short, the library quality was assessed with FastQC before alignment against the T. thermophila genome assembly SB210 and a custom built GTF file, which included gene annotations from T. thermophila (TGD, T_thermophila_June2014.gff3). Alignment was performed with STAR aligner version 2.7.3a (Dobin et al. 2013). Reads mapping to annotated features in the custom GTF file were counted with featureCounts (Liao et al. 2014). Initial CPM counts were calculated with DESeq2 (Love et al. 2014) in R (R Core Team 2022).

Further bioinformatic analysis

All further analysis was done with scripts developed in R (R Core Team 2022), incorporating ggplot2 for visualization (Wickham 2016) among other packages.

For proteome data, contaminants, reverse database hits, protein groups only identified by site, and protein groups with less than two peptides (at least one of them classified as unique) were removed. Additionally, only protein groups present in at least two out of four technical replicates were kept. Missing values were imputed by shifting a compressed beta distribution obtained from the LFQ intensity values to the limit of quantitation (between 0.2 and 2.5 percentile of the measured intensity distribution per sample). LFQ intensities were log₂ transformed, after which fold changes for individual comparisons of time points or strains could be calculated per protein; a Welch t-test was used to calculate P-values. The general protein enrichment threshold was set to a P-value lower than 0.05 and an absolute fold change higher than 1. All calculated values can be found in Supplemental Table S7.

For transcriptome data, transcripts which did not have any CPM value across the time points and treatments above the 25^th quantile of all CPM values (CPM < 1.673028) were removed. All CPM values were log₂ transformed. Differential regulation thresholds were set at L2FC > 1 or <−1, and adjusted P-value (FDR) < 0.05. All calculated values can be found in Supplemental Table S8.

Dynamicity of transcripts or proteins was calculated using the Gini ratio, as described before (Damgaard and Weiner 2000; Casas-Vila et al. 2017). Statistical testing of overlaps of dynamic genes was done with the R package SuperExactTest (Wang et al. 2015). Functional enrichment analysis was performed using Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al. 2023), Gene Ontology (The Gene Ontology Consortium 2021), and the ClusterProfiler R package (Yu et al. 2012; Wu et al. 2021) for statistical analysis. Terms for groups of enriched proteins were assessed for overrepresentation with a Fisher's exact test, against all terms found in our complete data set as background. The enrichment threshold was set to an adjusted (FDR) P-value < 0.05. Self-organizing map clustering was done with the help of the Kohonen package in R.