Shang-Yu Yu

Dr Gautier Schrodj

Dr Karine Mougin

Dr Joseph Denzer

Jean-Pierre Malval

Hsiao- Wen Zan

Dr Olivier Soppera
email: olivier.soppera@uha.fr

Arnaud Spangenberg
email: arnaud.spangenberg@uha.fr

S-Y Yu

Dr J-P Malval

Prof. H.-W Zan

Direct Laser Writing of Crystallized TiO2 and TiO2/Carbon Microstructures with Tunable Conductive Properties

Keywords: femtosecond direct laser writing, functional metal oxide, conductive micropattern, TiO2 / carbon nanocomposite, miniaturized pressure sensor

Metal oxides are an important class of materials for opto-electronic applications. In this context, developing simple and versatile processes for integrating these materials at the microscale and nanoscale has become increasingly important. One of the major remaining challenges is to control the microstructuration and electro-optical properties in a single step.

We show here that near-infrared femtosecond laser irradiation can be successfully used to prepare amorphous or crystalized TiO2 microstructures in a single step using a direct laser writing approach from a TiO2 precursor thin film doped with a suitable dye. When laser writing is conducted under a nitrogen atmosphere, simultaneous to the crosslinking of the Tioxide precursor, the graphitization of the organic species embedded in the initial film is observed. In this case, a carbon network is generated within the TiO2 matrix, which significantly increases the conductivity. Moreover, the TiO2/C nanocomposite exhibits piezoresistive behavior that was used in a pressure sensor device. Using this route, it is possible to use direct laser writing to fabricate microsized pressure sensors.

Metal oxides (MOs) have emerged as important functional materials for the design of devices that are used in a broad range of applications, such as health monitoring, electronic skin, displays or solar cells, among others. [1][START_REF] Ogale | Functional Metal Oxides: New Science and Novel Applications[END_REF][START_REF] Comini | [END_REF][4] Among all MOs, ZnO, TiO2, ZrO2, and HfO2, in particular, are attractive due to their excellent optical properties (transparency in the visible range and high refractive index) combined with electronic properties for use in conductors, semiconductors, and insulators. [5] Electronic properties are indeed closely linked to the nature of the metal and conditions of fabrication, which determine the structure of the final material at the atomic scale. These properties can be modified by doping with different metals, by adjusting the metal to oxygen stoichiometry or by generating an amorphous or crystal phase. [6-9] Today, integration of MOs as thin films, micro-patterns or nano-patterns by convenient and simple means remains a major challenge for the fabrication of electronic and optoelectronic devices. MOs are commonly deposited as thin films by sputtering. Though the sputtered MOs usually exhibit better electrical properties due to the accurate control of the composition and defect concentration at the atomic scale, solution-based processes have gained much attention during the last years due to their simplicity, cost-effectiveness, vacuum-free processes and high versatility.

In solution-based processes, thermal annealing is required to eliminate the solvent and organic ligands that are complexed on metal ions and to obtain material and phase compositions with suitable properties. This step is typically conducted at temperatures ranging from 300°C to 1000°C. Thus, such a fabrication strategy can be difficult to apply in the case of multistep integration processes involving different materials.

In this context, laser processing of MOs presents many advantages due to the control of the laser-matter interaction in space and time. [9] Laser curing allows a fast treatment that can be confined to a limited volume by focusing the laser beam. Such a laser-induced effect leads to a much lower energy consumption with respect to that of thermal annealing. For example, a pulsed laser in the ultra-violet (UV) or near-infrared (NIR) range has been successfully used for the crystal growth of amorphous MO films [10,11] or the hydrothermal growth of MO nanostructures. [12,13] Recently, laser processing has also been introduced for the direct writing of MO microstructures. For this purpose, UV or Deep-UV (DUV) wavelengths are usually used. At the molecular scale, the patterning proceeds according to the photo-decomposition of the metal precursors complexed by organic molecules. Patterning can be obtained at room temperature down to the nanoscale using direct laser writing (DLW), but usually, a thermal post-treatment is needed to obtain an MO with suitable optical or electrical properties. [6-8, 14, 15] For example, Yeh et al. have shown that laser irradiation of a zinc oxide (ZnO) sol-gel based photoresist using DUV lithography can produce nanostructures. [16] However, the resulting amorphous material exhibits a limited electrical performance. Nevertheless, after 1 hour of thermal annealing at 500°C, the field-effect mobility increased from 0.0075 cm 2 .V -1 to 0.034 cm 2 .V -1 , which was attributed both to the elimination of the organic residue inside the original photoresist and to the crystallization of ZnO. Kim et al. have also demonstrated that DUV irradiation combined with moderate heating (150°C) allows for suppression of thermal annealing post-treatment while preserving the final performance of the MO. [1] However, the feasibility of performing such a strategy for micropatterning has not been explored. In another study, TiO2 sol-gel precursors were directly excited by femtosecond laser pulses at 800 nm through a multiphoton process. The authors showed that the modification was strongly correlated to the laser fluence. However, no investigation of the resulting material was performed and thermal curing was needed to obtain the crystallized material. [17] In the present work, we propose a new strategy combining photothermal and photochemical activation of TiO2 sol-gel precursors using femtosecond NIR direct laser writing. TiO2 was chosen due to its potential as a transparent, highly refractive material with semi-conducting properties. TiO2 has shown great performance in opto-electrical applications (micro-optics, MEMS devices, dye-sensitized solar cells, electrochemical sensors, and photocatalysts). [18][19][20][21] For practical applications, a simple means of integrating TiO2 with a crystallinity control is needed. [22] We show that our strategy allows for the fabrication of micropatterned crystallized TiO2 in a single step under air and at room temperature. Activation of the TiO2 sol-gel precursor is efficiently triggered under a 2-photon absorption mechanism by a photoabsorber added in the metal oxide precursor solution. The laser induced chemical modification was investigated by coupling photostructuration experiments and Raman microspectroscopy. Remarkably, by adjusting the laser fluence and surrounding atmosphere, we also show the ability to tune the electrical properties via in situ generation of disordered carbon, leading to a novel and easy method for integrating TiO2/carbon composites in microdevices. As a proof of concept, electrical characterization of TiO2/carbon microstructures is performed, and an application as a miniaturized pressure microsensor is shown in the last section.

Figure 1 presents a schematic view of the fabrication process of microstructures using femtosecond DLW for titanium precursor solutions and their application as a pressure sensor.

The first step is the solution preparation. A commercial titanium precursor (titanium isopropoxide) was mixed with methacrylic acid in 1-propanol. Carboxylic acids are commonly used as stabilizers of transition metals since they can form strong bonds with metals. After 24 hr of stirring in an acid environment, a stable solution was obtained. The solution preparation conditions were chosen to limit the condensation reaction of the titanium precursor. The excess complexing agent associated with the presence of water was shown to lead to metal-oxo clusters (MOC) in the solution. [23][24][25] We assume that we formed such species in the early stages of the reaction. However, we were unable to directly characterize the size and structure of the MOC after spin-coating since partial condensation reactions may occur by reacting with atmospheric moisture. These conditions limit the size of the aggregates of metal-oxo clusters, which is very important for ensuring low surface roughness during the spin coating process and for limiting the line edge roughness after development. In the last step of the solution preparation, 4,4'-Bis(diethylamino)benzophenone (DABP), which is wellknown because of its photochemical properties in two-photon absorption [26,27] , was added to the solution. DABP serves as a photo-absorber for two-photon absorption at a wavelength of 800 nm. This compound was also chosen because of its good solubility, which is needed to obtain a homogeneous distribution of the photo-absorber inside the thin film.

After ageing for 24 hours, the solution was casted by spin-coating it on a glass slide. After a pre-baking process (100 °C, 10 s) to remove the solvent from the thin film, the sample was directly introduced into a chamber, enabling control of the atmosphere during the DLW fabrication process. With this setup, the fabrication process was carried out either under air or a pure nitrogen atmosphere. Irradiation was accomplished using a femtosecond laser with an adjustable irradiation time and power. After irradiation, the samples were developed over 20 s in cyclohexanone to remove the non-irradiated part, which revealed the microstructures. For electrical characterizations, the microstructures were connected to gold electrodes deposited by Au sputtering with shadow masks.

We first investigated the effect of femtosecond laser irradiation on the MOC thin film under various laser powers. In this first set of experiments, laser irradiation was conducted in air.

Figure 2a shows a SEM image of the line structures with different exposure powers. The formulation can be successfully patterned using femtosecond irradiation. For the chosen irradiation time (10 ms for each pixel), crosslinking is obtained for any laser power greater than 10.0 mW. It is observed that the line width increases along with the increase of the exposure power. A width of approximately 650 nm is obtained for the lowest power (10.0 mW). For the highest power (16.5 mW), the width increases to 3 m. The line edge roughness is limited and continuous lines were obtained for lower and higher powers. For intermediate powers, irregularities within the line pattern were systematically observed. Interestingly, such defects were obtained with conditions corresponding to the maximum release of organic functions from the thin film. Elimination of organic molecules could provoke this roughness.

For higher power, further densification and crystallization smooth the surface of the patterns.

