Increasing control over biomineralization in conodont evolution

Domain: earth_science
Published: 2024-06-20
Editorial outcome: minor_revision

Abstract
Vertebrates use the phosphate mineral apatite in their skeletons, which allowed them to develop tissues such as enamel, characterized by an outstanding combination of hardness and elasticity. It has been hypothesized that the evolution of the earliest vertebrate skeletal tissues, found in the teeth of the extinct group of conodonts, was driven by adaptation to dental function. We test this hypothesis quantitatively and demonstrate that the crystallographic order increased throughout the early evolution of conodont teeth in parallel with morphological adaptation to food processing. With the c-axes of apatite crystals oriented perpendicular to the functional feeding surfaces, the strongest resistance to uniaxial compressional stress is conferred along the long axes of denticles. Our results support increasing control over biomineralization in the first skeletonized vertebrates and allow us to test models of functional morphology and material properties across conodont dental diversity.

Introduction
Vertebrates are one of the few lineages of skeletonized animals that form their skeletal tissues out of the phosphate mineral apatite1 and the one that has achieved the highest disparity. Composite organic-apatite biomineral forming vertebrate skeletons range through internal bone scaffolds, body armour, teeth, and scales2,3. These versatile, remodelable tissues allowed a diversity of new feeding and protective structures, which triggered an evolutionary escalation of complex trophic interactions4. All other major, widespread animal groups build their hard tissues from softer minerals: calcite, aragonite, or robust organic polymers1. Vertebrates stand out in terms of the mineralogy of their skeletons and the disparity of structures used for feeding, locomotion, and defense. This disparity coincides with the outstanding evolutionary and ecological success of this group.

The earliest phosphatic skeletal tissues found in the evolutionary history of vertebrates are conodont teeth, which evolved in parallel to all other vertebrate teeth. Donoghue5 put forward the hypothesis that the biomineral structure and properties of conodont teeth evolved as an adaptation to dental function. This hypothesis is the basis of our approach here: hypotheses on the function of a skeletal tissue can be tested by predicting its material properties from its crystallographic structure. However, conodont skeletons are very heterogeneous, and in situ studies relating crystal orientations to individual functional surfaces have been very limited. We offer a quantitative test of the hypothesis by Donoghue5 by quantifying the degree of an organism’s control over the crystallographic order in conodonts.

Conodonts are paramount to understanding the origins of skeletal evolution thanks to their long evolutionary history and morphological diversity6,7. These marine vertebrates bore apparatuses composed of morphologically and functionally differentiated teeth (also referred to as elements). These elements consisted of composite materials formed by organic matter and hydroxyapatite (Ca10(PO4)6(OH)2)8–10. Such composite materials also form vertebrate bones, teeth, and mineralized scales. Conodont elements typically range from 0.1 mm to 2 mm in length, with rare giants reaching over 13 mm11. The morphology of their elements ranges from simple cone-like forms to molar-like structures with complex topography allowing accurate occlusal (mastication)12,13.

The oldest conodonts are represented by the paraphyletic grade of paraconodonts, which appeared in the Cambrian Period and were characterized by poorly mineralized phosphatic tissues rich in organic matter. Paraconodonts were succeeded by the more derived group, euconodonts, which dominated the lineage from the late Cambrian (ca. 490 Ma) until the extinction of conodonts at the end of the Triassic Period (ca. 200 Ma). Euconodonts were distinguished by hypermineralized crown tissues, which formed on top of less mineralized basal bodies with homologous structures to paraconodont elements. Crown tissues formed the functional part of the dental organs, which participated in the mechanical digestion of prey. Euconodont crown consists of two types of tissues: the white matter, unique to conodonts, and the enamel-like hyaline tissue7,14. The hyaline tissue is widely used by geochemists for paleoclimatic reconstructions thanks to its outstanding resistance to diagenesis15. This tissue grew by apposition throughout the life of the organism. Its complex developmental mechanism permitted conodonts to achieve enormous morphological diversity16,17. The known range of morphologies reflects adaptation to a broad range of trophic niches. For some conodont lineages, a trend towards developing complex surface topographies for stress dissipation and efficient food use has been reconstructed18,19. Other taxa evolved unparalleled sharpness interpreted as a mechanism of efficient food shearing and slicing20. Trophic diversification of conodonts is supported by distinct, species-specific geochemical values of proxies for trophic partitioning21.

In this study, we reconstruct the selective pressure that acted upon the earliest vertebrate biominerals and shaped early conodont evolution by examining the link between their ultrastructure and functional morphology.

Here we use Electron Backscatter Diffraction (EBSD) to characterize the crystallographic preferred orientation, in the following referred to as crystallographic texture, of conodont hyaline tissue across a range of morphologies and adaptations to food processing. Based on the textural development through conodont evolution, we infer structural adaptations using criteria linking crystallographic texture with resistance to compressional stress. The c-axis of an apatite crystal is the most resistant to uniaxial compressional stress22 and, thus, the most resilient to breakage when uniaxial compression is applied (Fig. 1). Consequently, a stronger c-axis alignment would disperse stress throughout the tooth better and allow for a stronger bite and reduced risk of breakage23,24.

Results
Crystallographic texture
The orientation of c-axes are mostly parallel to the long axes of the denticles and cusps. In the earliest diverging taxon, Pro. muelleri, the a-axes <11–20> show a girdle distribution resulting from gradual rotation around the c-axis, visible as a color gradient in Fig. 2. A similar distribution but with a lower dispersion around the c-axis was found in the element of the more derived coniform conodont, Pan. equicostatus (Fig. 2b). The range of rotation around the c-axis decreases systematically within each specimen (see further for the quantification), leading to a constrained distribution with the a-axes of crystals showing very little dispersion in W. excavata (Fig. 2c), T. gracilis (Fig. 2e) and the most derived taxon, Palmatolepis sp. (Fig. 2f). In B. cf. aculeatus (Fig. 2d), the crystals formed three distinct clusters (of the a-axis) in terms of their angles of rotation around the c-axis.

In biominerals, grains are areas with similar crystal orientations and appear in the same color on a crystallographic orientation map, even though they may consist of multiple nanocrystals36. Grain sizes can be used to characterize a material, and they impact its properties. Here, however, grain sizes have not been calculated, as this operation requires grain reconstruction following algorithms that rely on constant misorientation angles between adjacent grains. We demonstrate (Supplementary Fig. 2) that the misorientations vary between areas within one tissue depending on their position with respect to functional directions and the compressional and shear stresses acting in these directions. Thus, standard approaches to grain analysis, which rely on constant misorientation angle thresholds, may not be suitable for complex biominerals, where these thresholds changed during evolution and even within one tissue type. Furthermore, in abiotic crystalline materials, the long axes of grains commonly correspond to the crystallographic c-axis of the mineral. This is not always the case in biological tissues37. Although grains have not been computed using an algorithm, elongated clusters of points with similar crystallographic orientations can be distinguished here, especially in Pro. muelleri. In this specimen, in areas where such clusters are visible, their long axes do not correspond to crystallographic c-axes <0001>.

Spread of orientations
A decreasing tendency in the spread of orientations is visible in pole figures from the least derived taxon Pro. muelleri (Fig. 2a) to the most derived ones, T. gracilis and Palmatolepis sp. Pro. muelleri stands out with the lowest TI values ranging from 4 to 48 (mean 12, Fig. 3a). All five taxa with fully developed crown tissues overlap (Fig. 3), but the coniform Pan. equicostatus occupies a much smaller range of TI values (17 to 92, mean 29), compared to “complex” conodonts, B. cf. aculeatus (24 to 71, mean 36) and W. excavata (9 to 67, mean 35). Taxa with elements forming wide platforms, considered highly adapted to mastication18, have the highest TI values of all: T. gracilis (38 to 75, mean 57) and Palmatolepis sp. (13 to 92, mean 68).

