Charlotte Yvette Henderson
email: c.y.henderson@uc.pt

Valentina Mariotti

Frédéric Santos

Sébastien Villotte

Cynthia A Wilczak

The new Coimbra method for recording entheseal changes and the effect of age-at-death

Keywords: Asymmetry, Activity-patterns, Ageing, Degeneration Asymétrie, Pattern d'activités, Vieillissement, Dégénération

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

The most common method for identifying specific types of activities remains the observation and recording of entheseal changes (ECs) [START_REF] Lieverse | Upperlimb musculoskeletal stress markers among middle Holocene foragers of Siberia's Cis-Baikal region[END_REF][START_REF] Lieverse | Lowerlimb activity in the Cis-Baikal: Entheseal changes among middle Holocene Siberian foragers[END_REF][START_REF] Jurmain | Bioarchaeology's Holy Grail: The reconstruction of activity[END_REF]. These are changes from the normal biological appearance of the attachment of tendons and ligaments to the bone [START_REF] Villotte | In search of consensus: terminology for entheseal changes (EC)[END_REF]. They are frequently interpreted as indicators of repetitive loading of specific muscles, or muscle groups, and are therefore considered of good specificity to particular movements, which are then interpreted as tasks [START_REF] Jurmain | Bioarchaeology's Holy Grail: The reconstruction of activity[END_REF][START_REF] Hawkey | Activity-induced musculoskeletalstress markers (MSM) and subsistence strategy changes among ancient Hudson Bay Eskimos[END_REF][START_REF] Peterson | The Natufian hunting conundrum: Spears,atlatls, or bows? Musculoskeletal and armature evidence[END_REF][START_REF] Weiss | Muscle markers revisited: Activity pattern reconstruction with controls in a central California Amerind population[END_REF].

The enthesis itself is a poorly understood region of the muscle-bone complex, with limited current clinical research in this area. The research that exists focuses on changes associated with diseases, e.g. the seronegative spondyloarthropathies, which include diseases such as ankylosing spondylitis [START_REF] Benjamin | The skeletal attachment of tendons -Tendon "entheses[END_REF][START_REF] Benjamin | The anatomical basis for disease localisation in seronegative spondyloarthropathy at entheses and related sites[END_REF][START_REF] Benjamin | The "enthesis organ" concept: Why enthesopathies may not present as focal insertional disorders[END_REF][START_REF] Mcgonagle | Imaging the joint and enthesis: Insightsinto pathogenesis of psoriatic arthritis[END_REF]. While the development of the enthesis and its normal appearance are well described in the literature [START_REF] Benjamin | The skeletal attachment of tendons -Tendon "entheses[END_REF][START_REF] Villotte | Understanding entheseal changes: Definition and life course changes[END_REF], there has been little research on how the enthesis relates to its neighbouring structures, for example the interplay of tendon cross-sectional growth and enthesis size, or, tendon injury and response within the enthesis. The latter is of particular concern, as the inferences for using entheses as indicators of repetitive stress are reliant on a model based upon a direct interplay between tendon loading and enthesis damage.

What is clear from anthropological research is that ECs are associated with increasing age [START_REF] Alves Cardoso | Enthesopathy formationin the humerus: Data from known age-at-death and known occupation skeletal collections[END_REF][START_REF] Alves Cardoso | The categorisation ofoccupation in identified skeletal collections: a source of bias?[END_REF][START_REF] Cunha | What can bones tell about labourand occupation: the analysis of skeletal markers of occupational stress in the Identified Skeletal Collection of the Anthropological Museum of the University of Coimbra (preliminary results)[END_REF][START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF][START_REF] Henderson | Accounting for multiple effects andthe problem of small sample sizes in osteology: A case study focussing on entheseal changes[END_REF][START_REF] Mariotti | Enthesopathies -Proposal of a standardized scoring method and applications[END_REF][START_REF] Mariotti | The study of entheses: Proposal of a standardised scoring method for twentythree entheses of the postcranial skeleton[END_REF][START_REF] Mariotti | Musculoskeletalstress markers (MSM): methodological reflections[END_REF][START_REF] Michopoulou | Evaluating theefficiency of different recording protocols for entheseal changes in regards to expressing activity patterns using archival data and cross-sectional geometric properties[END_REF][START_REF] Milella | The effectof age, sex, and physical activity on entheseal morphology in a contemporary Italian skeletal collection[END_REF]. This has been shown for most recording methods when tested on identified skeletal collections (ibid.). It has also been used to support the argument that the changes represent cumulative repetitive stress, but this fails to take into account normal degeneration caused by increasing age, which is a cause of tendon tears [START_REF] Hashimoto | Pathologic evidence of degeneration as a primary cause of rotator cuff tear[END_REF]. Cumulative one-off trauma may also play a role, often in association with underuse because regular use promotes healthy tendon morphology [START_REF] Hess | Achilles tendon rupture: A review of etiology,population, anatomy, risk factors, and injury prevention[END_REF]. The concept of underuse injury is beginning to be discussed in sports medicine where injuries often associated with high loads are found in those not accustomed to activity [START_REF] Stovitz | Underuse" as a cause for musculoskeletal injuries: Is it time that we started reframing our message?[END_REF].

The new Coimbra method was developed from the original Coimbra method with the aim of creating a standardised scoring system for ECs while taking into account the normal biology and scoring a broad range of features seen at fibrocartilaginous entheses [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF][START_REF] Henderson | The new"Coimbra Method": A biologically appropriate method for recording specific features of fibrocartilaginous entheseal changes[END_REF][START_REF] Henderson | Commentary: Anupdate to the new Coimbra method for recording entheseal changes[END_REF]. These entheses consist of four tissues of differing mechanical properties that mediate the transfer of stress from the tendon to the bone [START_REF] Benjamin | The skeletal attachment of tendons -Tendon "entheses[END_REF]. The method separates entheses into two zones based on the anatomy of the enthesis: i.e. zone 2 the most fibrocartilaginous part and zone 1, which is more fibrous running along the margin of the enthesis at the most obtuse angle of tendon attachment. Alongside this, it focuses on recording the types of changes seen macroscopically, "called features" that are textural change, bone formation, erosions, fine porosity, macroporosity and cavitations. The scoring consists of three scores (as well as a score for unobservable) for most features except textural change that is only scored as present or absent as it is hard to identify unless it covers a relatively large surface area.

The original method was tested on a sample of 31 individuals with similar occupations curated as part of the SIMON collection in Geneva [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF]. It demonstrated first that left and right sides typically showed similar changes and that most of the changes were associated with ageing; although the small sample size and limited variability meant that the effects could not be fully interpreted. This paper repeats that approach using the new Coimbra method [START_REF] Henderson | The new"Coimbra Method": A biologically appropriate method for recording specific features of fibrocartilaginous entheseal changes[END_REF] on a larger sample of male skeletons, all with the same occupation. The aims are to test the relationship between ECs and age for each feature and determine bilateral asymmetry for ECs. To test the overall relationship between EC features and age a random forest approach was used, creating a prediction of age that could be compared to the real age, enabling the age effect to be determined.

Materials and methods

The Coimbra identified skeletal collection contains skeletal remains and documentary evidence of age-at-death, sex, occupation and cause of death for over 500 individuals who died in the early 20th century [START_REF] Rocha | Les collections ostéologiques humaines identifiées du Musée Anthropologique de l'Université de Coimbra[END_REF]. Diverse occupations are represented, but the predominant male occupation is 'trabalhador', meaning labourer. These individuals are likely to have predominantly worked the land as farm labourers undertaking seasonal work [START_REF] Cardoso | A portrait of gender in two 19th and 20thcentury Portuguese populations: A palaeopathological perspective[END_REF]. To fully capture age effects, male individuals aged from 16 upwards were included in this study. The minimum age of 16 was chosen because by this age all entheses examined in this study are fully developed and observable.

Several diseases, e.g. seronegative spondyloarthropathies and diffuse idiopathic skeletal hyperostosis, are widely described in the clinical literature as causing changes to the entheses [START_REF] Benjamin | The skeletal attachment of tendons -Tendon "entheses[END_REF][START_REF] Benjamin | The anatomical basis for disease localisation in seronegative spondyloarthropathy at entheses and related sites[END_REF][START_REF] Benjamin | The "enthesis organ" concept: Why enthesopathies may not present as focal insertional disorders[END_REF][START_REF] Mcgonagle | Imaging the joint and enthesis: Insightsinto pathogenesis of psoriatic arthritis[END_REF]. To exclude individuals whose ECs may be pathological, all individuals with more than two vertebrae ankylosed by new osseous tissue along the line of the anterior longitudinal ligament were excluded from the analysis [START_REF] Henderson | When hard work is disease: The interpretation of enthesopathies[END_REF]. Individuals with ligamentous ankylosis of the sacroiliac joint with at least two vertebrae ankylosed were also excluded (ibid.).

