Marc Le Vée

Amelie Moreau

Elodie Jouan

Claire Denizot

Yannick Parmentier

Irset Olivier Fardel
email: olivier.fardel@univ-rennes1.fr

Amélie Moreau

Inhibition of canalicular and sinusoidal taurocholate efflux by cholestatic drugs in human hepatoma HepaRG cells

Keywords: Bile salt export pump, cholestatic drug, efflux, HepaRG cells, taurocholate

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Alteration of hepatic bile acid transport is recognized as a major mechanism of drug-induced liver injury (DILI) [START_REF] Chatterjee | Drug-induced Cholestasis: Mechanisms, Models, and Markers[END_REF][START_REF] Fernández-Murga | Advances in drug-induced cholestasis: Clinical perspectives, potential mechanisms and in vitro systems[END_REF][START_REF] Mosedale | Drug-induced liver injury: Advances in mechanistic understanding that will inform risk management[END_REF]. It notably corresponds to the inhibition of the canalicular ATP-binding cassette (ABC) transporter bile salt export pump (BSEP/ABCB11) and results in bile acid retention into hepatocytes and cholestasis [START_REF] Dawson | In vitro inhibition of the bile salt export pump correlates with risk of cholestatic drug-induced liver injury in humans[END_REF][START_REF] Izat | Hepatic transporter-mediated pharmacokinetic drug-drug interactions: Recent studies and regulatory recommendations[END_REF]. Beside canalicular efflux, sinusoidal secretion of bile acids, mediated notably by the ABC transporters multidrug resistance-associated protein (MRP) 3 (ABCC3) and MRP4 (ABCC4), contributes to the regulation of bile acid accumulation into hepatocytes [START_REF] Keppler | Cholestasis and the role of basolateral efflux pumps[END_REF][START_REF] Köck | A perspective on efflux transport proteins in the liver[END_REF].

Inhibition of this sinusoidal bile acid efflux can consequently be considered as an additional cause of DILI [START_REF] Köck | Risk factors for development of cholestatic drug-induced liver injury: inhibition of hepatic basolateral bile acid transporters multidrug resistanceassociated proteins 3 and 4[END_REF]. Potential interactions of new molecular entities with BSEP and MRP3/MRP4 have therefore likely to be considered during their pharmaceutical development, in order to identify and discard candidate drugs that may cause DILIs [START_REF] Hillgren | Emerging transporters of clinical importance: an update from the International Transporter Consortium[END_REF]Kenna et al., 2018;[START_REF] Morgan | Interference with bile salt export pump function is a susceptibility factor for human liver injury in drug development[END_REF][START_REF] Morgan | A multifactorial approach to hepatobiliary transporter assessment enables improved therapeutic compound development[END_REF]. For such a purpose, various in vitro assays are available, such as sandwich-cultured hepatocytes-based efflux assays and BSEP-, MRP3-or MRP4-related vesicular assays [START_REF] Cheng | In vitro model systems to investigate bile salt export pump (BSEP) activity and drug interactions: A review[END_REF]Kenna et al., 2018;[START_REF] Kostrubsky | Evaluation of hepatotoxic potential of drugs by inhibition of bile-acid transport in cultured primary human hepatocytes and intact rats[END_REF][START_REF] Pedersen | Early identification of clinically relevant drug interactions with the human bile salt export pump (BSEP/ABCB11)[END_REF][START_REF] Van Brantegem | Vesicle-and Hepatocyte-Based Assays for Identification of Drug Candidates Inhibiting BSEP Function[END_REF][START_REF] Wolf | Use of cassette dosing in sandwich-cultured rat and human hepatocytes to identify drugs that inhibit bile acid transport[END_REF].

The use of cultured hepatocytes, but not that of vesicles, allows an integrated approach, with determination of drug effects towards both canalicular and sinusoidal efflux of bile acids such as taurocholate (TC) [START_REF] Susukida | Basal efflux of bile acids contributes to drug-induced bile acid-dependent hepatocyte toxicity in rat sandwich-cultured hepatocytes[END_REF]. In this context, owing to the differences between species with respect to hepatobiliary disposition of TC [START_REF] Yang | Species differences in hepatobiliary disposition of taurocholic acid in human and rat sandwich-cultured hepatocytes: implications for drug-induced liver injury[END_REF], primary human hepatocytes have likely to be preferentially considered. The poor and scarce availability of plateable human hepatocytes, together with the variability of their quality (depending on their production process and origin) and/or their drug detoxifying protein expression, may however limit their use, thus underlining the interest of surrogates for them, such as human hepatoma

Accepted Article

