Benjamin Lemasson

Thomas Christen

Xavier Tizon

Régine Farion

Nadège Fondraz

Peggy Provent

Christoph Segebarth

Emmanuel L Barbier

Philippe Genne

Olivier Duchamp

Msc Lemasson Benjamin

MSc Christen Thomas

PhD Tizon Xavier

Fondraz Nadège

PhD Provent Peggy

PhD Segebarth Christoph

PhD Genne Philippe

Duchamp Olivier

PhD Remy Chantal
email: chantal.remy@ujf-grenoble.fr

Assessment of multiparametric MRI in a

Keywords: Apparent diffusion coefficient of water (ADC), blood volume fraction (BVf), brain tumor, Contrast enhancement MRI, therapies, vessel size index (VSI) ADC, Apparent Diffusion Coefficient of water, BVf, Blood Volume fraction, CE, Contrast Enhancement, VSI, Vessel Size Index, RecA-1, Rat Endothelial Cell Antigen-1

Introduction

Glioblastomas are the most common subtypes of rapidly growing primary brain tumors in adults. Being among the most angiogenic human tumors, they are characterized by remarkable proliferative vascular components and are often necrotic and invasive [START_REF] Jain | Angiogenesis in brain tumours[END_REF]. Their aggressiveness is mainly due to their ability to stimulate the formation of new blood vessels [START_REF] Folkman | Role of angiogenesis in tumor growth and metastasis[END_REF].

Despite active efforts in therapeutics such as surgery, radiotherapy and chemotherapy, gliomas patients still have a poor prognosis and a high recurrence rate. Since 1980, the median survival of patients with glioblastoma has not increased and remains around 11 months after diagnosis [START_REF] Burnet | High grade glioma: imaging combined with pathological grade defines management and predicts prognosis[END_REF]. The lack of efficiency in conventional therapeutic approaches is related to drug delivery challenges inherent to chemoresistance. New therapeutic strategies directed against tumor vasculature or preventing angiogenesis have been developed with exciting preliminary results [START_REF] Tabatabai | Primetime for antiangiogenic therapy[END_REF][START_REF] Dietrich | Emerging antiangiogenic treatments for gliomas -efficacy and safety issues[END_REF]. Antiangiogenic drugs are likely to play a key role in the treatment of malignant glioma, mostly in association with other molecularly targeted compounds or with cytotoxic molecules [START_REF] Burnet | High grade glioma: imaging combined with pathological grade defines management and predicts prognosis[END_REF][START_REF] Nghiemphu | Bevacizumab and chemotherapy for recurrent glioblastoma: a single-institution experience[END_REF].

In the clinic, the current standard for assessing the efficacy of an anti-tumor therapy is the "Response Evaluation Criteria in Solid Tumor" (RECIST) which is mainly based on measuring tumor volume from radiographic images [START_REF] Husband | Evaluation of the response to treatment of solid tumours -a consensus statement of the International Cancer Imaging Society[END_REF]. Given the observed delay between treatment onset and its effect on tumor size (several weeks or months), an early biomarker of response to therapy would allow rapid treatment adaptation or guide treatment combination, on an individual basis. These biomarkers can arise from biochemical assays (e.g. tissular or circulant markers) or from imaging modalities [START_REF] Gagner | Angiogenesis in gliomas: imaging and experimental therapeutics[END_REF][START_REF] Jain | Biomarkers of response and resistance to antiangiogenic therapy[END_REF]. Among all imaging modalities, magnetic resonance imaging (MRI) appears as a good candidate for monitoring the effect of anti-tumor therapy. MRI is routinely used to estimate tumor volume using anatomical images, but also tumor associated edema and tumor cellularity, using the apparent diffusion coefficient (ADC) of water [START_REF] Hamstra | Evaluation of the functional diffusion map as an early biomarker of time-to-progression and overall survival in high-grade glioma[END_REF][START_REF] Ross | Evaluation of cancer therapy using diffusion magnetic resonance imaging[END_REF][START_REF] Cha | Neuroimaging in neuro-oncology[END_REF]. Moreover, numerous microvascular parameters that appear well-suited to assess the effect of antiangiogenic drugs can be estimated using MRI, such as blood volume, blood flow, microvessel diameter, or vessel wall permeability [START_REF] Barrett | MRI of tumor angiogenesis[END_REF][START_REF] Pike | Highresolution longitudinal assessment of flow and permeability in mouse glioma vasculature: Sequential small molecule and SPIO dynamic contrast agent MRI[END_REF][START_REF] Beaumont | Characterization of tumor angiogenesis in rat brain using iron-based vessel size index MRI in combination with gadolinium-based dynamic contrastenhanced MRI[END_REF]. A recent example shows that the extent of vascular normalization measured by MRI would be predictive of the outcome of anti-Vascular Endothelial Growth Factor (anti-VEGF) therapy in glioblastomas [START_REF] Sorensen | A "vascular normalization index" as potential mechanistic biomarker to predict survival after a single dose of cediranib in recurrent glioblastoma patients[END_REF].

Several studies have monitored one or two of these parameters to follow the effect of an anti-tumor therapy in addition to RECIST [START_REF] Flaherty | Pilot study of DCE-MRI to predict progression-free survival with sorafenib therapy in renal cell carcinoma[END_REF][START_REF] Zwick | Assessment of vascular remodeling under antiangiogenic therapy using DCE-MRI and vessel size imaging[END_REF][START_REF] Varallyay | Dynamic MRI using iron oxide nanoparticles to assess early vascular effects of antiangiogenic versus corticosteroid treatment in a glioma model[END_REF]. Very few therapeutic follow-ups have used an MRI protocol combining more than three parameters [START_REF] Batchelor | AZD2171, a pan-VEGF receptor tyrosine kinase inhibitor, normalizes tumor vasculature and alleviates edema in glioblastoma patients[END_REF].

The aim of this study was to analyze various MR parameters into one protocol to assess the impact of 2 different treatments on glioma model and to evaluate whether these parameters are relevant, i.e. sensitive to a given treatment.

To address this issue, we used a combined MRI / histology protocol to monitor the effect of a cytotoxic therapy (Carmustin / BCNU; an alkylating agent) and of an antiangiogenic therapy (Sorafenib; Nexavar ® ; a multikinase inhibitor) a glioma model obtained by orthotopically xenografting a human glioblastoma cell line (U-87 MG) in Nude rats. BCNU has been chosen because it is one of the standards of care for gliomas. Sorafenib is considered as an antiangiogenic molecule and has been found efficient for several solid tumors [START_REF] Wilhelm | Discovery and development of sorafenib: a multikinase inhibitor for treating cancer[END_REF]. Regarding brain tumors, it was shown that Sorafenib induces growth arrest on a medulloblastoma [START_REF] Yang | Sorafenib inhibits signal transducer and activator of transcription 3 signaling associated with growth arrest and apoptosis of medulloblastomas[END_REF] and a human glioblatoma [START_REF] Siegelin | Sorafenib exerts anti-glioma activity in vitro and in vivo[END_REF] xenografted grafted in mice. Vascular effects were detected in the last case [START_REF] Siegelin | Sorafenib exerts anti-glioma activity in vitro and in vivo[END_REF].

Material and methods

Glioma cell line and anti-tumor drugs

The U-87 MG human glioblastoma cell line was purchased from the American Type Culture Collection (Manassas, VA). Cells were cultured in RPMI1640 medium containing 10% fetal calf serum and 2 mM glutamine.

Sorafenib (Nexavar ® ) was purchased from Bayer Corporation (West Haven, CT, U.S.A) and prepared as a suspension in a mixture of 5% dimethyl sulfoxide (DMSO), 5% tween-20 and 90% NaCl 0.9%.

Carmustine (1,3-bis(2-chloroethyl)-1-nitroso-urea, BCNU) was purchased from Sigma-Aldrich (St. Louis, MO, USA) and dissolved in NaCl 0.9% containing 30% of ethanol.

In vivo intracerebral U-87 MG glioma model and treatments

This study was approved by the local committee for animal care and use (approval 081). Experiments were performed under permits n° A380521, A3851610004 and B3851610003 from the French Ministry of Agriculture.

Experimental design of this study is presented in figure 1a.