The fibrocartilaginous entheses included in this analysis were the infra-and supraspinatus insertion, subscapularis insertion, common extensor origin, common flexor origin (all located on the humerus) and biceps brachii insertion (on the radius). All ECs were recorded using the new Coimbra method [START_REF] Henderson | The new"Coimbra Method": A biologically appropriate method for recording specific features of fibrocartilaginous entheseal changes[END_REF], although the biceps brachii insertion was recorded using the footprint published in the first version of the method [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF]. The infra-and supraspinatus insertions were scored as one insertion because the fibres are known to merge in some areas close to the enthesis [START_REF] Minagawa | Humeral attachmentof the supraspinatus and infraspinatus tendons: an anatomic study[END_REF]. These entheses were all recorded by one observer (CH) with interand intraobserver error for this new method previously reported [START_REF] Wilczak | Trainingand interobserver reliability in qualitative scoring of skeletal samples[END_REF].

There are a total of six features recorded in the new Coimbra method [START_REF] Henderson | The new"Coimbra Method": A biologically appropriate method for recording specific features of fibrocartilaginous entheseal changes[END_REF]. Two features, bone formation [BF(Z1)] and erosion [ER(Z1)], are scored in zone 1. In zone 2, these features (labelled Z2) plus textural change (TC), fine porosity (FPO), macroporosity (MPO) and cavitations are scored. All, apart from textural change, are scored from 0 (absence of change) to 2 (maximal expression), while textural change is only scored as absence (0) or presence [START_REF] Lieverse | Upperlimb musculoskeletal stress markers among middle Holocene foragers of Siberia's Cis-Baikal region[END_REF]. The scoring has been simplified and standardised to three scores since the original method's publication [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF] and textural change, which was scored as a type of bone formation, is now scored separately.

The variability of scores for each enthesis and feature is described, as it has been found that some features are more common at some entheses than others. Asymmetry of scores by enthesis and feature was calculated by counting the number of individuals with equal scores on both sides; with right side scores higher and with left side scores higher [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF]. Frequencies for these results were also reported to enable comparison.

Ordinal regression was used to test the effect of age for each enthesis and feature within a single occupation group, as was done in the 2013 paper [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF] to reduce potential biomechanical effects. The skeletons (N = 51) meeting the above criteria had an age range from 16 to 66 (mean = 41.29, standard deviation = 13.94). The Shapiro-Wilk test of normality showed a normal distribution of age-at-death W = 0.96, P = 0.97. To study the effect of age, means with standard errors were plotted for each enthesis, alongside descriptive statistics. Where there was sufficient variability, ordinal regression was performed using a cumulative link model [START_REF] Christensen | Ordinal: regression models for ordinaldata[END_REF] in R version 3.2.2. A log-log link was used, as in the previous study, because lower scores were more probable than higher ones [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF]. Ordinal regression was performed for each feature, enthesis and side. Nagelkerke's pseudo-R-square was calculated to determine the effect of age using an R script [START_REF] Mangiafico | An R companion for the handbook of biological statistics, version 1.09[END_REF]. Ordinal regression was only performed where there was variability of scores based on the presence of at least five occurrences of changes (score 1 or 2).

Previous analyses using this recording method have only tested each enthesis and feature individually for the effect of age. While it was demonstrated that each feature relates to age differently [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF], no test has yet been performed using the variability from all entheses and all features simultaneously to study the impact of age. Our new approach is to create a prediction model of age, which can then be compared with true age to determine its effect. This could not be achieved with the labourers alone due to the relatively small sample size and its limited variability. For this study all males (N = 260) without signs of diseases associated with ECs (see above) were included. The age range was the same as for the ordinal regression study (mean = 44.8, standard deviation = 17.9). Random forests are a combination of regression trees and a bootstrap approach. Regression trees aim to partition the space into small homogeneous regions [START_REF] Breiman | Classification andregression trees[END_REF]. The trees are designed to build subsets with low within-variance and high between-variance, thereby partitioning the data according to an algorithm of rules. Thus, the algorithm will define some homogeneous regions of the feature space where the individuals have similar ages based on their overall EC expression. For random forests, several hundred regression trees are built, each of them using only a given portion of individuals and variables [36]. The overall prediction is obtained by taking the mean of the predictions of all trees. Random forests can handle ordinal predictors and have their own algorithm for missing value imputation. For this statistical analysis, the R package random forest was used [START_REF] Liaw | Classification and regression by random forest[END_REF]. The age of each individual was predicted using leave one out cross-validation. To determine the effect of age, the predicted age is plotted against the actual age such that if age is a major aetiology of these changes, the predicted age and the actual age should overlap, showing strong concordance.

Results

The overall distribution of changes (Table 1) demonstrates the low variability of scores with very few of the highest scores present. Only bone formation in zone 1 [BF(Z1)] shows a consistent presence of score 2 with the majority of individuals having scores of zero. No doubt this is also the reason most of the scores are equal between left and right sides (Table 2). Where there is asymmetry, the right side tends to have higher scores. It is noteworthy that the left side scores are more commonly higher than the right for all biceps brachii features.

The effect of age typically showed an increase in mean age with higher scores for features, but this was not the case for FPO or TC where higher scores often occurred in younger individuals (Fig. 1). The limited variability of scores meant that ordinal regression could not be performed for 3). Fine porosity also had an association with age for two of the entheses, the left infra-and supraspinatus and the right common flexor origin, with a maximum pseudo-Rsquare of 32%. The only other statistically significant association was for erosions in zone 2 [ER(Z2)] in the left infraand supraspinatus insertion with a pseudo-R-squared of only 21%. Using the whole data set, random forests support the ordinal regression by a minimal concordance between the true age and that predicted by the random forests (Fig. 2).

Discussion

Anthropological research has focussed on two aspects of EC presence: identifying occupations (activity-levels or tasks) [1-2,6-7,38] and ageing [START_REF] Alves Cardoso | Enthesopathy formationin the humerus: Data from known age-at-death and known occupation skeletal collections[END_REF][START_REF] Alves Cardoso | The categorisation ofoccupation in identified skeletal collections: a source of bias?[END_REF][START_REF] Cunha | What can bones tell about labourand occupation: the analysis of skeletal markers of occupational stress in the Identified Skeletal Collection of the Anthropological Museum of the University of Coimbra (preliminary results)[END_REF][START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF][START_REF] Henderson | Accounting for multiple effects andthe problem of small sample sizes in osteology: A case study focussing on entheseal changes[END_REF][START_REF] Michopoulou | Evaluating theefficiency of different recording protocols for entheseal changes in regards to expressing activity patterns using archival data and cross-sectional geometric properties[END_REF][START_REF] Milella | The effectof age, sex, and physical activity on entheseal morphology in a contemporary Italian skeletal collection[END_REF]. To determine whether it is possible to identify specific tasks or movements or even broader patterns of activity, e.g. heavy manual or nonmanual occupations, based on EC expression requires the understanding of their aetiology. From an anthropological perspective, it is easy to construct studies to take into consideration biological sex and age, and these studies have shown that this is necessary due to the increase in EC presence with age [START_REF] Alves Cardoso | Enthesopathy formationin the humerus: Data from known age-at-death and known occupation skeletal collections[END_REF][START_REF] Alves Cardoso | The categorisation ofoccupation in identified skeletal collections: a source of bias?[END_REF][START_REF] Cunha | What can bones tell about labourand occupation: the analysis of skeletal markers of occupational stress in the Identified Skeletal Collection of the Anthropological Museum of the University of Coimbra (preliminary results)[END_REF][START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF][START_REF] Henderson | Accounting for multiple effects andthe problem of small sample sizes in osteology: A case study focussing on entheseal changes[END_REF][START_REF] Michopoulou | Evaluating theefficiency of different recording protocols for entheseal changes in regards to expressing activity patterns using archival data and cross-sectional geometric properties[END_REF][START_REF] Milella | The effectof age, sex, and physical activity on entheseal morphology in a contemporary Italian skeletal collection[END_REF]. Previous studies of other methods have shown an increase in EC presence in fibrocartilaginous entheses with age-at-death, particularly in individuals aged over 50 [START_REF] Villotte | Enthesopathies asoccupational stress markers: Evidence from the upper limb[END_REF], while others have shown an increase with age that is more pronounced for mineralised tissue forming changes, with less of an effect on osteolytic changes [START_REF] Milella | The effectof age, sex, and physical activity on entheseal morphology in a contemporary Italian skeletal collection[END_REF]. The test of the Coimbra method to identify the effect of age on a variety of skeletal changes in 2013 also showed that the bone formation feature in either zone was most obviously affected by age [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF]. This updated study has found similar results with bone formation most consistently associated with age. Limited variability for many of the other scores made the effect of age hard to trace. It is clear that each enthesis is affected differently or has different typical expressions of ECs, making an overall pattern hard to determine. Fine porosity and textural change, in contrast to bone formation, showed higher scores in the lower ages. However, this was not completely consistent across all entheses. This is very similar to the original paper, which did not score textural change separately [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF][START_REF] Henderson | The new"Coimbra Method": A biologically appropriate method for recording specific features of fibrocartilaginous entheseal changes[END_REF]. Textural change was, in the original 2013 paper, scored as part of bone formation, and this update shows the importance of distinguishing between these, as the effect of age is very different for both features [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF].