The initial thickness of the film was set to 500 nm by adjusting the spin-coating parameters (from ellipsometry measurements). At low power (10 mW), AFM revealed a height of the lines of 350 nm. At highest power (16.5 mW), the thickness decreased to 80 nm. Such values correspond to a vertical shrinkage of 30 % and 85 % respectively, that corresponds to elimination of organic molecules and densification by condensation and crystallization. These structures were analyzed using Raman spectroscopy to monitor the chemical modification induced by the laser (Figure 2b). The first spectrum is related to the film deposited without any laser treatment. Peaks at 1420, 1594 and 1638 cm -1 are visible. From previous studies on metal oxo clusters stabilized by MAA, [14,28] these peaks were assigned to the MAA linked to Ti. More specifically, these peaks were assigned to s(COO -), as(COO -) and (C=C), respectively. [29] For the low power exposure (10.0 mW), a significant decrease of these peaks was observed, which was attributed to the cleavage of the metal to ligand bond and the release of MAA. Note that this power corresponds to the minimum value to obtain sufficient crosslinking for obtaining stable microstructures on the substrate. With the increased exposure power from 13.2 mW to 16.5 mW, new peaks appeared. The peaks at 435 and 610 cm -1 indicated the presence of the rutile phase (R mark in Figure 2b). [30,31] Two other peaks at 1360 cm -1 and 1597 cm -1 were attributed to the disordered carbon peaks at the D band and G band, respectively. [32,33] Usually, the D band is positioned at 1325 cm -1 . However, a shift of approximately 35 cm -1 was found in the Raman spectrum, which was attributed to the formation of the Ti-O-C bond. [34] It was also observed that the disordered carbon peaks had significantly decreased with the increased exposure power, which was explained by the oxidation of the carbon species under high laser fluency in the presence of oxygen from air. This point will be further discussed below.

The evolution of the Raman spectra with the laser power follows a similar trend to a thermal treatment with an increasing temperature. To verify this point, thin films with the same compositions and thicknesses were prepared in another experiment and were heated in a furnace at different temperatures from 200 °C to 800 °C (1 hour). Figure 2c shows the Raman spectra of the Ti-MAA thin film annealed at different temperatures. After annealing at 200°C, a significant decrease of the peaks at 1420, 1594 and 1638 cm -1 was observed. As expected, this evolution corresponds to the crosslinking of the Ti-oxo clusters due to thermally induced condensation. This peak evolution is thus similar to that recorded at 10.0 mW, which indicates that laser irradiation under these conditions has the equivalent effect as thermal annealing at 200 °C (1 hour). It was then evaluated whether the decrease of Raman peaks could not be assigned only based on the evaporation of free ligand. For this evaluation, temperature programmed desorption-mass spectroscopy (TPD-MS) was performed on the Ti-MAA thin film. As shown in Figure S1a, no strong desorption peak was found below 300 ˚C except for a small amount of water molecules. A group of high desorption peaks at approximately 400 ˚C was attributed to the thermal decomposition of organic compounds into gaseous substances. [7] The above results indicated that the decrease of the Raman peak at 200 ˚C was attributed to the condensation process of Ti-MAA. The above results also showed that the excess MAA that was not bonded to the metal was eliminated during spin-coating. With the increase of the annealing temperature to 500 °C and 800 °C , the peaks corresponding to the anatase phase (A mark) and those of the rutile phase (R mark) were found, as expected.

Microstructures were prepared by two photon structuration (TPS) with the Ti-MAA precursor to evaluate the structural performance and demonstrate the applicability of the DLW approach.

The crystallized TiO2 microstructures (Figure 2d) and hybrid organic/TiO2 microstructures (Figure 2e) were fabricated at a high exposure power (16.5 mW) and low exposure power (10 mW), respectively. From the SEM and AFM images, high reproducibility and low surface roughness were observed, highlighting the potential to combine this method with TPS for creating arbitrary structures with tunable compositions by simply adjusting the fabrication parameters.

These results demonstrate that the thin film formed from the titanium precursor can be successfully patterned at the microscale and that the final material composition can be tuned by changing the laser power. Depending on this parameter, the hybrid material, amorphous TiO2 or crystallized TiO2 microstructures that are eventually doped with carbon, can be obtained in a single step at room temperature without any further post-treatment step.

At this stage, the presence of carbon did not lead to conduction, as shown in Figures S1b andS1c, even when using the photonic parameters that led to the maximum concentration of carbon from the Raman data. This result reveals that the amount of carbon is not sufficient to create a continuous conductive network. As stated before, we assume that C is generated by the graphitization of the organic compounds embedded in the thin film due to laser irradiation.

Graphitization of the organic compounds is achieved at an intermediate laser power, which does not result in more C but in the oxidation of CO2 because of the presence of oxygen in the surrounding atmosphere.

To verify this assumption, laser fabrication was investigated in an oxygen free environment.

The samples were fabricated with the same exposure conditions as those in Figure 2 but in a nitrogen atmosphere by purging the chamber with a continuous N2-gas flow to reduce the oxidation of carbon. Raman analysis on these samples was performed, and the results are shown in Figures 3a and3b. Figure 3a, which include an additional sample prepared only with thermal annealing at 900 ˚C for 1 hour in a nitrogen atmosphere for comparison. First, we observed that laser microstructuring is still possible under nitrogen. Significant differences in the material structures are evidenced by Raman spectroscopy, especially for the highest powers. The main difference is the large increase of the carbon content, as revealed by the clearly distinguishable D and G band peaks that are present under all conditions, including at the highest power. These two peaks are also present in the sample cured at 900 °C under nitrogen. Interestingly, the peaks of crystalline TiO2 could barely be identified, unlike in laser irradiation in air. According to a previous study by Jin et al., the formation of a considerable amount of carbon is able to constrain the crystallization of TiO2. [START_REF] Jin | Fullerenes[END_REF] Another hypothesis relies on the lack of oxygen atoms in the thin film to reach the stoichiometric ratio. To further understand the effect of two-photon irradiation on the graphitization degree (defined as the D/G band ratio) [36][37][38] , we normalized the carbon peaks with the intensity of the G band, as shown in Figure 3b. Along with the increase of the exposure power, the relative intensity of the D band decreases. This result indicates a greater degree of graphitization for the highest laser power. We assume that as the laser power is increased, the graphitization is favored.

Interestingly, the graphitization degree generated by laser irradiation at a high power (0.64) is better than that of thermal annealing at 900 °C (0.87). To further investigate the structure of the carbon phase, electrical measurements of the TiO2/carbon composite lines were performed using a two-probe system. As shown in Figures S2a andS2b, the lines that were fabricated at a low exposure power, from 10 mW to 11.5 mW, revealed no significant conductivity. When the exposure power was increased to 12 mW, the current level was increased, meaning that a continuous carbon network was created within the TiO2 matrix. The typical Ohm's characteristics were obtained. The resistivity of the lines were evaluated using Ohm's law (Figure 3c). [39] The resistance was calculated from the slope of the I-V curve, and the aspect ratio was determined according to the AFM image and gap between the two electrodes (100 m). A significant decrease of the resistivity was obtained when increasing the exposure power. This behavior is in agreement with the trend of the D/G ratio revealed by the Raman spectra, thus indicating a relation between the nature of the carbon phase and the resistivity of the material. Nevertheless, due to the insulating host TiO2 matrix, the lowest resistivity of the nanocomposite line exhibited a slightly higher resistance than that of the pure disordered carbon. [40][41][START_REF] Ju-Seung | The electrical and electrochemical properties of poly(p-phenylene)-based carbon for lithium rechargeable batteries[END_REF] Conversely, it should be mentioned that another mechanism may account for the significant increase of the conductivity of the material after laser treatment under a nitrogen atmosphere.

The partial reduction of TiO2 into TiO2-x may also occur in the presence of C. This behavior was described by Lu et al.. [START_REF] Lu | [END_REF] The introduction of oxygen deficiencies with the resultant formation of TiO2-x is known to increase the electrical conductivity. [44] Figure 3d shows SEM images of the TiO2/carbon composite structures. The process described here enables the use of DLW to generate well-defined micropatterns of TiO2 with the inclusion of graphitized carbon leading to conductive properties. It has to be mentioned that the resolution of the structures is slightly degraded compared to that at a lower power.

The typical structure width was approximately 4 m, which was probably due to the lateral heat spreading after laser irradiation.

Room temperature one-step photopatterning is thus achieved with the in situ formation of the conductive carbon phase. Moreover, this conductive phase combined with the insulating matrix can potentially exhibit piezoresistive properties and can thus be used as a pressure sensor. To demonstrate this concept, 10 lines connected with two gold electrodes (Figure 4a) were prepared following the previously described method. An exposure power of 16.5 mW in a nitrogen atmosphere results in a TiO2/carbon composite. A controlled pressure was applied by a mechanical performance testing system. A small glass coverslip was added on the top of the microstructure (Figure S3a). Because the coverslip was wider than the microstructures, it ensured that a homogeneous force was applied and also avoided the lateral movement of the needle, which might cause scratches on the surface and further damage the structure. Several tests on the pressure sensor were then performed to characterize the piezoresistive response.