Male RH-rnu/rnu Nude rats (n=60, weight: 180-220g, Harlan Sprague Dawley, Indianapolis, USA) were prepared for xenograft by whole body irradiation with a 60 Co !-source (7 Gy), 24-48 h prior to human tumor cell inoculation. Animals were anesthetized with an intraperitoneal injection of Ketamine (75 mg/kg, Ketamine500 ® , Centravet, France) and Xylazine (5 mg/kg, Rompun ® , Centravet, France) in 0.9% NaCl solution and immobilized in a stereotactic frame (David Kopf Instrument, Germany). Five !l of U-87 MG cell suspension in serum-free RPMI1640 medium containing 10 5 tumor cells were inoculated in the right caudate nucleus of the rat brain by stereotactic injection through a 1-mm burr hole at 2.5 mm lateral to the bregma and at a depth of 4.5 mm from the dura. The injection was performed slowly over 15 minutes, and the needle was withdrawn over another 5 minutes. The burr hole was filled with bone wax to prevent extracranial extension of the tumor.

Eleven days after tumor implantation (D11), anatomical imaging (T2-weighted anatomical images) was performed to evaluate the volume of each tumor. Rats were then randomized and stratified in 3 groups (n=20 per group) with the similar average tumor volume per group (mean tumor volume across all rats at D11, 3.8 ± 2.2 mm 3 ).

Treatment began 14 days after tumor implantation (Fig. 1a; D14 (T0) ):

-BCNU group received two intravenous injections of BCNU (10 mg.kg -1 ) on the day of treatment onset and 13 days after (D14 (T0) and D27 (T13) ).

-SORA group received daily oral administration of Sorafenib (100 mg.kg -1 ) between treatment onset and the 14 th day after treatment start (D14 (T0) to D28 (T14) ).

-Untreated group received no treatment.

In vivo MRI experiments for tumor follow-up

Experiments were performed on a horizontal 2.35 Tesla (40 cm diameter) magnet equipped with actively shielded gradient coils (Magnex Scientific Ltd., Oxford, UK) and interfaced to a SMIS console (SMIS Ltd, Guildford, UK). MR acquisitions were performed under anesthesia: 5% isoflurane for induction and 2%

for maintenance in 70% air / 30% oxygen. Rectal temperature was maintained at 37.0 ± 0.5°C throughout the experiments. After anesthesia, a catheter was inserted into the tail vein allowing injection of contrast agent without removing the rat from the magnet.

Tumor volume, ADC, BVf, VSI and CE were mapped in the tumor and contralateral hemisphere one day before and 1, 4 and 14 days after treatment onset (Fig. 1; D13 (T-1) , D15 (T1) , D18 (T4) and D28 (T14) , respectively). The last time point was chosen to correspond to 3 days before the median survival time of untreated Nude rats bearing U-87 MG glioblastoma (31 days; based on the survival study -see Fig. 2a). The total duration of one MRI session was 1h30 per animal. MRI sequences were realized in the following order (Fig. 1b):

-Anatomical T2-weighted images were acquired using a spin-echo MRI sequence (TR/TE = 2000/80 ms, 19 slices with FOV = 30x30 mm 2 , matrix = 256x256 and voxel size = 117x117x1000 µm 3 ). Imaging time was 4 min and 17 s.