This raises an important point: what is the aetiology of these features? Enthesophytes are known to be associated with traction [START_REF] Abreu | Plantar calcanealenthesophytes: New observations regarding sites of origin based on radiographic, MR imaging, anatomic, and paleopathologic analysis[END_REF][START_REF] Thurston | Bone spurs: mechanism of production of different shapes based on observations in Dupuytren's diathesis[END_REF], although others have found an association with compression [START_REF] Kullar | A study of calcanealenthesophytes (spurs) in Indian population[END_REF][START_REF] Menz | Plantar calcaneal spurs in older people: longitudinal traction or vertical compression?[END_REF]. Research also demonstrates that enthesophytes could be caused by wear and tear or they have the effect of changing loading patterns thus minimising stress-related damage [START_REF] Kumai | Heel spur formation and the subcalcaneal enthesis of the plantar fascia[END_REF]. Bone formation in zone 2 has also been examined histologically with woven bone as well as bone formed through endochondral ossification identified [START_REF] Benjamin | Microdamageand altered vascularity at the enthesis-bone interface provides an anatomic explanation for bone involvement in the HLA-B27-associated spondylarthritides and allied disorders[END_REF]. The tidemark itself, the junction between the soft nonmineralised fibrocartilage and its mineralised form, has been found to be duplicated as it is in osteoarthritic joint cartilage, indicating that similar degenerative processes are at play [START_REF] Shaw | Structure-function relationships of entheses in relation to mechanical load and exercise[END_REF]. While this latter change cannot be identified visually in skeletal remains, it is worth bearing in mind when considering the presence of other mineralised tissue formations.

Pores and fissures have been identified histologically and some of these, those larger than one millimetre in diameter, have sometimes been found to be cysts [START_REF] Benjamin | Microdamageand altered vascularity at the enthesis-bone interface provides an anatomic explanation for bone involvement in the HLA-B27-associated spondylarthritides and allied disorders[END_REF]. In others blood vessels (sometimes with accompanying nerves) penetrate from the bone marrow to the edge of the mineralised fibrocartilage providing a link between the bone marrow and the enthesis such that new tissue can be formed and damage can be healed (ibid.). These are probably what anthropologists identify as porosity. In some areas of the enthesis, there is no layer of cortical bone, just mineralised fibrocartilage on top of trabecular bone (ibid.), which is likely to be identified as porosity or as an erosion in our recording system. Only elderly cadaver specimens are available to the anatomist, so finding these changes in younger individuals in skeletal remains may help to determine whether these are a response to wear and tear or are present earlier. Large fissures and macropores are highly visible on the rotator cuff entheses in adolescents during the development of the humeral head, and it is possible that some of the pores and erosions are left over from this period, but this requires further research.

Skeletal research, therefore, has a large role to play in identifying patterns of enthesis feature presence, which should improve our understanding of the causes of the changes described. What is clear from this specific piece of research is that there is a clear effect of increasing change presence with age for all changes, except fine porosity and textural change. However, it is unclear what this 'age effect' is. Ageing has numerous components and biomechanically can be considered in terms of normal age-related tissue degeneration, overuse (cumulative repetitive stress), cumulative one-off trauma or underuse injuries. The impact of these is harder to infer from skeletal remains than from in vivo studies. However, if normal age-related tissue degeneration were the cause, then left and right sides would be expected to show degeneration at the same rates. This study of asymmetry showed that most entheses and features had the same scores, indicating that the processes were occurring at the same rates in both sides (although it should be noted that this is also likely to be an effect of large numbers of zero scores). Where there were differences, there was a clear right side bias for all entheses, except the biceps brachii insertion. A clear side difference may indicate the Comparing the results of this paper to the original method effects of use, but further research with a larger sample size paper shows some stark differences. Only three enthesis feais needed to explore which use effects are more likely to be tures demonstrated higher scores on the left side in the origicausing the changes, due to the small variations in this study nal paper; these were erosions at the biceps brachii insertion and the limited activity patterns. (both ER(Z1) and ER(Z2)) and FPO at the common extensor Fig. 2 True age (in years) plotted against the random forest predicted age (in years). Line indicates concordance between the true and predicted age / Nuage de point suivant l'âge réel (en années) et l'âge prédit (en années) par les forêts aléatoires. La droite illustre une concordance parfaite origin [START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF]. While some of this may be due to a difference in the method, it is also probable that population differences and sampling effects play a role. However, biomechanical differences may also explain these results and further exploration of asymmetry is needed. In terms of the effect of age, similar patterns are seen for all features (excluding textural change which was not identified separately in the original method). The effect of age was then, as now, found to be present but was clearly not the only effect nor could the underlying cause of the effect of ageing be identified. Both studies suffer from the limited sample size available for one occupational group. Future studies are needed using larger samples with the same occupations to better characterise asymmetry and the effect of age. To study the cause of the age-related effects, large samples with multiple occupations are required, and improved methods to characterise and categorise occupations are also needed to enable identification of overuse and underuse effects.

Conclusions

The aim of this paper is to present the results of asymmetry and the effect of age on entheses and their features in a single occupation category. The study shows that feature expression varies by entheses with some features occurring only rarely at some entheses, at least in this sample.

Asymmetry is rare, with predominantly higher scores on the right side where it does occur. The effect of age is minimal overall. Where there is an effect of age it tends to be towards an increased score with age, except for fine porosity and textural change, which often have a higher score in younger individuals and may indicate a developmental origin for these changes. While there is an effect of age, it is unclear whether this is from normal age-related tissue degeneration, cumulative repetitive movement (overuse), cumulative one-off trauma or underuse injuries. Further research with larger samples of skeletal remains is needed to understand these effects. Longitudinal in vivo studies are also required to improve our understanding of the relationship between tendon degeneration and ECs.

  4 31 100.0 25 73.5 28 82.4 33 97.1 26 76.5 34 100.0 100.

Table 1

 1 Descriptive statistics showing the variability of enthesis scores / Statistiques descriptives illustrant la variabilité des scores The feature with the most variable scores, bone formation both BF(Z1) and BF(Z2), had the most consistent statistical significance (at 95%) with age, but the pseudo-R-squared showed that the effect of age was typically minimal, at most explaining 44% of the variability Table2Asymmetry scores for entheses by feature. Bold indicates a large number of higher scores on the left side / Score d'asymétrie pour chaque enthèse, par type de changement. Une dominance à gauche est indiquée en gras

	Enthesis	Side	Enthesis score	BF(Z1)	ER(Z1)	BF(Z2)	FPO	MPO	ER(Z2)	CA	TC
	Infra-	Left	NA	17		17	8	8	8	8	8	8
	an supraspinatus		0	31		33	33	19	37	26	42	42
			1	3	1	9	21	6	16	1	1
			2	0	0	1	3	0	1	0	0
		Right	NA	16		16	9	9	10	10	9	9
			0	31	32	30	19	38	20	41	41
			1	3	3	10	22	2	17	1	1
			2	1	0	2	1	1	4	0	0
	Subscapularis	Left	NA	4		4	4	4	4	4	4	4
			0	28	41	16	24	39	27	45	44
			1	12	6	26	23	8	19	2	3
			2	6	0	5	0	0	1	0	0
		Right	NA	9		9	5	5	5	5	5	5
			0	23	35	17	20	42	27	41	46
			1	13	6	24	25	4	18	5	0
			2	6	1	5	1	0	1	0	0
	Common extensor Left NA 15 15 12 12 12 12 12 12 origin 0 30 36 35 37 39 33 39 39
			1	5	0	4	2	0	5	0	0
			2	1	0	0	0	0	1	0	0
		Right	NA	17		17	9	9	9	9	9	9
			0	23	34	29	37	41	34	42	42
			1	6	0	11	2	1	8	0	0
			2	5	0	2	3	0	0	0	0
	Common flexor Left NA 18 18 13 13 13 13 13 13 origin 0 30 33 30 35 38 37 38 38
			1	3	0	6	3	0	1	0	0
			2	0	0	2	0	0	0	0	0
		Right	NA	20		20	15	14	14	14	14	14
			0	25	31	27	32	37	33	37	36
			1	4	0	7	5	0	4	0	1
			2	2	0	2	0	0	0	0	0
	Biceps brachii	Left	NA	6		6	2	2	2	2	2	2
			0	25	43	30	34	44	44	48	31
			1	10	2	15	11	4	5	1	18
			2	10	0	4	4	0	0	0	0
		Right	NA	6		6	7	7	7	7	7	7
			0	26	45	33	31	43	42	44	34
			1	9	0	10	9	1	0	0	9
			2	10	0	1	4	0	2	0	1
	Z1 refers to zone 1, Z2 to zone 2; BF: bone formation; ER: erosions; FPO: fine porosity; MPO: macroporosity; CA: cavitations; TC:
	textural change. For definitions, see Henderson et al. [26] / Z1 correspond à zone 1, Z2 à zone 2 ; BF : formation osseuse ; ER : érosions

; FPO : porosité fine ; MPO : macroporosité ; CA : géodes ; TC : changement de texture. Pour les définitions, voir Henderson et al.

[START_REF] Henderson | The new"Coimbra Method": A biologically appropriate method for recording specific features of fibrocartilaginous entheseal changes[END_REF]

all entheses and features.