In heterogeneous materials such as biogenic apatite, measures of spread depend on the size of the analyzed area (Supplementary Fig. 3). To eliminate sampling bias, random subsets were drawn from the measured areas and, for each subset, ODF and TI were calculated (Fig. 4). The smaller range of values for Pan. equicostatus cannot be explained by the smaller area of the element exclusively, as subsampling from other euconodont taxa yielded higher TI values for comparable areas (Fig. 4).

Because TI is quite sensitive to the area (or number of EBSD data points), we evaluated the robustness of our conclusions by repeating the analysis using another measure of the strength of preferred orientation proposed in the literature, the M-index38,39. It is calculated from the misorientation angles between uncorrelated grain pairs. Additionally, we evaluated the sharpness of the texture in the directions of c- and a-axes using the pole figure texture index pfJ31 (Fig. 4, Supplementary Fig. 4, and Supplementary Table 1). TI followed by pfJ showed the best sensitivity to differences in texture between our samples and allowed resolving the differences best.

Evaluation of potential diagenetic influence
Peak positions of the ν-PO4 3− band obtained through Raman analyzes overlap between studied conodonts (Fig. 5a). They range from 963.67 to 964.26 cm−1 for the band position and from 4.28 to 5.67 cm−1 for the full width at half maximum (FWHM; Fig. 5b)28. Palmatolepis sp. is the only exception, with lower values for both band position (mean 963.88 cm−1) and FWHM (mean 4.49 cm−1, both n = 5).

Diagenetic alteration such as uptake of elements post-mortem has been reported to proceed on the outer surface of the elements40. Therefore, we compared also the CAI values and distance from the center in each specimen against the position of the ν-PO4 3− band. No relationship between these variables and ν-PO4 3− band position was found (Supplementary Fig. 5).

Zhang et al.15, noted that higher values for the band position (lower redshifts) correspond to samples with little diagenetic alteration. Thus, if the specimens studied here had experienced different degrees of diagenesis, a high spread of the band position would be expected. Figure 5 shows the relationships between the peak center at maximum intensity (PCMI) and FWHM of the ν-PO4 3− band in our specimens compared with those reported by McMillan and Golding41, Rantitsch et al.42 and Zhang et al.15, which were the only published datasets available at this point.

In all datasets, the FWHM of the peak decreased with an increasing position of the band center (Fig. 5). The datasets by both McMillan and Golding41 and Rantitsch et al.42 showed high variability in the position and widths of the ν-PO4 3− band. Model II regression was fitted to the relationship for these two datasets, producing slopes equal to −0.92 (R2 = 0.63, N = 41) and −1.30 (R2 = 0.90, N = 153), respectively (Fig. 6). Both slope and intercept differed significantly between these two datasets at a confidence level of 0.05. In contrast, our dataset and that by Zhang et al.15 showed strong clustering close (Zhang et al.) or directly on the regression line fitted to the values measured by Rantitsch et al.42. Compared to the other studies, positions of the ν-PO4 3− band measured for our specimens showed consistently high values with a mean of 964.00 (median 963.99, standard deviation 0.13, n = 30).

Methods
Specimens
All specimens are stored in the collections of Utrecht University (contact person E. Jarochowska), except for Tripodellus gracilis, which is stored in the collections of the Institute of Paleobiology, Polish Academy of Sciences. Elements prior to sectioning are shown in Supplementary Fig. 1. For anatomical notation see Purnell et al.68 and for the phylogenetic position – Donoghue et al.55.

Proconodontus muelleri. A coniform aequaliform element (for anatomical notation see Müller, 1973) from the upper part of the Windfall Formation, at a section in Ninemile Canyon, Antelope Range, Eureka County, Nevada, USA (Supplementary Fig. 1A). It is from U.S. Geological Survey Cambrian-Ordovician sample locality USGS 12122-CO and is being stored under the accession number D000.000.389.576. The sample was collected on public land, and no permit is required for collection or export. The section was originally measured by J. Repetski, M. E. Taylor, and R. J. Ross, Jr. (all of USGS). The sample was collected by J. D. Loch (Central Missouri Univ.) and J. F. Taylor (Indiana University of Pennsylvania) and processed by J. Repetski. The conodont element was extracted from lime grainstone by processing in approx. 10% buffered acetic acid, sieved to 170 mesh (U.S. standard sieve size); the residue was not processed through heavy liquids. Age: late Cambrian, Stage 10; upper part of the Eoconodontus Biozone69.

Panderodus equicostatus. A coniform truncatiform element (for anatomical notation see Sansom et al., 1994) from Homerian (middle Silurian) Ternava Formation at Vrublivtsy, Ukraine (48°3703.90′′ N; 26°4602.63′′ E, Supplementary Fig. 1B). According to Ukrainian law, natural objects of value for the state’s cultural and natural heritage, such as rare fossils, are prohibited from export, whereas common rock fragments and fossils are not included in this definition and can be collected on public land and exported without a permit. The sampled rock fragment was labeled V-19.25, for details of age and the locality, see ref. 70. The rock sample was processed at GeoZentrum Nordbayern (Erlangen, Germany): it was dissolved in 7% buffered acetic acid according to the method of Jeppsson & Anehus (1995)71, wet sieved, and the 63 µm – 1 mm fraction was separated in sodium polytungstate. Accession number D000.000.389.574.

“Complex conodonts”: Bispathodus cf. aculeatus. A P1 element of a polygnathid species from the scree in the natural outcrop of Famennian (Upper Devonian) limestone at Köstenhof, north-west of the hamlet of Elbersreuth in Frankenwald, southern Germany (Tragelehn & Hartenfels72; Supplementary Fig. 1C). It was collected on public land. No permits are required for collecting and exporting common fossils from Germany73. The rock sample was processed at GeoZentrum Nordbayern (Erlangen, Germany): it was dissolved in 7% buffered acetic acid according to the method of Jeppsson & Anehus (1995)71, wet sieved, and conodonts were picked from the 63 µm – 2 mm fraction. Accession number D000.000.389.572.

“Complex conodonts”: Wurmiella excavata. An M element of an ozarkodinin species collected from the Ludfordian (upper Silurian) of Gotland, sample no. G14-19OB (Supplementary Fig. 1D). The rock piece from which the conodont element was extracted was picked up by Oskar Bremer from loose beach blocks at Barshageudd 2, located at the southernmost tip of Gotland (Swedish grid coordinates WGS84: 56°; 54ʹ 20.0” N, 18° 11ʹ 16.0” E). The collected material were loose scree blocks and the collection site is located on public land, thus no permit was required. The blocks were dissolved in buffered acetic acid at the Department of Geology at Lund University, according to the technique of Jeppsson & Anehus (1995)71. The phosphatic remains were separated using heavy liquid separation. For details of age and the locality, see Bremer et al.74. Accession number D000.000.389.573.

“Platform bearing complex conodonts”: Tripodellus gracilis. A P1 element extracted from an upper Famennian sample collected at the Kowala Quarry in the Holy Cross Mountains (central Poland). The sample number Ko-338 was located just above the Kowala Shale, which forms a distinct marker horizon in the quarry. It was collected in a private quarry with permission. For details of age and the locality, see Ziegler and Sandberg75. The rock fragment was processed at the Institute of Paleobiology, Polish Academy of Sciences, in Warsaw: it was dissolved in 10% acetic acid, and the insoluble residue was enriched using an electromagnetic separator. The specimen (Supplementary Fig. 1E) is stored in the collections of the Institute of Paleobiology, Polish Academy of Sciences, under the accession number ZPAL C16/3129.

“Platform bearing complex conodonts”: Palmatolepis sp. A P1 element of a palmatolepid species (Supplementary Fig. 1F) extracted from a limestone sample collected from the base of the Upper Kellwasser Horizon (Frasnian) at the now disused Schmidt Quarry located on public land outside of the hamlet of Braunau in the western part of the Rheinisches Schiefergebirge (central Germany). No permits are required for collecting and exporting common fossils from Germany73. The rock sample was processed at GeoZentrum Nordbayern (Erlangen, Germany): it was dissolved in 7% buffered acetic acid according to the method of Jeppsson & Anehus (1995)71, wet sieved, and conodonts were picked from the 63 µm – 2 mm fraction. For details of age and the locality, see Sandberg et al.76. Accession number D000.000.389.575.