HepaRG cells [START_REF] Andersson | The HepaRG cell line: a unique in vitro tool for understanding drug metabolism and toxicology in human[END_REF][START_REF] G12 Lübberstedt | HepaRG human hepatic cell line utility as a surrogate for primary human hepatocytes in drug metabolism assessment in vitro[END_REF]. Indeed, HepaRG cells exhibit high levels of drug detoxifying proteins, including drug transporters [START_REF] Hammer | Cross-species analysis of hepatic cytochrome P450 and transport protein expression[END_REF]. In particular, they display expression of bile acid transporters, including BSEP, MRP3 and MRP4, as well as functional bile canaliculi (Bachour-El [START_REF] Azzi | Comparative Localization and Functional Activity of the Main Hepatobiliary Transporters in HepaRG Cells and Primary Human Hepatocytes[END_REF][START_REF] Vee | Polarized expression of drug transporters in differentiated human hepatoma HepaRG cells[END_REF]. They have consequently permitted to characterize the toxic effects of cholestatic drugs, especially towards bile canaliculi dynamics [START_REF] Burbank | Early Alterations of Bile Canaliculi Dynamics and the Rho Kinase/Myosin Light Chain Kinase Pathway Are Characteristics of Drug-Induced Intrahepatic Cholestasis[END_REF]. In the present study, we demonstrate using a biphasic TC loading/efflux assay that HepaRG cells can be used for simultaneously investigating the effects of cholestatic drugs on canalicular and sinusoidal efflux of TC. This illustrates thus the potential relevance of HepaRG cells for in vitro assays aimed at predicting DILIs due to impairment of hepatobiliary transport of bile acids.

Material and Methods

Chemicals

Bosentan, cyclosporin A, flumazenil, glibenclamide, ketoconazole, probenecid, rifampicin, salicylic acid and troglitazone were provided by MedChemExpress (Sollentuna, Sweden). The cholestatic or non-cholestatic character of these drugs as well as their known inhibitory effects towards BSEP and MRP3/MRP4 activities are reported in Table S1. Radiolabelled [ 3 H(G)]-TC (specific activity = 15.5 Ci/mmol) was from Perkin-Elmer (Villebon-sur-Yvette, France).

Cell culture

Highly differentiated human HepaRG cells from passages 13 to 17 were cultured as already described (Le [START_REF] Vee | Polarized expression of drug transporters in differentiated human hepatoma HepaRG cells[END_REF]. Briefly, cells, plated at a density of 2 × 10 5 cells/cm 2 , were first grown in Williams' E medium supplemented with 10% (vol/vol), 100 IU/ml penicillin, 100 μg/ml streptomycin, 5 μg/ml insulin, 2 mM glutamine, and 5 × 10 -5 M hydrocortisone hemisuccinate for 2 weeks. Cells were next cultured for an additional 2 weeks in the same medium supplemented with 2% (vol/vol) dimethyl sulfoxide (DMSO), to get a full

Accepted Article

differentiation status of the cells [START_REF] Gripon | Infection of a human hepatoma cell line by hepatitis B virus[END_REF]. This differentiated status of DMSOtreated HepaRG cell cultures was routinely checked by phase-contrast microscopy analysis, demonstrating the presence of hepatocyte-like islands and the formation of bile canaliculi, as previously reported (Le [START_REF] Vee | Polarized expression of drug transporters in differentiated human hepatoma HepaRG cells[END_REF].

Biphasic TC loading/efflux assay

HepaRG cells were first incubated with 12.9 nM [ 3 H]-TC in transport assay buffer comprising 136 mM NaCl, 5.3 mM KCl, 1.1mM KHPO4, 0.8 mM MgSO4, 1.8 mM CaCl2, 10 mM HEPES, 11 mM D-glucose and adjusted to pH=7.4 for 30 min at 37°C (TC loading phase). After washing twice with ice-cold phosphate-buffered saline, cells were re-incubated at 37°C for up to 60 min in the absence of TC (TC efflux phase), in the transport medium described above (Calciumcontaining medium) or with the same buffer, except that CaCl2 and MgSO4 were withdrawn and 1 mM EGTA was added (Calcium-free medium). Incubation with this calcium-free buffer promotes disruption of tight junctions and opening of bile canaliculi networks [START_REF] Annaert | Pglycoprotein-mediated in vitro biliary excretion in sandwich-cultured rat hepatocytes[END_REF][START_REF] Liu | Biliary excretion in primary rat hepatocytes cultured in a collagen-sandwich configuration[END_REF] and thus allows to distinguish between cellular and canalicular compartments (Kenna et al., 2018;[START_REF] Kostrubsky | Evaluation of hepatotoxic potential of drugs by inhibition of bile-acid transport in cultured primary human hepatocytes and intact rats[END_REF]. Chemicals tested towards canalicular and sinusoidal TC transport were added in the assay buffer only during the efflux phase. Cells were next washing twice with ice-cold phosphate-buffered saline. Intracellular levels of TC were finally measured by scintillation counting and were normalised to protein content, determined by the Bradford's method [START_REF] Bradford | A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding[END_REF]. They were expressed as pmol TC/mg protein or as % of initial TC loading. Biliary excretion index (BEI) and sinusoidal excretion index (SEI) were calculated using the equations (A) and (B): with BEI Control = BEI in the absence of drug, BEI Drug = BEI in the presence of drug during the efflux phase, SEI Control = SEI in the absence of drug and SEI Drug = SEI in the presence of drug during the efflux phase.

𝐵𝐸𝐼 (%) = [𝑇𝐶 𝐸𝑓𝑓𝑙𝑢𝑥/+𝐶𝑎 2+ ] -[𝑇𝐶 𝐸𝑓𝑓𝑙𝑢𝑥/-𝐶𝑎 2+ ] [𝑇𝐶 𝐸𝑓𝑓𝑙𝑢𝑥/+𝐶𝑎 2+ ] × 100 (𝐴) 𝑆𝐸𝐼 (%) = [𝑇𝐶 𝐿𝑜𝑎𝑑𝑖𝑛𝑔 ] -[𝑇𝐶 𝐸𝑓𝑓𝑙𝑢𝑥/+𝐶𝑎 2+ ] [𝑇𝐶 𝐿𝑜𝑎𝑑𝑖𝑛𝑔 ] × 100 (𝐵)

Statistical analysis

Experimental data were usually expressed as means ± SD. They were statistically analysed through analysis of variance (ANOVA) followed by the Dunnett's or the Tukey's post-hoc test,

or by Spearman's rank correlation, using the Prism 8.4.3 software (GraphPad Software, San Diego, CA, USA). The criterion of significance was p < 0.05.

Results and Discussion

To specifically and simultaneously analyse the effects of drugs towards canalicular and sinusoidal efflux of TC in HepaRG cells, we developed a biphasic functional assay in the present study, with a first step of TC loading (without drug) and a second step of TC efflux (with or without drug). Indeed, our preliminary assays demonstrated that at least the two prototypical cholestatic drugs cyclosporin A and troglitazone (Table S1), used either at 10 or

Accepted Article

100 µM, markedly reduced TC loading in HepaRG cells (Figure S1). This likely reflects the inhibitory effects of these compounds towards the sodium-taurocholate co-transporting polypeptide (NTCP/SLC10A1) [START_REF] Grosser | Substrate Specificities and Inhibition Pattern of the Solute Carrier Family 10 Members NTCP, ASBT and SOAT[END_REF], responsible for the entry of bile acids at the sinusoidal pole of hepatocytes [START_REF] Kullak-Ublick | Enterohepatic bile salt transporters in normal physiology and liver disease[END_REF]. Initial loading of cells without drugs was therefore a pre-requisite to reach high intracellular levels of TC before efflux and to avoid false or misleading outcomes that may occur if the tested compound interferes with TC uptake [START_REF] Cheng | In vitro model systems to investigate bile salt export pump (BSEP) activity and drug interactions: A review[END_REF].

In order to optimize the efflux time of TC from HepaRG cells maintained either in the presence or absence of calcium, i.e., with or without functional bile canaliculi, we next determined TC retention kinetic over a 60-min efflux period in our functional assay. As shown in Figure 1A, cellular retention of TC was falling rapidly in the absence of calcium, to 11.4%

of initial loading after a 20 min efflux period. Cellular retention of TC also decreased with efflux time in the presence of calcium, but less quickly; it thus corresponded to 40.2% of initial loading after a 20 min efflux period (Figure 1A). With respect to SEI and BEI values, reflecting sinusoidal and canalicular secretion of TC, respectively, they increased with efflux time, to reach the high values of 60%-80% after a 20 min efflux time (Figure 1B and Figure 1C). This 20-min efflux time was consequently retained for further experiments investigating the effects of drugs toward both sinusoidal and canalicular efflux of TC.

BEI and SEI values in the absence or the presence of various drugs, used at 10 µM and 100 µM for most of them, are shown in Figure S2A and Figure S2B and corresponding inhibitions of canalicular and sinusoidal TC efflux, i.e., BEI and SEI inhibition %, are indicated in Table 1. The markedly increased retention of TC in response to the reference cholestatic agent troglitazone in HepaRG cells, in the presence or absence of calcium, as well as the lack of effect of the non-cholestatic compound salicylic acid, are additionally shown in Figure 1D.

Most of the cholestatic compounds included in the study, i.e., bosentan, cyclosporin A,

Accepted Article

glibenclamide, ketoconazole and troglitazone, were found to inhibit both canalicular and sinusoidal efflux of TC when used at 100 µM, whereas, by contrast, the non-cholestatic compounds flumazenil and salicylic acid failed to impair TC secretion (Table 1). This absence of TC leakage in response to cholestatic drugs used at 100 µM likely discards any major unspecific toxic effect of these chemicals; the fact that treatment by drugs, in the absence or presence of calcium, failed to alter cellular protein concentrations (Figure S3) also argues in favour of this conclusion. The dual repression of canalicular and sinusoidal TC efflux by cholestatic compounds indicates that these chemicals may act through inhibiting both BSEP and MRP3/MRP4, as already suggested using vesicular assays [START_REF] Köck | Risk factors for development of cholestatic drug-induced liver injury: inhibition of hepatic basolateral bile acid transporters multidrug resistanceassociated proteins 3 and 4[END_REF][START_REF] Morgan | A multifactorial approach to hepatobiliary transporter assessment enables improved therapeutic compound development[END_REF]. The fact that the percentages of BEI inhibition by drugs used either at 10 or 100 µM were significantly correlated with those of SEI inhibition (Figure S4) fully supports this hypothesis. Nevertheless, accurate comparison of IC50 values of cholestatic compounds towards BSEP-, MRP3 and MRP4-mediated transport of TC remains to be determined and would deserve further studies. Moreover, the concentrations required to inhibit canalicular or sinusoidal efflux of TC appeared to differ for some cholestatic drugs; indeed, glibenclamide used at 10 µM inhibited sinusoidal efflux of TC, but not canalicular one, whereas troglitazone used at 10 µM blocked canalicular, but not sinusoidal, efflux of TC (Table 1). This confirms that canalicular and sinusoidal efflux of TC are mediated by different transporters, that may exhibit differential sensitivities to inhibitors. This assertion may also be supported by the fact that probenecid, used at a 2 mM concentration known to block BSEP activity (Table S1) [START_REF] Dawson | In vitro inhibition of the bile salt export pump correlates with risk of cholestatic drug-induced liver injury in humans[END_REF], inhibited only canalicular efflux of TC in HepaRG cells (Table 1).

Rifampicin was also found to differently impair the two types of TC efflux in HepaRG cells;

thus, this compound, used either at 10 or 100 µM, markedly inhibited (by at least 50%) sinusoidal efflux of TC from HepaRG cells, but exerted a less potent inhibitory effect towards canalicular efflux (around 25% of inhibition), which moreover failed to reach a significant level

Accepted Article

(Table 1). Such a lack of major inhibition of canalicular TC efflux from HepaRG cells may be rather surprising, because rifampicin is usually considered as a potent inhibitor of BSEP, as notably demonstrated by vesicular assays [START_REF] Dawson | In vitro inhibition of the bile salt export pump correlates with risk of cholestatic drug-induced liver injury in humans[END_REF][START_REF] Stieger | Drug-and estrogen-induced cholestasis through inhibition of the hepatocellular bile salt export pump (Bsep) of rat liver[END_REF]. The nature of the assays used for investigating BSEP inhibition (vesicular assays versus cells-based assays) may be hypothesized to contribute to this discrepancy. In this context, it is noteworthy that our HepaRG cells-based assay is likely more physiological than vesicles-based counterparts, and notably allows an holistic approach, taking into account both sinusoidal and canalicular secretion of TC. By this way, interplay of BSEP and other bile acid efflux transporters such as MRP3/MRP4, and possibly also the canalicular transporter MRP2 (ABCC2) and the sinusoidal organic solute transporter (OST) α/β (SLC51A/SLC51B), which may additionally be implicated in cholestasis (Kenna et al., 2018;[START_REF] Malinen | Organic solute transporter OSTα/β is overexpressed in nonalcoholic steatohepatitis and modulated by drugs associated with liver injury[END_REF], are considered. Moreover, TC is used as substrate for sinusoidal efflux from HepaRG cells, whereas by contrast, the MRP3/MRP4 probe substrate in vesicular assays is usually an estradiol derivative [START_REF] Morgan | A multifactorial approach to hepatobiliary transporter assessment enables improved therapeutic compound development[END_REF], therefore less relevant for bile acid transport issues. Besides, the possible role of drug metabolites in bile salt transporter inhibition can likely be explored in HepaRG cells, which display high expression of drug metabolizing enzymes [START_REF] Andersson | The HepaRG cell line: a unique in vitro tool for understanding drug metabolism and toxicology in human[END_REF]; this may allow to take into account the transporter-enzyme interplay, relevant for hepatic processing of drugs [START_REF] Matsunaga | Sandwich-Cultured Hepatocytes for Mechanistic Understanding of Hepatic Disposition of Parent Drugs and Metabolites by Transporter-Enzyme Interplay[END_REF]. The use of calcium-free conditions, which corresponds to a nonphysiological situation, may however impair the exposure of canalicular transporters to drugs in HepaRG cells, which may constitute a limitation of the experimental system.

The present study fully supports the use of HepaRG cells for predicting bile acid transport inhibition by drugs and, beyond DILIs, as already claimed [START_REF] Burbank | Early Alterations of Bile Canaliculi Dynamics and the Rho Kinase/Myosin Light Chain Kinase Pathway Are Characteristics of Drug-Induced Intrahepatic Cholestasis[END_REF][START_REF] De Bruijn | Hepatic bile acid synthesis and secretion: Comparison of in vitro methods[END_REF]. This also confirms that HepaRG cells may serve as surrogates for human hepatocytes, as proposed for drug metabolism, transport and toxicity studies [START_REF] Andersson | The HepaRG cell line: a unique in vitro tool for understanding drug metabolism and toxicology in human[END_REF][START_REF] Guillouzo | The human hepatoma HepaRG cells: a highly differentiated model for studies of liver metabolism and toxicity of xenobiotics[END_REF][START_REF] Vee | Functional expression of sinusoidal and canalicular hepatic drug transporters in the differentiated human hepatoma HepaRG cell line[END_REF][START_REF] G12 Lübberstedt | HepaRG human hepatic cell line utility as a surrogate for primary human hepatocytes in drug metabolism assessment in vitro[END_REF], including high-

Accepted Article

throughput ones [START_REF] Ooka | Application of In Vitro Metabolism Activation in High-Throughput Screening[END_REF]. It should however be kept in mind that in vitro prediction of DILIs due to transporter inhibition remains rather complex [START_REF] Chan | Measures of BSEP Inhibition In Vitro Are Not Useful Predictors of DILI[END_REF][START_REF] Garzel | A Change in Bile Flow: Looking Beyond Transporter Inhibition in the Development of Drug-induced Cholestasis[END_REF][START_REF] Kenna | Do In Vitro Assays Predict Drug Candidate Idiosyncratic Drug-Induced Liver Injury Risk?[END_REF][START_REF] Kotsampasakou | Predicting drug-induced liver injury: The importance of data curation[END_REF]. In particular, cutoff values for BSEP IC50 values with respect to "concerning levels" of BSEP inhibition have to be harmonized [START_REF] Hafey | A Two-Tiered In Vitro Approach to De-Risk Drug Candidates for Potential Bile Salt Export Pump Inhibition Liabilities in Drug Discovery[END_REF]. In addition, the unbound drug concentration in hepatocytes at the site of interaction with BSEP or MRP3/MRP4 remains a key-parameter to consider for predicting in vivo inhibition of hepatic bile acid efflux from in vitro data (Kenna et al., 2018;[START_REF] Riede | Assessing the Risk of Drug-Induced Cholestasis Using Unbound Intrahepatic Concentrations[END_REF]. Hepatic partition coefficient (Kp) is also a notable point to take into account for improving hepatic exposure to cholestatic drugs [START_REF] Martin | Early Drug-Induced Liver Injury (DILI) Risk Screening: "Free", as good as it gets[END_REF]. In this context, direct measurement and modeling of intracellular unbound drug concentrations in HepaRG cells may help to evaluate in vivo inhibition of canalicular and/or sinusoidal secretion of bile acids, as recently described for assessment of acetaminophen toxicity through integration of in vitro data from cultured HepaRG cells and physiologically-based pharmacokinetic modeling (Zhang et al., 2020). Finally, it is noteworthy that regulation of hepatic transporter expression, including that of BSEP, can be analysed in HepaRG cells (Le [START_REF] Vee | Functional expression of sinusoidal and canalicular hepatic drug transporters in the differentiated human hepatoma HepaRG cell line[END_REF], which are therefore also suitable for studying drug-induced repression of BSEP, a possible cause of cholestasis recently reported for metformin [START_REF] Garzel | Metformin Disrupts Bile Acid Efflux by Repressing Bile Salt Export Pump Expression[END_REF].

Conclusion

Human hepatoma HepaRG cells were demonstrated to exhibit canalicular and sinusoidal efflux of TC through a biphasic functional assay consisting of a first step of TC loading and a second step of TC efflux, in the presence or absence of functional bile canaliculi. These canalicular and sinusoidal secretions of TC were found to be inhibitable by cholestatic compounds, thus underlining the interest of HepaRG cell cultures for the in vitro prediction of bile acid secretion inhibition by drugs and, beyond, of DILIs.

  Accepted Articlewith [TC Efflux/+Ca 2+] = TC concentration after the efflux phase in the presence of Ca 2+ , i.e., in cells + bile canaliculi after efflux, [TC Efflux/-Ca 2+] = TC concentration after the efflux phase in the absence of calcium, i.e., in cells after efflux, and [TC Loading ] = TC accumulation after the loading phase in the presence of calcium, i.e., in cells + bile canaliculi after loading.

	Effects of drugs towards canalicular or sinusoidal TC efflux, corresponding to % of BEI or SEI
	inhibition, were finally determined according to the equations (C) and (D):
	𝐶𝑎𝑛𝑎𝑙𝑖𝑐𝑢𝑙𝑎𝑟 𝑇𝐶 𝑒𝑓𝑓𝑙𝑢𝑥 𝑖𝑛ℎ𝑖𝑏𝑖𝑡𝑖𝑜𝑛 % = 𝐵𝐸𝐼 𝑖𝑛ℎ𝑖𝑏𝑖𝑡𝑖𝑜𝑛 % =	𝐵𝐸𝐼 𝐶𝑜𝑛𝑡𝑟𝑜𝑙 -𝐵𝐸𝐼 𝐷𝑟𝑢𝑔 𝐵𝐸𝐼 𝐶𝑜𝑛𝑡𝑟𝑜𝑙	× 100 (𝐶)
	𝑆𝑖𝑛𝑢𝑠𝑜𝑖𝑑𝑎𝑙 𝑇𝐶 𝑒𝑓𝑓𝑙𝑢𝑥 𝑖𝑛ℎ𝑖𝑏𝑖𝑡𝑖𝑜𝑛 % = 𝑆𝐸𝐼 𝑖𝑛ℎ𝑖𝑏𝑖𝑡𝑖𝑜𝑛 % =	𝑆𝐸𝐼 𝐶𝑜𝑛𝑡𝑟𝑜𝑙 -𝑆𝐸𝐼 𝐷𝑟𝑢𝑔 𝑆𝐸𝐼 𝐶𝑜𝑛𝑡𝑟𝑜𝑙	× 100 (𝐷)

Conflict of interest

The authors have no conflict of interest to disclose. (A) HepaRG, initially loaded with 12.9 nM [ 3 H] TC for 30 min, were re-incubated in TC-free medium in the presence or absence of calcium for various times (from 5 min to 60 min). TC retention was next determined by scintillation counting. Data are expressed as percentages of initial TC loading and are the means ± SD of at least 3 independent experiments. (B) TC sinusoidal excretion index (SEI) and (C) biliary excretion index (BEI) values were determined for various times of TC efflux (from 5 to 60 min), as described in Materials and Methods; they are the means ± SD of at least 3 independent experiments. (D) HepaRG cells, initially loaded with 12.9 nM [ 3 H] TC for 30 min, were re-incubated for 20 min in TC-free medium in the presence or absence of calcium, and without (control) or with troglitazone or salicylic acid, used at 10 µM or 100 µM. TC retention was finally determined by scintillation counting. Data are expressed as percentages of initial TC loading. They are the means ± SD of 3 independent experiments and have been analysed through ANOVA followed by the Tukey's post-hoc test (***, p<0.001 when compared to corresponding controls, i.e., control in the presence of calcium or control in the absence of calcium). 81.0 ± 9.5*** 49.3 ± 6.5*** a data are expressed as means ± SD of three independent experiments and have been analysed though ANOVA and Dunnett's post-hoc test (**, p<0.01 and ***, p<0.001 when compared to control cells not exposed to drugs).

Accepted Article

Legend to figure

Accepted Article


Niels J Dingemanse

Anne R Hertel

Raphael Royauté

Moving away from repeatability: a comment on Stuber et al. 1

Frankfurt am Main, Germany 6

Behavioural ecologists increasingly acknowledge that among-and within-individual processes are 7 important; this notion also reached movement ecology [START_REF] Hertel | A guide for studying among-individual 62 behavioral variation from movement data in the wild[END_REF]. [START_REF] Stuber | Spatial Personalities: a meta-analysis of consistent individual 69 differences in spatial behavior[END_REF] show 8 that movement behaviours are individually repeatable but caution that "repeatability itself does not 9 indicate the mechanisms by which consistent differences arise". We here point out that comparing 10 repeatabilities is, in fact, a comparison of apples and oranges. This is because repeatability (𝑅) is a 11 proportion that varies as a function of multiple variance components: 12

𝑅 = 𝑉 𝐼 𝑉 𝑃 = 𝑉 𝐼 𝑉 𝐼 +𝑉 𝐸 +𝑉 𝑒 (Eqn. 1) 13
In its simplest form, the phenotypic variance (𝑉 𝑃 ) is the sum of among-individual variance (𝑉 𝐼 ), within-14 individual variance ( 𝑉 𝐸 ), and measurement error ( 𝑉 𝑒 ). 𝑅 is the proportion of 𝑉 𝑃 explained by 𝑉 𝐼 . 15 Differences in 𝑅 therefore do not necessarily imply equivalent differences in 𝑉 𝐼 (the statistical definition 16 of 'consistent differences'; [START_REF] Dingemanse | Behavioural reaction norms: animal personality 59 meets individual plasticity[END_REF]). We give an example (with exaggerated made-up 17 numbers). Location fixes form the input of most movement metrics. Spatial measurement errors 18 associated with location fixes are smaller when comparing, for example, satellite-based (GPS-tags) with 19 telemetry-based (VHF-transmitters) fixes. Even when 𝑉 𝐼 (set to 0.3) and 𝑉 𝐸 (set to 0.6) do not vary with 20 methodology, repeatability will thus be greater for satellite-based (𝑅 𝐴 ) vs. telemetry-based (𝑅 𝐵 ) datasets: In Eqn. 2, we set 𝑉 𝑒 𝐴 to 0.1 and 𝑉 𝑒 𝐵 to 1.1. Measurement error can also vary with the environment, e.g. 24 precision is typically higher for GPS-fixes acquired from open versus forested areas. From this, we can 25 learn three lessons. 26 First, effects of fixed/random effects estimated from repeatability meta-analyses result both from 27 biology (𝑉 𝐼 , 𝑉 𝐸 ) and methodology (𝑉 𝑒 ). Drawing unambiguous conclusions is hard. Second, acquiring 28 unbiased estimates necessitate sampling designed to estimate measurement error. This requires not just 29 repeated measures over time spans so long that the animals might have moved but also repeated measures 30 over time spans so short they cannot have moved. Only then can 𝑉 𝐼 , 𝑉 𝐸 and 𝑉 𝑒 be estimated 31 simultaneously, and 𝑅 re-defined in biological terms (cf. [START_REF] Falconer | Introduction to Quantitative Genetics[END_REF]: 32

𝑅 = 𝑉 𝐼 𝑉 𝐼 +𝑉 𝐸 (Eqn. 3) 33
In our worked example, 𝑅 𝐴 would then correctly equate 𝑅 𝐵 . A third lesson is that there is important 34 biology hidden in repeated measures datasets. While repeatability is often discussed as a proxy for 𝑉 𝐼 , 𝑉 𝐸 35 is also of key interest. 𝑉 𝐸 results from learning, habituation, and reversible plasticity-all important 36 proximate drivers of movement behaviour. 37 Building upon Stuber et al., we recommend a re-focus from repeatability toward variances. This 38 will reveal processes underpinning movement behavior. For example, high within-individual variability 39 suggests a major role for reversible plasticity. This mechanism makes populations able to respond to 40 environmental change within generations. High among-individual variability instead presents greater 41 opportunity for adaptive evolution among generations. Evolutionary biologists have made similar calls 42 to move away from proportions. This is because variances-but not proportions-are comparable across 43 samples, at least when traits are measured in units with natural nulls, and their variances expressed as 44

Senckenberg Biodiversity and Climate Research Centre (SBiK-F), Senckenberganlage 25, 60325

of coefficients of among-individual variation (𝐶𝑉 𝐼 )[START_REF] Holtmann | Metabolic rates, and not hormone levels, are a likely mediator 64 of between-individual differences in behaviour: a meta-analysis[END_REF] in meta-analytical studies of

individual variation (e.g., Royauté and Dochtermann, 2021):


Raphaël Teboul

Michalina Grabias

Jessica Zucman-Rossi

Eric Letouzé
email: eric.letouze@inserm.fr

Discovering cryptic splice mutations in cancers via a deep neural network framework

Somatic mutations can disrupt splicing regulatory elements and have dramatic effects on cancer genes, yet the functional consequences of mutations located in extended splice regions is difficult to predict. Here, we use a deep neural network (SpliceAI) to characterize the landscape of splice-altering mutations in cancer. In our in-house series of 401 liver cancers, SpliceAI uncovers 1,244 cryptic splice mutations, located outside essential splice sites, that validate at a high rate (66%) in matched RNAseq data. We then extend the analysis to a large pan-cancer cohort of 17,714 tumors, revealing >100,000 cryptic splice mutations. Taking into account these mutations increases the power of driver gene discovery, revealing 126 new candidate driver genes. It also reveals new driver mutations in known cancer genes, doubling the frequency of splice alterations in tumor suppressor genes. Mutational signature analysis suggests mutational processes that could give rise preferentially to splice mutations in each cancer type, with an enrichment of signatures related to clock-like processes and DNA repair deficiency. Altogether, this work sheds light on the causes and impact of cryptic splice mutations in cancer, and highlights the power of deep learning approaches to better annotate the functional consequences of mutations in oncology.

INTRODUCTION

RNA splicing is a key process in biology where non-coding regions of RNA (introns) are removed and exons joined together, to generate the mature messenger RNA (mRNA) that will be translated [START_REF] Black | Mechanisms of alternative pre-messenger RNA splicing[END_REF]. In eukaryotes, alternative splicing occurs in ~95% of multi-exonic genes and generates a diversity of protein isoforms from a single gene [START_REF] Pan | Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing[END_REF]. This process provides evolutionary flexibility and may have contributed to the development of multicellular organisms [START_REF] Romero | Alternative splicing in concert with protein intrinsic disorder enables increased functional diversity in multicellular organisms[END_REF]. In humans, splicing involves a series of reactions catalyzed by a large ribonucleoprotein complex called the spliceosome, which recognizes key regulatory sequences on the pre-mRNA. The two bases at each extremity of the intron constitute the canonical donor (GU) and acceptor (AG) splice sites. Other cis-regulatory elements include the polypyrimidine tract and the branch point located around 30 nucleotides upstream the acceptor site, as well as exonic and intronic splicing enhancers and silencers [START_REF] Wang | Splicing regulation: From a parts list of regulatory elements to an integrated splicing code[END_REF]. Yet, these sequences are less evolutionary conserved than the essential splice sites. RNA splicing is thus carefully regulated, and alterations of the splice machinery components or cis-regulatory sequences underlie many diseases [START_REF] Sterne-Weiler | Exon identity crisis: disease-causing mutations that disrupt the splicing code[END_REF][START_REF] Xiong | The human splicing code reveals new insights into the genetic determinants of disease[END_REF][START_REF] Scotti | RNA mis-splicing in disease[END_REF]. In cancer, core spliceosome components including SF3B1, SRSF2 or U2AF1 are recurrently altered in both hematological and solid tumors [START_REF] Anczuków | Splicing-factor alterations in cancers[END_REF][START_REF] Sveen | Aberrant RNA splicing in cancer; expression changes and driver mutations of splicing factor genes[END_REF]. In addition, splice mutations are a frequent mechanism leading to the inactivation of tumor suppressor genes, and have been reported in 13% of genes in the Cancer Gene Census [START_REF] Sveen | Aberrant RNA splicing in cancer; expression changes and driver mutations of splicing factor genes[END_REF][START_REF] Sondka | The COSMIC Cancer Gene Census: describing genetic dysfunction across all human cancers[END_REF]. Mutations affecting essential splice sites lead to various degrees of intron retention or exon skipping and are thus annotated as damaging in cancer genomics projects [START_REF] Calabrese | Genomic basis for RNA alterations in cancer[END_REF]. In contrast, the functional consequences of mutations in extended splice regions is difficult to predict because other regulatory elements largely involve degenerate sequence features [START_REF] Witten | Understanding splicing regulation through RNA splicing maps[END_REF]. Several computational methods have been developed to automatically detect splice sites in DNA sequences and to predict the splicing impact of mutations. One of the earliest tools for detecting splice sites in eukaryotic mRNA was a decision tree via maximal dependence decomposition [START_REF] Burge | Prediction of complete gene structures in human genomic DNA[END_REF]. This idea was enhanced in GeneSplicer by Markov models that capture additional dependencies among neighboring bases in splice regions [START_REF] Pertea | GeneSplicer: a new computational method for splice site prediction[END_REF]. Other methods were introduced to model splice sites with Bayes networks [START_REF] Cai | Modeling splice sites with Bayes networks[END_REF] or pairwise correlations [START_REF] Arita | Modeling splicing sites with pairwise correlations[END_REF]. Yeo and Burge proposed the MaxEntScan framework modeling sequence motifs based on the maximum entropy principle [START_REF] Yeo | Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals[END_REF]. This tool was shown to outperform previous probabilistic models and is used by variant annotation tools like the Variant Effect Predictor (VEP) [START_REF] Shamsani | A plugin for the Ensembl Variant Effect Predictor that uses MaxEntScan to predict variant spliceogenicity[END_REF]. Neural networks were introduced to identify splice sites in NNSplice [START_REF] Reese | Improved splice site detection in Genie[END_REF]. More recently, a 32-layer deep residual neural network (SpliceAI) was developed to predict whether each position in a pre-mRNA transcript is a splice donor, acceptor or neither [START_REF] Jaganathan | Predicting Splicing from Primary Sequence with Deep Learning[END_REF]. In contrast to previous methods that only consider tens of bases around exon-intron boundaries, SpliceAI uses 10,000 nucleotide sequences as input, allowing to capture long-range determinants of splicing. This tool largely outperformed previous methods in terms of both specificity and validation rate [START_REF] Jaganathan | Predicting Splicing from Primary Sequence with Deep Learning[END_REF] and was successfully used to unravel pathogenic cryptic splice mutations in patients with rare genetic diseases.

Here, we use a large series of liver tumors with matched whole genome/exome sequencing and RNA sequencing data to benchmark SpliceAI in a cancer setting. We then apply this tool to uncover cryptic splice mutations in a large pan-cancer series totaling >17,000 tumors. We show that SpliceAI outperforms previous widely used tools, and reveals numerous cryptic splice mutations that enhance the landscape of cancer driver genes and mutations.

MATERIAL AND METHODS

Datasets

We used an in-house series of 401 liver cancers with matched RNA-seq and DNA variants identified by whole exome (WES, n=275) or whole genome (WGS, n=126) sequencing to benchmark the performance of SpliceAI. Initial bioinformatic analyses (alignment and variant calling) were already performed in previous studies [START_REF] Schulze | Exome sequencing of hepatocellular carcinomas identifies new mutational signatures and potential therapeutic targets[END_REF][START_REF] Letouzé | Mutational signatures reveal the dynamic interplay of risk factors and cellular processes during liver tumorigenesis[END_REF][START_REF] Bayard | Cyclin A2/E1 activation defines a hepatocellular carcinoma subclass with a rearrangement signature of replication stress[END_REF][START_REF] Hirsch | BAP1 mutations define a homogeneous subgroup of hepatocellular carcinoma with fibrolamellar-like features and activated PKA[END_REF][START_REF] Péneau | Hepatitis B virus integrations promote local and distant oncogenic driver alterations in hepatocellular carcinoma[END_REF][START_REF] Hirsch | Integrated Genomic Analysis Identifies Driver Genes and Cisplatin-Resistant Progenitor Phenotype in Pediatric Liver Cancer[END_REF][START_REF] Molina | Bi-allelic hydroxymethylbilane synthase inactivation defines a homogenous clinico-molecular subtype of hepatocellular carcinoma[END_REF] and we directly used processed data in this study. We then analyzed a large pan-cancer series comprising 17,714 tumor samples across 25 cancer types from the ICGC-PCAWG (n=2,699 WGS) [START_REF] Campbell | Pan-cancer analysis of whole genomes[END_REF], Hartwig Medical Foundation (n=4,818 WGS) [START_REF] Priestley | Pan-cancer whole-genome analyses of metastatic solid tumours[END_REF] and TCGA (n=10,197 WES) [START_REF] Hoadley | Cell-of-Origin Patterns Dominate the Molecular Classification of 10,000 Tumors from 33 Types of Cancer[END_REF] projects. Somatic mutations were downloaded directly from the online repositories of each project.

Predicting splice-affecting mutations with SpliceAI

We used SpliceAI (v1.3.1) to predict the splicing impact of almost 300 million somatic mutations. This computation-intensive task was performed on the Jean Zay supercomputer, the world's 10th most powerful supercomputer, with a computing power reaching 28 petaflops. Computations were heavily performed on a GPU V100 allocation which allowed to shrink the computation time down to 0.003 seconds. For comparison, the run time on a personal computer's CPU was 1.2s per variant, so the analysis of 300 million mutations would take 4,167 days. SpliceAI returns 4 delta scores (DS) for each type of splicing abnormality (donor loss/gain, acceptor loss/gain). We retrieved for each mutation the strongest DS, the corresponding splicing consequence and its position relative to the mutation.

Mutations with a DS ≥ 0.5 were considered to have a significant impact on splicing, as suggested in the original publication [START_REF] Jaganathan | Predicting Splicing from Primary Sequence with Deep Learning[END_REF].

Validation of splice-affecting mutations using matched RNA-seq data

Let 𝑝 be the position of a predicted abnormal splicing event. Let 𝑅 # be the set of reads crossing position 𝑝 in a control sample with a MapQ score ≥ 20 and 𝑅 $ be the set of reads crossing position p in the mutated sample of interest with a MapQ score ≥ 20. We define 𝑀 & and 𝑀 ' as a partition of 𝑅 $ such as 𝑅 $ = 𝑀 & ∪ 𝑀 ' . We define 𝐶 & and 𝐶 ' as a partition of 𝑅 # such as 𝑅 # = 𝐶 & ∪ 𝐶 ' . 𝑀 ' (respectively 𝑀 & ) represents the ensemble of abnormal (respectively normal) reads in the mutated sample of interest. 𝐶 ' (respectively 𝐶 & ) represents the ensemble of abnormal (respectively normal) reads in the control samples. We calculated the Relative Usage of the Novel Junction (RUNJ, Figure 1C) as previously described [START_REF] Jaganathan | Predicting Splicing from Primary Sequence with Deep Learning[END_REF]:

𝑅𝑈𝑁𝐽 = #𝑀 ' #𝑀 ' + #𝑀 & - #𝐶 ' #𝐶 ' + #𝐶 &
Two custom control panels were used for each mutation: a control PON (Panel of Normal) composed of 5 non-tumor samples from the same sequencing series as the sample of interest, and a control POT (Panel of Tumors) composed of 5 tumor samples from the same sequencing series and without a somatic mutation at the same locus. The POT allows to exclude tumor-specific effects not directly linked with the mutation of interest (e.g. altered gene expression of isoform usage) that could lead to abnormal RNA-seq reads in tumors. RUNJ scores were calculated in two ways. First, the RUNJdirect score considers only abnormal reads that support the most direct consequence of the splicing alteration (e.g. reads passing trough the essential donor for donor loss events). However, the observed consequence on the RNA can be different from the expected one. For example, a donor loss could favor usage of a cryptic donor site upstream within the exon, leading to a truncated exon or even full exon skipping. In that case, we would not detect RNA-seq reads continuing from the exon to the intron, but rather abnormal junctions linking the cryptic donor to the next acceptor site. To account for these situations, we thus introduced a RUNJall score considering all abnormal reads at the considered exon-intron boundary (Supplementary Figure S1). 𝑅𝑈𝑁𝐽 1234#5,789 and 𝑅𝑈𝑁𝐽 '::,789 were calculated considering the Panel of Normal as controls. 𝑅𝑈𝑁𝐽 1234#5,78< and 𝑅𝑈𝑁𝐽 '::,78< were calculated considering the Panel of Tumors as controls.

Finally, a prediction was validated if and only if all of the following hypotheses are verified:

-#𝑀 ' ≥ 2 -𝑅𝑈𝑁𝐽 1234#5,789 ≥ 0.01 or 𝑅𝑈𝑁𝐽 '::,789 ≥ 0.01 -𝑅𝑈𝑁𝐽 1234#5,78< ≥ 0.01 or 𝑅𝑈𝑁𝐽 '::,78< ≥ 0.01 or not applicable.

To facilitate the validation of large numbers of variants, we created the PyRNA library that allows retrieving relevant reads from RNA-seq bam files and calculating RUNJ scores: https://github.com/FunGeST/pyRNA.

Benchmarking the performance of SpliceAI against GeneSplicer and MaxEntScan

We used the two widely used tools GeneSplicer [START_REF] Pertea | GeneSplicer: a new computational method for splice site prediction[END_REF] and MaxEntScan [START_REF] Yeo | Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals[END_REF] to benchmark the performance of SpliceAI in predicting splice-affecting mutations. We downloaded the MaxEntScan and GeneSplicer softwares from http://genes.mit.edu/burgelab/maxent/download/ and http://www.cs.jhu.edu/genomics/GeneSplicer/ respectively. Both tools were run with default parameters to predict the impact on splicing of the 1,704,461 somatic mutations in th LICA-FR series.

For comparison, we applied thresholds for GeneSplicer and MaxEntScan scores in order to have the same number of positive predictions as SpliceAi for three DS thresholds (0.2, 0.5, 0.8). We validated the predictions of the 3 tools using matched RNA-seq data as described above. We computed the proportion validated predictions and the area under the precision-recall curve (PR-AUC) for each tool.

Driver gene analysis

We used MutSigCV to identify genes with significantly more mutations than expected by chance considering mutation categories and genomic covariates, as previously described [START_REF] Lawrence | Mutational heterogeneity in cancer and the search for new cancer-associated genes[END_REF]. For each tumor series, we first run MutSigCV with classical categories of coding mutations (missense, nonsense, essential splice, frameshift, inframe indels or synonymous). We then re-run the analysis after converting to splice variants the cryptic splice mutations uncovered by SpliceAI (DS ≥ 0.5). We then compared MutSigCV q-values obtained with both settings to identify new genes reaching significance when adding cryptic splice variants. We also calculated the mutation frequency of each gene in the Cancer Gene Census [START_REF] Sondka | The COSMIC Cancer Gene Census: describing genetic dysfunction across all human cancers[END_REF] considering only classical categories of coding mutations or adding cryptic splice sites.

Mutational signature analysis

We used Palimpsest [START_REF] Shinde | Palimpsest: an R package for studying mutational and structural variant signatures along clonal evolution in cancer[END_REF] to quantify the contribution of known COSMIC signatures in 2,699 tumors of the ICGC PCAWG data set. To ensure that only relevant signatures were included, we analyzed each cancer type separately and we included only signatures that were shown to be present in that cancer type [START_REF] Alexandrov | The repertoire of mutational signatures in human cancer[END_REF]. We then used Palimpsest to quantify the probability 𝑝 $,B2C of each mutation m being due to each mutational process with signature sig, as previously described [START_REF] Letouzé | Mutational signatures reveal the dynamic interplay of risk factors and cellular processes during liver tumorigenesis[END_REF]. Finally, we compared the contribution of signatures to splice mutations and other mutations in each series. The contribution of signature sig to a set of mutations Sm was estimated as the sum:

D 𝑃(𝑚, 𝑠) $J5∈L M
To identify mutational processes preferentially giving rise to splice mutations, we compared the distribution of probabilities 𝑃(𝑚, 𝑠) for each signatures s between splice mutations and others using Wilcoxon rank sum tests.

RESULTS

SpliceAI outperforms traditional tools for predicting cryptic splice mutations in cancer

We used an in-house series of 401 liver cancers [START_REF] Schulze | Exome sequencing of hepatocellular carcinomas identifies new mutational signatures and potential therapeutic targets[END_REF][START_REF] Letouzé | Mutational signatures reveal the dynamic interplay of risk factors and cellular processes during liver tumorigenesis[END_REF][START_REF] Bayard | Cyclin A2/E1 activation defines a hepatocellular carcinoma subclass with a rearrangement signature of replication stress[END_REF][START_REF] Hirsch | BAP1 mutations define a homogeneous subgroup of hepatocellular carcinoma with fibrolamellar-like features and activated PKA[END_REF][START_REF] Péneau | Hepatitis B virus integrations promote local and distant oncogenic driver alterations in hepatocellular carcinoma[END_REF][START_REF] Hirsch | Integrated Genomic Analysis Identifies Driver Genes and Cisplatin-Resistant Progenitor Phenotype in Pediatric Liver Cancer[END_REF][START_REF] Molina | Bi-allelic hydroxymethylbilane synthase inactivation defines a homogenous clinico-molecular subtype of hepatocellular carcinoma[END_REF] with whole exome (n=275) or whole genome sequencing (n=126) and matched RNA-seq data to predict and validate cryptic splice mutations using SpliceAI, MaxEntScan and GeneSplicer (Figure 1A). A total of 1,704,461 somatic mutations were analyzed with the three tools. SpliceAI returned predictions for 432,773 mutations, covering most intragenic mutations. By contrast, GeneSplicer returned predictions for only 13,589 mutations, and MaxEntScan 6,317 mutations, mostly located in splice regions. SpliceAI provides for each mutation a delta score (DS) representing the probability that a mutation creates (positive DS) or disrupts (negative DS) a splice site. 2,056 mutations were predicted to alter splicing with an absolute DS ≥ 0.5 (threshold recommended by SpliceAI developers), including 649 acceptor loss (AL), 753 donor loss (DL), 263 acceptor gain (AG) and 391 donor gain (DG) events (Supplementary Table S1 and Supplementary Figure S2). Mutations with an absolute DS ≥ 0.9 were strongly enriched in essential splice sites (67%).

By contrast, 60% of mutations with an absolute DS ≥ 0.5 were located outside canonical splice sites (Figure 1B). For comparison, we applied thresholds on GeneSplicer and MaxEntScan scores so that the number of mutations predicted to impact splicing was equivalent for the 3 tools. We then used matched RNA-seq data to validate the predictions. For each mutation predicted to affect splicing, we identified reads consistent with abnormal vs. normal splicing. We then calculated the ratio of abnormal reads over the total number of reads (UNJ score), which we normalized against UNJ scores in a series of normal and non-mutated tumor samples to obtain the RUNJ scores (Figure 1C). Mutations with a RUNJ ≥ 0.01 and at least 2 reads supporting the abnormal splicing events were considered validated. 357 mutations were predicted to alter splicing by all the methods and showed the best validation rate (90%, Figure 1D). SpliceAI and MaxEntScan had more consistent predictions, with 1,159 mutations predicted to impact splicing by both tools but not GeneSplicer. Among tool-specific mutations, SpliceAIonly mutations showed the best validation rate (64 % vs. 41% for MaxEntScan & 31% for GeneSplicer).

The precision-recall curve shows that the best overall performance was achieved with SpliceAI, with a PR-AUC=0.54 (Figure 1E). Two examples of cryptic splice mutations affecting liver cancer driver genes uncovered by SpliceAI are shown in Figure 1. These include a TP53 mutation initially annotated as synonymous that actually disrupts the adjacent splice donor site, leading to various degrees of intron retention (Figure 1F). Another interesting example is a missense mutation in the tumor suppressor gene ARID1A that creates an acceptor gain 3 bases downstream. RNA-seq reads show that this cryptic acceptor is used in 12% of the transcripts, leading to a dual consequence of the mutation: missense or frameshift depending on the transcripts (Figure 1G). Altogether, these data indicate that SpliceAI detects splice-impacting mutations with better precision and recall than previous tools, and can reveal new cryptic splice mutations in cancer driver genes.

Pan-cancer analysis of cryptic splice mutations

We next applied SpliceAI to uncover cryptic splice mutations in a large compendium of 17,714 tumor samples across 25 cancer types from the ICGC-PCAWG (n=2,699 WGS, mostly primary tumors (28)), Hartwig Medical Foundation (n=4,818 WGS of metastatic cancers [START_REF] Priestley | Pan-cancer whole-genome analyses of metastatic solid tumours[END_REF]) and TCGA (n=10,197 WES, mostly primary tumors (30)) projects (Figure 2A). To analyze the 277,415,611 somatic mutations, we used the Jean Zay supercomputer V100 quadri GPU partition, which took 0.003 second per variant to compute a prediction. We identified 203,889 mutations predicted to impact splicing with a DS ≥ 0.5, with a predominance of donor/acceptor losses (71.5%) over donor/acceptor gains (28.5%, Figure 2B). Of these, 43.5% corresponded to essential splice sites whereas 56.5% were cryptic splice mutations located in extended splice regions (13%), introns (14%) or exons (26%, Figure 2C,D). We used matched RNA-seq to validate 3,000 randomly selected mutations predicted to impact splicing in the TCGA series, and we obtained a validation rate of 67%, consistent with the results in our in-house series (Supplementary Table S2). In addition to the two first intronic bases corresponding to essential splice sites, donor loss mutations frequently involved the last exonic base (21%) and the 5th intronic base (9%), indicating an important role of these positions in regulating the recognition of donor sites (Figure 2E). Interestingly, donor gain mutations also frequently involved the same positions (0 and +5, 9%, Figure 2F). Thus, by altering these bases that are important for the regulation of the native donor sites, these mutations favor the usage of alternative cryptic donors. Cryptic acceptor loss mutations were mostly located between the 3rd and 15th intronic bases (17%), corresponding to the polypyrimidine tract. Acceptor gains were quite evenly distributed between positions -15 and +3 relative to the closest essential acceptor site.

Cryptic splice mutations enhance the landscape of driver genes and mutations

We next investigated whether considering cryptic splice mutations would uncover new cancer driver genes and new damaging mutations in known cancer genes. To address this point, we used MutSigCV (31) to identify genes with significantly recurrent mutations in the HMF and TCGA series, considering classical VEP annotations or adding cryptic splice mutations. MutSigCV identified a median of 13.5 significantly mutated genes (q < 0.05) per series (range 1-121) with classical annotations. Including cryptic splice mutations added a total of 126 significant genes (median=3.5 per series, range=1-16, Figure 3A, and Supplementary Figure S3). Extra significant genes were enriched in known cancer genes (n=32, 25%, P=1.1e-18, Supplementary Table S3), validating the biological relevance of the new findings. For example, ARID1A, a known driver gene in cholangiocarcinoma, was not significant in the relatively small TCGA_CHOL series (n=36) but became significant when adding cryptic splice mutations. This was also the case for SMAD4 in pancreatic cancer, TBX3 in colorectal cancer, or TSC2 in liver cancer. Adding cryptic splice mutations also revealed extra-significant genes known to be drivers in other cancer types. For example, we identified in the HMF breast metastasis series significantly recurrent mutations of CDKN1A (known driver in bladder cancer), MAX (known in paragangliomas, endometrioid and colon carcinomas) and CDK12 (known in ovarian cancer). Finally, 94 of the genes that became significant when adding cryptic splice mutations were not in the Cancer Gene Census and may correspond to new driver genes (Supplementary Table S4). These included GPS2, negative regulator of RAS-and MAPK-mediated signal [START_REF] Spain | Two human cDNAs, including a homolog of Arabidopsis FUS6 (COP11), suppress G-proteinand mitogen-activated protein kinase-mediated signal transduction in yeast and mammalian cells[END_REF], recurrently mutated in prostate cancer; BRD7, member of the SWI/SNF chromatin remodeling complex, recurrently altered in urothelial cancer; GPAM, associated with steatosis and liver damage [START_REF] Jamialahmadi | Exome-Wide Association Study on Alanine Aminotransferase Identifies Sequence Variants in the GPAM and APOE Associated With Fatty Liver Disease[END_REF], recurrently mutated in liver cancer or FANCM, known susceptibility gene for breast cancer [START_REF] Kiiski | Exome sequencing identifies FANCM as a susceptibility gene for triple-negative breast cancer[END_REF] but identified here as a somatic driver.

In addition to finding new driver genes, cryptic splice mutations may reveal additional driver mutations in known cancer genes, incorrectly annotated as harmless variants otherwise (e.g. intronic or synonymous). To quantify this, we estimated the mutation frequency of the 695 genes from the Cancer Gene Census present in our large pan-cancer series, considering only classical coding mutations (missense, nonsense, essential splice, frameshift and in-frame indels) or adding cryptic splice mutations (Supplementary Table S5). On average, cryptic splice mutations increased splice mutation frequencies of driver genes by 100%, and particularly in tumor suppressor genes (Figure 3C). For example, TP53 mutation frequency increased from 3.4% to 4.7% when considering cryptic splice mutations. In total, we identified 7,043 cryptic splice mutations in cancer genes versus 7,475 essential splice mutations, thus doubling the prevalence of splice mutations in those genes. Considering cryptic splice mutations thus improves the power of driver gene discovery, and increases the mutation frequencies of driver genes by rescuing variants misclassified as harmless.

Mutational processes at the origin of splice mutations

We next investigated the mutational processes at the origin of splice mutations. To that aim, we extracted the contribution of COSMIC mutational signatures [START_REF] Alexandrov | The repertoire of mutational signatures in human cancer[END_REF] to each of the 2,699 tumors from the ICGC-PCAWG project. We then used Palimpsest [START_REF] Shinde | Palimpsest: an R package for studying mutational and structural variant signatures along clonal evolution in cancer[END_REF] to estimate the probability of each splice mutation being due to each mutational process. By summing these probabilities over splice mutations (essential splice + cryptic splice mutations identified by SpliceAI) and other mutations, we could estimate the contribution of signatures to both mutation sets in each cancer type (Figure 4A). The contribution of signatures to splice and other mutations were globally consistent and highly dependent on the tumor type. For example, splice mutations were mostly attributed to single base substitution signature SBS4 (tobacco-related) in lung, wheras they were mostly attributed to SBS7a/b (UV-related) in melanoma.

Yet, we identified signatures with significantly higher probabilities of giving rise to splice than other mutations (Supplementary Table S6). Interestingly, clock-like signatures SBS1 (due to the spontaneous deamination of mCpGs) and SBS5 (unknown etiology, age-related) were significantly enriched in splice mutations of several cancer types (P < 6e-3, Wilcoxon rank-sum test), including esophageal, stomach and colorectal cancers (Figure 4B). Several signatures associated with DNA repair defects were also enriched in splice mutations (Figure 4C). This was the case for the mismatch repair deficiency signatures SBS26 (uterus) and SBS44 (stomach, colorectal), as well as the homologous repair deficiency signature SBS3 (breast, ovary). This may indicate that the preferential sequence contexts of these mutational processes overlap with splicing regulatory sequences. In contrast, other signatures were significantly less prevalent in splice mutations, including SBS17b signature (associated with fluoropyrimidine-based chemotherapy) in esophageal (P=4.6e-32) and stomach cancer (P=1.4e-10). Altogether, these data indicate that the mutational processes at the origin of splice mutations are cancer-dependent, and that some processes (clock-like, DNA repair deficiencies) may be more prone to generating splice mutations.

DISCUSSION

Our comparison of SpliceAI, GeneSplicer and MaxEntScan in a cancer consistent with a previous benchmark in a context of germline variants [START_REF] Jaganathan | Predicting Splicing from Primary Sequence with Deep Learning[END_REF], showing a higher sensitivity and specificity of SpliceAI. Of note, SpliceAI includes long sequences as input (10 kb around the mutation), which allows the tool to account for long-range sequence determinants. This may explain at least in part the higher performance of SpliceAI. In addition, SpliceAI returns predictions for mutations that can be located deeply in the introns, while GeneSplicer and MaxEntScan are restricted to mutations close to the essential splice site. Another approach, MiSplice, was recently proposed and applied to large pancancer data sets [START_REF] Cao | Discovery of driver non-coding splice-site-creating mutations in cancer[END_REF]. While this approach unraveled and validated many cryptic splice sites, it requires as input both DNA and RNA sequencing data. By contrast, SpliceAI requires only the DNA sequence and can thus be used to annotate any variant list identified by WES/WGS. SpliceAI identified 115,197 cryptic splice mutations in our large pan-cancer cohort, 34% are synonymous or intronic variants hence classified as low impact variants by classical variant annotation tools. We showed that rescuing these cryptic splice mutations increases the power of driver gene discovery with MutSigCV, and doubles the frequency of splice alterations in known tumor suppressor genes. Cryptic splice variants may explain part of the 5% of cancer cases without any driver alteration identified [START_REF] Campbell | Pan-cancer analysis of whole genomes[END_REF]. In addition, 28% of cryptic splice mutations are annotated as missense variants. These mutations may have dual consequences, leading to missense and/or splice alteration depending on the RNA molecule, exemplified by the ARID1A mutation displayed in Figure 1G.

Splice mutations may occur in different sequence contexts and can be generated by distinct mutational processes depending on the cancer type. Yet, we identified signatures enriched in splice mutations, related to clock-like processes, mismatch repair and homologous repair deficiencies. These associations may result from an overrepresentation of the sequence contexts in which these processes are particularly active in splice regulatory elements, which needs to be validated with experimental data.

The contribution of DNA repair signatures to splice mutations may lead to a higher risk of acquiring splice mutations in patients carrying defects in DNA repair enzymes. This study highlights the power of deep learning to uncover previously overlooked regulatory mutations in cancer sequencing data, and emphasizes the importance of better annotating cryptic splice mutations in research series and clinical sequencing data.

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2. Pan-cancer analysis of cryptic splice mutations with SpliceAI. (A) Number, origin and cancer type of samples analyzed. (B) Proportion of splice altering mutations identified in each series classified by type of alteration (DL: donor loss, AL: acceptor loss, DG: donor gain, AG: acceptor gain). (C) Proportion of initial mutation categories for mutations predicted to impact splicing by SpliceAI in each cohort (DS > 0.5). (D) Overall proportion of initial mutation categories for the 203,889 mutations predicted to impact splicing with a DS > 0.5 in the whole data set. (E) Position of donor and acceptor

  by SpliceAI relative to the closest essential splice acceptor and donor. (F) Position of donor and acceptor gain mutations identified by SpliceAI relative to the closest essential splice acceptor and donor.

Figure 3 .

 3 Figure 3. Cryptic splice mutations enhance the landscape of driver genes and mutations. (A) Number of significant genes identified by MutSigCV algorithm when considering only classical coding mutations (in blue) or when adding cryptic splice mutations revealed by SpliceAI (in yellow). (B) Comparison of MutSigCV q-values when considering only classical coding mutations or when adding cryptic splice mutations. Four series are shown as examples. Results for all series are shown in Supplementary Figure S2. (C) Frequency of splice mutations of Cancer Gene Census genes considering only essential or both essential and cryptic splice mutations. Tumor suppressor genes and oncogenes are distinguished with a color code.

FUNDING

This work was performed using HPC resources from GENCI-IDRIS (Grant 2021-AD011011412R1). This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute); EL supported by the Région Pays de la Loire (Connect Talent program); Nantes Métropole and the SIRIC ILIAD (INCa-DGOS-Inserm-12558).

AVAILABILITY

LICA-FR data used in this manuscript are hosted in the European Genome-Phenome Archive (https://ega-archive.org). We downloaded public cancer genomic data from the ICGC-PCAWG (https://dcc.icgc.org/pcawg), TCGA (https://gdc.cancer.gov/about-data/publications/pancanatlas) and HMF (https://www.hartwigmedicalfoundation.nl/en/data/database/) websites. The pyRNA library that we developed to mine RNA-seq data in Python is available on GitHub (https://github.com/FunGeST/pyRNA) Zenodo (DOI: 10.5281/zenodo.7660932).

SUPPLEMENTARY DATA

Supplementary data available at NAR online.

ACKNOWLEDGEMENT

We thank members of the thesis committee of Raphaël Teboul (Tatiana Popova, Andrei Zinovyev, Josh Waterfall and Stéphane Minvielle) for fruitful discussions and advice.

CONFLICT OF INTEREST

This authors declare that they have no conflict of interest.