Table 3

 3 Results of the ordinal regression / Résultats de la régression ordinale

	Enthesis	Side	Feature	Estimate Std. Error z value	Pr(>z)	Nagelkerke	2013 paper
								pseudo R	Nagelkerke
								squared	pseudo R
									squared
	Infra-	Left	BFZ1	na	na	na	na	na	-
	and supraspinatus		ERZ1	na	na	na	na	na	-
			BFZ2	0.038	0.025	1.542	0.123	0.082	-
			ERZ2	0.054	0.021	2.518	0.012 *	0.205	-
			FPO	-0.038	0.018	-2.120	0.034 *	0.125	-
			MPO	0.027	0.031	0.866	0.386	0.033	-
			CA	na	na	na	na	na	-
			TC	na	na	na	na	na	-
		Right	BFZ1	na	na	na	na	na	-
			ERZ1	na	na	na	na	na	-

Table 3

 3 (suite) NA indicates lack of variability,indicates that the enthesis or feature was not tested, bold indicates a pseudo R-squared score > 0.29 and *indicates level of significance. The pseudo R-squared for this study is presented next to that from the 2013 paper[START_REF] Henderson | Recordingspecific entheseal changes of fibrocartilaginous entheses: initial tests using the Coimbra method[END_REF]. Abbreviations:

	Enthesis	Side	Feature	Estimate Std.	z value	Pr(>z)	Nagelkerke	2013 paper
						Error				pseudo R	Nagelkerke
										squared	pseudo R
											squared
				MPO	na	na	na		na	na	na
				CA	na	na	na		na	na	na
				TC	na	na	na		na	na	-
	Common flexor Left	BFZ1	na	na	na	na	na	-origin ERZ1	na	na
	na	na	na	-
				BFZ2	0.071	0.033	2.188	0.029 *	0.195	-
				ERZ2	na	na	na		na	na	-
				FPO	na	na	na		na	na	-
				MPO	na	na	na		na	na	-
				CA	na	na	na		na	na	-
				TC	na	na	na		na	na	-
		Right	BFZ1	0.053	0.033	1.610	0.107	0.121	-
				ERZ1	na	na	na		na	na	-
				BFZ2	0.052	0.028	1.846	0.065	0.133	-
				ERZ2	na	na	na		na	na	-
				FPO	-0.107	0.046	-2.329	0.020 *	0.322	-
				MPO	na	na	na		na	na	-
				CA	na	na	na		na	na	-
				TC	na	na	na		na	na	-
	Biceps brachii	Left		BFZ1	0.047	0.018	2.651	0.008 **	0.190	na
				ERZ1	na	na	na		na	na	na
				BFZ2	0.038	0.018	2.161	0.031 * 0.119	0.04 ERZ2	0.040
				0.036	1.108	0.268	0.056	na
				FPO	-0.004	0.018	-0.223	0.823	0.001	0.10
				MPO	na	na	na		na	na	na
				CA	na	na	na		na	na	na
				TC	-0.033	0.019	-1.788	0.074	0.090	-
		Right	BFZ1	0.051	0.020	2.592	0.010 **	0.186	0.15
				ERZ1	na	na	na		na	na	na
				BFZ2	0.080	0.030	2.662	0.008 **	0.270	0.04
				ERZ2	na	na	na		na	na	na
				FPO	0.013	0.020	0.641	0.522	0.012	0.08
				MPO	na	na	na		na	na	na
				CA	na	na	na		na	na	na

Table 1 /

 1 NA indique une absence de variabilité, -indique que cette enthèse ou ce changement n'a pas été testé, les valeurs en gras indiquent un score de pseudo R carré supérieur à 0,29 et * indique une valeur p inférieure à 0,05. Abréviations : voir tableaux 1

Acknowledgements This paper is based on a presentation at the Working your Fingers to the Bone Conference, held in July 2016 in Coimbra, Portugal. The authors would like to express their thanks to the Department of Life Sciences, University of Coimbra, for providing access to the identified skeletal collection. Charlotte Henderson is funded by Fundação para a Ciência e a Tecnologia grant number SFRH/BPD/82559/2011.

Author's contribution

All authors contributed to the research design of the study. Henderson collected the data and undertook all statistics excluding the random forest tests, which were run by Santos. All authors reviewed and critiqued the drafting, contributing to the final version. The authors have no conflicts of interest.

Conflict of interest: None.


Diying Huang
email: dyhuang@nigpas.ac.cn

Fu

Yanzhe

Lian Xinneng

André Nel
email: anel@mnhn.fr

Sinaspidoneura magnifica nov. gen., nov. sp., first Chinese Caloneurodea (Insecta: Archaeorthoptera)

Keywords: Aspidoneuridae, paleobiogeography, Laurussia, Southern China, nov. gen., nov. sp., Middle Permian extinction

We describe the first Chinese Caloneurodea, Sinaspidoneura magnifica nov. gen., nov. sp. from the middle Permian Yinping Formation, in the small family Aspidoneuridae, previously known from two genera and species, one from the latest Carboniferous of France and another from the late early Permian of North America. This discovery shows that this order was more widespread during the Middle Permian than previously supposed, under a great variety of palaeoclimates. This clade is still unknown in the Late Permian, and possibly became extinct because of the crisis of biodiversity that happened at the end of the Middle Permian.

Introduction

The small extinct order Caloneurodea [START_REF] Martynov | Permian fossil insects from the Arkhangelsk district. Part 5. The family Euthygrammidae and its relationships, with the description of a new genus and family from Chekarda[END_REF] is currently only known between the late Carboniferous and the end of the middle Permian, in Western Europe, North America, and Russian Federation. With 22 described genera, it is one of the smallest order of Paleozoic insect orders. Its position was uncertain till [START_REF] Béthoux | The extinct order Caloneurodea (Insecta: Pterygota: Panorthoptera): wing venation, systematics and phylogenetic relationships[END_REF] discovered that their pattern of wing venation was that of the superorder Archaeorthoptera [START_REF] Béthoux | Venational pattern and revision of Orthoptera sensu n. and sister groups. Phylogeny of Palaeozoic and Mesozoic Orthoptera sensu n[END_REF], the total group of the recent Orthoptera Olivier, 1789 and some other extinct clades such as the Titanoptera [START_REF] Sharov | Filogeniya ortopteroidnykh nasekomykh. Trudy Paleontologicheskogo Instituta[END_REF]. The Caloneurodea are well-characterized by a series of characters of the forewing venation, the most obvious being the veins 'MP+CuA+CuPaα and CuPaβ close, parallel and straight'. Very little is known about their body structures, young stages and biology. They could adapt to the important climatic changes and environments that happened between the Carboniferous and the middle Permian. They were widespread in the Laurussia but still unknown elsewhere during the Middle Permian.

For China, Hong and Wang (1976: 82-83, pl. 51, fig. 5-5a) described Sinogramma reticularis in the caloneurodean Pleisiogrammatidae [START_REF] Carpenter | The Lower Permian insects of Kansas. Part 9. The orders Neuroptera, Raphidiodea, Caloneurodea and Prothortoptera (Prosbnisidae), with additional Protodonata and Megasecoptera[END_REF] from the Lower Cretaceous Guyang Formation (Inner Mongolia). It is in fact the cubito-anal area of a Dictyoptera [START_REF] Leach | The zoological miscellany; being descriptions of new interesting animals[END_REF]. Hong (1984: 147-148, fig. 22, pl. 63, fig. 3) described another early Cretaceous taxon (from Hebei) he attributed to the Caloneurodea, in the family Mesogrammatidae [START_REF] Hong | Tracheata, Insecta[END_REF], for the genus and species Mesogramma divaricata [START_REF] Hong | Tracheata, Insecta[END_REF] these names being later replaced by Hebeigrammidae, and Hebeigramma divaricata [START_REF] Hong | Hebeigramma nom. nov., a new name for Mesogramma Hong, 1984 (Caloneurodea) from the Lower Cretaceous of Hebei Province, China[END_REF]. This fossil has not the main synapomorphy of the Caloneurodea, viz. veins 'MP+CuA+CuPaα and CuPaβ close, parallel and straight'. After the original photograph of the holotype, the general shape of venation rather fits with a Phasmatodea Jacobson and Bianchi, 1902, but a revision of the type material would be necessary to precise its affinities.

Here we report the first discovery of a Chinese Caloneurodea from the middle Permian Yinping Formation near the Houdong Village, Chaohu City, Anhui Province, Southern China, corresponding at that time to a group of islands in the eastern margin of the paleo-Thetys, under a climate very different from those of the northern part of the Permian Pangea.

Material and methods

The fossil was preserved in the black shale of the Yinping Formation at the Paomaling Section near the Houdong Village, Chaohu City, Anhui Province, Southern China (detail map see Szwedo and Huang, 2019, fig. 2). Its age was assigned within the late part of Capitanian, Middle Permian (Zhang et al., 2018), or considered around the boundary of middle-late Permian [START_REF] Yao | Astronomical forcing of a Middle Permian chert sequence in Chaohu, South China[END_REF]. Our fossil was collected from the lower part of the formation that would be more convincing of the latest part of the middle Permian age.

The fossil wing was prepared using a sharp knife. Photographs were taken using a digital camera attached to a Zeiss Discovery V16 microscope. Line drawings were drafted with Adobe Illustrator CC 2018 graphic software.