Preparation and EBSD data acquisition
Electron backscatter diffraction (EBSD) is a scanning electron microscopy (SEM) based analytical technique that can be used to quantify crystal properties such as size, shape, and orientation in situ in direct comparison to the functional surface of the element. Owing to the bioapatite’s high sensitivity to damage by the electron beam, previous studies faced beam damage to such an extent that crystal orientations could not be analyzed47,77. Thanks to improved sample preparation78, the creation of an extremely flat surface, and access to highly sensitive CMOS-EBSD detectors79, we were able to overcome this limitation. The samples were suspended in epoxy resin and polished using a succession of 6, 3 and 1 µm diamond suspension (for full protocol see Shirley et al.78). They were mounted on SEM stubs and coated with 3.5 nm of carbon coating. Backscatter electron images can be accessed in Supplementary Data 1 and source EBSD files are provided by Shirley et al.80.

The initial data (Bispathodus cf. aculeatus) were collected using a Hitachi SU70 FEG SEM. Pro. muelleri, Pan. equicostatus and W. excavata EBSD data was acquired at the Department of Werkstoffwissenschaften WW1 of Friedrich-Alexander-Universität Erlangen-Nürnberg on a Helios NanoLab 600i DualBeam. T. gracilis, Palmatolepis sp. maps were collected at the Electron Microscopy Centre of Utrecht University using a Zeiss Gemini 450. All SEMs were fitted with an Oxford Instrument Symmetry EBSD detector using the AZtec software. This detector combines an optimized CMOS sensor with fiber optics, delivering a significant sensitivity increase compared to previous-generation EBSD detectors. This sensitivity, when combined with meticulous sample preparation, allowed us to lower the electron dose required for each diffraction pattern, thus minimizing tissue damage by the electron beam. Acquisition parameters were adjusted individually to maximize indexing rates (Supplementary Table 2). Indexing rates were between 57% and 87% for sample areas shown in Fig. 2, which is low compared to abiotic minerals, but very high for beam-sensitive materials. The dominant phase detected was apatite with lattice parameters listed in Supplementary Table 2 with an admixture of hydroxyapatite. Hydroxyapatite was converted to apatite for subsequent analysis.

Raman spectrometry
All six sectioned conodonts were analyzed with Raman spectrometry using a WiTec alpha300 Raman microscope equipped with a 532 nm laser and a grating of 1800 grooves per mm−1. The system was calibrated using a white light prior to use, and instruments' spectral shifts were monitored using a silicon wafer standard. A 50× long working distance lens from Zeiss was used for the analysis, which provided a lateral spot size of 1 µm and approximately 2 µm depth penetration. The laser power was lowered to 2 mW to shield the samples from destructive effects of the Raman laser. Five Raman spectra were acquired per specimen beginning at the centre of each denticle moving towards the outer edge (i.e. occlusal surface) in 5 µm steps. Each spectrum was acquired for 10 seconds and averaged over 10 acquisitions to optimize the signal to noise ratio. Peak center at maximum intensity (PCMI) and full width at half maximum (FWHM) of the ν-PO4 3− band in the Raman spectrum of apatite were calculated using a pseudo-Voigt function in Project FIVE 5.2 software by WiTec.

Raman spectra, band positions, and widths (PCMI and FWHM, respectively)28 were extracted from the Project FIVE 5.2 software and plotted using R Software 4.3.0 (Fig. 4). The ν-PO4 3− band parameters were then compared with those published in other studies15,41,42. Ranged major axis regression of FWHM over PCMI was calculated for each study using the package lmodel239,81.

Quantification of crystal texture
EBSD data was processed in the MTEX 5.9.0 toolbox for MATLAB82–84 using code provided by Shirley et al.39. Once imported, acquired areas were cropped to focus on the functional surfaces. The quality of the data obtained was visually assessed with band contrast images (Supplementary Fig. 4). No grain reconstruction, filtering or smoothing has been applied to data shown in Figs. 2, 3, or 7. Two maps of crystallographic orientations were constructed for each specimen. The X, Y, Z axes refer to directions in the specimen (Supplementary Fig. 5).

The X axis is defined by the long axis of the denticle. The Y axis points to the north and the Z axis into the plane of view. The colors of the maps represent the crystallographic orientations at each pixel of the map, relative to the X and Z axes, respectively. This is in contrast to the pole figures, which are also included in Fig. 2, which shows the distribution of crystallographic directions in the specimen coordinate system. The contoured upper hemisphere pole figures were created using the calcDensity function in MTEX using the la Vallee Poussin kernel and a half width parameter of 4° (Fig. 2).

Calculation of the Texture Index: equations
ODF f(g) = dV / (V × dg) (Eq. 1)

Orientation distribution function (ODF), a function on the orientation space that associates each orientation (g) to the volume percentage (V) of individual measurements in an area with this orientation32. V is the sample volume, dV is the volume of all grains i with the orientation g in the orientation range (angular element) dg = 1 sinΦ dφ dΦ dφ , f(g) dg = 1 and φ, Φ and φ are Euler angles.

TI = ∫ ODF f(g)^2 dg (Eq. 2)

Texture Index, a measure of deviation of individual recorded orientations (g) from the orientation distribution function (ODF) calculated for the entire studied area. For ODF f(g) and dg, see Eq. 1.

Statistical evaluation of the Texture Index between taxa
The null hypothesis that the TI values for all taxa originate from the same distribution was tested using the non-parametric Kruskal-Wallis test as implemented in R Software85 using the code provided by Shirley et al.39. Pairwise tests between samples, reported in Supplementary Table 3, were carried out using the Wilcoxon test in the same software environment, with the Bonferroni correction applied. Confidence intervals reported in Supplementary Table 4 were calculated with the Bonferroni correction, i.e. the significance level of 5% was divided by 6 (number of samples), resulting in a significance level for each pairwise comparison equal to 0.83% and a corresponding confidence interval of 99.17%.

The effect of measured area on Texture Index
Sizes of conodont elements and other vertebrate dental organs span multiple orders of magnitude, reflecting their range of body sizes and sizes of consumed prey. To investigate whether areas of functional surfaces might confound TI as a measure of the degree of structural control, we investigated the relationship between TI and the area from which it is calculated. For each sample (the entire scanned area, which might have been larger than the functional areas shown in Fig. 2), between 184 and 332 squares had been subsampled. For each square, ODF, TI, M-index34 and pfJ31 were calculated. To test the null hypothesis of no relationship (ρ = 0) we calculated Pearson’s correlation coefficient (ρ) between TI and the area. As comparison, well-documented biominerals, EBSD datasets of aragonitic nautiloid nacre and oyster prismatic layer were subject to the same procedure (Supplementary Table 1). All tests show a negative correlation between the TI and area, significant at α = 0.05. It shows that the larger an examined area is, the lower the TI will be (Supplementary Fig. 3c). This indicates that the size of larger specimens, which in our dataset tend to have more complex morphology (i.e. Bispathodus and Wurmiella are larger than the coniform taxa), does not explain higher TI values.

Comparison of metrics for quantifying crystal texture
There are several methods to quantify crystal texture from a material, the majority of which are derived from the Orientation Distribution Function (ODF), a mathematical representation of the distribution of crystallographic orientations within a polycrystalline material. TI and M-index are both measures of the degree of crystallographic texture in a material. However, they are calculated using different methods and are based on different mathematical representations of the texture. Texture Index (TI, also referred to as the J-index33) quantifies the maximum intensity of the (ODF), providing information about the dominant orientation(s) in the material. TI is commonly applied to geological materials to measure fabric strength34. This metric, however, has the disadvantage of requiring a substantial number of discrete orientations, which make it impractical to obtain a meaningful measure of fabric strength and its difficulty to interpret86. The M-index, on the other hand, characterizes the spread or distribution of crystallographic orientations, indicating how broad or diffuse the texture is (Supplementary Fig. 4).