As shown in Figure 4b, the output current signals under different forces from 0 to 1 N display a linear and dynamic response, indicating an instantaneous recovery of the pressure sensor to the loading-unloading of the pressure. To further confirm the reliability, 30 cycles of dynamic forces with different values ranging from 0.2 to 2.4 N were applied to the sensor. The responses of the pressure sensor to the applied forces for several cycles are shown in Figure 4c. A good repeatability of the sensing response in the linear and saturation regimes was observed, demonstrating a stable sensing system without degradation. Nevertheless, the stiffness of the microstructure is also important for the pressure sensor due to the possibility of having a huge force be accidentally applied. Therefore, we investigated the stiffness of the pressure sensor by gradually applying a force up to 3 N corresponding to two times the maximum response force in the linear regime, as shown in Figure 4d. A stable response plateau in the saturation regime and an accurate recovery curve after unloading were observed.

These results indicated the strong stiffness of our microstructure, which might be attributed to the good mechanical properties of the metal oxide network within our hybrid material structure.

The results show the good dynamic responses and reliability of the pressure sensor with no hysteresis. However, because of the use of the coverslip, the precise value of the applied pressure and sensitivity could not be defined due to the unknown contact area. Therefore, another experiment performed with an atomic force microscope was conducted to define the value of the applied pressure, as shown in Figure 4e. The samples were prepared in a single line structure with a width of 5.5 m. A special large plateau AFM tip with a circular flat top that was 10 m in diameter and a calibrated stiffness was used to apply the pressure. Thus, the contact area was defined as shown in Figure S3b. The tip gradually approached the sensor, and five load/unload cycles were performed, as shown in Figure 4f. It is noted that during the AFM approach, according to the default parameters of the AFM system, the cantilever will keep approaching until reaching a certain position, and then, after a small withdrawal, five load/unload cycles will be started. The pressure (P) can be calculated with the following equation: P = (VDK)/A, where V is voltage of the AFM deflection signal; D and K denote the deflection and spring constant of cantilever, respectively; and A is the contact area. The applied force was evaluated to be in the range of 23.9 Pa. The sensing of this range of pressure shows the sensitivity at the microscale from a single line. A repeatable response was observed. The sensitivity could be calculated as 0.186 kPa -1 based on the definition of S = (∆I/I0)/P. This value is comparable to the device responses fabricated by other means [45][46][47][48] , which shows the great potential of this new fabrication technique.

A possible pressure sensing mechanism of TiO2/C composite was proposed. At the nanoscale, as schematically illustrated in Figure 4g, the pressure applied on the microstructure leads to the compression of the structure. As a result, a diminution of the conductive pathway is expected, thereby leading to a wider conductive network. [49] Interestingly, while approaching the tip, a negative response was observed slightly before the contact point. This response could be explained by the attractive Van der Waals forces that are applied at a long distance in the approach step, thus causing a snap of the AFM tip at the surface. [50,51] It demonstrates the sensitivity of the pressure sensor to the attractive forces with a high sensitivity since the force was evaluated at 150 pN. This set of experiments demonstrates that the material prepared by femtosecond DLW can be used for monitoring forces with good sensitivity and reliability. The fabrication process allows for the preparation of miniaturized pressure sensors at well-defined positions or in high density arrays without thermal curing, which can be useful for applications in artificial skin or cell mechanical response investigations.

To conclude, a new strategy to produce TiO2 or TiO2/carbon microstructures using the direct laser writing process was proposed. By changing the laser power and surrounding atmosphere, it is possible to tune the final material composition and confer conductive and piezoresistive properties to the material. The on demand design of complex microstructures is possible due to the laser writing approach. This process greatly simplifies the integration of these materials for the future fabrication of complex devices. The use of femtosecond laser irradiation associated with a two-photon absorption pathway allows the microscale resolution to be reached. The minimum line width was 650 nm. With experimental conditions needed to obtain crystallized TiO2 and TiO2/C, typical line widths were 3 m. Pressure sensors were chosen to illustrate the potential applications of such an approach. Other applications for energy storage, solar cells or photocatalysis can also be developed from these results.

Experimental Section

Preparation of Solutions: One milliliter of titanium(IV) tetraisopropoxide (C12H28O4Ti, Sigma-Aldrich) and 2 mL of methacrylic acid (99%, Sigma-Aldrich) were mixed together and stirred for 5 min. Then, 2 mL of 1-propanol (anhydrous grade, 99.7%, Sigma-Aldrich) was added as a solvent and stirred for another 10 min. After that, 0.9 mL of acidified water (HCL, 0.37 M) was introduced into the solution. After 24 h of vigorous stirring, 1 wt% (to titanium precursor) of DABD (4,4'-Bis(diethylamino)benzophenone, Sigma-Aldrich) was finally added as a photo-absorber before use. The typical formulation composition was thus Ti:MAA:H2O:solvent = 1:7:15:8.

Two-photon fabrication of TiO2 structure:

The solution was filtered through a 0.2 m pore diameter filter and spin-coated on a glass substrate to obtain a xerogel film with a thickness of approximately 500 nm. A pre-baking step was used between spin-coating and laser irradiation (100°C, 10 sec) to remove the solvent from the thin film. The sample was directly transferred to a special chamber holder to control the atmosphere during fabrication. In particular, the chamber could be filled with dry nitrogen. For the microfabrication, we used a commercial set-up purchased from Microlight based on 3D direct laser writing technology. Thanks to this no mask technology, arbitrary and complex 2D and 3D shapes can be easely fabricated with 0.2 to 3 microns resolution in the x,y plane. This 3D microprinting machine include a standard microscope (Zeiss Observer D1), a piezo stage (PI Mars 300 µm 3 ) a sub-ns microlaser (TeemPhotonics) and a dedicated software to optimize fabrication paths and to control the machine. However, in this study, to favor the two-photon absorption process from the chosen dye, we have replaced the sub-ns source by a femtoseond laser source (Chameleon Ultra II, 140 fs @ 800 nm). The laser beam was introduced via a 40X objective (NA: 0.65), and the exposure time was set to 10 ms with different exposure powers from 10 mW to 16.5 mW. The beam diameter was evaluated to 750 nm. The local dose is thus comprised between 0.10 to 0.16 mJ and the surface energy between 22.6 and 37.3 kJ/cm 2 . After the laser exposure, the sample was developed by immersion in cyclohexanone for 20 s and dried by N2 purging.

Micro Raman Spectroscopy Analysis: The preparation of samples for thermal annealing was performed by spin-coating the solution onto a quartz substrate and annealing it in a furnace.

To monitor the laser induced modifications, the DLW system was used to fabricate a square structure with dimensions of 30x30 m 2 . Then, micro Raman spectroscopy analysis was conducted using a Spectrometer Labrum HR with a 532 nm excitation laser and CCD detector.

The Raman signals were acquired in a backscattering configuration with a 100X objective.

Characterization Methods: An SEM (Quanta 400 FEL) and an AFM (Picoplus 5500 Agilent) in resonant mode were used to characterize the samples in this study. For SEM characterization, a thin layer of Au around 10 nm was sputtered on samples.

Temperature Programming Desorption -Mass Spectrometry (TPD-MS) Analysis:

To investigate the degradation of the spin-coated thin film, a homemade temperatureprogrammed desorption coupled to a mass spectrometry setup was used. For the detailed description of the TPD-MS, one can refer to our recently published paper. [15] Fabrication and Electrical measurement of TiO2/Carbon structure: Two gold electrodes with a distance of 100 m were first sputtered onto a clean glass substrate using a shadow mask.

Then, the solution was spin-coated onto the glass substrate. Using DLW process, 10 lines with a width of approximately 5.5 m were fabricated between the two electrodes. The electrical performance was measured using a digital source meter (Keithley 2636B).

Performance test of TiO2/Carbon Pressure Sensor:

To apply a controlled force to the substrate, a mechanical performance testing system (Lhomargy DY 34) was used. The electrical response of the sensor was recorded with a digital source meter (Keithley 2636B).

The samples were prepared in the same way as those for the electrical measurements. A needle was used to apply the force in a confined area, and a small glass substrate was placed beneath pressure sensor to ensure that a homogeneous force was applied and to protect the structure. The sensor was operated with a voltage of 3 V to record the changes of the current.

For the microscale pressure tests, a force AFM (Nanosurf C3000) was used to apply pressure with a specially designed flat top AFM tip (Nanosensors SD-PL-NCH-SPL). The diameter of this flat top tip was 10 m. The pressure (P) was defined as P = (VDK)/A, where V is the deflection voltage of the AFM deflection signal; D and K denote the deflectivity and spring constant of the cantilever, respectively; and A is the contact area. In the task, five repetitions with a deflection voltage of 0.001 V were performed. By knowing the D and K values of the cantilever (125 nm/V and 10.5 N/m, respectively), a pressure of 23.86 Pa was calculated.

Figure 1 .

 1 Figure 1. Schematic of the preparation of the solution and the one-step micro pressure sensor fabrication process: Titanium isopropoxide is first chelated with the MAA ligand and a DABP photoabsorber is added. Then, the solution is spin-coated on a glass substrate, irradiated by an NIR femtosecond laser focused beam in a controlled atmosphere and developed in cyclohexanone. Finally, the structure is connected to an electrical circuit, and the electrical properties under a controlled pressure can be monitored.

Figure 2 .

 2 Figure 2. Laser curing in air of the TiO2 precursor material. a) SEM image of the line structure fabricated with different powers from 10.0 mW to 16.5 mW. b) Raman spectra of the microstructures generated by direct laser writing with different irradiation powers. c) Raman spectra of the spin coated thin film after thermal annealing at different temperatures (1 hour). R and A stand for "Rutile" and "Anatase" in the Raman spectra, respectively. d) Finger type microstructure fabricated at 16.5 mW. e) AFM image of the IS2M logo fabricated under a low exposure power (10.0 mW).