We follow the wing venation terminology of [START_REF] Béthoux | Venational pattern and revision of Orthoptera sensu n. and sister groups. Phylogeny of Palaeozoic and Mesozoic Orthoptera sensu n[END_REF] and the classification of the Caloneurodea of [START_REF] Béthoux | The extinct order Caloneurodea (Insecta: Pterygota: Panorthoptera): wing venation, systematics and phylogenetic relationships[END_REF]. Abbreviations of wing venation: ScP subcostal posterior; RA radius anterior; RP radius posterior; MA media anterior; MP media posterior; CuA cubitus anterior; CuP cubitus posterior, AA anal anterior. urn:lsid:zoobank.org:pub:E83918BF-B83E-4381-84AF-2087540FBEE7

Systematic paleontology

Superorder Archaeorthoptera [START_REF] Béthoux | Venational pattern and revision of Orthoptera sensu n. and sister groups. Phylogeny of Palaeozoic and Mesozoic Orthoptera sensu n[END_REF]Subclade Panorthoptera Crampton, 1928(sensu Béthoux and Nel 2002) Order Caloneurodea [START_REF] Martynov | Permian fossil insects from the Arkhangelsk district. Part 5. The family Euthygrammidae and its relationships, with the description of a new genus and family from Chekarda[END_REF] Family Aspidoneuridae [START_REF] Carpenter | Studies on Carboniferous insects of Commentry, France. Part 3. The Caloneurodea[END_REF] Type genus. Aspidoneura [START_REF] Carpenter | The Lower Permian insects of Kansas. Part 9. The orders Neuroptera, Raphidiodea, Caloneurodea and Prothortoptera (Prosbnisidae), with additional Protodonata and Megasecoptera[END_REF]. Other genera. Homaloptila Handlirsch, 1919, Sinaspidoneura gen. nov. Genus Sinaspidoneura nov. urn:lsid:zoobank.org:act:99DD11E0-CE59-4E1C-9E22-D1A323BDE239 Derivation of the name: Named after 'Sina', Latin name for China, and Aspidoneura. Gender feminine.

Type species: Sinaspidoneura magnifica nov. sp., here designed and by monotypy.

Diagnosis: Forewing ScP extending to wing apex; anterior branch of MA strongly sigmoidal; three branches of RP; base of RP very close to that of M+CuA, just distal to base of MA; vein AA2 absent or strongly reduced; wing basally petiolate with a reduced anal area.

Sinaspidoneura magnifica nov. gen., nov. sp. Paratype with same venation (Fig. 3).

Discussion

Sinaspidoneura nov. gen. falls in the Caloneurodea because of the following characters [START_REF] Béthoux | The extinct order Caloneurodea (Insecta: Pterygota: Panorthoptera): wing venation, systematics and phylogenetic relationships[END_REF] Sinaspidoneura nov. gen. differs from all the described Caloneurodea, except Aspidoneura flexa [START_REF] Carpenter | The Lower Permian insects of Kansas. Part 9. The orders Neuroptera, Raphidiodea, Caloneurodea and Prothortoptera (Prosbnisidae), with additional Protodonata and Megasecoptera[END_REF] and Homaloptila similis [START_REF] Meunier | Nouveaux insectes du houiller de Commentry[END_REF], in the conjunction of the following characters: 'ScP extending to wing apex' and 'anterior branch of MA strongly sigmoidal' [START_REF] Zalessky | Sur de nouveaux insectes permiens[END_REF][START_REF] Carpenter | The Lower Permian insects of Kansas. Part 9. The orders Neuroptera, Raphidiodea, Caloneurodea and Prothortoptera (Prosbnisidae), with additional Protodonata and Megasecoptera[END_REF][START_REF] Carpenter | Fossil insects from New Mexico[END_REF][START_REF] Béthoux | The extinct order Caloneurodea (Insecta: Pterygota: Panorthoptera): wing venation, systematics and phylogenetic relationships[END_REF]. Some taxa, viz. [START_REF] Carpenter | Carboniferous insects from Pennsylvania in the Carnegie Museum and the Museum of Comparative Zoology[END_REF] and Isadistica issada, have a long ScP nearly reaching wing apex, but they differ from Sinaspidoneura nov. gen. in the presence of four branches of RP or more, and the anterior branch of MA parallel to the posterior branch, with a narrower area between the two than in Sinaspidoneura nov. gen. [START_REF] Carpenter | Carboniferous insects from Pennsylvania in the Carnegie Museum and the Museum of Comparative Zoology[END_REF][START_REF] Rasnitsyn | 2.1. New Fossil Insects (Insecta: Caloneurida, Hypoperlida, Palaeomanteida, Jurinida) from the Middle and Upper Permian of European Russia[END_REF]. Pleisiogramma medialis [START_REF] Carpenter | The Lower Permian insects of Kansas. Part 9. The orders Neuroptera, Raphidiodea, Caloneurodea and Prothortoptera (Prosbnisidae), with additional Protodonata and Megasecoptera[END_REF] has also a very long ScP, but a simple MA and four branches of RP. Synomaloptila longipennis Martynov, 1938 has a long ScP, but four branches of RP and an anterior branch of MA strongly curved but not sigmoidal [START_REF] Martynov | Permian fossil insects from the Arkhangelsk district. Part 5. The family Euthygrammidae and its relationships, with the description of a new genus and family from Chekarda[END_REF].

Caloneurella carbonaria

Sinaspidoneura nov. gen. differs from Aspidoneura flexa in the less sigmoidal anterior branch of MA, the presence of three branches of RP, instead of two, and the base of RP very close to that of M+CuA, just distal to base of MA (Béthoux et al., 2004, fig. 15). Homaloptila similis has four branches of RP, but also a more pronounced sigmoid in the anterior branch of MA and the base of RP well distad of the base of MA (Béthoux et al., 2004, fig. 16).

Aspidoneura and Homaloptila were classified together in the same family Aspidoneuridae with the following diagnosis: 'Sc extending to about the apex of the wing ; CuA and CuP very close together and nearly parallel; MP forked broadly and unevenly, the anterior branch (MP1+2) arching strongly away from MP3+4 ; 3 anal veins' (Carpenter, 1961, p. 151), that is, with the current terminology: ScP extending to wing apex; MP+CuA+CuPaα and CuPaβ very close together and nearly parallel (a character present in nearly all the Caloneurodea); anterior branch of MA strongly away from its posterior branch; three veins CuPb, AA1 and AA2 present.

Sinaspidoneura nov. gen. has all these characters except the presence of a vein AA2.

More precisely, if AA2 is present, it should be a very short vein in the basally petiolate wing.

The strongly sigmoidal anterior branch of MA is a potential synapomorphy supporting the family Aspidoneuridae. The polarization of the character 'ScP extending to wing apex' is more uncertain. Thus we propose to place Sinaspidoneura nov. gen. in this family.

Conclusion

Homaloptila is known from the latest Carboniferous of Commentry in France, while Aspidoneura is known from the early Permian of Elmo (Wellington Formation, Oklahoma, USA). These two outcrops correspond to very different palaeoenvironments, viz. an equatorial hot and humid swamp-forest for Commentry while Elmo is supposed to have been 'a humid spot in a regional environment of more or less pronounced and long-continued aridity' (Carpenter, 1930, p. 72). Tasch (1964) interpreted the layers of Elmo with insects as limnic deposits; Zeller (1968: 50) indicated that the 'Wellington is predominantly shale with minor amounts of limestone and dolomite, siltstone, and gypsum and anhydrite', with 'marine and brackish-and fresh-water deposits'; [START_REF] Hall | Depositional facies and diagenesis of the Carlton Member (sic) (Kansas) and the Midco Member (Oklahoma) of the Wellington Formation (Sumner Group, Leonardian) (abs.)[END_REF] suggested a 'deposition in a marginalmarine environment'; and Beckemeyer (2000, p. 3) considered that the region was a 'coastal plain at the edge of an inland sea'. Thus the Aspidoneuridae could adapt to a more contrasted climate and environment between the Latest Carboniferous and the early Permian.

The Yinping Formation at Chaohu area, from where Sinaspidoneura nov. gen. comes, was a group of islands in the paleo-Thetys ocean, separated from the Pangea (and Laurussia) during the Middle Permian [START_REF] Campi | The Permian -A time of major evolutions and revolutions in the history of life[END_REF]Zhang et al., 2019, fig. 1). It was at the level of equator and probably under a much more humid and equally warm than the paleo-regions where the other Permian taxa of the Caloneurodea were found [START_REF] Chumakov | Climate during Permian-Triassic Biosphere Reorganizations. Article 1: Climate of the Early Permian[END_REF][START_REF] Wang | The Carboniferous -Permian transition on the North China microcontinent -Oceanic climate in the tropics[END_REF]. The very broad distribution of these insects during the Permian suggests that they could adapt to very different climates and biotas. Apparently the Caloneurodea are unrecorded in the late Permian. The younger ones are late Capitanian. Maybe they became extinct during the middle Permian crisis of diversity, maybe in relation to the 'major changeover in plant taxa' that happened at the end of the Guadalupian [START_REF] Rampino | The end-Guadalupian (259.8 Ma) biodiversity crisis: the sixth major mass extinction?[END_REF].