-ADC was computed from 3 diffusion weighted spin-echo images (in X, Y, and Z directions) with b = 900 s.mm -2 and a reference image (b " 0 s.mm -2 ) (TR/TE = 2000/80 ms, 7 slices with FOV = 30x30 mm 2 , matrix = 128x128 and voxel size = 234x234x1000 µm 3 ). Imaging time was 17 min.

-BVf and VSI were estimated using a multi gradient echo and spin-echo MRI sequence (MGESE) (TR = 6 s, 7 evenly spaced gradient-echoes = for r 1 and r 2 respectively (data provided by Guerbet).

-CE was assessed using T1-weighted anatomical images acquired using a 3D-MDEFT sequence (TI = 605 ms, alpha = 22°, TR/TE = 15/5 ms, FOV 40x40x30 mm 3 , matrix = 128x128x96 and voxel size = 313x313x313 µm 3 ). Imaging time was 4 min and 15 s. To limit the duration of the MR session per animal, a simple subtraction approach was preferred to a dynamic contrast enhanced MRI approach.

Thus, one 3D-MDEFT sequence was performed before and 3 min after administration of P846 via the catheter in about 20 sec (50 µmoles/kg body weight; flush with 200µL of saline). P846 (Guerbet, Roissy, France) is a Gd-based contrast agent with a molecular weight of 3.5 kDa and a plasmatic half-life in rats of 51.1 ± 1.1 min at the concentration used (data provided by Guerbet). The relaxivities of P846, measured in water at 2.35T and 37°C are 28 and 39 s -1 .mM -1 for r 1 and r 2 respectively (data provided by Guerbet). In this study we used high molecular weight contrast agent (P846: 3.5 kD), instead of a routinely used clinical agent (Gd-DOTA or Gd-DTPA: around 0.5 kD), to increase CE sensitivity to different vessel maturity [START_REF] Jacquier | MR imaging assessment of the kinetics of P846, a new gadolinium-based MR contrast medium, in ischemically injured myocardium[END_REF].

For each group, 4 rats, called the "longitudinal" subgroup were imaged at each time point (D13 (T-1) , D15 (T1) , D18 (T4) and D28 (T14) ). After the last MRI measurement, animals were euthanized for subsequent ex vivo experiments. Per group, the sixteen remaining rats were imaged once and then euthanized at the end of imaging session to enable ex-vivo studies (4 animals/time point/group). They form the "single-time" subgroup (Fig. 1a). Due to experimental and technical problems (death of animals due to anesthesia, data unusable because of problems during acquisitions) only 7 out of 8 rats per group yielded acceptable data at D13 (T-1) , D15 (T1) , and D18 (T4) and only 4 out of 8 rats per group at D28 (T14) . Details of data usable per group (untreated, BCNU and SORA) and subgroup (single-time and longitudinal) are presented in table 1.

Ex vivo experiments

At the end of the MRI experiment, the animals were humanely sacrified (Fig. 1a). The brain was quickly removed, frozen in -40°C isopentane and stored at -80°C until processing. Brains were sliced at -20°C with a cryostat (10 µm thick sections).

First, anatomic changes associated with BCNU and Sorafenib treatment were analyzed using hematoxylin and eosin (HE) staining.

Secondly, vascular components were analyzed by immunohistochemistry of the vessel basal lamina with collagen IV labeling and of endothelial cells using rat endothelial cell antigen-1 labeling (RecA-1). These two stainings assess the impact of anti-tumor treatments on vessel integrity. Slices were rehydrated in phosphatebuffered saline (PBS, 0.01M) and fixed in 4% paraformaldehyde. After saturation in PBS-Tween 0.01%-BSA 3% for 30 min at room temperature, primary antibodies (a

goat

Survival Study

The effects of both treatments on the mean survival of Nude rats bearing a human glioblastoma xenografts were assessed on 6 additional rats per group (untreated, BCNU, and SORA groups) implanted and treated as described above.

Animals were sacrificed at the onset of neurological signs (according to Redgate et al. ( 25)). Increased lifespan (ILS) was estimated for each treated group compared (in percent). ILS was calculated using the mean survival time (MST) of each group as following: ILS (%) = (MST untreated -MST treated )/ MST untreated

Control MRI experiment: evaluation of USPIO extravasation

The MR method used to estimate BVf and VSI uses as a starting hypothesis that the contrast agent remains intravascular [START_REF] Valable | Assessment of blood volume, vessel size, and the expression of angiogenic factors in two rat glioma models: a longitudinal in vivo and ex vivo study[END_REF]. To estimate the amount of USPIO that extravasates during BVf and VSI measurements, dynamic contrast-enhanced MRI (DCE-MRI: T 1 -based) was performed on another set of tumor bearing rats (n=12, in addition to those already mentioned). For this 'control' experiment, we used a dedicated MRI protocol (Fig. 1c). At each time point, the 4 animals of each group Since we neglected the plasmatic clearance of USPIO in our analysis, we overestimated the T 1 effects and thus obtained an upper limit of the extravascular iron concentration.

Quantitative histology: Sections corresponding to MRI images (up to three microscopic fields per ROI and animal) were digitized using a CCD camera (Olympus, Rungis, France). Collagen IV images were thresholded and vascular parameters (mean vessel density, fractional vascular surface area, mean vessel radius) were obtained using ImageJ software (Rasband,W.S., ImageJ) as described by Valable et al. [START_REF] Valable | Assessment of blood volume, vessel size, and the expression of angiogenic factors in two rat glioma models: a longitudinal in vivo and ex vivo study[END_REF]. The number of vessels stained by Collagen IV and RecA has been account manually onto the immunohistologies images and then the percentage of vessels doubly stained was calculated. We also realized a manual quantification of the cell density from the HE sections (cells/mm 2 ).

Statistical analysis

Paired t tests were used for comparing tumor versus contralateral striatum inside the same group. Comparisons between each group were subjected to factorial analysis of variance with a Protected Least Significant Difference (PLSD) Fisher's test as post hoc test (intergroup analysis). p<0.05: *; p<0.01: ** and p<0.001: ***.

Survival was visualized using Kaplan-Meier curves and differences between groups assessed with a log-rank test. All statistics were performed with SPSS (SPSS Inc, Chicago, Ill, USA). Results are presented as mean ± standard deviation and are given by group, by time point and by region of interest.

Results

Accuracy of the BVf and VSI measurements

Longitudinal MRI follow-up

In a longitudinal study, repeated anesthesia and contrast agent injections might affect the physiological status of the animal, the rate of tumor growth and the accuracy of the MRI measurements. To check this, we compared for each group the "single-time" and the "longitudinal" subgroups at each time point and we observed no difference in body weight, tumor size, ADC, BVf, VSI and in CE between. In this experiment, a longitudinal study (repetition of MRI protocol) has no effect on the assessed parameters. Consequently, for the 3 groups, values from the "single-time"

and "longitudinal" subgroups were pooled at each time-point.

In vivo anti-tumor activity of BCNU and Sorafenib treatments

MRI data in contralateral rat brain hemisphere

No significant difference over time and across the 3 groups in ADC, BVf and VSI was detected contralaterally (720±85 !m$.s -1 , 3.3±0.5% and 5.6±1.1!m for ADC, BVf and VSI mean across time points and groups respectively). Consequently, the mean contralateral value across groups is represented at each time point, for the sake of clarity (Fig. 3; dotted line). However, intragroup comparisons between tumor and contralateral tissues were performed using the contralateral values of each group.

U-87 MG tumor development (untreated group)

Median survival time of untreated U-87 MG glioblastoma cells bearing Nude rats was 31 days (Fig. 2a). The volume of the untreated U-87 MG tumors increased rapidly between D13 (T-1) and D28 (T14) (4.5 ±1.4 to 117.1 ± 22.9 mm 3 respectively; Fig. 2b). ADC, BVf and VSI in U-87 MG tumors were significantly higher than in the contralateral striatum and remained stable with time (tumor: 851 ± 43 !m$.s -1 , 4.4 ± 0.3 % and 7.9 ± 1.2 !m; contralateral striatum: 720 ± 85 !m$.s -1 , 3.3 ± 0.5 % and 5.5 ± 0.3 !m for ADC, BVf and VSI respectively, mean value across the four time points, p<0.001; Fig. 3a-c). In the untreated group, glioma vessels were permeable to P846 during follow-up (CE: 190 ± 47 %, mean across all time points; Fig. 3d).

Effect of BCNU on U-87 MG glioma (BCNU group)

BCNU treatment slightly increased (but not significantly) the median survival of U-87 MG glioma bearing rats (+16%; Fig. 2a). BCNU treatment strongly inhibited tumor growth compared to untreated group (5.7 ± 4.9 vs. 117.1 ± 22.9 mm 3 , for the BCNU and untreated groups respectively, at D28 (T14) , p<0.001; Fig. 2b). Intratumoral ADC was comparable between the BCNU and untreated groups at all time points except at D28 (T14) when ADC in the BCNU group became significantly larger than in the untreated group (1048 ± 10 vs. 794 ± 49 !m$.s -1 , p<0.001; Fig. 3a). VSI, BVf and CE in the BCNU group did not differ from those measured in the untreated group (Figs. 3b-d).

Effect of Sorafenib on U-87 MG glioma (SORA group)

Sorafenib treatment induced a small increase (but not significantly) of the median survival time of the U87-MG glioma bearing rats as compared to the untreated group (+23%; Fig. 2a). Sorafenib treatment significantly inhibited tumor growth as compared to untreated group (28.7 ± 11.1 vs. 117.1 ± 22.9 mm 3 , respectively, at D28 (T14) , p<0.001; Fig. 2b). We also observe, at D28 (T14) , that the tumoral volume into the SORA group were statistically greater than that of the BCNU group (28.7 ± 11.1 vs. 5.7 ± 4.9 mm 3 , respectively, at D28 (T14) , p<0.001; Fig. 2b). At D28 (T14) , ADC in the SORA group was larger than in the untreated group (976 ± 37 and 794 ± 49 !m$.s -1 , respectively, p<0.01; Fig. 3a). while it remained stable in the untreated group (D15 (T1) : 5.2 ± 0.9 to D28 (T14) :

All vascular parameters

5.5 ± 0.8 %; Fig. 3c). In contrast to what happens in untreated and BCNU group, in the SORA group the extravasation of P846 was significantly reduced as early as 1 day after Sorafenib treatment start (CE: 116 ± 65% vs 184 ± 42% for the SORA and untreated groups respectively at D15 (T1) ) and was not detectable at D28 (T14) (Fig. 3d).

HE staining

From D13 (T-1) to D18 (T4) , U-87 MG tumors in the untreated, BCNU and SORA groups presented very similar anatomical histological structures on HE stained sections. Cell density was higher in the tumor core with some tumor cells aligned along vessels, forming thin bundles. Tumors were surrounded by a ring of highly edematous tissue. Neither necrosis nor pseudo-cyst was observed (Fig. 4a). At D28 (T14) , in the SORA and BCNU groups tumor cell density was reduced compared to previous time points and edematous areas appeared inside the tumor tissue (Fig. 4a).

Vascular staining

To assess variation in vasculature properties, collagen IV and RecA labelings were performed (Fig. 4). For all groups and all time points, vessels in contralateral striatum exhibited similar diameter, density and vascular surface area (4.4 ± 0.2 !m, 337 ± 37 vessels per mm 2 and 3.76 ± 0.34 % vascular surface area; mean across all groups and all time points; Fig. 4b-d). All vessels were also stained by Collagen IV and RecA (Fig. 4a). Before treatments, in all groups, tumor vessel diameters were larger than in the contralateral striatum (5.2 ± 0.6; 5.4 ± 0.6 and 6 ± 0.3 !m versus 4.2 ± 0.2 !m for untreated, BCNU and SORA groups versus mean of contralateral striatum across all group, respectively; p<0.05; Fig. 4b) and vascular surface area was higher than in contralateral striatum (9.1 ± 1.4; 6.2 ± 1.3 and 7.7 ± 0.3 % versus 3.7 ± 0.5 % for untreated, BCNU and SORA groups versus mean of contralateral striatum across all group, respectively; p<0.05; Fig. 4c). Vascular parameters were stable with time in the untreated and BCNU groups (Fig. 4b-d).

At D28 (T14) , in the SORA group, tumor vessel density was lower compared to untreated group (141 ± 21 versus 301 ± 60 vessels per mm$ respectively; p<0.01) and vascular area was similar to that in contralateral striatum (3.7 ± 1.4 versus 3.7 ± 0.3 %, respectively; p = 0.45; Fig. 4c), but vessel diameter was higher than in contralateral striatum (5.8 ± 0.7 versus 4.2 ± 0.1 vessels per mm$, respectively; p<0.05; Fig. 4b). Finally, only a fraction (63 %; Fig. 4a) of the tumor vessels stained by Collagen IV antibody was also stained by RecA antibody in SORA group as compared to untreated and BCNU groups (100 %; Fig. 4a).

Discussion

The As a first prerequisite of this study, we determined whether our extensive MR protocol could be repeated to monitor therapy response on brain tumor. First, we compared for each group the "single-time" and the "longitudinal" subgroups. In each "longitudinal" subgroup repeated anesthesia, repeated injection of contrast agent, and treatment (BCNU or Sorafenib) had no detectable effect on body weight (data not shown), tumor volume, or on MR-determined parameters as compared to the "single-time" subgroup. The second prerequisite was to assess whether the USPIO extravasated in this tumor model within the time frame of the MR session. A separate experiment, performed on 4 animals per group followed at each time point (D13 (T-1) , D15 (T1) , D18 (T4) and D28 (T14) ), showed that the extravascular concentration of USPIO in the tumor was at most 0.3% of the plasmatic concentration 20 min after injection.

This result validates the main hypothesis of the VSI and BVf measure, considering the USPIO as an intravascular contrast agent. Moreover, it showed that the complete imaging protocol could fit within approximatively 90 minutes, which allows designing follow-up experiments with a sufficient number of animals per group. All these results indicate that our experimental protocol was appropriate to monitor the effect of antitumoral therapies in our U-87 MG tumor model.

The goal of this study was to assess the sensitivity of multiparametric MRI to monitor the anti-tumor activity of BCNU and Sorafenib.

ADC in tumor was increased by both treatments compared to the untreated group 14 days after treatment onset thus showing the slowest evolution under both treatments. This increase is consistent with the cytotoxic effect of BCNU, and the ADC is also known to be affected by changes in cell density (or cellularity) [START_REF] Koh | Diffusion-weighted MRI in the body: applications and challenges in oncology[END_REF] and by the occurrence of vasogenic edema. HE staining showed that cell density was smaller in the tumors of treated groups than in the untreated group 14 days after treatment onset. HE staining also showed intratumoral edema in treated groups.

Even if the mechanisms involved in ADC changes are complex, these histological findings could explain the increase in intratumoral ADC in BCNU and SORA groups observed 14 days after treatment onset.

Tumor microvasculature properties assessed by both MRI and histology were not changed by the BCNU treatment. Conversely, Sorafenib induced significant modifications in tumor microvasculature. Four days after antiangiogenic treatment onset, the mean intratumoral BVf was markedly reduced compared to the untreated group and was close to its normal value. This result suggests a tendency towards normalization of tumor vasculature: under antiangiogenic treatment, the tumor vasculature becomes more similar to healthy tissue. However, at the same time point (D18 (T4) ), tumoral VSI was increased by Sorafenib, with values higher than that observed in animals of untreated group. The quantitative immunohistochemistry of tumor vessels confirmed these observations on vessel diameters and also indicated that the vessel density had decreased, as found by Siegelin et al. [START_REF] Siegelin | Sorafenib exerts anti-glioma activity in vitro and in vivo[END_REF]. The decrease in vessel density observed in immunohistochemistry images is consistent with the decrease in BVf and the concomitant increase in vessel observed by MRI.

In addition to these morphological changes, 1 day after the start of Sorafenib treatment (D15 (T1) ), immunohistochemistry showed some collagen IV-positive vessels without RecA labeling. These observations indicate that some tumor vessels were formed only with basal lamina. This discrepancy between both labelings was even more pronounced after 4 and 14 days of treatment. This has been previously reported in several studies [START_REF] Mancuso | Rapid vascular regrowth in tumors after reversal of VEGF inhibition[END_REF][START_REF] Inai | Inhibition of vascular endothelial growth factor (VEGF) signaling in cancer causes loss of endothelial fenestrations, regression of tumor vessels, and appearance of basement membrane ghosts[END_REF]. Inai et al. observed collagen IV-positive vessels without endothelial marker in spontaneous pancreatic islet tumors treated with a VEGFRs inhibitor.

We also observed a decrease in tumor contrast enhancement as soon as 1 day after antiangiogenic treatment onset. This reduction was more pronounced 4 days after treatment onset and no P846 contrast agent extravasation was detectable 14 days after treatment onset. BVf and VSI estimates using an intravascular contrast agent (USPIO) injection show perfusion from at least a fraction of tumoral vessels under antiangiogenic treatment. These different BBB permeabilities between contrast agents with different sizes have already been described in the study by Turetschek et al. [START_REF] Turetschek | Tumor microvascular changes in antiangiogenic treatment: assessment by magnetic resonance contrast media of different molecular weights[END_REF]. Together these results could indicate that the absence of contrast enhancement of tumors after 14 days of Sorafenib treatment is due to modifications of vessel permeability and not due to a lack of tumor perfusion. Modifications of vessel integrity observed in our study are consistent with results of Flaherty et al. who observed, by MRI, a decrease in vessel permeability/perfusion (K trans ) after 3-12 weeks of Sorafenib treatment in patients bearing a renal cell carcinoma [START_REF] Flaherty | Pilot study of DCE-MRI to predict progression-free survival with sorafenib therapy in renal cell carcinoma[END_REF].

In our study, BVf and CE parameters indicated a normalization of tumoral vessels under antiangiogenic treatment, which has already been described in several studies [START_REF] Batchelor | AZD2171, a pan-VEGF receptor tyrosine kinase inhibitor, normalizes tumor vasculature and alleviates edema in glioblastoma patients[END_REF][START_REF] Jain | Normalization of tumor vasculature: an emerging concept in antiangiogenic therapy[END_REF]. In addition to these two MR parameters, estimates of VSI indicate that tumor vessels had abnormal size compared to vessels in healthy tissue. A normalization of vascular permeability/perfusion associated to an increase in vessel diameters under antiangiogenic treatment has already been described [START_REF] Nakamura | KRN951, a highly potent inhibitor of vascular endothelial growth factor receptor tyrosine kinases, has antitumor activities and affects functional vascular properties[END_REF]. These results highlight the necessity to measure VSI, BVf and CE parameters together to monitor the different aspects of tumor microvasculature under antiangiogenic

treatment. An important step in future studies will be to demonstrate that all the MR parameters used to evaluate the impact of anti-tumoral treatments are independent of each other. After this demonstration we could relate, using multivariate statistical methods, MR parameters with animal survival or tumor progression as proposed by Sorensen et al. with 2 MR parameters plus one biological measurement associated with patient outcome [START_REF] Sorensen | A "vascular normalization index" as potential mechanistic biomarker to predict survival after a single dose of cediranib in recurrent glioblastoma patients[END_REF].

Glioblastomas are tumors that are highly resistant to many types of treatment.

Monotherapies with antiangiogenic or cytotoxic drugs, like Sorafenib or BCNU respectively, demonstrated moderate efficacy in the U-87 MG glioblatoma. However, despite their low effectiveness on the overall survival of animals, MRI biomarkers reveal that drugs effectively reach their target and eventually induce both morphological and physiological changes for which drugs were designed. This could be helpful for pharmacological treatment optimization or combination.

The multiparametric MRI approach described in this study allowed the measurement of various characteristics of cellular and vascular remodeling that occurs due to therapy. To further complete this characterization, additional MRI developments can be foreseen or refined analysis could be evaluated. In our study we use a CE-MRI method to evaluate vessel wall integrity. However, this parameter depends on blood flow, vessel wall permeability and diffusion of contrast agent into the extravascular compartment. An appropriate data model [START_REF] Tofts | Quantitative analysis of dynamic Gd-DTPA enhancement in breast tumors using a permeability model[END_REF][START_REF] Tofts | Measurement of the blood-brain barrier permeability and leakage space using dynamic MR imaging. 1. Fundamental concepts[END_REF] could be used to quantify the change in vascular permeability and in blood flow. This requires a protocol optimization to allow measurements of VSI, BVf (using USPIO) and permeability/perfusion (using Gd-based contrast agent) in the same MRI session [START_REF] Beaumont | Characterization of tumor angiogenesis in rat brain using iron-based vessel 30 size index MRI in combination with gadolinium-based dynamic contrastenhanced MRI[END_REF]. Blood flow could also be obtained using Dynamic Susceptibility Contrast MRI or Arterial Spin Labeling [START_REF] Barbier | Methodology of brain perfusion imaging[END_REF]. Using a quantitative analysis of T 2 * decay [START_REF] He | Quantitative BOLD: mapping of human cerebral deoxygenated blood volume and oxygen extraction fraction: default state[END_REF], one could estimate the oxygen extraction fraction in the tumor, a possible marker of tumor hypoxia.

In conclusion, this study indicates that numerous changes occur in the tumor microvasculature and cellularity after antiangiogenic or cytotoxic therapy. This study also indicates that the sensitivity of the studied MR parameters varies for a given treatment and is also treatment dependent. Moreover, these changes appear complex and may not be summarized by a single physiological measurement. MRI appears as an imaging modality that can follow, non-invasively, several microvascular and cellular parameters with a strong potential for the evaluation of new drugs on brain tumors. In addition, these parameters could provide methods to better describe and understand the mechanisms of action of antiangiogenic or antivascular therapies, allowing the rational design of optimized combinations with chemotherapy and/or radiotherapy.

  Before and three time points after treatment onset, the following MRI parameters were estimated: tumor volume, water ADC, blood volume fraction (BVf), vessel diameter (vessel size index, VSI), and vessel wall integrity (contrast enhancement, CE). The effect of treatments on each of these MRI parameters was assessed and

  antibody against collagen IV, Southern Biotech, ref 1340-01, 1/2000; a mouse antibody against RecA-1, AbD Serotec ref. MCA9070R, 1/200) were incubated overnight at 4°C in PBS-Tween 0.01%-BSA 1%. Secondary antibodies were Alexa 546-linked Donkey anti-goat IgG (Invitrogen, A11030, 1/200) and Alexa 488 linked Donkey anti-mouse IgG (Invitrogen, A21202, 1/200).

(

  SORA, BCNU and untreated groups) were imaged with multiple T 1w spin-echo images (TR/TE = 600/17ms, 11 slices with FOV = 30x30 mm 2 , matrix 128x128 and voxel size = 234x234x1000 µm 3 , 1 min 17 s/image; 17 repetitions; total acquisition time = 21 min 49s). After one repetition (i.e. 1 min 17s after imaging start), Sinerem ® was administered via the catheter in about 20 sec (200 µmoles of iron/kg body weight; flush with 200 µL of saline).Data analysisDetermination of tumor volume (Fig.1b): Tumor volume was obtained by manually delineating the tumors on anatomical T 2 -weighted MR images from adjacent slices containing the lesions, counting voxels within the tumor boundaries, and scaling with the voxel volume.ADC, BVf and VSI (Fig.1b): Data were processed as described in Valable et al.[START_REF] Valable | Assessment of blood volume, vessel size, and the expression of angiogenic factors in two rat glioma models: a longitudinal in vivo and ex vivo study[END_REF] using in-house developed software in the Matlab ® environment. To avoid partial volume effects that occur in slices located on the edge of the tumor (healthy and tumoral tissue are mixed), regions of interest (ROIs) were manually delineated over the tumor on the 3 T 2 -weighted slices containing the largest tumor area. Contralateral ROIs were defined in the contralateral striatum on the same slices. ROIs were transferred to ADC, BVf and VSI maps. Then, within each ROI and each map, voxels for which the analysis could not be performed were identified (e.g. voxels with nonconverging fit and voxels with values outside the range of validity of the method (ADC>3500 !m 2 .s -1 ; BVf>17%, VSI>50 !m) and excluded from the analysis. The total number of voxels, the number of excluded voxels, and the mean and standard deviation over the 3 analyzed slices for a rat were computed for each ROI and each parameter. Rejected voxels per ROI represented 4.5% for tumor and 0.5% contralateral striatum of the total number of voxels, for all groups. CE: Vessel integrity was assessed using the T 1w images acquired before (T 1w(before) ) and 3 min after (T 1w(after) ) injection of P846 as describe by equation 1. The mean and standard deviation were computed for each group using the same ROI described above. equation 1 SI : Signal Intensity USPIO extravasation (Fig. 1c): Since DCE-MRI images were acquired with a TE = 17 ms, a decrease in signal (due to T 2 * reduction, intravascular USPIO) was observed just after injection of USPIO, after which, in case of extravasation, signal increased (T 1 reduction, extravascular USPIO). The relative increase in signal between the first image acquired after USPIO injection and the last image of the DCE-MRI protocol was assigned to T 1 effect. Eventually, this signal enhancement was converted to extravascular iron concentration as described by Valable et al. (26).

  assessed in this study (VSI, BVf and CE) were modified by Sorafenib. While tumoral VSI in the SORA and untreated groups were similar before treatment (D13 (T-1) ), tumoral VSI in the SORA group became significantly larger than in the untreated group as early as D18 (T4) and up to D28 (T14) (D18 (T4) : 7.2 ± 1.8 vs 5.8 ± 1.8 !m; p<0.001 and D28 (T14) : 12.4 ± 1.7 vs 8.7 ± 1.7 µm; p<0.01, respectively; Fig 3b). In the SORA group, the tumoral BVf decreased from treatment onset (D15 (T1) : 5.0 ± 0.8 to D28 (T14) : 2.6 ± 1 %; Fig. 3c)

  effect of two different treatments (BCNU, an alkylating agent and Sorafenib, a multikinase inhibitor) on an orthotopic human glioblastoma xenografted in Nude rats (U-87 MG cell line) was investigated by MRI and by immunohistology. MRI highlighted an important inhibition of tumor growth induced by both Sorafenib and BCNU treatments. As previously demonstrated, ADC, BVf and VSI values were higher in the tumor tissue than in the contralateral striatum. Sorafenib treatment induced significant differences first on tumor contrast enhancement (D13 (T-1) ), then on BVf and VSI (D18 (T4) ), and finally on ADC and tumor size (D28 (T14) ) compared to control animals. In the BCNU group, only tumor size and ADC were modified 14 days after treatment onset.

LegendsFigure 1

 1 Figure 1. a) Experimental protocol. At D0 60 rats were orthotopically implanted with 10 5 U-87 MG glioma cells. Eleven days later (D11) animals were randomized and

Figure 2 .

 2 Figure 2. a) Percent of animal survival as function of time. Open squares show the

Figure 3 .

 3 Figure 3. In vivo MRI estimated parameters. (a) Evolution of ADC, (b) microvessel

Figure 4 .

 4 Figure 4. Ex-vivo study. a) Histological and immunohistological images at D28 (T14) . Example of hematoxylin/eosine (HE) staining representative of contralateral and

  Using our model, accurate BVf and VSI measurements require that the ironbased contrast agent remains intravascular. At each time point, USPIO extravasation data were not statistically different across all groups (data were thus pooled). In contralateral tissue, iron concentration remained below the detection sensitivity of our DCE-MRI approach (#1!M)[START_REF] Valable | Assessment of blood volume, vessel size, and the expression of angiogenic factors in two rat glioma models: a longitudinal in vivo and ex vivo study[END_REF]. Nineteen minutes after USPIO injection (corresponding to the end of MR data acquisition in BVf/VSI measurement), the average iron concentration in the tumor across all groups and time points was 2.6 ± 0.4 !M. This concentration corresponds to 0.3% of the plasma iron

concentration (3420 ± 970 !M, data reported by Valable et al. (

26

)).

Table 1

 1

	Group	Subgroup	D13 (T-1)	D15 (T1)	D18 (T4) D28 (T14)
	Untreated	Longitudinal	4	4	4	3
		Single-time	3	3	3	1
		Total	7	7	7	4
	BCNU	Longitudinal	4	4	3	2
		Single-time	3	3	4	2
		Total	7	7	7	4
	SORA	Longitudinal	4	4	3	2
		Single-time	3	3	4	2
		Total	7	7	7	4

Table 1 :

 1 Data available for each group (untreated, BCNU and SORA),

Acknowledgements

B.L. benefited from a doctoral grant from 'Agence national de la Recherche et de la technologie'. We thank Guerbet for providing Sinerem ® and P846, the 'Institut National du Cancer', the 'Région Rhône-Alpes' and the 'Cancéropôle Lyon Auvergne Rhône-Alpes' for financial support.


Claudio Ottoni
email: claudio.ottoni@med.kuleuven.be

François-Xavier Ricaut

Nancy Vanderheyden

Nicolas Brucato

Marc Waelkens

Ronny Decorte

Mitochondrial analysis of a Byzantine population reveals the differential impact of multiple historical events in South Anatolia

Keywords: Ancient DNA, mtDNA, Anatolia, Byzantine

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Anatolia, the Asian part of Turkey (Asia Minor) is at the crossroad between the Balkans, the Near East and the Caucasus, and has acted during major population expansions as a bridge between Europe, Africa and Asia. Presence of modern humans in Anatolia is attested since the upper Palaeolithic [START_REF] Kuhn | Paleolithic Archeology in Turkey[END_REF] . This region was also suggested to be a refuge during the Last Glacial Maximum 2 , at the end of which an improving climate allowed the coastal hunter-gatherers to venture inland, where they established temporary camps in search of game and flint [START_REF] Vermeersch | Late Palaeolithic at the Dereköy Karain Cave[END_REF] . Anatolia is also considered as an area through which the 'Neolithic revolution' spread into Europe [START_REF] Ammerman | The Neolithic transition and genetics of populations in Europe[END_REF][START_REF] Renfrew | Archaeogenetics: Towards a population prehistory of Europe[END_REF][START_REF] Richards | Tracing European founder lineages in the Near Eastern mtDNA pool[END_REF][START_REF] Bellwood | Examining the Farming/Language Dispersal Hypothesis[END_REF][START_REF] Zeder | Domestication and early agriculture in the Mediterranean Basin: Origins, diffusion, and impact[END_REF] . However, the major cultural changes in Anatolia occurred from the Bronze Age, when various civilization and empires succeeded one another over time and increasing interactions in the eastern Mediterranean resulted in the development of trading networks, military campaigns and human colonizations, particularly during the Mycenaean, the Greco-Hellenistic, Roman Republican and Imperial, and later the Islamic periods [START_REF] Edwards | The Middle East and the Aegean Region, c.1800-1380 BC[END_REF][START_REF] Sahoglu | The Anatolian trade network and the Izmir region during the early bronze age[END_REF][START_REF] Waelkens | The Late Antique to Early Byzantine city in Southwest Anatolia. Sagalassos and its territory: a case study[END_REF][START_REF] King | Differential Y-chromosome Anatolian influences on the Greek and Cretan Neolithic[END_REF][START_REF] Cruciani | Tracing past human male movements in northern/eastern Africa and western Eurasia: new clues from Y-chromosomal haplogroups E-M78 and J-M12[END_REF] .

From the early Middle Bronze Age onwards city-states and federations of the Hittites and the Luwians, two originally ethnically and linguistically related groups from Transcaucasia, started to infiltrate Central to Eastern and Western to Southern Anatolia respectively [START_REF] Waelkens | Sagalassos and Pisidia during the Late Bronze Age[END_REF][START_REF] Oettinger | Indo-germanische Sprachträger lebten schon im 3. Jahrtausend v. Chr. In Kleinasien. Die Ausbildung der anatolischen Sprachen: Die Hethiter und ihr Reich[END_REF] . Shortly afterward, from the 16 th millennium BC, Minoan and Mycenaean colonization took place on the Anatolian west and southwest coast.

The Early Iron Age in Central and West Anatolia is characterized by the kingdom of the Phrygians (10 th to early 7 th century BC), followed after its collapse by the Lydians.

From 546 until 334/333 BC Anatolia was incorporated into the Achaemenid Persian Empire, and later into that of Alexander the Great. Then, from the early 2 nd century BC onwards Pisidia was part of the realm of the Romans [START_REF] Waelkens | Romanization in the East: A case study. Sagalassos and Pisidia (SW Turkey)[END_REF][START_REF] Waelkens | Ein Blick von der Ferne: Seleukiden und Attaliden in Pisidien[END_REF] . It has been estimated that in the Late Roman Period, the size of the Anatolian population reached 12 million individuals [START_REF] Russell | Late Ancient and Medieval Population[END_REF] .

In the mid of the 6 th century AD the late Roman civilisation and Empire, to which Anatolia belonged, became part of the Byzantine culture and Empire, Greek in language and ideology, but Christian in religion (6 th -15 th century AD). Under Byzantine rule the intermittent Arabic invasions took place, before the Seljuk Turks conquered Anatolia in the late 11 th century AD. After their crushing by the Mongols, followed by the Beylik period (Turkish principalities), the Anatolian peninsula was unified again into the Ottoman Empire, from the 14 th to the early 20 th century AD, when it became the major part of the current Republic of Turkey [START_REF] Edwards | The Middle East and the Aegean Region, c.1800-1380 BC[END_REF][START_REF] Sahoglu | The Anatolian trade network and the Izmir region during the early bronze age[END_REF][START_REF] Mitchell | Anatolia:land, men, and gods in Asia Minor I, The Celts in Anatolia and the Impact of Roman Rule New York[END_REF] . Evidently, the Anatolian population history is the result of complex processes involving different population groups, the impact and contribution of which to the genetic pool of the Anatolian population is still relatively unknown. Molecular genetic analyses of modern human populations provided important clues to reconstruct the past migratory trajectories in the eastern Mediterranean and Anatolia 2,20-24 . Nevertheless, these results are limited by the fact that various demographic events could have obscured over time the evolutionary history inferred from modern populations. In this context, genetic analysis of ancient populations appears to be a promising way to investigate the origin of human populations and to validate hypotheses based on modern genetic data and other fields of study (e.g. historical and archaeological data).

In the present study, we describe the mitochondrial DNA (mtDNA) variation of a Byzantine population (11 th -13 th century) from the archaeological site of Sagalassos, in Southwestern Anatolia. The importance of this population is due to both (i) a continuous human occupation of its territory from the later 6 th millennium BP up to the 13 th century AD [START_REF] Waelkens | Man and environment in the territory of Sagalassos, a classical city in SW Turkey[END_REF] , and (ii) its very well known historical and archaeological background resulting from 25 years of research (http://www.sagalassos.be) which revealed the different influences that Sagalassos had undergone since the Bronze Age through its affiliation to successive empires/civilizations in this region (see Ricaut & Waelkens 26 for a more detailed review of the historical context of Sagalassos).

The purpose of this study is to estimate the maternal genetic signature potentially left in the Byzantine population from Sagalassos by the historic events impacting the Sagalassos region and South-western Anatolia since the Bronze Age, and whether these events may have been associated to large population movements, to more restricted elite dominance, or only to cultural diffusion. More particularly, on the basis of the available historical and archaeological records, our objective is to test whether differential contributions to Sagalassos population were made by populations from the Mediterranean area (Balkan and Italian peninsula), the Caucasus, the Near/Middle East and Central Asia.

MATERIALS AND METHODS

Samples

The Classical to early Byzantine archaeological site of Sagalassos is located in the Taurus Mountains of Southwestern Turkey at an altitude of around 1,500 m (Figure 1).

The excavations of this site revealed a city with monumental architecture and several necropoles, which were given up during the 13 th century AD [START_REF] Waelkens | The Late Antique to Early Byzantine city in Southwest Anatolia. Sagalassos and its territory: a case study[END_REF] .

The present work analyzed 85 individuals excavated from two different locations: 57 skeletons were recovered from pits in the context of a slope composed of earthquake debris caused by a seismic event in the late 6 th century AD, flanking the area of the 'Lower Agora', while 28 individuals were found in a neighbouring churchyard, within the adjoining 'Apollo Klarios Shrine' [START_REF] Waelkens | The Late Antique to Early Byzantine city in Southwest Anatolia. Sagalassos and its territory: a case study[END_REF] . The archaeological context pointed to lowsocial status burials which followed the Christian tradition (see Supplementary text).

All graves have been dated to the 11 th -13 th century cal AD by AMS carbon-14 dating of human bones (Analytic Inc. http://www.radiocarbon.com) and/or from the stratigraphical sequencing coupled with the general historical reconstruction of the site.

DNA extractions were made from more than 169 bone and tooth samples from the 85 skeletons.

Contamination avoidance

Genetic analyses were performed in the ancient DNA (aDNA) facilities of the Laboratory of Forensic Genetics and Molecular Archaeology in Leuven. Standard contamination precautions as previously described in literature were employed for the aDNA analysis [START_REF] Cooper | Ancient DNA: do it right or not at all[END_REF][START_REF] Hofreiter | Ancient DNA[END_REF][START_REF] Pääbo | Genetic analyses from ancient DNA[END_REF] , amongst which: pre-and post-PCR procedures carried out in physically separated laboratories; restricted and controlled access to the pre-PCR laboratory; mock-extraction controls in the DNA extractions and amplification reactions; independent replication of data; analysis of associated fauna remains. More details about contamination avoidance procedures are reported as Supplementary text.

DNA extraction and PCR amplifications

After bone and tooth samples decontamination and pulverization, aliquots of 0.5-1 g powder were digested in a Proteinase-K lysis buffer and DNA was extracted through silica-based spin columns [START_REF] Yang | Improved DNA extraction from ancient bones using silica-based spin columns[END_REF] . Analysis of the first and the second hypervariable segments (HVS-I and HVS-II) of the mtDNA control region was accomplished by amplifying and direct sequencing respectively five and two overlapping fragments (ranging from 109 to 166 bp in size), in order to read respectively 357 bp, from nucleotide position (np) 16009 to 16365, and 217 bp, from np 49 to 265. Furthermore, a selection of samples was used for cloning of the PCR products. This allowed us to confirm the results obtained by direct sequencing and to have an indication about post-mortem DNA damage in the samples (Supplementary text). Eight diagnostic markers in the mtDNA coding region were selected and analysed by means of multiplex amplification and minisequencing reactions in order to confirm and in some instances to resolve the haplogroup classification. List of primers used and details about the amplification, sequencing and minisequencing reactions are reported in Supplementary table 1 and Supplementary text.

Data analysis

On the basis of the HVS-I/HVS-II haplotype and the SNPs in the mtDNA coding region each individual was assigned to a haplogroup (hg) according to the latest mtDNA phylogeny [START_REF] Van Oven | Updated comprehensive phylogenetic tree of global human mitochondrial DNA variation[END_REF] (http://www.phylotree.org/). Comparative surveys with mtDNA data of up to 157 African and Eurasian populations reported in literature (Supplementary tables 2,3,4, Supplementary text) were done in order to detect the genetic signature left by the populations considered descendant of those which, from historical documents and archaeological evidence, may have contributed culturally and/or demographically to the ethnogenesis of the Sagalassos population since the Bronze Age [START_REF] Waelkens | The Late Antique to Early Byzantine city in Southwest Anatolia. Sagalassos and its territory: a case study[END_REF][START_REF] Waelkens | Romanization in the East: A case study. Sagalassos and Pisidia (SW Turkey)[END_REF][START_REF]The Chora of Sagalassos: The evolution of the settlement pattern on the territory of Sagalassos (Pisidia, SW Turkey) from prehistory until recent times[END_REF] . Search of haplotypes shared and analysis of Fst genetic distances based on the HVS-I and HVS-II variation were carried out in order to infer relationship between Sagalassos and the comparative populations (ARLEQUIN 3.11 software). The statistical significance of Fst values was estimated by permutation analysis, using 10,000 permutations. The pattern of genetic differentiation was visualized by MDS (XLstat v.7.5.2) and geographical interpolation method (Surfer v.8.0, Golden Software Inc 2002). Furthermore, a Principal Component Analysis (PCA) was carried out on the basis of haplogroup frequencies from 80 Eurasian populations present in literature (Statistica v.8, StatSoft Inc). See Supplementary text for detailed method.

RESULTS

Even when applying strict precautions, contamination from modern human DNA cannot yet be totally ruled out [START_REF] Pääbo | Genetic analyses from ancient DNA[END_REF] , and a strong logical chain of evidence is required to authenticate aDNA results [START_REF] Gilbert | Assessing ancient DNA studies[END_REF] . In the present study, on the basis of the procedures followed and the criteria used to assess the reliability of results, we could exclude contamination with a high level of confidence and attest to the authenticity of our aDNA data (Supplementary text).

Genetic diversity and kinship

HVS-I and HVS-II sequences were determined in 53 out of the total 85 initial individuals (62%). Analysis of SNPs in the coding region of the mtDNA was done in a selection of samples and haplogroup affiliation could be ascertained in all the individuals (Supplementary table 5). Ten haplotypes were shared between two or more individuals. On the basis of the burial location in the cemetery we assumed a putative maternal relationship in two couples of individuals buried in relatively close graves (SA2005 AK41a, SA2005 AK41b2 and SA2005 AK41b, SA2005 AK41c; see

Supplementary text).

After excluding the related individuals, the final Sagalassos dataset encompassed 51 individuals, revealing 39 different haplotypes. This accounts for the fairly high haplotype diversity observed in Sagalassos (H=0.989 ± 0.006).

All Sagalassos sequences fall into a set of haplogroups within macrohaplogroups N (X, W and N1b) and R (R0a, H, V, HV, U, K, J, T) that are characteristic of West Eurasians. None of the sequences from Sagalassos belong to the East/South Asian macrohaplogroup M, and sub-Saharan haplogroups (L1, L2, and L3A).

Population comparisons

To infer phylogenetic relationships between the Sagalassos population and 76 Eurasian and African populations we performed a non-metric multidimensional scaling (MDS) analysis based on the Fst distance matrix calculated from the HVS-I/HVS-II haplotypes (Supplementary table 6). The results of the MDS (figure 2) suggest that geography can explain the grouping of populations on the two-dimensional plot. Sagalassos is located very close to the samples from Turkey, Crimea and Iran [START_REF] Calafell | From Asia to Europe: mitochondrial DNA sequence variability in Bulgarians and Turks[END_REF][START_REF] Comas | Admixture, migrations, and dispersals in Central Asia: evidence from maternal DNA lineages[END_REF][START_REF] Derenko | Phylogeographic analysis of mitochondrial DNA in northern Asian populations[END_REF] which present a nonsignificant differentiation with Sagalassos (Fst<0.01) (Supplementary table 8) and their position in the middle of the plot reflects a geographic bridge between Asia and Europe.

Regarding Asian populations, a genetic proximity in the MDS plot of some Uzbek samples (Fst<0,01) to Sagalassos has been observed.

The maternal gene pool of Sagalassos was then compared to a higher number of populations (i.e. 157 populations) but from a shorter mtDNA fragment (HVS-I), and the Fst values obtained (Supplementary table 7) have been graphically plotted on a map of Eurasia and North Africa (Figure 3). Results are consistent with those from the HVS-I/HVS-II comparisons, showing the genetic proximity of Sagalassos with its with Central Asian samples, one with Tajiks (the haplotype AK15a, assigned to haplogroup N1b) and one with Uygurs (AK60, assigned to haplogroup T2b) were detected.

DISCUSSION

The high rate of successful genetic analyses (i.e. 62% of the total individuals sampled) indicates the good state of preservation of the aDNA in the Byzantine samples from Sagalassos. This is likely related to the high altitude of the archaeological site (i.e., approximately 1,500 m), which provides suitable thermal conditions for the preservation of DNA even at the low latitudes of South-West Turkey. Among the samples resulting in unsuccessful analyses, a large fraction (91%) is from the Lower Agora site. This site was excavated 10 to 15 years ago, and differently from Apollo Klarios, skeletons from the Lower Agora were unearthed without the necessary precautions for avoiding contamination, this confirming how crucial the postexcavations procedures are for aDNA investigations [START_REF] Pruvost | Freshly excavated fossil bones are best for amplification of ancient DNA[END_REF] .

Within the context described in the introduction, the mtDNA pool of the Byzantine population of Sagalassos brings new elements to estimate the genetic signature associated to the historic/demographic events which took place in this region since the Bronze Age. The results of our study point to a West-Eurasian composition of the mtDNA gene pool of the 11 th -13 th -century Byzantine population. This genetic pattern is in agreement with results from previous studies investigating Anatolian/Turkish populations history and obtained from different biological markers, e.g. proteins, Alu sequences, mtDNA, Y chromosome, and morphological data [START_REF] Cinnioglu | Excavating Y-chromosome haplotype strata in Anatolia[END_REF][START_REF] Benedetto | DNA diversity and population admixture in Anatolia[END_REF][START_REF] Berkman | Alu insertion polymorphisms and an assessment of the genetic contribution of Central Asia to Anatolia with respect to the Balkans[END_REF][START_REF] Ricaut | Cranial discrete traits in a Byzantine population and eastern Mediterranean population movements[END_REF][START_REF] Calafell | From Asia to Europe: mitochondrial DNA sequence variability in Bulgarians and Turks[END_REF][START_REF] Brega | Study of 15 protein polymorphisms in a sample of the Turkish population[END_REF][START_REF] Tambets | The topology of the maternal lineages of the Anatolian and Trans-Caucasus populations and the peopling of Europe: Some preliminary considerations[END_REF] .

More particularly, a genetic affinity of Sagalassos with Anatolian and Balkan populations has been revealed by our study (Figure 2, Figure 3, Supplementary tables 8 and 10), suggesting a relatively similar background of the maternal gene pool in these populations. The affinity with current Anatolian populations, suggests that the same historic and demographic events that shaped the mtDNA pool of Sagalassos, might have left a genetic signature in the modern Turkish populations.

The genetic proximity with Balkan populations may probably be the result of the long term interactions between Anatolia and Balkan regions since at the latest the Early Iron Age, as attested by linguistic, archaeological, mythological and biological studies 9- 10,26,32 . This might be a genetic signature of the Phrygians originating from the Balkans, who filled the political hiatus in Central Anatolia left by the Hittites from the 10 th century BC onwards. Yet, it is uncertain whether domination of the Phrygians ever extended as far as Sagalassos, although Phrygian influences are present in its territory.

Alternatively, the Balkan genetic signature might date back to the 3 rd -2 nd century BC, when Macedonian Seleucid dynasty ruling over Pisidia enlisted many mercenaries in their armies, this leading to mixed marriages. Particularly in Sagalassos a whole series of Seleucid colonies populated by Macedonian veterans and Jews were created [START_REF] Kosmetatou | The 'Macedonian Shields of Sagalassos[END_REF] .

Our analyses also detected an affinity between Sagalassos and populations from Iran and the Italian peninsula. Within the historical context described above, this pattern may result from historical events which facilitated the installation of Persian and Roman settlers. When ruling over Sagalassos, the Persian Achaemenids (546-333/332 BC) enlisted many Pisidians in their armies, but also acquired estates near the city of Sagalassos and even must have had settlers in the city itself. As for the Italian affinity, this can easily be explained by the fact that Augustus, after his confiscation of the Kingdom of Galatia, to which Sagalassos had belonged from 39 to 25 BC, established at least eight colonies of Roman veterans from Italy within a radius of a 100 km around Sagalassos [START_REF] Waelkens | Romanization in the East: A case study. Sagalassos and Pisidia (SW Turkey)[END_REF] .

A Levantine component has been also detected in Sagalassos, which is likely related to the high frequencies of haplogroups U3 (5.66%) and N1 (11.32%) (Supplementary table 9, Figures S1 andS2). This result may be explained by an intense trade system that was established with the Levantine populations from Archaic, over Hellenistic and Roman Imperial to Early Byzantine time [START_REF] Arndt | Roman trade relationships at Sagalassos (Turkey) elucidated by ancient DNA of fish remains[END_REF] , when Southwest Anatolia formed a kind of cultural koinè in its material culture with Western Cyprus, Pamphylia and the Levant itself. This is further supported by the proximity of the sample from Cyprus in the PCA (figure S1).

As an alternative scenario, it should be mentioned that the introduction of Levantine mtDNAs may result from the allocation of Pisidian mercenaries by the Macedonian dynasty of the Ptolemies in the Levantine coasts, where they may have married local women with whom they eventually returned to their home cities [START_REF] Bracke | Pisidia in Hellenistic Times[END_REF] .

It is worth noting that as far as the U3 lineages are concerned, the hypothesis of a Caucasian origin cannot be totally ruled out [START_REF] Richards | Tracing European founder lineages in the Near Eastern mtDNA pool[END_REF][START_REF] Bermisheva | Phylogeografic analysis of mitochondrial DNA Nogays: the[END_REF] . The presence of this haplogroup in the Sagalassos territory might represent an ancient genetic legacy of the Indo-European tribes from Trans-Caucasus which occupied Anatolia at the latest around the beginning of the 2 nd millennium BC. In particular, the Hittites and the Luwians settled in Central and Southern Anatolia, whereby the territory of Sagalassos became part of a Luwian state [START_REF] Waelkens | Sagalassos and Pisidia during the Late Bronze Age[END_REF] .

Overall, the analyses that we performed point to the absence of a significant Central Asian genetic signature in Sagalassos, at least at the maternal side. With this regard, the apparent affinity detected with some Central Asian populations (e.g. the Uzbeks in figure 2, Turkmen and Kurds in figure S1, Tajiks and Uygurs in Supplementary table 10) likely result from the high frequencies of West Eurasian specific haplogroup in these samples (up to 42% in Uygurs and Uzbeks 47 ). As far as evidenced by the archaeological data, presence of the Seljuk Turks, a nomadic tribe from Central Asia which invaded Anatolia in the late 11 th century AD, is not attested in Sagalassos before the first half of the 13 th century AD, when they built a hamam and according to Arabic sources also a caravansary in the neighbour site of Ağlasun 48 . Even if it is not excluded that they had already occupied part of the former territory of Sagalassos before, the 11 th to 13 th century AD cemetery which is analyzed here shows a clearly homogeneous Christian population. Remarkably, a 12 th -13 th century AD fortified site on the southern part of Sagalassos was destroyed and erased by the Seljuks in the course of the 13 th century AD, making a definitive end to the occupation of the site and showing that the relations between the two population groups seem to have been hostile.

In conclusion, despite the limitations of investigating a single maternally inherited locus, and the impossibility to perform statistical methods to test genealogical relationships between past and present populations (within the framework of the Sagalassos project, research is restricted to the period preceding the capture of Constantinople by the Ottoman Turks in 1453 AD), our study provides new clues to pinpoint the impact of historical events, and correspondingly of the related population groups, on the maternal gene pool of the Byzantine population from Sagalassos.

The analysis of the Fst genetic distances and the shared haplotypes highlighted a genetic proximity of Sagalassos with populations from the Balkan and Italian peninsulas, and the Levant. Notably, no signature from Central Asia has been detected, while an ancient genetic legacy of the Indo-European tribes from Trans-Caucasus which occupied Anatolia at the beginning of the 2 nd millennium BC cannot be excluded. Whether this relatively good agreement between female history (mtDNA) and historical/archaeological records is also found in the paternal lineages, will have to be explored through investigation of Y chromosome markers.

ACKNOWLEDGMENTS

We wish to thank Maarten Larmuseau for critical reading and comments on earlier drafts of the present manuscript. This research was supported by the Belgian Programme on Interuniversity Poles of Attraction (IAP 6/22) and the Research Fund of the K.U.Leuven (BOF-GOA 07/02), next to projects G.0286.00 and G.0421.06 of the Fund for Scientific Research, Flanders (FWO). M. Waelkens is one of the people who received a 'Methusalem grant' from the Flemish Ministry for Science Policy.

The mtDNA sequences reported in this paper have been deposited in GenBank (Accession Numbers: HQ007946-HQ008051).

neighbouring populations, particularly some Eastern Mediterranean groups (Turkey/Anatolia), the Balkans (mainland Greece, Bulgaria, Croatia) and Crimea 20,34,37- 39 .

Overall, it is remarkable that the Fst values point to a proximity of Sagalassos with the samples from Turkey, Crimea, Iran, Italy (Campania, Puglia), Cyprus and the Balkans (Bulgaria, Croatia, Greece), which show the lowest non-significant Fst values (Fst<0.005, Supplementary Table 8).

Haplogroup composition in Sagalassos was investigated by means of the PCA (Supplementary table 9). In the two-dimensional plot (Figure S1), which explains 32% of the total variation, Sagalassos is located in a marginal position, close to a cluster of populations from the Caucasus (the Georgian sample being the closest to Sagalassos), the East Mediterranean (Cyprus and Rhodes), and the samples from the Near/Middle East. In this pool of populations, contribution of haplogroups U3, N1, U1, and HV to the first and second PC was found to be particularly high (Figure S2).

Twenty HVS-I/HVS-II haplotypes from Sagalassos were found to match modern sequences (Supplementary table 10 and Supplementary text). In agreement with the results of the Fst based analyses, we detected relevant haplotype matches with the populations that historically may have contributed to the make-up of the Sagalassos mtDNA pool, particularly the Balkan area (e.g. Northern Greece, Macedonia, Bosnia and Bulgaria), Italy and Iran. Of particular note is the presence of haplotype LA35 exclusively in Macedonia (frequency 0.5%), and AK41a and LA20 exclusively in Italy (both at 0.25%), as well as the high frequency of AK15a in Iran (15%). Shared haplotypes were also found in the Near/Middle East (e.g. AK62 and AK56-69_2, 2.33% in Saudi Arabia), in particular LA54 found only in Jordan (1%). Finally, two matches

CONFLICT OF INTEREST

TITLES AND LEGENDS TO FIGURES

SUPPLEMENTARY INFORMATION Supplementary tables (Ottoni_et_al_Supplementary_tables.xls):

-Supplementary table 1. List of primers used in this study and product size.

-Supplementary table 2. References of the populations compiled in the HVS-I database.

-Supplementary table 3. References of the populations compiled in the HVS-I/HVS-II database.

-Supplementary table 4. References of the populations used in the PCA analysis -Supplementary table 5. Polymorphic positions observed in the HVS-I, HVS-II and the coding region of the mtDNA, and haplogroups affiliation in 53 ancient individuals from Sagalassos (11th -13th centuries AD) -Supplementary table 6. Matrix of Fst genetic distances between the 77 populations included in the HVS-I/HVS-II database.

-Supplementary table