Figure 3 .

 3 Figure 3. Characteristics of the TiO2 microstructure generated in a N2 atmosphere. a) Raman spectrum of the TiO2 microstructures with different fabrication powers and a spin coated thin film annealed in 900 ˚C 1 h with a vertical shift of 0.3 unit. b) Zoomed in spectra in the range of 1100-2100 cm -1 and normalized by the carbon G bond. c) Resistivity and D/G bond ratio plot versus the exposure power of the TiO2 microstructures. d) SEM image of the arbitrary TiO2/carbon microstructure generated with 16 mW.

Figure 4 .

 4 Figure 4. Reliability, sensitivity, durability and mechanism of the TiO2/carbon composite pressure sensors. a) Optical microscopy image of the TiO2/carbon microstructure pressure sensor device. b) The change in the current under an applied force ranging from 0.2 N to 1.1

Figure

  Figure S2. a) and b) Linear and log scale of the I-V curve of TiO2/carbon line with different exposure powers from 10.0 mW to 16.5 mW, respectively.

Acknowledgements

The authors would like to thank Patrick Lamielle for designing gaz control chamber, and Simon Gree for helpfull discussion regarding the Raman microspectroscopy. This work was supported by funding from Agence Nationale pour la Recherche (Project PHOTOMOC, ANR-14-CE26-0039 & Project 2PhotonInsight, ANR-16-CE08-0020), Région Grand Est (MIPPI 4D) and Ministry of Science and Technology, Taiwan (Project 104-2923-E-009-001-MY3).

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.

Supporting Information


Mario Modesto-Mata

Rebeca García-González

Yuliet Quintino

Cecilia García-Campos

Marina Martínez De Pinillos

Laura Martín-Francés

María Martinón-Torres

Yann Heuzé

Eudald Carbonell

Juan Luis Arsuaga

Yuliet Quintino

Cecilia García- Campos

Juan Luis Arsuaga

M Christopher Dean

José María Bermúdez De Castro

Early and Middle Pleistocene hominins from Atapuerca (Spain) show differences in dental developmental patterns

à 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

This study examines relative dental development in two important samples of Lower and Middle Pleistocene humans, from Gran Dolina and Sima de los Huesos (Atapuerca, Burgos, Spain). Previous work showed that these samples presented on average, enamel formation time 27% faster than in modern humans [START_REF] Modesto-Mata | Short and long period growth markers of enamel formation distinguish European Pleistocene hominins[END_REF]. While this finding may suggest a faster overall biological clock in these earlier hominins, the sequence of tooth formation is also relevant to reconstruct life history of these hominins.

The variation in the sequence of dental maturation has been investigated between pairs of teeth (Tompkins, 1996 a;[START_REF] Tompkins | Relative dental development of Upper Pleistocene hominids compared to human population variation[END_REF] or taking into account the chronological ages [START_REF] Smith | Dental development in Australopithecus and early Homo[END_REF]. However, the study of the entire dentition rather than single or pairs of teeth provides a more complete picture of the dental development of an individual. One approach to investigating the sequence of dental maturation is to adopt a probabilistic approach using all teeth. [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF] proposed an analytical approach using Bayesian statistics to quantify variability in sequences of key events during tooth mineralization, independent of chronological age. This method considers the tooth within the mineralization sequence as a series of dependent units, growing within a developmental and hierarchical module. A statistical probability is determined for the likelihood of a particular dental mineralization sequence (DMS) to be present within known modern human variation. The modular system of dental development [START_REF] Raff | The Shape of Life[END_REF] presents the following features: 1) dental development presents an autonomous, genetically discrete organization, 2) the modular system of dental development is formed by hierarchical units, 3) dental modules have a physical place within the developing system, 4) the degree of connectivity between dental modules occurs at different levels, and 5) dental development displays temporal transformations. Overall, the variability in the patterns or DMSs is relevant in human evolution, as these processes allow to assess possible foundation for morphological changes [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF]. This approach has been recently applied to Neandertals (Bayle et al., 2009 a;[START_REF] Quam | The Neandertals of northeastern Iberia: New remains from the Cova del Gegant (Sitges, Barcelona)[END_REF], to an Upper Paleolithic specimen from La Madeleine (Bayle et al., 2009 b), and to the Lagar Velho individual [START_REF] Bayle | Dental maturational sequence and dental tissue proportions in the early Upper Paleolithic child from Abrigo do Lagar Velho, Portugal[END_REF], a possible Neandertal-H. sapiens hybrid [START_REF] Duarte | The early Upper Paleolithic human skeleton from the Abrigo do Lagar Velho (Portugal) and modern human emergence in Iberia[END_REF]. Remarkably, both Neandertals and the Lagar Velho specimen showed a DMS that is not known to exist in modern human variation, whereas the La Madeleine individual presented a DMS that falls within the known variation of modern H. sapiens. Moreover, the pattern of dental development of the Lagar Velho specimen is only present within the Neandertal variation. Overall, Neandertals show a proportionally advanced stage of mineralization of the first molar in comparison to their proportionally delayed stages of mineralization of the incisors. Interestingly, this statistical approach has rarely been applied to hominins older than Neandertals. Therefore, the Atapuerca hominins may give some clues on the changes of the pattern of dental development on earlier stages of human evolution.

The two archaeological sites from Atapuerca included in this study are Gran Dolina and Sima de los Huesos. Level TD6.2 of the Gran Dolina site (~0.9-0.8 Ma) contains more than 160 human fossil remains representing at least 8 individuals attributed to the species Homo antecessor [START_REF] Moreno | New radiometric dates on the lowest stratigraphical section (TD1 to TD6) of Gran Dolina site (Atapuerca, Spain)[END_REF][START_REF] Bermúdez-De-Castro | Homo antecessor: The state of the art eighteen years later[END_REF][START_REF] Duval | The first direct ESR dating of a hominin tooth from Atapuerca Gran Dolina TD-6 (Spain) supports the antiquity of Homo antecessor[END_REF]. Sima de los Huesos (SH) (~0.43 Ma) contains more than 7000 human fossils attributed to at least 29 individuals (Bermúdez de [START_REF] Bermúdez De Castro | The Sima de los Huesos Middle Pleistocene hominin site[END_REF] whose genetic and morphological data strongly suggest they belong in the Neandertal clade [START_REF] Arsuaga | Neandertal roots: Cranial and chronological evidence from Sima de los Huesos[END_REF].

Previous studies attempted to evaluate the patterns of dental maturation from these two fossil populations by treating tooth types as independent units. The information provided by the study of three hominins from the Gran Dolina site [START_REF] Bermúdez De Castro | A modern human pattern of dental development in Lower Pleistocene hominids from Atapuerca-TD6 (Spain)[END_REF][START_REF] Bermúdez De Castro | New immature hominin fossil from European Lower Pleistocene shows the earliest evidence of a modern human dental development pattern[END_REF] indicated relatively advanced calcification in the M3 compared to the M1. Hominin XVIII from the Sima de los Huesos site (Fig. 1) displayed a relatively delayed development of the lower and upper canines, a more advanced development of the lower second molars, and notably advanced development in the upper and lower third molars [START_REF] Bermúdez De Castro | Pattern of dental development in Hominid XVIII from the Middle Pleistocene Atapuerca-Sima de los Huesos site (Spain)[END_REF]. Bermúdez de Castro and [START_REF] Bermúdez De Castro | Pattern of dental development in Hominid XVIII from the Middle Pleistocene Atapuerca-Sima de los Huesos site (Spain)[END_REF] analyzed dental development with a non-Bayesian frame. However, this specimen was studied under this Bayesian analysis [START_REF] Bayle | Analyses quantitatives par imagerie à haute résolution des séquences de maturation dentaire et des proportions des tissus des dents déciduales chez les Néanderthaliens et les Hommes modernes[END_REF] using the modern human samples in [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF], showing that there were not differences in regard to modern humans. In our study we apply the Bayesian statistical analysis for the complete permanent dentition of hominin XVIII, including M3s.

Here, we present a Bayesian analysis of the relative dental development of five hominins from Gran Dolina and Sima de los Huesos. We aim to test the null hypothesis that DMSs of Atapuerca hominins are found within the range of modern humans, although they may present a more advanced molar development than our species. This study complements the one published recently on the methods of analyzing the absolute timing and pattern of lateral enamel formation [START_REF] Modesto-Mata | Short and long period growth markers of enamel formation distinguish European Pleistocene hominins[END_REF]. The results are compared with published data on Neandertals, to identify similarities or differences in the DMSs that might shed light about potential taxonomic signals.

Materials and Methods

Materials

The specimens analysed in this paper include two individuals from Sima de los Huesos-SH (XVIII and XXV) and three individuals of H. antecessor (H1, H3 and H11) from Gran Dolina-TD6.2 (Supplementary Table 1 andSupplementary Table 2, respectively). The criteria to determine that isolated teeth, either from Gran Dolina or Sima de los Huesos specimens, belonged to specific individuals is detailed in Bermúdez de [START_REF] Bermúdez De Castro | The Sima de los Huesos Middle Pleistocene hominin site[END_REF].

The complete permanent dentition (32 teeth) and the four deciduous second molars are present in hominin XVIII (Fig. 1). Hominin XXV includes a complete permanent lower dentition (16 teeth) and two deciduous second molars (Supplementary Fig. 1). All of the teeth from SH were isolated finds. The dentitions of the hominins from Gran Dolina are incomplete. H1 is composed of 16 upper and lower teeth (Supplementary Fig. 2), whereas H3 is composed of 9 upper teeth (Supplementary Fig. 3) .

The mineralization stages of the H. antecessor and Sima de los Huesos teeth were obtained by analyzing high-quality photographs of the fossils, along with 3D volume renderings from micro-CT scans. The developmental stages of the permanent dentition of fossil hominins were scored from A to H [START_REF] Demirjian | A new system of dental age assessment[END_REF], while deciduous dentition was scored from A to H2 [START_REF] Liversidge | Variation in crown and root formation and eruption of human deciduous teeth[END_REF]. The teeth and mineralization stages of the TD6.2 hominins and the SH specimens are presented in Table 1. Tooth mineralization stages of fossil hominins were identified independently by three co-authors (RG-G, YQ, MM-M). Discrepancies appeared in less than 4% of the teeth, and they were present over the latest stages of development of the root. In order to increase the utility of these data for other researchers who may want to reanalyze by using the [START_REF] Moorrees | Age variation of formation stages for ten permanent teeth[END_REF] scoring system, we also provide these scores in Supplementary Table 3.

When establishing the mineralization stages of the TD6.2 hominin H11 (Bermúdez de [START_REF] Bermúdez De Castro | New immature hominin fossil from European Lower Pleistocene shows the earliest evidence of a modern human dental development pattern[END_REF], a new undocumented tooth situated in its alveolus has been identified by using the micro-CT. This new tooth in the TD6.2 fossil hypodigm is the lower right fourth premolar of the hominin H11 mandible. This tooth is only represented by the initial formation stage of the buccal cusp enamel (Supplementary Fig. 4). Its uncompleted mesiodistal diameter measures ~3.4 mm.

Modern Human Reference Samples

Because our fossil sample includes individuals with only permanent teeth (only mandibular, only maxillary, or both), and other individuals with a mixed dentition (Table 1), it was necessary to use different comparative samples to perform different comparisons.

Burgos mandibular I: permanent mandibular dentition without M3

The first comparative sample is called "Burgos mandibular I" which was drawn from three different subsets and consists of DMSs of mandibular teeth except M3. The third molar was excluded of this sample due to the tempo of development of this tooth, as it is highly variable among modern populations with a high frequency of agenesis. Thus, in a first attempt, we did not include the developmental stage of M3 in the DMSs in order to avoid confounding findings due to this effect.

The first subset is composed of 415 cross-sectional standardized orthopantomographs of Spanish children between 4 and 16 years old. These orthopantomographs were chosen from patients attending to different dental clinics to diagnostic and treatment. The inclusion criteria were the following: a) availability of panoramic radiographs with high clarity and good contrast, b) no systemic diseases or craniofacial anomalities, c) normal dental conditions (e.g. no hypodontia, gross pathology and missing mandibular permanent teeth except third molar) and d) no previous orthodontic treatment. The developmental stages of the different teeth were assesed following the system developed by [START_REF] Demirjian | A new system of dental age assessment[END_REF]. They were first scored by YQ and thus, independently validated by RG-G.