Figs

  Figs 1, 2, 3 rn:lsid:zoobank.org:act:DA69A668-67E2-4D5C-B604-3843FB61002C Derivation of the name: Named after the wonderful state of preservation of the holotype. Holotype: NIGP173204 (a complete wing, with part and counterpart, perfectly preserved); paratype NIGP173730, stored at the Nanjing Institute of Geology and Palaeontology, Chinese Academy of Sciences, Nanjing, China. Type horizon and locality: Late Capitanian, middle Permian; Yinping Formation, Houdong Village, Chaohu City, Anhui Province, Southern China. Diagnosis: As for the genus. Description: (Based on holotype specimen) preserved part of wing 19.5 mm long, wing ca. 20.0 mm long, 4.1 mm wide; ScP reaching anterior wing margin at wing apex; RA and RP diverging ca. 3.8 mm distal of wing base; RA simple; first posterior branch from RP 8.6 mm

  : in forewing, MA with only two branches; MP+CuA+CuPaα simple; fusion of CuPaα with M+CuA; MP+CuA+CuPaα (convex) and CuPaβ (concave) close, parallel and straight (main synapomorphy); MP running fused with CuA+CuPaα at length. It has also the following additional characters, present in a majority of Caloneurodea: RP posteriorly pectinate when branched; anterior branch of MA strongly sigmoidal; marked convexity of crossveins; absence of anterior branches of RA. Béthoux et al. (2004) gave a list of the species of the Caloneurodea, to which it is necessary to add the following taxa: Ligogramma wichita Beckemeyer, 2009 (early Permian, Okhlahoma, USA); Lusitaneura covensis Loureiro et al., 2010 (latest Carboniferous, Portugal); Isadistica issada Rasnitsyn et Aristov, 2013 (Isady, Vologda, Russia; Severodvinian); Aviogramma gracilis Prokop et al., 2014 (Avion, France Early Late Carboniferous);Gallogramma galadrieli[START_REF] Garrouste | Insects in the Red Middle Permian of Southern France: first Protanisoptera (Odonatoptera) and new Caloneurodea (Panorthoptera), with biostratigraphical implications[END_REF] Paleuthygramma cf. acuta Carpenter, 1943 (Lodève, France, middle Permian).

Fig. 1 .

 1 Fig. 1. Sinaspidoneura magnifica gen. et sp. nov., holotype NIGP173204, photographs of forewing. A. Part (NIGP173204a). B. Enlargement of A showing details of wing base. C. Enlargement of A showing details of wing apex. D. Counterpart (NIGP173204b). Scale bars: 2 mm (A,D), 1 mm (B,C).

Fig. 2 .

 2 Fig. 2. Sinaspidoneura magnifica gen. et sp. nov., line drawing, small black arrow: base of MA. Scale bar: 2 mm.

Fig. 3 .

 3 Fig. 3. Sinaspidoneura magnifica gen. et sp. nov., paratype NIGP173730, photographs of forewing. A. Part. B. Counterpart. Scale bars: 1 mm.

Acknowledgements

We sincerely thank Pedro Correia and an anonymous referee for the useful comments on the first version of the paper. We are grateful to Jian Gao and Chenyang Cai for field work helps. This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB18000000 and XDB26000000), and the National Natural Science Foundation of China (41925008 and 41688103).


Natarajan Vasudevan

Eric Wimmer

Elvina Barré

Daniel Cortés-Borda

Mireia Rodriguez- Zubiri

François-Xavier Felpin
email: fx.felpin@univ-nantes.fr

Mireia Rodriguez-Zubiri

Xavier Felpin

Direct C-H arylation of indole-3-acetic acid derivatives enabled by an autonomous self-optimizing flow reactor

Keywords: Self-Optimizing flow reactor, Flow chemistry, Direct C-H arylation, Autonomous synthesis, Indole

Described herein is a continuous-flow strategy for the palladium-catalyzed direct C-H arylation of indole-3acetic acid derivatives with arenediazonium salts. A fully autonomous self-optimizing flow platform was used to efficiently optimize the coupling reaction in a 3-dimensional space. The flow methodology developed is experimentally simple, mild, broad in scope and safer than traditional batch approaches. Our continuous-flow approach is particularly convenient to prepare precursors of pharmaceutically relevant compounds.

Introduction

The 2-aryl-3-carbonylmethyl indole skeleton is a ubiquitous moiety in compounds of pharmaceutical relevance (Figure 1). For instance, paullone 1 and kenpaullone 2 exhibit potent biological activities as ATP-competitive inhibitors of cyclin-dependent kinases (CDKs), glycogen-synthase kinases-3 (GSKs 3) and mitochondrial malate dehydrogenase (mMDH). [1] FS-554 3 is a potent trypanothione synthetase inhibitor with promising activity against parasites. [2] Compound 4 is the most active member of a large family of farnesyl pyrophosphatecompetitive inhibitors of farnesyl protein transferase. [3] FGIN-1-27 5 is a well-known potent and selective ligand for mitochondrial diazodiazepine receptors, producing anxiolytic effects. [4] Because these compounds are frequently used by biologists as biochemical tools and references in drug discovery, they are occasionally commercially available from specialized sellers but only on milligram scale and at a very prohibitive cost. Therefore, the development of methodologies allowing a straightforward and efficient access to 2-aryl-3-carbonylmethyl indole derivatives is of huge interest for biologists and medicinal chemists.

The traditional approach to access the 2arylindole-3-acetic acid motif involves the arylation of indole acetic acid derivatives through Suzuki-Miyaura or Stille cross-coupling reactions. [5,6] While versatile and efficient, these strategies suffer from the necessity to pre-functionalize the indole moiety at the C2 position with either a halogen atom or a metal (Sn or B). A more step-and atom-efficient approach consists of the direct C-H arylation of indole acetic acids with aryl halides using palladium-catalyzed coupling reactions. [7] While generally high yielding, this batch approach suffers from extensive reaction times and high temperatures. Our long standing experience in diazonium chemistry associated with pioneering works from Correia, [8] Noël [9] and Fairlamb [10] describing the arylations of indoles with arenediazonium salts in batch, prompted us to propose our own flow-based strategy. Indeed, with the objective of developing an efficient arylation strategy of indole acetic acid derivatives, we acknowledged that the use of arenediazonium salts as highly reactive aryl halides surrogates would allow mild experimental conditions and fast reaction rates. However, considering the high reactivity and hazardous behavior of arenediazonium salts, we envisioned the development of a flow-based strategy to mitigate the safety concerns. In this contribution, we describe our efforts toward the development of an efficient, safe and operationally simple palladiumcatalyzed direct C-H arylation strategy of indole acetic acid derivatives, assisted by an autonomous self-optimizing flow reactor. The methodology shows a broad substrate scope and is applied to the preparation of a key intermediate of farnesyl pyrophosphate-competitive inhibitors of farnesyl protein transferase 4.

Results and Discussion

The arylation of indole 6 with tolyldiazonium tosylate 7a, in the presence of Pd(OAc) 2 as catalyst, was selected as benchmark transformation to be optimized for the development of a continuous flow strategy (Scheme 1). As our objective was to develop a simple, robust and environmentally acceptable methodology, we excluded the use of additives such as ligand, base, oxidant and stabilizing agent (ammonium salt). Traditionally, the optimizations of chemical transformations conducted either in batch or in flow, are mainly based on the (often biased) know-how of chemists through one-variable-at-a-time optimizations. Whilst this popular optimization method is easily set up, it is generally poorly efficient for multi-dimensional optimizations as it neglects interactions between variables. The use of design-ofexperiment (DoE) methods partially addresses these drawbacks, but the high number of experiments required to locate a satisfactory optimum associated to the absence of feedback are strong handicaps to develop efficient automations. The development and optimization of continuous flow processes can be greatly speeded-up with the use of autonomous selfoptimizing flow reactors. [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] These powerful automated devices combine flow reactors with process control instrumentations, in-line/online analyses and optimization algorithms to assist the decision-making process of chemists. [27][28][29][30][31][32][33][34] Our lab also recently contributed to this emerging field of research with the development of a reconfigurable flow platform for the screening of reagents and the autonomous optimization of chemical transformations including natural product synthesis. [35][36][37] Based on our experience in autonomous selfoptimizing reactors, we designed an automated twostream flow setup integrating an online HPLC unit for the benchmark arylation of indole 6a with tolyldiazonium tosylate 7a under palladium catalysis (Figure 2). Before starting an optimization campaign with the flow platform, a careful assessment of the solubility profile of all known chemicals involved in the reaction is an important pre-requisite as any precipitate formed during the process would inevitably clog the flow reactor. This is particularly true when using arenediazonium tosylates which are poorly soluble in many organic solvents and generate inorganic salts. Arenediazonium tosylates are often ideally handled in MeOH or EtOH in which they are generally well soluble and stable over hours. [38] By contrast, Pd(OAc) 2 is rapidly reduced by MeOH and EtOH into Pd nanoparticles which precipitate within minutes. Considering these features, we used a solution of diazonium 7a in MeOH in the first stream while in the second stream indole 6a, Pd(OAc) 2 and naphthalene 9 (used as internal standard for HPLC analysis) were solubilized in DMF. Therefore, a solution of diazonium salt 7a was continuously pumped (0.1 M) and meet, in a T-shaped mixer, a DMF solution of indole 6a (0.1 M), Pd(OAc) 2 (5 × 10 -3 M) and naphthalene 9 (0.1 M). The reaction occurred in a stainless steel reactor coil (0.75 mm i.d., 10 mL) placed in an oven. The outlet of the reactor was connected to an automatic 2-way 6-port switch valve which injected an aliquot of the crude mixture into the HPLC line. An automated data processing allowed to integrate the chromatogram and to calculate the reaction yield. The value was automatically sent to the algorithm which proposed a new set of experimental conditions and the process control software modified the variables accordingly. The machine autonomously continues this iterative process until it reaches an optimum. The reaction yield was optimized in a 3-dimensional space where three variables, i.e., equivalent of 7a, residence time and temperature, were considered in the range of 1-2 equiv, 15-90 min, and 25-60 °C, respectively. The temperature was limited to 60 °C to preserve the diazonium salt integrity. The initial experiment X 0 was fixed at 1 equiv of 7a, 15 min of residence time and 25 °C with d value of 0.2 equiv, 15 min and 7 °C, respectively. The optimization progressed rapidly until reaching a maximum reaction yield of 62% after 17 experiments at 1.2 equiv of 7a, 65 min of residence time and 48 °C (Figure 3a, blue curve). The algorithm reached a stopping criteria in experiment 20 after 5 consecutive simplexes without improvement (see Table S1 in Supporting Information for details on experimental conditions). A close examination of the crude mixture showed the presence of a significant amount of toluene which likely arose from a rapid proto-dediazotization of 7a. We suspected that this unwanted behavior was, at least partially, favored in DMF which has been reported to promote protodediazotization of arenediazonium salts. [39] We, therefore, reconsidered the selection of solvents and reasoned that the volume of DMF, which is required to solubilize Pd species, should be significantly reduced through the use of a co-solvent. After having evaluated various combinations of DMF/co-solvent in different ratios, we selected a mixture of DMF/EtOAc in a 1/4 ratio. This choice was justified by the fact that tolyldiazonium tosylate 7a is stable in EtOAc (no degradation observed after 30 min at room temperature) which is, in turn, considered as a sustainable solvent. The DMF/EtOAc ratio (1/4, respectively) was adjusted in order to reduce the volume of DMF to a minimum while keeping a clear solution with a complete dissolution of Pd species. We re-started the optimization process with this new solvent system while conserving the same search space and starting parameters. A spectacular effect was observed on the initial simplex since the reaction yields, ranging from 5 to 14% in MeOH/DMF, increased from 19 to 41% in MeOH/DMF/EtOAc (Figure 3b, red curve). After the initial simplex, the algorithm smoothly progressed by increasing the value of all variables until reaching an optimum in the 15 th experiment at 44 °C, 74 min of residence time and 1.6 equivalents of 7a, corresponding to 90% HPLC yield (88% isolated yield, see Figure S1 in supporting information for HPLC chromatogram). The algorithm was unable to improve the reaction yield in the subsequent 5 simplexes and, therefore, reached a stopping criterion after experiment 28 (see Table S2 in Supporting Information for details on experimental conditions).