This has the potential to provide more meaningful information on the crystal texture, in particular with datasets with limited orientation data. Another method is Pole Figure T-Index (pfJ), a measure of crystallographic texture in a material utilizing the same method as the TI, specifically based on data obtained from pole figures. This allows crystal texture to be quantified relative to an individual crystal axis.

We conducted texture analysis on our six samples using a combination of TI and M-Index measurements (Supplementary Fig. 4). TI exhibited a consistent trend across all samples, showing a clear increase in values indicative of a stronger and more pronounced texture with evolution. This suggests a progressive shift towards a preferred orientation within the denticle tissue. The M-Index showed a similar increase, but it was less pronounced compared to the changes observed in TI. Several factors may contribute to this difference. One possibility is that TI is sensitive to textures with high maximum intensities. The M-Index, which predominantly characterizes the spread of orientations, may show a more moderate response to changes in maximum intensity. A similar trend was seen in the pfJ for the a-axis and c-axis, where the differences are smaller than in the M-Index. It appears from our results that, in the context of textures analyzed here, TI has more sensitivity to subtle changes to texture differences, allowing for a higher degree of separation between samples.

We tested the influence of the sample area on these texture indices. We found that TI exhibited a consistent decrease as the area of measurement increased. This trend suggests that TI is sensitive to the measured area. In contrast, the M-index presented a more stable response to changes in the sampled area. It maintained a relatively constant value in the same material regardless of the spatial extent. This may indicate that the M-index may be less influenced by variations in the area of measurement (Supplementary Fig. 4).

Both indices investigated appear to support our hypothesis on increasing crystallographic organization. Here we predominately discuss TI due to its apparent higher sensitivity to more subtle changes in crystallographic texture (Supplementary Fig. 4).