The second subset was derived from the data included in the Electronic Encyclopedia on Maxillo-Facial, Dental and Skeletal Development [START_REF] Demirjian | Dental development, CD-ROM[END_REF]. These data, whose developmental stages were scored by the same author, come from a longitudinal study of French-Canadian orthopantomographs of children in Montreal, conducted in the 1960s and 1970s. A total of 494 radiographs were utilized from this collection. The inclusion criteria for these radiographs was that at least one tooth was developing.

The third subset was selected from a sample of 75 mandibles of subadult individuals from a medieval archaeological population excavated in the Dominican Monastery of San Pablo [START_REF] Ge | Memoria de la actuación arqueológica en el antiguo monasterio de San Pablo[END_REF] that are now housed at the Laboratory of Human Evolution of the University of Burgos. These specimens were scanned using a YXLON Compact (YXLON International X-Ray GmbH, Hamburg, Germany) industrial multislice computedtomography (CT) scanner, located at the University of Burgos. The mandibles were aligned along the long axis of the right mandible corpus. Scanner energy was set at 160KV and 4 mA and the field of view was between 111.1 and 187.5 mm to encompass the entire mandibles. The Mimics™ (Materialise, Belgium) software program was used to visualize the CT images and make the virtual reconstructions. In this sample, the inclusion criteria for these mandibles were that they preserved seven permanent teeth (from I1 to M2) and that at least one of them was developing. In this way, the final sample is composed of 32 mandibles. In this sample, developmental stages of each tooth were scored by RG-G.

Burgos mandibular II: permanent mandibular dentition with M3

The second comparative sample is called "Burgos mandibular II" and consist of 462 crosssectional standardized orthopantomographs of Spanish children in which at least M3 was developing. All the previous 415 cross-sectional orthopantomographs from the first subset of "Burgos mandibular I" are included in this comparative sample. This sample was used to explore the effect of the calcification status of this tooth in the differences and analogies between DMSs of fossil and modern humans. As in the case of "Burgos mandibular I," the developmental stages of the different teeth were first scored by YQ and thus, independently validated by RG-G.

Burgos mandibular III: permanent and deciduous teeth

The third comparative sample is denoted as "Burgos mandibular III" and consists of 24 CTscans from San Pablo collection in which both permanent and deciduous teeth were preserved. The total sample consist of 24 individuals.

Burgos maxillar

The fourth comparative sample, called "Burgos maxillar" consist of 380 orthopantomographs of maxillary teeth chosen from patients attending to different dental clinics to diagnostic and treatment. The inclusion criteria and assessment of the developmental stages of each tooth were performed in the same way that "Burgos mandibular I." This sample was used in comparison of DMSs of fossil specimens with upper dentition preserved. More especifically, it was used to compare with hominin XVIII from Sima de los Huesos and hominins H1 and H3 from Gran Dolina-TD6.2.

Bordeaux sample

All the above samples consist of European or European origin populations. Studies comparing African-derived versus European-derived populations have shown differences in the relative calcification of several teeth (Tompkins, 1996 a). Thus, differences between DMSs of fossil specimens and our reference sample, do not necessarily imply that these DMSs are unique of fossil specimens. Taking this into account, in those cases in which we detected differences between DMSs of TD6 and SH individuals and our reference sample, we performed an additional comparison of these fossils with other modern human sample. To do that, we used a sample from the University of Bordeaux that consists of 2387 children (1346 girls and 1041 boys) aged from 2 to 16 years. Their geographic origins include Southern France, Iran and the Ivory Coast (see [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF] for further details). This reference sample is based on cross-sectional standardized panoramic radiographs of teeth, and children selected were clinically free of anomalies in tooth number, size or shape. Only lower teeth were scored (from I1 to M2), thus, we have only been able to make comparisons between this sample and Sima de los Huesos specimens with lower dentition. Despite these modern human samples are relatively large, we acknowledge that more samples from larger geographic areas and temporal frames should be taken into consideration to have an overall perspective of modern human variation.

Bayesian Statistical Approach

The comparison of the DMSs in the H. antecessor sample, the Sima de los Huesos hominins, and the samples of modern humans was performed using a Bayesian statistical approach [START_REF] Braga | Non-adult dental age assessment: Correspondence analysis and linear regression versus Bayesian predictions[END_REF][START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF], which produces a probability indicating the likelihood that the developmental pattern of a fossil individual may be found within the variation of the modern human population. The underlying hypothesis is that any DMS represents a developmental module. This DMS is composed of hierarchical units, which show varying degrees of interaction. The teeth are, therefore considered as statistically dependent units in the Bayes's rule of conditional probability. As stated previously, the Bayesian approach integrates the concept of modularity, as opposed to tooth-by-tooth analysis [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF].

The interactions between teeth can be measured by deconstructing the original DMS into two subsets with no shared elements. Following this approach, each DMS was separated into a finite number of combinations that correspond to the conditional probabilities of observing one or more teeth at a certain developmental stage (subset_1) conditioned on the attained developmental status of the remaining teeth (subset_2) [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF]. This probability can be expressed as follows:

P(subset_1 | subset_2) = P(subset_2 |



Pprior(subset_1) is the probability that the DMS corresponding to the subset_1 may be found within the references samples.  Pprior(subset≠1) is the probability that the DMS is different to that of the subset_1 may be found within the references samples.

As can be seen in this formula, the posterior probabilities depends on the prior probabilities. These prior probabilities can be calculated as the relative frequencies of this subset_1 in our reference sample. In this way, we were assuming that our reference samples priors are representative of the fossil individuals. However, the references samples were constructed by "availability sampling." It implies that our reference samples exhibit biased distributions of the different dental maturity stages. Thus, posterior probabilities derived this way will tend to over or underestimate in the fossil individuals. In these cases, the best option is to choose priors by using external knowledge independent of the data [START_REF] Couvreur | Insights into the Influence of Priors in Posterior Mapping of Discrete Morphological Characters: A Case Study in Annonaceae[END_REF]. Thus, we have two options to estimate priors: either use demographic data or to assume an unbiased and uniform frequency distribution of the different developmental stages of each teeth [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF]. In this case, we selected the second option, since the posterior probabilities of the different DMSs of fossil specimens will be estimated independent of the distribution of the dental maturity stages in the reference samples used in this study [START_REF] Konigsberg | Estimation of age structure in anthropological demography[END_REF].

Another question related to the Bayesian analysis is that there may be no observations for "subset 2" that are the same as for the fossils. This, in turn, give division by zero. In these cases the particular posterior probabilities have been dropped. For this reasons, we depicted in each case the number of posterior probabilities that could have been calculated.

In addition, we have also calculated the number of times in which the exact fossil sequence appears in the reference sample.

Probabilities were classified based on the thresholds of p<0.25 and p>0.75. While values of 0.25 and 0.75 does not represent an absolute cutoff in a continuous probability distribution ranging from 0 to 1, in Bayesian approaches they represent a threshold to assess the likelihood of different events [START_REF] Braga | Quantifying variation in human dental development sequences: An EVO-DEVO perspective[END_REF]. In this way, probabilities lower than 0.25 indicate very unlike events, probabilities greater than 0.75 represent very likely events and probabilities between 0.25 and 0.75 are more likely associated to random events.

Results

DMSs from H. antecessor from Gran Dolina and the Sima de los Huesos hominins are displayed in Table 1. The H. antecessor individual H11 is the least mature of the TD6.2 sample. Hominin H1 is the most mature individual within the H. antecessor hypodigm, and hominin H3 has an intermediate state of dental maturity. The individuals XVIII and XXV from Sima de los Huesos have completely formed incisors and first molars, with premolars and second and third molars that are still developing. Hominin XVIII is slightly less mature than hominin XXV, as the former's developing teeth are in an earlier mineralization stage.

In fact, the canine roots of the XVIII hominin were still forming when this individual died.

Dental development of H. antecessor

Table 2 shows the distribution of probabilities that the DMSs of these fossil specimens falls within our reference sample from Burgos. The hominin H11 of H. antecessor displays the highest probabilities that its DMS falls within the modern human range, when deciduous teeth are both incorporated or excluded, with values above 0.95 (Fig. 2). In the case of H1, which has both lower and upper dentition preserved, the probabilities, when the M3 is excluded, are between 0.73 and 0.83 respectively. However, when the M3 is incorporated in the analysis, the probability of the upper dentition decreases to 0.60, although the probability of the lower dentition is barely altered. The probability of the upper dentition sequence of the H3 specimen when the M3 is excluded is 0.85. However, when the M3 is included, the probability is 0, as this specific pattern of DMS is not present in the modern human sample.

Within our reference sample the most frequent dental maturational sequence is that showed by H1 lower dentition (excluding M3), and the less frequent is that showed by H3 upper dentition (including M3). They were found 69 and zero times, respectively (7.44% and 0% with respect to the total DMSs, respectively).

Dental development of Sima de los Huesos hominins

Regarding Sima de los Huesos individuals, within our reference sample the maturational sequences of hominin XVIII upper and lower dentition are found 1 and 11 times respectively (0.21% and 1.19% with respect to the total DMSs, respectively) Table 2. For the upper dentition this sequence is found in one male aged 8 years. In the case of the mandibular dentition, we found this DMS in six females aged between seven and nine years and in five males aged between seven and ten years. Interestingly, none of these modern individuals show at the same time the maxillar and mandibular maturational sequence than hominin XVIII. The DMS of hominin XXV is found one time in our reference simple, in one male aged 11 years.

The probabilities of the upper and lower dentitions of the Sima de los Huesos hominins XVIII and XXV (Fig. 3) when the M3 is included are 0, for the same reason stated above. Interestingly, when M3 is removed, the DMS of SH specimens are compared with Burgos's sample, most of probabilities are lower than 0.25. Thus, we can assume that the DMSs showed by hominins XVIII and XXV are unlikely to occur is in this sample.

Nonetheless, when DMSs of SH individuals are compared with Bordeaux's sample, most of probabilities are higher or equal 0.75 (Fig. 3), namely, these sequences are likely in this sample. These differences may be related to inter-populational differences in dental development in Homo sapiens. Based on this, we cannot discard that hominins from SH present a pattern of development similar to that of Homo sapiens.

However, in some combinations the probabilities showed by hominins from SH are lower than 0.25 relative to both modern samples. In the case of hominin XVIII these combinations correspond to the following conditional probabilities: I1I2M1 if CP3P4M2 and CP3P4M2 if I1I2M1. For the hominin XXV these conditional probabilities are:

CM2 if I1I2P3P4M1, P3P4 if I1 I2CM1M2, I1I2CM1M2 if P3P4 and I1I2P3P4M1 if CM2.
Therefore, in these sequences the developmental status of one or some teeth deviate significantly from those observed in our two reference samples. We did not observe probabilities lower than 0.25 in the combinations corresponding to the developmental status of each of these teeth taken separately versus the others. Thus, we concluded that none of these teeth separately show an abnormal developmental status in these hominins from SH.

In the case of hominin XVIII, the more informative DMS is I1I2M1 if CP3P4M2, since this specimen showed I1, I2 and M1 completely formed. When we fixed in our "Burgos mandibular I" sample the sequence C=F, P3=E, P4=E and M2=E, the developmental stages of I2 and M2 are more delayed. More precisely, a 64% of individulas showed the I2 and M2, respectively, in a developmental stage F or G. This could point to a relative advancement of development of I2 and M2 in this SH specimen.

In the case of hominin XXV, the analysis of sequences CM2 if I1I2P3P4M1 and I1I2P3P4M1 if CM2, point to a relative advance in the developmental status of canine and second molar. If we fixed in our "Burgos mandibular I" sample the developmental stages of P3 and P4 as F, a 88% of individuals showed the canine and second mandibular in stages of development more delayed than XXV hominin. The developmental stages of I1, I2 and M1 are equal in our comparative sample then in hominin XXV.

In the other way, the analysis of sequences P3P4 if I1 I2CM1M2 and I1I2CM1M2 if P3P4 show that P3 and P4 are relatively more advanced in our comparative sample. Fixing the developmental stages of I1, I2, C, M1 and M2 a 82% of individulas within our comparative sample showed a developmental stage H or G.

Thus, it seems that SH hominins could be characterized by a relative advancement of development of M1 and M2 and/or a relative delay in the development of both premolars.

Discussion and Conclusion

This study represents the first attempt to apply a Bayesian statistical approach to quantify DMS variation in fossil populations older than Neandertals. Here, we used this method on two extinct hominin populations from the Early and Middle Pleistocene sites at Atapuerca (Spain).

In H. antecessor, the DMS of hominin H1 has a high range of probabilities of belonging to a modern human population. However, the DMS of this hominin shows fully developed incisors, canines, premolars and first molars, all at stage H (Table 1). This prevents the evaluation of differences between anterior and posterior teeth during development, as well as the comparison of the developmental pattern with that of modern humans. The high probabilities compared to modern humans are due principally to the relative development of the second and third molars, which are still forming. It is also very likely that the DMS of the hominin H11 of the same population, which has the permanent dentition from the I1 to the M1, including the two deciduous molars, is within modern human variation. However, as M2 and M3 are not preserved in this specimen, the comparison between anteriorposterior dichotomy also remains incomplete. Finally, H. antecessor H3 does present anterior and posterior teeth that are still forming. In this case, the probabilities vary from 0 to 0.85 depending on the inclusion or exclusion of the M3, respectively. In this regard, H. antecessor H3 suffered from a unilateral impactation of the M2 as a result of the ectopic position of the developing M3, likely due to the lack of space in the maxilla [START_REF] Martín-Francés | Ectopic maxillary third molar in Early Pleistocene Homo antecessor from Atapuerca-Gran Dolina site (Burgos, Spain)[END_REF]. However, it has been documented that there exist significant differences in root growth between unilateral impacted and non-impacted mandibular M2s coming from modern Israelis and Chinese-Americans [START_REF] Shapira | Mandibular second molar impaction. Part I: Genetic traits and characteristics[END_REF]. By expanding these results to their maxillary counterparts, our analysis should be taken with some caution. Anyhow, the relative dental development of the M3 in H. antecessor appears to be developmentally advanced in relation to the M1 and M2 when they are compared with modern humans.

Overall, H. antecessor dental development follows modern human patterns when looking at only anterior or posterior teeth separately, and when M3 is excluded. When both regions are compared, the M3 is advanced in its mineralization in respect to modern humans. However the relative development of I1 and M1 fit within the range of modern humans, which corresponds with previous observations (Bermúdez [START_REF] Bermúdez De Castro | A modern human pattern of dental development in Lower Pleistocene hominids from Atapuerca-TD6 (Spain)[END_REF][START_REF] Bermúdez De Castro | New immature hominin fossil from European Lower Pleistocene shows the earliest evidence of a modern human dental development pattern[END_REF].

In respect to the hominins from the Sima de los Huesos site, different probabilities are produced, depending on the modern sample used for comparison, and on the presence of M3 in the analyses. The DMSs do not correspond to either modern human sample when the M3 is included, indicating that this tooth is advanced in its relative development when compared to modern humans. When the M3 is removed, the probabilities vary depending on which reference sample of modern humans is used for comparison.

The inclusion or exclusion of third molars in the Bayesian analyses offered different probabilities of H. antecessor and Sima de los Huesos hominins to be grouped with modern humans. It is known that mandibular third molars are highly variable in their timing of maturation in modern humans [START_REF] Liversidge | Timing of human mandibular third molar formation[END_REF], and they also have the highest frequency of polymorphism, malposition, impactation and agenesis [START_REF] Nanda | Agenesis of the third molar in man[END_REF][START_REF] Garn | Third molar agenesis and size reductions of the remaining teeth[END_REF][START_REF] Anderson | Evolutionary dental changes[END_REF][START_REF] De Castro | Third molar agenesis in human prehistoric populations of the Canary Islands[END_REF]. Moreover, statistical differences within modern humans were found in the delay of M3 formation in Caucasian children from both London and Cape Town in comparison to black South African children [START_REF] Liversidge | Timing of human mandibular third molar formation[END_REF]. This delay is also present in a French-Canadian population compared to black South Africans, in both the second and third molars (Tompkins, 1996 a).

In this regard, the delay of molar formation in non-black Africans compared to black Africans might be responsible for the different probabilities observed in the Sima de los Huesos lower dentition when their mineralization sequences are compared with modern human variation. The sample from the University of Bordeaux includes individuals of black sub-Saharan provenance (Ivory Coast), whereas the sample from Univerity of Burgos does not include any. This may be the reason for the apparent similarity between the Burgos sample and the Sima de los Huesos individuals, and the apparently advanced lower molar developmental sequences of the ancient individuals in comparison to the Bordeaux sample.

As the Burgos sample was used when calculating the probabilities in H. antecessor and the hominins from Sima de los Huesos when the M3 is removed, the fact that the Early Pleistocene population has higher probabilities than the Middle Pleistocene indicates that the M2 is relatively advanced in its development in the Sima de los Huesos individuals, as it was previously postulated for the XVIII specimen (Bermúdez de Castro and [START_REF] Bermúdez De Castro | Pattern of dental development in Hominid XVIII from the Middle Pleistocene Atapuerca-Sima de los Huesos site (Spain)[END_REF].

It is remarkable the pattern of development between upper and lower third molars. For instance, it has been observed in a sample of white Americans that the formation of maxillary M3s was slightly advanced in comparison to the mandibular M3s [START_REF] Mincer | The ABFO study of third molar development and its use as an estimator of chronological age[END_REF]. However, this pattern is not shared in either the Sima de los Huesos hominins or the sample of H. antecessor. Both M3s of hominin H1 of the Gran Dolina-TD6.2 site are at the same stage of mineralization (C), whereas the maxillary M3 of the hominin XVIII from Sima de los Huesos is delayed in its formation in comparison to the mandibular M3 (B and C, respectively). This fact could be key to differentiate European Lower and Middle Pleistocene hominins among them and also in respect with modern humans.

Leaving aside the analysis of third molars on their own, the lower first and second molars of both the H. antecessor H1 and the hominin XXV of Sima de los Huesos are at stages H and G, respectively (Table 1). Differences between the hominins from these sites emerge when the premolars and third molars are considered. Premolars in the hominin XXV from Sima de los Huesos are developmentally delayed in comparison to hominin H1 of H. antecessor (stages F-F and H-H, respectively). The third molar, however, is developmentally advanced in XXV in comparison to hominin H1 (stages D and C, respectively). Caucasian French-Canadian, Native American and black African modern human samples do not display differences in their relative premolar calcification across these three groups (Tompkins, 1996 a), so differences observed in both Pleistocene hominins could be explained by a taxonomic signal. Taking into account third premolars, it was stated that P3 in Neandertals was delayed in comparison to modern humans [START_REF] Tompkins | Relative dental development of Upper Pleistocene hominids compared to human population variation[END_REF], whereas in Sima de los Huesos hominins both premolars appear to be delayed.

Neandertals are key to be used as comparison to H. antecessor and Sima de los Huesos, not only because they share the same geographical area, but also because typical Neandertal features were identified in both human populations at Atapuerca Bermúdez- [START_REF] Bermúdez-De-Castro | Homo antecessor: The state of the art eighteen years later[END_REF]. Following the Bayesian statistical approach, Neandertals display probability values that exclude them from belonging to modern humans, as shown in the Roc de Marsal (Bayle et al., 2009 c) and the Cova del Gegant [START_REF] Quam | The Neandertals of northeastern Iberia: New remains from the Cova del Gegant (Sitges, Barcelona)[END_REF] specimens. The Cova del Gegant individual is dated in 52.3 ± 2.3 ka [START_REF] Daura | Stratigraphic context and direct dating of the Neandertal mandible from Cova del Gegant (Sitges, Barcelona)[END_REF], whereas the Roc de Marsal Neandertal is located most probably between 60 and 70 ka [START_REF] Guérin | Multi-method (TL and OSL), multi-material (quartz and flint) dating of the Mousterian site of Roc de Marsal (Dordogne, France): Correlating Neanderthal occupations with the climatic variability of MIS 5-3[END_REF]. In both cases, the Bayesian probability is zero, which means that their DMSs are not present in the modern human reference samples employed. In addition, the Roc de Marsal child has a relatively advanced stage of mineralization of the first molar in respect to the comparatively delayed maturation levels of its incisors (Bayle et al., 2009 c). Bearing in mind that the relative development of I1 and M1 in H. antecessor, and likely that of the Sima de los Huesos hominins, is within modern human variation, this asynchrony in the Neandertals I1/M1 relative development could be interpreted as exclusive to H. neanderthalensis. Although this observation was made only on the Roc de Marsal Neandertal, a claim for caution is made, as this feature could be key to differentiating Upper Pleistocene Neandertals from their ancestors in the European Middle Pleistocene. Interestingly, the Neandertal specimen Spy VI, represented only by four deciduous teeth (lower i1, i2, c; and upper i1), shows a Bayesian probability above 0.75 in respect to modern humans [START_REF] Crevecoeur | The Spy VI child: A newly discovered Neandertal infant[END_REF]. However, this high probability is expected as it is only based on deciduous incisors and one canine, therefore, only anterior teeth and not anterior vs posterior dentition. Finally, The Lagar Velho 1 child, a potential hybrid specimen between Neandertals and anatomically modern humans, with an age of ~24.5 ka B.P. [START_REF] Duarte | The early Upper Paleolithic human skeleton from the Abrigo do Lagar Velho (Portugal) and modern human emergence in Iberia[END_REF], also presents a dental maturation pattern not represented in the modern human variation [START_REF] Bayle | Dental maturational sequence and dental tissue proportions in the early Upper Paleolithic child from Abrigo do Lagar Velho, Portugal[END_REF].

It is remarkable that the relative relationships between anterior and posterior teeth were found to be responsible for differences in extinct hominins [START_REF] Tompkins | Relative dental development of Upper Pleistocene hominids compared to human population variation[END_REF]. This could explain why Neandertals that preserve both types of teeth and the Lagar Velho specimen have a pattern of dental development not present in modern humans, whereas Spy VI, only represented by anterior dentition, shows higher probabilities. In contrast, the Upper Paleolithic child from La Madeleine (LM4), a fully anatomically modern human child with an age of 10,190 ± 100 years [START_REF] Gambier | Datation de vestiges humains présumés du Paléolithique supérieur par la méthode du Carbone 14 en spectrométrie de masse par accélérateur / Accelerator mass spectrometry radiocarbon dates of human remains from Upper Palaeolithic[END_REF], shows comparatively higher probabilities of belonging to modern humans (Bayle et al., 2009 b). In particular, 30% of LM4 DMSs probabilities are superior to 0.75, and 70% fall between 0.25 and 0.75.

El Sidrón J1 Neandertal preserves the permanent dentition from I1 to M2 and some deciduous teeth [START_REF] Rosas | The growth pattern of Neandertals, reconstructed from a juvenile skeleton from El Sidrón (Spain)[END_REF]. Only the M3 crypt is present so that its relative stage of mineralization with respect to the modern human sample employed in this study remains unknown. However, the DMS of this specimen is not found in any of the 10901 modern humans radiographs used to calculate probability density plots. From a non-Bayesian approach, probability density plots for mean age of transition entering each individual permanent tooth stage indicate that El Sidrón J1 fell within the modern human range [START_REF] Rosas | The growth pattern of Neandertals, reconstructed from a juvenile skeleton from El Sidrón (Spain)[END_REF]. However, premolar and molar plots show slightly older, advanced, ages than anterior dentition plots in respect to the modern human sample, although a high degree of overlap exists among both fields. [START_REF] Smith | Dental evidence for ontogenetic differences between modern humans and Neanderthals[END_REF], then M2 and M3 are just two stages apart and comparable to individual XXV from Sima de los Huesos. This both fits with the range of more advanced or delayed chronological ages determined for Neandertals at various stages of development [START_REF] Macchiarelli | How Neanderthal molar teeth grew[END_REF][START_REF] Smith | Dental evidence for ontogenetic differences between modern humans and Neanderthals[END_REF][START_REF] Rosas | The growth pattern of Neandertals, reconstructed from a juvenile skeleton from El Sidrón (Spain)[END_REF] but also underscores the likely developmental overlap between the TD6.2, SH and Neandertal hominins.

It is worth mentioning that the hominins from Sima de los Huesos and H. antecessor are older in their ages-at-death than Gegant-5, Roc de Marsal 1, Lagar Velho 1 or La Madeleine 4. Although the previous comparisons are directly observed and measured, we acknowledge that the comparisons were not made in similar ages-at-death.

In a broader perspective, dental development is not only assessed by approaching the relative development of the teeth, but also by taking into account the time of formation of their dental tissues. In this regard, lateral enamel formation times of the entire dentition in both H. antecessor and the Sima de los Huesos hominins was ~27% shorter than in modern humans [START_REF] Modesto-Mata | Short and long period growth markers of enamel formation distinguish European Pleistocene hominins[END_REF]. Despite the high probability of some dental mineralization sequences of H. antecessor and Sima de los Huesos in respect to modern human variation, the fact that both populations display an advanced molar development and more rapid enamel formation times, distinguishes them from H. sapiens.

These two different lines of evidence shed some light on the processes of growth and development in H. antecessor and the Sima de los Huesos hominins. They support a working hypothesis that both Pleistocene populations had a shorter period of general growth, ontogeny and skeletal maturation than modern human average, being placed most likely at the lower tail of modern human distribution today. Although more data are needed to further test this hypothesis, such as accurate estimations of the timing and rate of root formation, cuspal enamel formation times and ages at eruption in these hominins, some findings about skeletal development seem to confirm it. For instance, clavicular growth and development in H. antecessor is faster than in modern humans [START_REF] García-González | Étude analytique d'une clavicule complète de subadulte d'Homo antecessor (site de Gran Dolina, Sierra d'Atapuerca, Burgos, Espagne)[END_REF], a trait shared with H. neanderthalensis and H. ergaster. Some evidence points to a slower height growth rate in Neandertals during infancy and early childhood compared to modern humans, which would explain differences in adult height between these populations [START_REF] Martín-González | Differences between Neandertal and modern human infant and child growth models[END_REF]. However, this fact explains differences in the growth rate models but not in the overall timing of maturation. In modern humans, dental development and skeletal growth are moderately correlated and thus, individuals that are dentally advanced relative to their peers also tend to be skeletally advanced (Šešelj, 2013). Thus, if we can assume the same relationship between dental and skeletal development in H. antecessor, Sima de los Huesos and modern humans, these findings, along to those of this study, suggest an advanced development in H antecessor and Sima de los Huesos hominins, probably reaching adulthood at around 14-15 years of age.

In conclusion, both H. antecessor and the Sima de los Huesos hominins show advanced development of the M3 in comparison to modern H. sapiens. However, the Sima de los Huesos hominins appear also to show advanced development of the M2 in respect to both H. antecessor and modern humans. When anterior and posterior dentitions of H. antecessor are compared independently with modern humans, they present high probabilities of statistical variation. However, when both anterior and posterior teeth are compared synchronously, their probabilities decrease, indicating differences between the developing anterior and posterior dentitions between modern and extinct human populations. The advanced developmental sequence of M2 and M3, along with the more rapid enamel formation times, points to a shorter ontogeny in these hominins compared to modern humans.

Tables and Figures

Table 1: Mineralization stages of H. antecessor (TD6.2) and Sima de los Huesos (SH) teeth following stages defined by [START_REF] Demirjian | A new system of dental age assessment[END_REF] for permanent dentition and [START_REF] Liversidge | Variation in crown and root formation and eruption of human deciduous teeth[END_REF] for deciduous dentition**. These stages were employed in the Bayesian statistical approach to determine whether their sequences belong to H. sapiens. brk = broken.

Site Specimen Position I1 I2 C P3 P4 M1 M2 M3 dm1 dm2

  subset_1) * Pprior(subset_1) / [P(subset_2 | subset_1) * Pprior(subset_1)] + [P(subset_2 | subset≠1) * Pprior(subset≠1)] where:  P(subset_2 | subset_1) is the observed proportion of individuals in the reference samples showing the DMS corresponding to the subset_2 given the DMS corresponding to the subset_1.  P(subset_2 | subset≠1) is the observed proportion of individuals in the reference samples showing the DMS corresponding to the subset_2 given a DMS different to that of subset_1.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Buccal view of the complete dentition of the Sima de los Huesos hominin XVIII. Top row: upper dentition; bottom row: lower dentition. dm2s are placed above the crown of their respective P4s. Top left: upper right M3; bottom right: lower left M3. Scale bar = 1 cm.

Fig. 3 :

 3 Fig.3: Bayesian probabilities of the Sima de los Huesos (SH) dental sequences in respect to the DMS variation in modern humans. Two hominins are represented: XVIII and XXV. The M3 has not been included in the calculation of the Bayesian probabilities for both specimens. A) lower dentition of hominin XXV; B) lower dentition of hominin XVIII; C) upper dentition of hominin XVIII. Green bar: probabilities equal to or higher than 0.75; red bar: probabilities equal to or lower than 0.25. Red lines and dots: probabilities calculated by using the modern human reference sample from the University of Burgos; blue lines and dots in respect to the reference sample from University of Bordeaux.