With the optimized conditions in hand, the reaction scope was examined with various arenediazonium salts and indoles, using the continuous flow setup depicted in Figure 4. The first stream, equipped with a 5 mL injection loop 1 and loaded with a solution of arenediazonium salts in MeOH M, 82 µL/min), met in a T-shaped piece a second stream equipped with a 2 mL injection loop 2 and loaded with a solution of indole (0.1 M) and Pd(OAc) 2 (5 × 10 -3 M) in DMF/EtOAc (1/4, 51 µL/min). The merged streams (133 µL/min) entered in a stainless steel reactor coil (0.75 mm id, 10 mL, 74 min residence time) heated at 44 °C, and the resulting 2-arylindole-3-acetic acid derivatives 8a-8q were collected in vials. The arylation of indole-3-acetic acid alkyl esters was explored with a broad set of arenediazonium salts. While the optimization was conducted with indole-3-acetic acid ethyl ester 6a, the reaction scope was extended to the methyl ester derivative 6b-f to avoid any unwanted transesterification process with MeOH. The arylation of indoles 6a-b was found to be compatible with arenediazonium salts bearing alkyl substituents at either para-or ortho-positions (8a-e). The excellent yield obtained for the ortho, ortho'substituted indole 8c suggests that the transformation was poorly sensitive to steric hindrance. Switching an alkyl group for an aryl substituent on the diazonium salt, did not significantly affect the reaction outcome (compounds 8g). One of the salient advantages of using arenediazonium salts as super-electrophile is the possibility to introduce very electron-rich substituents such as methoxy groups which usually significantly deactivate aryl halides and compromise the C-H arylation efficiency. Another unique advantage of arenediazonium salts is the possibility to work without any base, resulting in a high chemoselectivity at the diazonium group in the presence of halogen positions and the compatibility of base-sensitive functional groups as esters. For instance the of 4-chloro-and 4bromophenyl diazonium tosylates provided indoles 8h-i with no evidence for a competing electrophilic reactivity of the halide atoms. Interestingly, we did not observe any significant difference in reactivity of 1-methylindole 6c and N-Boc indole 6d which behaved similarly than free NH-indoles (see 8l-n versus 8d, 8k and 8b, respectively). Decorating the indole skeleton with substituents (OMe and Br) significantly extended the reaction scope of the process (8o-q). Despite the lower yields observed with bromo-indoles, no evidence for a competing electrophilic reactivity of the bromine atom was observed and the drop in reactivity was instead attributed to a lower solubility of the starting indole 6e. The inertness of halogen atoms in the reaction conditions might allow further chemical functionalizations through standard transition-metal catalyzed reactions.

Scheme 2. Synthesis of the key intermediate 8r of FPTase inhibitor 4.

To further confirm the efficiency of the methodology developed herein and, if this should prove necessary, highlight the importance of such transformations, we applied our process to the synthesis of Methyl 2-(2-phenyl-1-methyl-1H-indol-3-yl)acetate 8r which was isolated with 81% yield (Scheme 2). A continuous use of the flow platform over a period of 10 hours allowed to isolate ca. 0.5 g of indole 8r (79%), showing the facile scalability of our methodology in flow. Indole 8r is a key intermediate in the synthesis of farnesyl pyrophosphate-competitive inhibitor of farnesyl protein transferase 4. [3]

Conclusion

In we developed a mild, base-free, ligandfree and safe continuous-flow protocol for the palladium-catalyzed direct C-H arylation of indole-3acetic acid derivatives with arenediazonium salts. Highlights of this work include (i) the development of a fully autonomous self-optimizing flow platform, (ii) the use of arenediazonium salts as superelectrophiles, and (iii) the synthesis of a key intermediate in the synthesis of a farnesyl pyrophosphate-competitive inhibitor of farnesyl protein transferase. The capability of our autonomous self-optimizing flow platform to simultaneously screen several continuous variables without a priori knowledge of the optimal objective function is an indisputable attractive feature to optimize chemical transformations in a short time frame and reliable experimental conditions. While arenediazonium salts are often considered as hazardous reagents, our continuous-flow approach mitigates the risk associated with these attractive electrophiles which allow high reaction rates and mild experimental conditions. We believe that this work is a representative contribution in the area of machineassisted chemistry which should become the standard of chemical laboratories 4.0.

Experimental Section

General Information. All commercially available chemicals were used as received unless otherwise noted. High-field 1 H and 13 C NMR spectra were recorded at 300 or 400 and 75 or 100 MHz, respectively. 1 H and 13 C NMR spectra were referenced to the residual signal of the internal deuterated solvent (CDCl 3 ) at 7.26 and 77.16 ppm, respectively, and coupling constants were measured in Hertz. The following abbreviations were used to explain the multiplicities: s = singlet, d = doublet, t = triplet, q = quartet, m = multiplet, br = broad. FT-IR spectra were recorded in the ATR mode. Wavelengths of maximum absorbance (ν max ) are quoted in wave numbers (cm -1 ). High resolution mass spectrometry (HRMS) was recorded on a microTOF spectrometer equipped with orthogonal electrospray interface (ESI). Reactions were monitored by thin layer chromatography (TLC) with pre-coated silica gel plates (60 F 254 ), and visualization was accomplished with either a UV lamp at 254 nm or by immersion in ethanolic solution of phosphomolybdic acid (PMA). Flash column chromatography was performed using silica gel 60 (40-63 μm). In order to remove any trace of waste (Pd particles, inorganic solts…) which progressively deposits on the wall of the reactor coil, the tubbing was washed every 100 hours of use with an aqueous solution of nitric acid (1M, 50 mL) followed by a thorough washing with water (100 mL) and CH 3 CN (100 mL).

Details of the experimental setup of the selfoptimizing reactor. HPLC pumps (JASCO PU2080 and JASCO PU4185) equipped with a RS-232 port were employed to flow the solution through the system. The reactor coil was heated with a heating plate (Heidolph, MR Hei-Connect) equipped with a RS-232 port. A 2-way 6port valve (VICI, Cheminert C2-3006D) equipped with a RS-232 port was used to inject an aliquot of the crude mixture within the on-line HPLC unit. The HPLC column (Agela Promosil C18, 3.5 mm × 150mm, 5 µm) outlet was connected to a UV detector (JASCO, UV 2075) equipped with a RS-232 port. The flow outlet was connected to a programmable fraction collector (Advantec, CHF 1225C). All units equipped with a RS-232 port were autonomously controlled with MATLAB ® through the use of communication protocols provided by the manufacturers.