Reconstruction of crystal grains
Conodont skeletal tissues form grains, also referred to as crystallites, which are visible in topography in broken or etched sections13,14 and in EBSD maps (Fig. 2). In EBSD datasets, algorithms based on misorientation angles are applied to calculate grain sizes. This allows inferring properties of the grains such as the area, grain orientation and length. However, grain reconstruction is done based on a threshold misorientation angle which, in most cases, is obtained from other analytical methods. For enamel and enamel-like tissues such as conodont hyaline tissue, to the best of our knowledge, no such threshold has been established. We carried out multiple grain reconstructions for conodont crown tissues using threshold values from 1° to 17° with a step of 4°. Supplementary Fig. 2 illustrates an example of calculation of grain boundaries compared to a band contrast map of a sample subset of Bispathodus cf. aculeatus. Band contrast is a measure of image quality derived from the Hough transformation. It describes the average intensity of Kikuchi bands in comparison to the overall intensity within the area of interest. The example in Supplementary Fig. 2 shows an end-product of grain reconstruction using the calcGrains() function of MTEX set to follow non convex outer boundaries based on a user defined misorientation angle. After the first round of grain reconstruction, grains smaller than 2 pixels were discarded and the function was applied again to eliminate spurious grains appearing as a result of misindexed pixels. Following this the data was then filled in Supplementary Fig. 2 using the fill() function of MTEX which attempt to fill grains with pixels of similar orientations. The resulting maps (Supplementary Fig. 2) may convey an impression of well resolved texture with a great variety of grain sizes and shapes. However, comparison with the band contrast image reveals that the inferred grain boundaries mostly correspond to non-indexed areas. Thus, obtained grain distribution would reflect the pattern of missing data and not the “true” grain shapes and sizes. No cleaning algorithms were used in our post processing workflow. We identify two major barriers to grain reconstruction in conodont tissues: (1) low indexing rate due to beam damage and challenging preparation, (2) internal heterogeneity of tissues, such as small misorientation angles between adjacent grains, typical for biominerals and often reflect the mechanism of formation.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a very interesting manuscript and I concur that, if and once revised appropriately, it may be\npublished in Nature Communications.\nThe first vertebrate skeleton is found in conodonts as dental elements made of apatite. Thanks to the\nfirst successful application of Electron Backscatter Diffraction to bioapatite, the authors show that the\ncrystallographic order varies from one sample to another. They argue that the crystallographic order\nincreases throughout the evolution of conodont teeth, in parallel with morphological adaptation to food\nprocessing, thereby suggesting increasing control over biomineralization throughout the evolution of\nthe first skeletonized vertebrates. This opens interesting perspectives, also for the application of this\ntechnique to other vertebrates.\nYet, in my opinion, the authors fail, for now, to justify convincingly their sampling choices: they\nanalyzed 6 taxa only and they argue they are representative of ‘increasing degrees of adaptation to\ndental function’. Besides the fact that I do not understand what ‘adaptation to dental function’ exactly\nmeans and how this was quantified here, there seems to be a risk of circular reasoning and/or of\n‘cherry picking’ the data that agrees with one’s hypothesis. To me the authors’ figure 2a does not\nnecessarily demonstrate a temporal trend of texture index but rather suggests that texture index\nincreases from coniform to blade to platform conodont elements. Since these three types of conodont\nelements can also co-occur, the authors should discuss the polarity of this sequence in conodont\nevolution (btw, does it really reflect increased ‘sophistication’?). I would also encourage them to\ndiscuss the possibility that these three geometries correspond to three (or more) ways to process food\nand whether their hypotheses about the occlusal axes are robust and whether considering only\nuniaxial compressive stresses is relevant here.\nI am looking forward to reading the revised version of this ms.\nBest regards,\nNicolas Goudemand\nother (mostly minor) comments:\nline 24: ‘throughout the evolution of conodont teeth’: not convinced that these 6 taxa represent the\nwhole story. Please nuance.\nline 57: ‘precise occlusal’: ‘precise’ does not seem to be the right term here\nline 65: ‘digestion’: replace with ‘processing’?\nline 71: ‘for many conodont lineages’: replace ‘many’ with ‘some’ ?\nline 84: ‘load’: only when that load is uniaxial compression\nline 93: precise what is ‘V’\nline 104: what do you mean with ‘increasing degrees of adaptation to dental function’? Coniform\nelements may be as adapted to their particular function as blade or platform elements are to theirs…\nline 114: how do you define those axes? What if the mastication movement includes rotations?\nline 122: ‘orientations’\nlines 120 and 124: what are those brackets?\nlines 146: what if one analyzes a more ‘derived’ (e.g. Triassic) coniform element and a less ‘derived’\nplatform element?\nline175: ‘least specialized feeding in Proconodontus’ : How do you know that feeding in Proconodontus\nis less specialized than in Palmatolepis?\nline 199: ‘however’: ?? unclear logical relationship\nline 212: ‘sophistication’: not demonstrated\nline 233: too strong a statement\nline 236: ‘selectively’: selectivity?", "reviewer_2": "This is a significant and nicely illustrated conodont manuscript describing the variation of degree of\ncrystallographic order in selected conodont material. The approach is innovative. My main criticism is\nthat expanding conclusions to all groups of early vertebrates is hazardous.\nIn addition, the introductory part is overestimating the role of vertebrates.\n- Line 17: \"outstanding diversity of skeletal tissue\": what about calcareous shelled skeletal organisms?\nThey did the same (also in the microworld) or even better.\n- Line 18: \"extreme hardness and elasticity\": siliceous shells are even harder. Furthermore, elasticity\ndoes not reflect simple mineralization, but strongly depends on the presence/quantity of organic\nmaterial.\n- Line 20: \"teeth\" or \"tooth-like\" elements as reported in line 52?\n- Lines 35-37: \"Vertebrates are ... the one that has achieved the highest diversity and disparity\":\nWhat about insects? As regards disparity, it is a delicate topic: disparity is not only morphology, but\ndisparity can be any phenotypic aspect including behaviour and ecology.\n- Line 45: \"first vertebrates\": the manuscript deals only with conodonts (six specimens ranging in age\nfrom the Cambrian to the Devonian). In that time frame many other vertebrates appeared in the sea\nand were going to jump on the land. Why wasn't that material investigated?\nIn general:\n- What about the role of diagenesis on your material?\n- You base your conclusions on a very few specimens (six): how confident can you be on such a\nlimited set?", "reviewer_3": "1. Short summary of opinion\nThis MS was of interest because EBSD analysis from apatite biominerals are very rare and is in its\ninfancy, so this study has a good potential to be a significant study. At the same time, unfortunately,\nthe robustness of data (“sufficiently strong evidence”) is not as good as what I expect for a usual Nat.\nComm paper.\n2. Key results\nShirley et al. analyzed six conodont taxa (from “primitive” to “derived” ones) with EBSD, which is a\nsuperb tool for crystallographic data of (bio)minerals. The authors succeeded in getting the orientation\nmaps and pole figures, which were interpreted in the context of conodont evolution. To my knowledge,\nearlier studies tried EBSD analysis to conodont, and some were even successful at getting a good map\nbut getting the data from several taxa is the major accomplishment of this study. I myself once tried\nEBSD analysis to bones (bioapatite) and was disappointed by the absence of any Kikuchi bands, so I\nfully understand that what the authors accomplished herein is a very good contribution.\n3. Validity\nWith that being said, I have concerns for the validity of the research. Please note that I do not mean\nthat it is not a good study design. I just considered that Nat. Comm. is a very selective journal, so I\nalso used stringent criteria, which I expect for a paper published in Nat. Comm.\n3.1 Here, the authors did not consider any taphonomic status of each specimen. I felt that the authors\ntreated the specimens as if they were neontological specimens, which are free from taphonomic\nconsideration. For example, Proconodontus muelleri is Cambrian in age, and Palmatolepis sp. is\nDevonian in age, so there is nearly 100 Myr time gap between the two specimens. So one may expect\nthat ‘biogenic’ minerals in Cambrian Proconodontus muelleri is comparatively less well-preserved than\nthose of Devonian Palmatolepis sp. and it may influence the difference between the crystallographic\nalignments of two specimens. But there is no consideration like this in the current study.\nBesides, even if ages of two specimen are similar (e.g. Devonian Tripodellus and Palmatolepis), spatial\nsetting can affect the preservation. For example, when a formation of conodont A was deeply buried or\nunder the influence of significant tectonic pressures and heat, whereas conodont B was buried in a\n‘peaceful formation’, which was relatively free from pressure and heat, the preservation status of\nconodont A and B may not be the same.\nIn my opinion, neglecting the taphonomic status of specimens can significantly undermine the validity\nof this study. Please note that consideration of taphonomic status is an important topic for many\npaleontological studies (e.g., see a great example of Casella et al. 2018, Palaeo3,\nhttps://doi.org/10.1016/j.palaeo.2018.03.011).\n3.2 Although I fully understand getting EBSD maps from conodonts is the very hard task, I still have\nconcern for the small sample size (i.e., 6). There are many families in conodont\n(https://en.wikipedia.org/wiki/Conodont), so I expect that the number of genera is larger than that of\nfamilies. So ‘six’ may be an insufficient number to represent the pattern in conodont evolution and it\nmay weaken the robustness of the conclusion. In other words, a future study that uses much more\ndiverse conodont taxa may report a different conclusion, and the current study is vulnerable to this\ntype of challenge. In addition, if I understand correctly, just a single map was acquired from each\nspecimen? In my experience, the crystallographic orientation (i.e. pole figure) is pretty local\n(especially a-axis direction), so if a different part of the same specimen is analyzed with EBSD, the\nresults (e.g. pole figure) may be different. In other words, just a single map per a specimen may not\nguarantee the solid reproducibility, which is required for Nat. Comm. papers.\nAgain, please do not get me wrong. I truly think ‘six’ is already a great accomplishment, but I just feel\nthat it is somewhat weaker than the ‘level of support’ criterion of Nat. Comm.\n4. Significance\nWhen it comes to significance, I think it’s A+. Because this research proved that spatiotemporally\ndiverse conodonts can be analyzed with EBSD, many conodont researchers will adopt EBSD for their\narsenal. In addition, many other vertebrates have teeth (e.g. fish, Felidae, T. rex, the authors, the\nreviewer [though some are replaced with amalgam], etc), so the accomplishments of this study will\ninfluence many fields of science such as paleontology, odontology, zoology, etc.\nI just wish that the contribution of this study and that of Atakul-Özdemir et al. (2021, R. Soc. Open\nSci.) should have been clearly explained in the main text to guarantee the significance of this study.\n5. Data and methodology\n5.1 About the methodology: The panel B of Extended Data Figure 2 used a ‘cleaning’ to enhance the\nmapping. However, compared to the panel A, the cleaned maps in the panel B look ‘too clear’ and it\nmay mean that the authors used too radical cleaning method for the panel B.\nIf this cleaning method was used in Figure 1, I would say that Figure 1 may contain pretty much\n‘false-positive’ data, which are attributable to the radical cleaning. If Figure 1 used ‘raw data’, which\nwas not influenced by any cleaning, please dismiss my current comment (and I believe the authors did\nnot apply ‘radical cleaning’ to Figure 1).\nAnyway, please provide this sort of information and all other steps applied to present Figure 1 to\nguarantee the full reproducibility. Currently, I felt that the method section is not detailed enough, so\nsome steps could not be evaluated with enough confidence.\n5.2 About the quality of presentation: I think there are some room for improvements. Please see\nbelow “11. Additional comments”.\n6. Analytical approach\nThis issue was explained in the ‘3’ and ‘5’ above.\n7. Suggested improvements\nThis part is directly related to ‘3. Validity’ above.\n7.1 Taphonomic consideration: I strongly recommend that Casella et al. (2018) be a helpful literature\nthat explains why considering taphonomic status is important for biomineral studies. There, one can\nsee that taphonomic influences degrade the calcite in brachiopods and it can make ‘fake news’ from\nthe old past. Admittedly, calcite is soft mineral as authors pointed out, so bioapatite fossils like\nconodont may not be severely affected by diagenesis (against my concerns). However, at least\ninvestigating the taphonomic status of the materials used in this study may be able to defend the\nconcerns of mine and potential future readers.\nI suggest that investigating the Raman spectra of conodont may yield useful information at least for\nthermal alteration. Please see McMillan and Golding (2019, Palaeo3) for details. If it turns out that\nCambrian materials are heavily thermally influenced while younger Devonian materials (the ones used\nin Figure 1) are less affected by thermal alteration, I think this factor may contribute to ‘false-positive’\npattern (highly aligned Devonian signals vs weakly aligned Cambrian signals) in Figure 1. So I wish\nthat the authors eliminate this sort of possibility completely with additional Raman data.\nAt the same time, I know Raman spectroscopy is not always available, and learning new skills used in\nMcMillan and Golding (2019) can be a harsh request as a revision. So I think somewhat more practical\nsuggestion (though I prefer Raman data) would be to investigate the classic conodont Color Alteration\nIndex of the materials used in this study. If CAI shows that older materials are not always more\nheavily thermally altered, I believe the pattern reported in this study can be more safely understood\nas the biological signal as authors claimed (against my “devil’s advocate” concern).\n7.2 Ideally, I believe more wider sampling of taxa or getting a few more EBSD maps from each taxon\n(to guarantee the reproducibility of Figure 1) will definitely improve the robustness of the conclusion.\nBut I know how hard it is to prepare new EBSD specimens and running EBSD analysis. So I will leave\nit to the editor’s discretion.\n8. Clarity and context\nThe text is generally clear and explanation of context at Introduction is excellent (another A+). But\nsome parts are hard to follow. Please see my further comments in “11. Additional comments”. Also, I\nbelieve the title tries to encompass too wider implications. In my opinion, even if the title is simply\nconstrained to conodont instead of wider ‘early vertebrate’, this study is equally awesome, so please\nconsider this issue during the revision.\n9. References\nI am not a conodont specialist, so I am not qualified to officially say something about the quality of\nbibliography. That said, many great journals and renowned paleontologists’ name in the bibliography\nmade me believe that the literature cited are appropriate.\n10. My expertise\nWhen it comes to statistical analysis, I am not good at it, so please follow the opinion of other\nreviewers. About the validity of textural analysis and technical issues of EBSD analysis, such as\nOrientation Distribution Function and Texture Index, other reviewer(s) will provide more solid opinions\nthan I can.\nAlso, because I am not a conodont expertise, I am not qualified to provide critical comments to\nDiscussion. I hope other conodont experts can do this.\n11. Additional comments (about writing)\n11.1 Title\nIn my opinion, it needs to be more focused on the core contents. When I read the title, I anticipated\nthat this study used wide variety of vertebrates fossils (especially diverse Agnatha of Paleozoic) or at\nleast analyzed much more diverse taxa of conodont. Please do not get me wrong. I strongly believe\nthat getting the crystallographic data from conodonts is a very important contribution, but current title\nmay not be representative for the contents of the paper because all results were obtained from the\nconodonts, and someone may think that just one clade (=conodont) is not enough to represent ‘early\nvertebrate’. Although I know that a concise title is preferred for a journal like Nat. Comm., I still\nbelieve that the best title should be a direct indicator for the contents.\n11.2 Abstract\nIt is a well written summary. I point out two things: 1. There are four citations, but to my knowledge,\nNat. Comm. does not allow citations in abstracts; 2. It is composed of more than 150 words, which is\na length limit for abstracts of Nat. Comm.\nI feel the current abstract is more like a that of Nature, so please revise the abstract to meet the\nguideline of Nat. Comm.\n11.3 Introduction\nThe background information is very well presented and justifies the research. However, the\narrangement needs improvements. The findings of other research and the contribution of current\nresearch is somewhat jumbled. I feel clear separation of ‘what others did (i.e., literature review)’ and\n‘what we did in this research (i.e., Here we, blah blah)’ will make a better Introduction.\nThe authors also mentioned c-axis of apatite and it is an essential terminology for the study. My small\nconcern is that very few readers will know what ‘c-axis’ means. I would recommend that the authors\nprepare a simple supplementary figure to show a unit cell apatite and its axes direction. I believe it\nwill help readers to follow your thoughts.\nPlease see also my minor comments in the annotated PDF.\n11.4 Methods (Preparation and data acquisition)\nIf space permits, please briefly explain how this study succeeded in getting clear Kikuchi bands from\nthe polished surface of conodont unlike earlier research, which were not that successful. Although\nShirley et al.10 was cited, readers may want to know the difference without visiting Shirley et al.10.\nPlease see also my minor comments in the annotated PDF.\n11.4 Extended Data\n(1) Figure 1. Although SEM images are very good at showing the morphological details, I recommend\nstereoscopic microscope images (e.g. Panels D, E) because they show color. I think color image is\nespecially important (please see my “7.1 Taphonomic consideration” above). If possible, please\nprovide.\n(2) Figure 2(a). It is very hard to find any difference among the different grain boundary thresholds. It\nis understandable because just one degree of difference is very subtle, but I am not sure then what\nangle is the most suitable for the match between ‘human eye’ and grain boundary threshold. I\nanticipated the answer from the authors somewhere in Lines 441–454, but there was no answer (also\nin Lines 135–138). It may confuse readers.\nFigure 2(b). Because there is no detailed information for the “standard grain smoothing and filter\n(here half-quadratic)” (Line 446), I cannot provide clear feedback. However, if ‘Filled grains’ maps are\nthe results of ‘cleaning (=extrapolation of pixels)’ of the original data presented in Panel (a), I think it\nis problematic because the ‘cleaning’ step was too aggressive and radical. As authors pointed out, the\n‘boundaries’ in the cleaned maps may not show the precise positions of real boundaries, grain shapes,\nand sizes. Most importantly, orientation information may be affected by this aggressive, anthropogenic\nfactor. Did the authors used the “standard grain smoothing and filter (here half-quadratic)” cleaning\nmethod to make Figure 1 (I hope not)?\n11.6 Results\n(1) Lines 118–122: Here, the description is very hard to follow. What are the ‘elongated apatite grains’? There are many shapes of grain in SIX EBSD maps, so it’s hard to know what I should see. In\naddition, due to the diverse shapes of apatite, I am pretty skeptical that the description is generally\napplicable to the maps in Figure 1.\n(2) Figure 1: Although the maps are most probably ‘inverse pole figure maps’, no legend was provided\nso it is hard to fully understand the meaning of colors in the maps. In addition, there is no explanation\nfor the numbers above the pole figures. I believe they are MUD (multiple of uniform density) and if it\nis correct, please explain at the caption.\n11.7 Discussion\n(1) Lines 197–202: Similar to Lines 118–122, the authors keep mentioning ‘elongated grains’, but it is\nnot that clear to readers. Please clarify what are ‘elongated grains’ and get more generality for this\nobservation.\n(2) Lines 218–223: Perhaps due to the too long sentence, I cannot clearly understand what this part\ntries to explain. Could you please rephrase?", "reviewer_4": "Review to the manuscript titled “Increasing control over biomineralization in early vertebrates\".\nThe work documented by Shirley et al. suggested the evolution of the first mineralized vertebrate\ntissues in conodont by progressive adaptation to dental function. The indexing of crystallographic axes\nof apatite in conodont fossil by EBSD analysis was successful and the results well showed that c-axes\nof apatite consisting of hypermineralized crown tissues were aligned parallel/subparallel to the occlusal\n(biting) axes. It is interesting that the preferred crystallographic orientations of conodont’s crown\ntissues could be developed by evolving conodont feeding mechanisms. However, I feel that further\nconsideration of two issues is required to interpret the EBSD results.\n1. How did you rule out the influence of diagenesis on the development of preferred crystallographic\norientations of conodont fossils? You mentioned that “the hyaline tissue is widely used by geochemists\nthanks to its outstanding resistance to diagenesis (lines 67-68)”, but this does not guarantee the\nresistance of recrystallization. The related reference (15. Zhang et al., 2017) also does not suggest\nthat recrystallization of conodont apatite does not occur easily by diagenesis, but suggests the\nintegrated approach to recognize the nature and degree of diagenetic alteration in conodont specimen.\nIs there any other evidence that development of preferred crystallographic orientations of biomineral\ndid not be affected by diagenesis? In this study, the spatial distribution of conodont samples is varied,\nsuch as U.S.A., Ukraine, Germany, Sweden, Canada, and Belgium (maybe?). The time variation is also\nlarge, from upper Cambrian to upper Devonian. These may indicate that the environments of\ndiagenesis of samples should be different. That is, how did you rule out the effect of pressure and\ntemperature on the crystallographic texture?\nIn addition, you mentioned that “the discordance between the long axes of grains and the c-axes is\nunusual in abiotic minerals but common in biological tissues (line 121-122)”. What are the long axes of\nelongated apatite grains? From what I understand, the long axis of apatite is the direction of\nmaximum length of same colored grain in grain color maps within Figure 1. In Fig.1f, the elongated\nlong axes of apatite are likely to be perpendicular to the crystallographic c-axis of apatite. However,\nthis trend is rarely observed in other samples. In Fig.1a and 1b, the crown appears to consist of\nalmost a large single crystalline grain. Especially, in Fig.1d, the elongated long axes of apatite are\nparallel to the c-axes of apatite. In this case, is there a possibility of recrystallization in these samples?\n2. In this study, six hyaline tissues of conodont dental elements were analyzed for six taxa of\nconodonts. That is, only one section of hyaline tissue was selected for each taxon. In this case, each\nspecimen was understood to be representative of each taxon of conodonts. However, is the cross-section\nrepresentative sufficient to interpret the skeletal evolution of conodonts according to changes in\nfeeding function? Of course, this result is the first successful analysis of EBSD for bioapatite, but I\nthink that only one section for each taxon would not be enough to representative of one taxon of\nconodonts. Although Pro. muelleri and Pan. equicostatus have a single coniform element in each, the\nother four samples appear to have multiple hyaline tissues. Are there similar trends in EBSD results\nfor other hyaline tissues in each specimen?\nIn addition to these major points, I’d like to suggest some recommendations for minor issues.\n1. In Fig.1, what is the meaning of the reference frame (X- and Y-directions)? Do the black areas in\nschematic diagrams indicate the grain color maps in each specimen? From what I understand, you\nwant to show the crystallographic c-axes of apatite grains aligned perpendicular to the occlusal\nsurfaces (that is, aligned parallel to the occlusal (biting) axes). For this purpose, it is recommended to\nuse an inverse pole figure (IPF) map in terms of the occlusal axis.\nOf course, Fig.1 can show the c-axes aligned parallel to the occlusal surface, but an IPF or IPF map is\ngenerally used to show where the crystallographic axes of minerals are aligned in a specific direction\n(like in the reference 26. Kilian and Heilbronner, 2017). Because you processed the EBSD data using\nMTEX, a MATLAB toolbox, I believe that you can rotate the Euler angles to the specific direction. If X-\nand Y-directions are the sample coordinates in the SEM, you can make the X-direction aligned parallel\nto the occlusal surface, and Y-direction aligned perpendicular to the occlusal surface (that is, parallel\nto the occlusal axes). Then, when you plot the IPF map or IPF in terms of the Y-direction, it is easy to\nsee which axis is perpendicular to the occlusal surface (in this case, maybe c-axis).\nOtherwise, please indicate the meaning of the X- and Y- directions and match the directions between\nthe schematic diagram and grain color map in each specimen.\n2. You calculated the texture index (TI) to quantify the spread of crystallographic orientation of\nbiogenic apatite. You tried many re-calculations of ODF and TI for random subsets to eliminate the\nbias of the analyzed size of sample, but you excluded the calculation of grain size because of no\nexistence of known misorientation angle to discernible grain boundary in enamel and enamel-like\ntissues. Of course, it is not easy to define the reference angle of the grain boundary. However, in\nmaterials science (generally used for alloys), boundary with misorientation angle of ~15° is used to\ndefine a high-angle grain boundary (indicating the general grain boundary), and boundary with\nmisorientation angle less than 15° is considered as a low-angle grain boundary (indicating the\nsubgrain boundary). In the field of geology, this angle is considered to be about 10° in the rock sample.\nIn my opinion, the reference angle of the grain boundary for biogenic apatite would be similar to that of\nabiogenic (abiotic) minerals because of the similar mechanical properties between biogenic and abiogenic\nminerals. If you calculate the grain boundary of the sample, you can easily interpret the trend of TI\nvalues considering a grain size distribution.\nTI values of T. gracilis and Palmatolepis sp. are higher than those of Pro. muelleri and Pan.\nequicostatus because former two samples are composed almost large single crystalline texture. The\nnarrow TI value ranges of Pro. muelleri and Pan. equicostatus indicate the homogeneous grain size\ndistribution. I don’t know why grain size distribution is differ among the samples, but the crystallite of\nhyaline tissue appears to generally tend to become larger towards the later taxa of conodonts.\nIn addition, it seems that boundaries with misorientation angle less than 10° should be regarded as\nlow-angle grain boundaries (that is, subgrain boundary). In that case, the kernel angle misorientation\n(KAM) map or grain reference orientation deviation (GROD) angle map is more suitable to show low-\nangle grain boundaries in the biogenic apatite minerals. The KAM map shows the misorientation less\nthan the predefined threshold value calculated by the mean orientation between a point and its\nneighbors. The GROD map is generated based on the deviation between the mean orientation of a\nreference point and those of other points. You mentioned that “the misorientations vary between areas\nwithin one tissue depending on their position with respect to functional directions and the\ncompressional and shear stresses acting in these directions (lines 133-135)”. So, these maps can\neasily represent the intracrystalline misorientation within one tissue regardless of the existence of\nknown misorientation angle to discernible grain boundary. It may require some data processing, but I\nbelieve that you can plot these maps by using codes within MTEX.\n3. In lines 93-95, are reference numbers 28 and 29 correct for ODF? Please check and correct with\nappropriate references (maybe Bunge, H.-J., 1982. Texture Analysis in Materials Science:\nMathematical Models. Butterworths, London, UK).\nIn lines 100-101, a texture index (TI) also needs a citation. Please cite the appropriate\nreferences (probably Bunge (1982), or generally used as “J-index” by Mainprice and Silver (1993)).\nIf you want to just quantify the sharpness of pole figure of c-axes, pole figure J-index (pfJ-index) is\nalso considerable. This value is analytically defined in a similar manner of J-index\n(Michibayashi and Mainprice, 2004), but as far as I know, there are no code in MTEX unfortunately.\nIf you want to quantify the intracrystalline misorientation within a hyaline tissue, it is considerable to\nuse a “misorientation index (M-index)” calculated from the misorientation angles among the\nuncorrelated grain pairs (Skemer et al., 2005). Unlike J-index (texture index), M-index has an\nadvantage that it is hardly affected by the number of data, but also has a disadvantage of loss of\nspatial and textural information of the original dataset. I’m not sure whether you can get more\nreasonable result than J-index, but it is worth considering the M-index to quantify the crystallinity of\nhyaline tissue. However, be careful if you calculate the M-index by MTEX code. You have to do it\ndirectly using raw EBSD data, not using the result of ODF calculation.\nOverall, the results of this manuscript will be interesting to geologists, especially mineralogist, and I\nthink it may have the potential to provide important inspiration to paleontologists. However, several\nissues require clarification before it is considered by Nature Communications, at least for major issues.", "reviewer_author_response": "Response to comments by the reviewers\nReponses are marked in blue\nReviewer #1 (Remarks to the Author):\nThis is a very interesting manuscript and I concur that, if and once revised appropriately, it\nmay be published in Nature Communications.\nThe first vertebrate skeleton is found in conodonts as dental elements made of apatite.\nThanks to the first successful application of Electron Backscatter Diffraction to bioapatite, the\nauthors show that the crystallographic order varies from one sample to another. They argue\nthat the crystallographic order increases throughout the evolution of conodont teeth, in\nparallel with morphological adaptation to food processing, thereby suggesting increasing\ncontrol over biomineralization throughout the evolution of the first skeletonized vertebrates.\nThis opens interesting perspectives, also for the application of this technique to other\nvertebrates.\nYet, in my opinion, the authors fail, for now, to justify convincingly their sampling choices:\nthey analyzed 6 taxa only and they argue they are representative of ‘increasing degrees of\nadaptation to dental function’.\nThe number of specimens reflects the difficulties of this analytical methods. The six taxa is\nthe outcome of a 5-year project. This is a second study in which any EBSD data has been\nobtained and the only previous one (which corresponded to an entire postdoctoral project)\nexamined only one taxon and yielded data of limited use, because the datasets and\nprocessing methods were not provided. In our sampling, we focused on the earliest and on\nthe least derived taxa to capture the early stages of the evolution of biomineralization.\nTherefore the youngest taxa are represented poorly: they show smaller disparity than the\nearly conodonts, based on which we made the assumption that the highest ultrastructural\ndiversity should be sought in older taxa. The complete datasets and code are provided with\nthe manuscript, which is unique for a EBSD study in palaeobiology.\nBesides the fact that I do not understand what ‘adaptation to dental function’ exactly means\nand how this was quantified here, there seems to be a risk of circular reasoning and/or of\n‘cherry picking’ the data that agrees with one’s hypothesis.\nThe “adaptation to dental function” in this study follows the definition proposed by\nDonoghue (2001), as we directly test his hypothesis. This has been spelled out in the text.\nFurthermore, we specify that “adaptation to dental function” is not a general metric measured the same\nway in all organisms. It is related to how the tooth (or its individual surfaces) operates, e.g. if\nit’s a shearing, crushing or grinding tooth and in which directions it occludes. We specify in\nthe text that in case of conodont denticles or single cones, that adaptation is reflected in the\nstrongest resistance to uniaxial compressional stress.\nTo me the authors’ figure 2a does not necessarily demonstrate a temporal trend of texture\nindex but rather suggests that texture index increases from coniform to blade to platform\nconodont elements. Since these three types of conodont elements can also co-occur, the\nauthors should discuss the polarity of this sequence in conodont evolution (btw, does it really\nreflect increased ‘sophistication’?).\nThis seems to be a misunderstanding between a temporal succession and the degree of\nevolutionary change. A younger taxon does not need to be more derived than an older one.\n“Primitive” (preserving many ancestral characters) organisms can coexist for long periods\nwith their more derived relatives. Therefore we do not present a temporal trend and in that\nsense the reviewer is correct that this is not what Figure 3a (formerly 2a) shows. It shows a\nphylogenetic tree with evolutionary distance distinguishing the taxa. Derived conodonts are\nthose showing more specialized morphologies. Thus, more complex and more derived\nmorphologies can co-occur with simpler and less derived ones: they are still separated by a\nlong evolutionary distance.\nI would also encourage them to discuss the possibility that these three geometries\ncorrespond to three (or more) ways to process food and whether their hypotheses about the\nocclusal axes are robust and whether considering only uniaxial compressive stresses is\nrelevant here.\nWe distinguished between occlusion models of entire apparatuses and the occlusal axes of\nsingle denticles and cusps. The distinction is now emphasized in the text. The article focuses\non comparison of the single morphological unit (denticle or cusp) in a conodont element,\nwhich can be compared across the entire lineage.\nline 114: how do you define those axes? What if the mastication movement includes\nrotations?\nWe do consider rotations and the asymmetry of the movement in the results and discussion.\nBut we do not go into this detail in the introduction because it will be confusing for the reader\nbefore they see the figures.\nAll the corrections below have been applied:\nother (mostly minor) comments:\nline 24: ‘throughout the evolution of conodont teeth’: not convinced that these 6 taxa\nrepresent the whole story. Please nuance.\nline 57: ‘precise occlusal’: ‘precise’ does not seem to be the right term here\nline 71: ‘for many conodont lineages’: replace ‘many’ with ‘some’ ?\nline 84: ‘load’: only when that load is uniaxial compression\nline 93: precise what is ‘V’\nline 104: what do you mean with ‘increasing degrees of adaptation to dental function’?\nConiform elements may be as adapted to their particular function as blade or platform\nelements are to theirs…\nline 122: ‘orientations’\nline175: ‘least specialized feeding in Proconodontus’ : How do you know that feeding in\nProconodontus is less specialized than in Palmatolepis?\nline 199: ‘however’: ?? unclear logical relationship\nline 212: ‘sophistication’: not demonstrated\nline 233: too strong a statement\nlines 120 and 124: what are those brackets?\n- These brackets describe a set of directions or axes equivalent to [hkl] under the\nsymmetry group of the crystal; a direction or axis in the set. This is standard convention in\nmineralogy; we did not see the point of explaining it in the main text since it is not required to\nfollow the text for a person without a mineralogical background and anyone familiar with\nmineralogy will be familiar with this notation.\nlines 146: what if one analyzes a more ‘derived’ (e.g. Triassic) coniform element and a less\n‘derived’ platform element?\nYounger does not mean more derived, as explained earlier.\nline 65: ‘digestion’: replace with ‘processing’?\n“Digestion” is commonly used in the literature in this context so we decided not to follow this\nsuggestion.\nReviewer #2 (Remarks to the Author):\nThis is a significant and nicely illustrated conodont manuscript describing the variation of\ndegree of crystallographic order in selected conodont material. The approach is innovative.\nMy main criticism is that expanding conclusions to all groups of early vertebrates is\nhazardous.\nWe revised the text to assure that our conclusions focus on conodonts. Those mentions of\nvertebrates which explain the broader context of the study were left in the text. Although we\nfocused on conodonts as the earliest biomineralizing vertebrates, our study proposes a\nmethodology that applies to other vertebrate groups.\n- Line 17: \"outstanding diversity of skeletal tissue\": what about calcareous shelled skeletal\norganisms? They did the same (also in the microworld) or even better.\nThat’s across the many calcifying phyla. Any single phylum of calcifying organisms will show\na much smaller skeletal diversity. The sentence emphasizes diversity within one taxonomic\ngroup (vertebrates), not a competition between organisms with phosphatic vs. calcareous\nskeletons.\n- Line 18: \"extreme hardness and elasticity\": siliceous shells are even harder. Furthermore,\nelasticity does not reflect simple mineralization, but strongly depends on the\npresence/quantity of organic material.\nWe agree with the reviewer. The sentence has been corrected to “outstanding combination of\nhardness and elasticity”, which is unique to vertebrate tissues. Siliceous shells are harder, but more\nbrittle and have lower elasticity.\n- Line 20: \"teeth\" or \"tooth-like\" elements as reported in line 52?\nWe stick to “teeth” for readability. Conodont teeth are not homologous with the teeth of other\nvertebrates, but so are eyes of arthropods and vertebrates or wings of pterosaurs and\ninsects; yet they are still referred to as “eyes” and “wings”. We added the clarification\n“conodont teeth, which evolved in parallel to all other vertebrate teeth” at the beginning of\nthe introduction to make sure there is no misunderstanding.\n- Lines 35-37: \"Vertebrates are ... the one that has achieved the highest diversity and\ndisparity\": What about insects? As regards disparity, it is a delicate topic: disparity is not only\nmorphology, but disparity can be any phenotypic aspect including behaviour and ecology.\n- This was indeed not clear. We removed “disparity” as the term has several meanings and\nin terms of species richness vertebrates are not the most diverse group. However, we keep\n“disparity” because, compared to insects, vertebrates clearly dominate in terms of the range\nof environments (from deep marine active swimmers through arctic birds and humans\nlanding on the Moon) and of body sizes.\n- Line 45: \"first vertebrates\": the manuscript deals only with conodonts (six specimens\nranging in age from the Cambrian to the Devonian). In that time frame many other\nvertebrates appeared in the sea and were going to jump on the land. Why wasn't that\nmaterial investigated?\nWe tested the hypothesis by Donoghue (2001) about the FIRST vertebrate dental tissues.\nThe appearance of the teeth precedes the emergence from the sea to the land by hundreds\nof millions of years and there is no causal link between these events. In fact, many toothed\nvertebrate lineages never left the oceans. We focus on the origins of the first dental tissues\nand hope that the scientific community will endorse and adapt our methodology to study\nfurther questions on teeth evolution.\nIn general:\n- What about the role of diagenesis on your material?\nThe revised version has been expanded with a thorough Raman analysis and discussion of\npossible influence of diagenesis.\n- You base your conclusions on a very few specimens (six): how confident can you be on\nsuch a limited set?\n- All quantitative analyses (Texture Index, parameters of the 𝜈 -PO 3− band carried out\n1 4\nusing Raman spectroscopy) and all our conclusions are based on statistical tests carried out\nat the confidence level of 0.05. For example, the differences in TI between taxa are so big\nthat they are significant even for such a small number of samples.\nReviewer #3 (Remarks to the Author):\n3.1 Here, the authors did not consider any taphonomic status of each specimen. I felt that\nthe authors treated the specimens as if they were neontological specimens, which are free\nfrom taphonomic consideration. For example, Proconodontus muelleri is Cambrian in age,\nand Palmatolepis sp. is Devonian in age, so there is nearly 100 Myr time gap between the\ntwo specimens. So one may expect that ‘biogenic’ minerals in Cambrian Proconodontus\nmuelleri is comparatively less well-preserved than those of Devonian Palmatolepis sp. and it\nmay influence the difference between the crystallographic alignments of two specimens. But\nthere is no consideration like this in the current study.\nBesides, even if ages of two specimen are similar (e.g. Devonian Tripodellus and\nPalmatolepis), spatial setting can affect the preservation. For example, when a formation of\nconodont A was deeply buried or under the influence of significant tectonic pressures and\nheat, whereas conodont B was buried in a ‘peaceful formation’, which was relatively free\nfrom pressure and heat, the preservation status of conodont A and B may not be the same.\nIn my opinion, neglecting the taphonomic status of specimens can significantly undermine\nthe validity of this study. Please note that consideration of taphonomic status is an important\ntopic for many paleontological studies (e.g., see a great example of Casella et al. 2018,\nPalaeo3"}, "author_response": ""}]

Decision letter