Table 2 :

 2 Results of Bayesian analysis of the dental mineralization sequences of H. antecessor (TD6.2) and Sima de los Huesos (SH) hominins. Arc: arch (U, upper; L, lower); DMS: dental maturational sequence (I, DMS excluding M3; II, DMS sequence including M3; III, DMS including deciduous teeth). N: comparative sample size; NTC: number of combinations; NS: number of times the fossil DMSs are found within our modern human reference sample; NS (%): NS expressed in percentage with this formula (NS*100)/N.

	TD6 H1	Lower	-H H H H H	G	C	-	-
	TD6 H1	Upper	--H H -	H	G	C	-	-
	TD6 H3	Upper	-G F E E H	E	B	-	-
	TD6 H11	Lower	D D C B A E	-	-	H2	H1
	SH	XVIII	Lower	H H F E E H	E	C	-	brk
	SH	XVIII	Upper	H H F E E H	E	B	-	brk
	SH	XXV	Lower	H H H F F	H	G	D	-	brk

Acknowledgments

Funding for this study comes from the Spanish Ministry of Science/FEDER, grant numbers CGL2012-38434-C03-01/02/03, CGL2014-52611-C2-1-P, CGL2015-65387-C3-2-P, CGL2015-65387-C3-3-P, PGC2018-093925-B-C31, PGC2018-093925-B-C33, HAR2014-55131, HAR2017-86509-P), Acción Integrada Francia-España (HF2001-0115), Consejería de Educación de la Junta de Castilla y León (CEN074A12-2), The Leakey Foundation (through the personal support of Gordon Getty 2013 and Dub Crook 2014, 2015, 2016, 2018 to author MM-T). MMP and LMF are beneficiaries of postdoctoral grants by the Atapuerca Foundation and IdEx, respectively. CG-C, MM-M and MMP were funded by doctoral grants from the European Social Fund (BOCYL-D-02102014-10 and BOCYL-D-20052013-14) through the Consejería de Educación de Castilla y León (Spain) with the European Social Fund. Most of the data acquisition was performed at the Laboratory of Microscopy of the CENIEH-ICTS (Spain) in collaboration with CENIEH staff. Restoration and conservation work on the material was carried out by Pilar Fernández-Colón from the Conservation and Restoration Area of CENIEH-ICTS, Lucía López-Polín from IPHES and Maicu Ortega from the Centro de Evolución y Comportamiento Humanos UCM-ISCIII. J.M. Álverez-Tovar, M. Esteves Dias, A. Castillo Fernández and José Miguel Carretero provided the orthopantomographies of modern children. Marina Lozano, Ignacio Martínez, Cáceres Museum and Extremadura Government helped in the assessment of the scientific material. Without the persistent effort of the Atapuerca Research Team this work would have not been possible. Descriptions in the paper use the taxonomy provided by CRediT. We describe contributions to the paper using the taxonomy provided by CRediT. Conceptualization: MM-M, JMBC, MCD; Formal analysis: MM-M, RG-G, YQ; Investigation: MM-M, RG-G, YQ; Methodology: MM-M, RG-G, YQ; Resources: MCD, JMBC, EC, JLA; Software: MM-M; Supervision: JMBC, MCD; Visualization: MM-M; Writing-original draft: MM-M;

Writing -review & editing: MM-M, MCD, RG-G, YQ, YH, CG-C, MMP, LMF, MM-T, EC, JLA, JMBC. We are especially grateful to the referees and Associated Editor for their comments and suggestions, which have greatly improved the paper.

Data availability statement

The dental mineralization sequences (DMS) of the two Pleistocene populations of Atapuerca are available in the article. Modern human DMSs may be available from Rebeca García (Burgos sample), Yuliet Quintino (Burgos sample) and Yann Heuzé (Bordeaux sample) upon reasonable request. All other data to support the findings of this study are available in the supplementary material of this article.

Conflict of interest

The authors declare no conflicts of interest.

Supplementary Materials Supplementary Tables

Supplementary Table 1: Teeth from the two Sima de los Huesos hominins employed in the calculations of the Bayesian probabilities in respect to the DMS variability in modern humans. Hominins: XVIII and XXV. Position (U = upper; L = lower); side (L = left; R = right).

Supplementary Figures

Supplementary Fig.