Experimental procedure of the self-optimization for compound 8a. The experimental setup consisted of two streams as depicted in Figure 2. The first stream contained a solution of tolyldiazonium tosylate 7a (0.1 M) in MeOH while the second contained a solution of indole 6a (0.1 M), Pd(OAc) (5 × 10 -3 M) and naphthalene 9 (0.1 M) as the internal standard in DMF. Each solution was continuously pumped with two independent pumps at the required flow rate, and met in a stainless steel T-shaped piece (internal volume: 0.57 μL). The resulting mixture was introduced in a PEEK reactor coil (10 mL, 0.75 mm i.d.) heated at the desired temperature. The reactor outlet was connected to an automatic 2-way 6-port switch valve which injected 0.2 µL of the crude mixture in the HPLC unit while the remaining stream was collected in a fraction collector. A mixture of MeOH/H 2 O (60/40, v/v) was used as mobile phase for the HPLC analysis at a flow rate of 0.4 mL/min. A UV detector was connected to the outlet of the HPLC column (Agela Promosil C18, 3.5 mm × 150 mm, 5 µm) to follow the absorbance at a wavelength of 254 nm. Peak integration and yield calculation were under full MATLAB automation. The calculated yield was automatically sent to the algorithm which set new experimental conditions to the units via RS-232 ports. A 3-D optimization of the reaction yield was conducted using the temperature, residence time and stoichiometry as the variables. The initial experiment of the simplex was: 25 °C, 15 min and 1 equivalent of diazonium salt 7a. The d values were 7 °C, 15 min and 0.2 equivalent. The lower and upper boundaries of the research space were the following: 25-70 °C, 10-90 min and 1-2 equiv for the temperature, residence time and stoichiometry, respectively. The details of the 28 experiments carried out with the autonomous flow reactor are given in supporting information (Table S2). An optimum giving 90% HPLC yield was found in experiment 15 at 44 °C, 74 min of residence time and 1.60 equiv. of diazonium salt 7a. An analytical sample of ethyl 2-(2-(p-tolyl)-1H-indol-3-yl) acetate 8a was obtained by flash chromatography (10% AcOEt-petroleum ether) as a white solid (51.6 mg, 88%). mp 97 °C. IR (ATR) ν 3346, 1719, 764 cm -1 . 1 Purification by flash chromatography on silica gel (5% EtOAc-petroleum ether) gave 8c as a white solid (44.1 mg, 79%). mp 98 °C [Lit. [40] 102 °C]. IR (ATR) ν 3362, 1719, 1167, 694 cm -1 . 1 H NMR (CDCl 3 , 300 MHz) δ 8.03 (br s, 1H), 7.66 (dd, 1H, J = 0.9, 8.0 Hz), 7.42 -7.28 (m, 5H), 7.24 (dt, 1H, J = 1.4, 7.0 Hz), 7.18 (dt, 1H, J = 1.2, 7.1 Hz), 3.65 (s, 2H), 3.63 (s, 3H), 2.27 (s, 3H). 13 Purification by flash chromatography on silica gel (5% EtOAc-petroleum ether) gave 8d as a white solid (54.6 mg, 89%). mp 124 °C. IR (ATR) ν 3309, 1703, 1343, 738 cm -1 . 1 H NMR (CDCl 3 , 300 MHz) δ 8.10 (br s, 1H), 7.66 (dt, 1H, J = 0.6, 7.7 Hz), 7.58 (dt, 2H, J = 2.0, 8.3 Hz), 7.40-7.34 (m, 3H), 7.22 (dt, 1H, J = 1.4, 7.1 Hz), 7.16 (dt, 1H, J = 1.2, 7.1 Hz), 3.85 (s, 2H), 3.72 (s, 3H), 2.98 (sept, 1H, J = 6.9 Hz), 1.30 (d, 6H, J = 6.9 Hz). 13

  Figure 1. Pharmaceutically relevant 2-aryl-3-carbonylmethyl indoles.

Figure 2 .

 2 Figure 2. Autonomous self-optimizing flow reactor for the optimization of 8a. The red dashed lines indicate the units controlled by the computer via MATLAB.

Figure 3 .

 3 Figure 3. (a) Maximization of the yield of indole 8a. The blue curve corresponds to the MeOH/DMF solvent system while the red curve corresponds to the MeOH/DMF/EtOAc solvent system. (b) Representation of the threedimensional experimental conditions for the maximization

Figure 4 .

 4 Figure 4. Scope of the C2-arylation of indole-3-acetic acid derivatives.

Figure 4 .

 4 The first stream (flow rate: 82 µL/min), equipped with a PEEK injection loop (5 mL) loaded with a solution of arenediazonium tolylate in MeOH (0.1 M), meet in a stainless steel T-shaped piece (internal volume 0.57 μL) a second stream (flow rate: 51 µL/min) consisting of a solution of indoles (0.1 M) and Pd(OAc) 2 (5 × 10 -3 M) in DMF/EtOAc (1/4) loaded in a second PEEK loop (2 mL). The merged streams entered in a PEEK reactor coil (10 or 15 mL, 0.75 mm id, 44 °C) at a flow rate of 0.133 mL/min and the resulting indoles 8b-r were collected in vials. Volatiles were removed under reduced pressure and the crude indoles 8b-r were purified by flash chromatography on silica gel.Methyl 2-(2-(p-tolyl)-1H-indol-3-yl)acetate 8b.Purification by flash chromatography on silica gel (5% EtOAc-petroleum ether) gave 8b as a colorless oil (48 mg, 86%). IR (ATR) ν 3373, 1723, 1456, 740 cm -1 . 1 H NMR (CDCl 3 , 400 MHz) δ 8.14 (br s, 1H), 7.66 (d, 1H, J = 7.7 Hz), 7.54 (d, 2H, J = 8.1 Hz), 7.36 (d, 1H, J = 8.0 Hz), 7.30 (d, 2H, J = 8.2 Hz), 7.22 (dt, 1H, J = 1.1, 7.1 Hz), 7.17 (dt, 1H, J = 0.9, 7.0 Hz), 3.85 (s, 2H), 3.72 (s, 3H), 2.43 (s, 3H).13 C NMR (CDCl 3 , 100 MHz) δ 172.7, 138.0, 136.3, 135.6, 129.6, 129.5, 129.1, 128.1, 122.4, 120.0, 119.1, 110.8, 105.2, 52.0, 30.9, 21.3. HRMS (ESI) m/z [M + H] + Calcd for C 18 H 18 O 2 N 280.1338; Found 280.1336.

  H NMR (CDCl 3 , 300 MHz) δ 8.10 (brs, 1H), 7.67 (dd, 1H, J = 0.7, 7.6 Hz), 7.56 (d, 2H, J = 8.1 Hz), 7.37 (dm, 1H, J = 6.4 Hz), 7.32 (d, 2H, J = 7.8 Hz), 7.22 (dt, 1H, J = 1.3, 7.1 Hz), 7.15 (dt, 1H, J = 1.2, 7.1 Hz), 4.18 (q, 2H, J = 7.2 Hz), 3.82 (s, 2H), 2.42 (s, 3H), 1.26 (t, 3H, J = 7.1 Hz).13 C NMR (CDCl 3 , 75 MHz,) δ 172.3, 137.9, 136.3, 135.6, 130.0, 129.6, 129.5, 129.1, 128.1, 122.3, 119.9, 119.2, 115.0, 110.8, 105.3, 60.8, 31.2, 21.3, 14.2. HRMS (ESI) m/z [M + H] + Calcd for C 19 H 20 O 2 N 294.1494; Found 294.1493.

	General experimental setup for the synthesis of
	indoles 8b-r. The experimental setup consisted of two
	streams as depicted in

  C NMR (CDCl 3 , 75 MHz) δ 172.4, 137.8, 136.0, 135.5, 131.6, 131.1, 130.3, 128.8, 128.1, 125.7, 122.2, 119.8, 119.0, 110.8, 106.3, 51.8, 30.6, 19.9. HRMS (ESI) m/z [M + H] + Calcd for C 18 H 18 O 2 N 280.1338; Found 280.1341.

	Methyl	2-(2-(4-isopropylphenyl)-1H-indol-3-
	yl)acetate 8d.

  C NMR (CDCl 3 , 75 MHz) δ 172.8, 148.9, 136.3, 135.6, 129.8, 129.0, 128.2, 127.0, 122.4, 120.0, 119.1, 110.8, 105.1, 52.0, 33.9, 31.0, 23.9. HRMS (ESI) m/z [M + H] + Calcd for C 21 H 24 O 2 N 308.1651; Found 308.1663. Purification by flash chromatography on silica gel (5% EtOAc-petroleum ether) gave 8e as a white solid (54.6 mg, 85%). mp 138 °C. IR (ATR) ν 3316,1704,

Methyl 2-(2-(4-(tert-butyl)phenyl)-1H-indol-3yl)acetate 8e.

Acknowledgements

We gratefully acknowledge the University of Nantes, the ''Centre National de la Recherche Scientifique'' (CNRS) and the "Région des Pays de la Loire" in the framework of the SmartCat project for financial support. N. Vasudevan thanks the UBL for a grant. We acknowledge Julie Hémez (CEISAM, University of Nantes) for HRMS analyses and Paul De Bonfils (CEISAM, University of Nantes) for technical assistance with FTIR.