Hao Chen
email: hao.chen@inria.fr

Benoit Lagadec
email: benoit.lagadec@esifrance.net

Francois Bremond
email: francois.bremond@inria.fr

Enhancing Diversity in Teacher-Student Networks via Asymmetric branches for Unsupervised Person Re-identification

The objective of unsupervised person re-identification (Re-ID) is to learn discriminative features without laborintensive identity annotations. State-of-the-art unsupervised Re-ID methods assign pseudo labels to unlabeled images in the target domain and learn from these noisy pseudo labels. Recently introduced Mean Teacher Model is a promising way to mitigate the label noise. However, during the training, self-ensembled teacher-student networks quickly converge to a consensus which leads to a local minimum. We explore the possibility of using an asymmetric structure inside neural network to address this problem. First, asymmetric branches are proposed to extract features in different manners, which enhances the feature diversity in appearance signatures. Then, our proposed cross-branch supervision allows one branch to get supervision from the other branch, which transfers distinct knowledge and enhances the weight diversity between teacher and student networks. Extensive experiments show that our proposed method can significantly surpass the performance of previous work on both unsupervised domain adaptation and fully unsupervised Re-ID tasks.

Introduction

Person re-identification (Re-ID) targets at retrieving a person of interest across non-overlapping cameras. Since there are domain gaps resulting from illumination condition, camera property and view-point variation, a Re-ID model trained on a source domain usually shows a huge performance drop on other domains.

Unsupervised domain adaptation (UDA) targets at shifting the model trained from a source domain with identity annotation to a target domain via learning from unlabeled target images. In the real world, unlabeled images in a target domain can be easily recorded, which is almost labor-free. It is intuitive to use these images to adapt a pretrained Re-ID model to the desired domain. Fully unsupervised Re-ID further minimises the supervision by removing pre-training on the labelled source domain.

State-of-the-art UDA Person Re-ID methods [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF][START_REF] Yang | Asymmetric co-teaching for unsupervised cross domain person re-identification[END_REF] and unsupervised methods [START_REF] Lin | Unsupervised person re-identification via softened similarity learning[END_REF] assign pseudo labels to unlabeled target images. The generated pseudo labels are generally very noisy. The noise is mainly from several inevitable factors, such as the strong domain gaps and the imperfection of clustering. In this way, an unsupervised Re-ID problem is naturally transferred into Generating pseudo labels and Learning from noisy labels problems.

To generate pseudo labels, the most intuitive way is to use a clustering algorithm, which gives a good starting point for clustering based UDA Re-ID [START_REF] Zhang | Self-training with progressive augmentation for unsupervised cross-domain person re-identification[END_REF][START_REF] Fu | Self-similarity grouping: A simple unsupervised cross domain adaptation approach for person re-identification[END_REF]. Recently, Ge et al. [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] propose to add a Mean Teacher [START_REF] Tarvainen | Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results[END_REF] model as online soft pseudo label generator, which effectively reduces the error amplification during the training with noisy labels. In this paper, we also use both clustering-based hard labels and teacher-based soft labels in our baseline.

To handle noisy labels, one of the most popular approaches is to train paired networks so that each network helps to correct its peer, e.g., two-student networks in Coteaching [START_REF] Han | Co-teaching: Robust training of deep neural networks with extremely noisy labels[END_REF] and two-teacher-two-student networks in MMT [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF]. However, these paired models with identical structure are prone to converge to each other and get stuck in a local minimum. There are several attempts to alleviate this problem, such as Co-teaching+ [START_REF] Yu | How does disagreement help generalization against label corruption?[END_REF], ACT [START_REF] Yang | Asymmetric co-teaching for unsupervised cross domain person re-identification[END_REF] and MMT [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF]. These attempts of keeping divergence between paired models are mainly based on either different training sample selection [START_REF] Yu | How does disagreement help generalization against label corruption?[END_REF][START_REF] Yang | Asymmetric co-teaching for unsupervised cross domain person re-identification[END_REF] or different initialization and data augmentation [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF]. In this paper, we propose a strong alternative by designing asymmetric neural network structure in the Mean Teacher Model. We use two independent branches with different depth and global pooling methods as last layers of a neural network. Features extracted from both branches are concatenated as the appearance signature, which enhances the feature diversity in the appearance signature and allows to get better clustering-based hard labels. Each branch gets supervision from its peer branch of different structure, which enhances the divergence between paired teacher-student networks. Our proposed decoupling method does not rely on different source domain initializa-tions, which makes it more effective in the fully unsupervised scenario where the source domain is not available.

In summary, our contributions are:

1. We propose to enhance the feature diversity inside person Re-ID appearance signatures by splitting last layers of a backbone network into two asymmetric branches, which increases the quality of clusteringbased hard labels.

2. We propose a novel decoupling method where asymmetric branches get cross-branch supervision, which avoids weights in paired teacher-student networks converging to each other and increases the quality of teacher-based soft labels.

3.

Extensive experiments and ablation study are conducted to validate the effectiveness of each proposed component and the whole framework.

Related Work

Unsupervised domain adaptive Re-ID. Recent unsupervised cross-domain Re-ID methods can be roughly categorized into distribution alignment and pseudo label based adaptation. The objective of distribution alignment is to learn domain invariant features. Several attempts [START_REF] Wang | Transferable joint attribute-identity deep learning for unsupervised person re-identification[END_REF][START_REF] Lin | Multi-task mid-level feature alignment network for unsupervised cross-dataset person re-identification[END_REF] leverage semantic attributes to align the feature distribution in the latent space. However, these approaches strongly rely on extra attribute annotation, which requires extra labor. Another possibility is to align the feature distribution by transferring labeled source domain images into the style of target domain with generative adversarial networks [START_REF] Wei | Person transfer gan to bridge domain gap for person reidentification[END_REF][START_REF] Zhong | Generalizing a person retrieval model hetero-and homogeneously[END_REF][START_REF] Chen | Instanceguided context rendering for cross-domain person reidentification[END_REF]. Style transferred images are usually combined with pseudo label based adaptation to get a better performance. Pseudo label based adaptation is a more straightforward approach for unsupervised cross-domain Re-ID, which directly assigns pseudo labels to unlabelled target images and allows to fine-tune a pre-trained model in a supervised manner. Clustering algorithms are widely used in previous unsupervised cross-domain Re-ID methods. UDAP [START_REF] Song | Unsupervised domain adaptive re-identification: Theory and practice[END_REF] provides a good analysis on clustering based adaptation and use a k-reciprocal encoding [START_REF] Zhong | Reranking person re-identification with k-reciprocal encoding[END_REF] to improve the quality of clusters. PCB-PAST [START_REF] Zhang | Self-training with progressive augmentation for unsupervised cross-domain person re-identification[END_REF] simultaneously learns from a ranking-based and clustering-based triplet losses. SSG [START_REF] Fu | Self-similarity grouping: A simple unsupervised cross domain adaptation approach for person re-identification[END_REF] assigns clustering-based pseudo labels to both global and local features. To mitigate the clustering-based label noise, researchers borrow ideas from how unlabeled data is used in Semi-supervised learning and Learning from noisy labels. ECN [START_REF] Zhong | Invariance matters: Exemplar memory for domain adaptive person re-identification[END_REF] uses an exemplar memory to save averaged features to assign soft labels. ACT [START_REF] Yang | Asymmetric co-teaching for unsupervised cross domain person re-identification[END_REF] splits the training data into inliers/outliers to enhance the divergence of paired networks in Co-teaching [START_REF] Han | Co-teaching: Robust training of deep neural networks with extremely noisy labels[END_REF]. MMT [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] adopts two student and two Mean Teacher networks. Two students are initialized differently from source pre-training in order to enhance the divergence of paired teacher-student networks. Each mean teacher network provides soft labels to supervise peer student network. However, despite different initializations and different data augmentations used in peer networks, the decoupling is not encouraged enough during the training. We directly use asymmetric neural network structure inside teacher-student networks, which encourages the decoupling at all epochs.

Teacher-Student Network for Semi-Supervised Learning. Unsupervised domain adaptation can be regarded to some extent as Semi-Supervised Learning (SSL), since both of them utilize labeled data (source domain for UDA) and large amount of unlabeled data (target domain for UDA).

A teacher-student structure is commonly used in SSL. This structure allows student network to gradually exploit data with perturbations under consistency constraints. In Π model and Temporal ensembling [START_REF] Laine | Temporal ensembling for semisupervised learning[END_REF], the student learns from either samples forwarded twice with different noise or exponential moving averaged (EMA) predictions under consistency constraints. Instead of EMA predictions, Mean-teacher model [START_REF] Tarvainen | Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results[END_REF] uses directly the EMA weights from the student to supervise the student under a consistency constraint. Authors of Dual student [START_REF] Ke | Dual student: Breaking the limits of the teacher in semi-supervised learning[END_REF] point out that the Mean Teacher converging to student along with training (coupling problem) prevents the teacher-student from exploiting more meaningful information from data. Inspired by Deep Co-training [START_REF] Qiao | Deep co-training for semi-supervised image recognition[END_REF], they propose to train two independent students on stable samples which have same predictions and enough large feature difference. However, in unsupervised cross-domain Re-ID, labeled source domain and unlabeled target domain do not share the same identity classes, which makes traditional close-set SSL methods hard to use.

Fully unsupervised Re-ID. Recently, several fully unsupervised Re-ID methods are proposed to further minimize the supervision, which does not require any Re-ID annotation. A bottom-up clustering framework is proposed in BUC [START_REF] Lin | A bottom-up clustering approach to unsupervised person re-identification[END_REF], which trains a network based on the clusteringbased pseudo labels in an iterative way. [START_REF] Lin | Unsupervised person re-identification via softened similarity learning[END_REF] replaces clustering-based pseudo labels with similarity-based softened labels. Different to image-based unsupervised Re-ID, [START_REF] Wu | Tracklet self-supervised learning for unsupervised person reidentification[END_REF] learns tacklet information with clustering-based pseudo labels. MMT [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] can be transferred into an unsupervised method by removing the pre-training in source domain. However, without different source domain initializations, divergence between peer networks can not be enough encouraged in MMT. Instead of different source domain initializations, divergence is encouraged by asymmetric network structures, which is more suitable for fully unsupervised Re-ID.

Proposed Method

Overview

Given two datasets: one labeled source dataset D s and one unlabeled target dataset D t , the objective of UDA is to adapt a source pretrained model M pre to the target dataset with unlabeled target data. To achieve this goal, we propose a two-staged adaptation approach based on Mean Teacher Model. We focus on the coupling problem (teacher and student converge to each other) existing inside the original Mean Teacher. Asymmetric branches and cross-branch supervision are proposed in this paper to address this problem and to enhance the diversity in the network, which show great effectiveness for UDA Re-ID.

Asymmetric branches

A multi-branch structure is widely used in the fully supervised Re-ID methods, especially in global-local feature based methods [START_REF] Fu | Horizontal pyramid matching for person reidentification[END_REF][START_REF] Dai | Batch dropblock network for person re-identification and beyond[END_REF][START_REF] Hao Chen | Learning discriminative and generalizable representations by spatial-channel partition for person re-identification[END_REF]. Such structure keeps independence between branches, which makes features extracted from different branches diversified. In the unsupervised Re-ID, we conduct clustering on appearance signatures computed from person images to generate pseudo labels. The quality of appearance signatures can be improved by extracting distinct meaningful features from different branches. Thus, we duplicate last layers of a backbone network and make them different in the structure, which we call Asymmetric Branches.

Asymmetric branches are illustrated in Figure 1. For a ResNet-based [START_REF] He | Deep residual learning for image recognition[END_REF] backbone, the layer 4 is duplicated. The first branch is kept unchanged as the one used in the original backbone: 3 bottlenecks and global average pooling (GAP). The second branch is composed of 4 bottlenecks and global max pooling (GMP). The GAP perceives global information, while the GMP focuses on the most discriminative in-formation (most distinguishable identity information, such as a red bag or a yellow t-shirt). Asymmetric branches improve appearance signature quality by enhancing the feature diversity, which is validated by source pre-training performance boost in Table 3 as well as examples in Figure 5. They further improve the quality of pseudo labels during the adaptation, which is validated by target adaptation performance in Table 3.

Asymmetric Branched Mean Teaching

We call our proposed adaptation method Asymmetric Branched Mean Teaching (ABMT). Our proposed ABMT contains two stages: Source pre-training and Target adaptation.

Source domain supervised pre-training

In the first stage, we train a network in the fully supervised way on the source domain. Thanks to this stage, the model used for adaptation obtains a basic Re-ID capacity, which helps to alleviate pseudo label noise. Given a source sample x s i and its ground truth identity y i , the network (with weight θ) encodes x s i into average F a (x s i |θ) and max features F m (x s i |θ) and then gets two predictions P a (x s i |θ) and P m (x s i |θ). Cross-entropy L ce and batch hard triplet [START_REF] Hermans | In defense of the triplet loss for person re-identification[END_REF] L tri losses are used in this stage as shown in Figure 1.

The whole network is trained with a combination of both losses:

L scr =λ s ce L ce (P a (x s i |θ), y i ) + λ s ce L ce (P m (x s i |θ), y i )+ λ s tri L tri (F a (x s i |θ), y i ) + λ s tri L tri (F m (x s i |θ), y i ) (1)

Target domain unsupervised adaptation

The adaptation procedure is illustrated in Figure 2. It contains two components: Clustering-based hard label generation and Cross-branch teacher-based soft label training. After adaptation, only teacher network is used during the inference.

Clustering-based hard label generation. In previous UDA Re-ID methods, distance-based K-Means [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] and density-based clustering DBSCAN [START_REF] Yang | Asymmetric co-teaching for unsupervised cross domain person re-identification[END_REF][START_REF] Song | Unsupervised domain adaptive re-identification: Theory and practice[END_REF] are main approaches to generate pseudo labels. We follow the state-of-the-art DBSCAN based clustering method presented in [START_REF] Song | Unsupervised domain adaptive re-identification: Theory and practice[END_REF]. To adapt it to our proposed asymmetric branches, we concatenate the average and max features from asymmetric branches in the teacher network as appearance signatures. Images belonging to the same identity should have the same nearest neighbors in the feature space. Distance metric for DBSCAN are obtained by kreciprocal re-ranking encoding [START_REF] Zhong | Reranking person re-identification with k-reciprocal encoding[END_REF] between target domain and source domain samples. The density-based clustering generates unfixed cluster numbers at different epochs, which means old classifiers from the last epoch can not be reused after a new clustering. Thus, we simply create new classifiers depending on the number of clusters at the beginning of each epoch. We take normalized mean features of each cluster from the average branch to initialize the average branch classifiers and similarly normalized max features from max branch to initialize the max branch classifiers. We call these classifiers with flexible dimension "Dynamic Classifiers". With the help of these Dynamic Classifiers, the student is trained on cluster components (outliers are discarded) with cross-entropy loss:

L ce = - i (y i log(P m (x t i |θ))) - i (y i log(P a (x t i |θ))) (2
) where y i is the clustering based hard label and P a (x t i |θ) and P m (x t i |θ) are student predictions from both asymmetric branches.

Cross-branch teacher-based soft label training. Clustering algorithms generate hard pseudo labels whose confidences are 100%. Since Re-ID is a fine-grained recognition problem, people with similar clothes are not rare in the dataset. Hard pseudo labels of these similar samples can be extremely noisy. In this case, soft pseudo labels (confidences < 100%) are more reliable. Learning with both hard and soft pseudo labels can effectively alleviate label noise.

The Mean Teacher Model [START_REF] Tarvainen | Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results[END_REF] (teacher weights θ ) uses the EMA weights of the student model (student weights θ). The Mean Teacher Model shows strong capacity to handle label noise and avoids error amplification along with training. We define θ t at training step t as the EMA of successive weights:

θ t = θ t , if t = 0 αθ t-1 + (1 -α)θ t , otherwise (3)
where α is a smoothing coefficient that controls the selfensembling speed of the Mean Teacher. Despite these advantages of Mean Teacher, such selfensembling teacher-student networks (the teacher is formed by EMA weights of the student, and the student is supervised by the teacher) face the coupling problem. We use the Mean Teacher soft label generator as in [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] and address the coupling problem by cross-branch supervision. Each branch in the student is supervised by a teacher branch which has different structure. Weight diversity between the paired teacher-student can be better kept. Given one target domain sample x t i , the teacher (teacher weights θ ) encodes it into two feature vectors from two asymmetric branches, average features F a (x t i |θ ) and max features F m (x t i |θ ). The dynamic classifiers then transform these two feature vectors into two predictions respectively P a (x t i |θ ) and P m (x t i |θ ). Similarly, features of the student (student weights θ) are F a (x t i |θ) and F m (x t i |θ), while predictions are P a (x t i |θ) and P m (x t i |θ). The predictions from the teacher supervise those from the student with a soft cross-entropy loss [START_REF] Hinton | Distilling the knowledge in a neural network[END_REF] in a cross-branch manner, which can be formulated as

L a→m sce = - i (P a (x t i |θ ) log(P m (x t i |θ))) (4)
L m→a sce = - i (P m (x t i |θ ) log(P a (x t i |θ))) (5)
To further enhance the teacher-student networks' discriminative capacity, the features in the teacher supervise those of the student with a soft triplet loss [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF]:

L a→m stri = - i (T a (x t i |θ ) log(T m (x t i |θ))) (6)
L m→a stri = - i (T m (x t i |θ ) log(T a (x t i |θ))) (7)
where

T (x t i |θ) = exp( F (x t i |θ)-F (x t p |θ) 2 ) exp( F (x t i |θ)-F (x t p |θ) 2 )+exp( F (x t i |θ)-F (x t n |θ) 2 )
is the softmax triplet distance of the sample x t i , its hardest positive x t p and its hardest negative x t n in a mini-batch. By minimizing the soft triplet loss, the softmax triplet distance in a mini-batch from the student is encouraged to get as close as possible to the distance from the teacher. The positive and negative samples within a mini-batch are decided by clustering-based hard pseudo labels. It can effectively improve the UDA Re-ID performance. The teacher-student networks are trained end-to-end with Equation ( 2), ( 4), ( 5), ( 6), [START_REF] Fu | Horizontal pyramid matching for person reidentification[END_REF].

L target =λ t ce L ce + λ t sce (L a→m sce + L m→a sce ) + λ t stri (L a→m stri + L m→a stri ) (8)

Coupling Problem in Mean Teacher Based Methods

The Mean Teacher Baseline is illustrated in Figure 3 (a) where the student gets supervision from its own EMA weights. In the Mean Teacher Baseline, the student and the teacher quickly converge to each other (coupling problem), which prevents them from exploring more diversified information. Authors of MMT [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] propose to pre-train 2 student networks with different seeds. As illustrated in Figure 3 (b), two Mean Teacher networks are formed separately from two students, which alleviates the coupling problem. However, different initializations decouple both teacher peers only at first epochs. Without a diversity encouragement during the adaptation, both teachers still converge to each other along with training. In Figure 3 (c), our proposed asymmetric branches provide a diversity encouragement during the adaptation, which decouples both teacher peers at all epochs.

To validate our idea, we propose to measure Euclidean distance of appearance signature features between two teacher networks or two teacher branches. We extract feature vectors after global pooling on all images in the target training set. Then, we calculate the Euclidean distance between feature vectors of both teachers and sum up the distance of every image as the final feature distance. If the feature distance is large, we can say that both teacher peers extract diversified features. Otherwise, the teacher peers converge to each other. As we can see from the left curves in Figure 4, the feature distance between two teachers in MMT is large at the beginning, but it decreases and then stabilizes. Differently, the feature distance between two branches in our proposed method remains large during the training. Moreover, we visualize the Euclidean distance of appearance signature features on all target training samples between teacher and student networks in Figure 4 right curves. Our method can maintain a larger distance, which shows that it can better decouple teacher-student networks.

Experiments

Datasets and Evaluation Protocols

Our proposed adaptation method is evaluated on 3 Re-ID datasets: Market-1501, DukeMTMC-reID and MSMT17. Market-1501 [START_REF] Zheng | Scalable person re-identification: A benchmark[END_REF]

Implementation details

Hyper-parameters used in our proposed method are searched empirically from the Market→ Duke task and kept the same for the other tasks. To conduct fair comparison with state-of-the-arts, we use a ImageNet [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF] pre-trained ResNet-50 [START_REF] He | Deep residual learning for image recognition[END_REF] as our backbone network. The backbone can be extended to ResNet-based networks designed for cross domain tasks, e.g., IBN-ResNet-50 [START_REF] Pan | Two at once: Enhancing learning and generalization capacities via ibn-net[END_REF]. An Adam optimizer with a weight decay rate of 0.0005 is used to optimize our networks. Our networks are trained on 4 Nvidia 1080Ti GPUs under Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] framework. Detailed configurations are given in the following paragraphs. Stage2: Target domain unsupervised adaptation. For the clustering, we set the minimum cluster samples to 4 and the density radius r=0.002. Re-ranking parameters for calculating distances are kept the same as in [START_REF] Song | Unsupervised domain adaptive re-identification: Theory and practice[END_REF] for UDA setting. Re-ranking between source and target domain is not considered for fully unsupervised setting. The Mean Teacher network is initialized and updated in the way of Equation 3 with a smoothing coefficient α = 0.999. We set λ t ce = 0.5, λ t sce = 0.5 and λ t stri = 1 in Equation 8. The adaptation epoch E ada is set to 40. For each epoch, the networks are trained R ada = 400 iterations with a fixed learning rate 0.00035. For each iteration, 64 images of 16 clustering-based pseudo identities are resized to 256*128 and fed into networks with Random erasing [START_REF] Zhong | Random erasing data augmentation[END_REF] data augmentation.

Comparison with State-of-the-Art Methods

We compare our proposed methods with state-of-theart UDA methods in Table 1 for 4 cross-dataset Re-ID tasks: Market→ Duke, Duke → Market, Market → MSMT and Duke → MSMT. Post-processing techniques (e.g., Reranking [START_REF] Zhong | Reranking person re-identification with k-reciprocal encoding[END_REF]) are not used in the comparison. Our proposed method outperforms MMT [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] (cluster number is set to 500, 700 and 1500 respectively). We can also adjust the density radius in DBSCAN depending on target domain size to get a better performance, but we think it is hard to know the target domain size in the real world. With an IBN-ResNet50 [START_REF] Pan | Two at once: Enhancing learning and generalization capacities via ibn-net[END_REF] backbone, the performance on 4 tasks can be further improved. Examples of retrieved images are illustrated in Figure 5. Compared to MMT, embeddings from our proposed method contain more discriminative appearance information (e.g., shoulder bag in the first row), which are ro-bust to noisy information (e.g., pose variation in the second row, occlusion in the third row and background variation in the fourth row). This qualitative comparison confirms that appearance signatures of our proposed method are of improved quality.

We compare unsupervised Re-ID methods in Table 2. Since the Mean Teacher is designed for handling label noise, it is interesting to see the performance without source pre-training, which introduces more label noise during the adaptation. This setting corresponds to an unsupervised Re-ID. We use ImageNet pretained weights as initialization. Our proposed method outperforms previous unsupervised Re-ID by a large margin, which shows that ImageNet initialization can provide basic discriminative capacity for Re-ID.

MMT [START_REF] Ge | Mutual meanteaching: Pseudo label refinery for unsupervised domain adaptation on person re-identification[END_REF] is the first Mean Teacher based UDA Re-ID method. Authors of MMT propose to use 2 students and 2 teachers with different initialization and stochastic data augmentation to address the coupling problem. We also use Mean Teacher soft pseudo labels but propose a different decoupling solution. Features in asymmetric branches are always extracted in different manners during the adaptation. Compared to MMT, our proposed method has less parameters (approximately 10% less parameters and 20% less operations) but achieves better performance. Moreover, in the unsupervised scenario, we can not pre-train MMT with different seeds to obtain different Re-ID initializations. This decoupling strategy becomes inappropriate. Our decoupling strategy relies on structural asymmetry instead of different initializations, which is much more effective in the unsupervised scenario.

ACT [START_REF] Yang | Asymmetric co-teaching for unsupervised cross domain person re-identification[END_REF] uses 2 networks, in which each network learns from its peer. Input data are split into inliers and ouliers after DBSCAN. Then, the first network selects small entropy inliers to train the second network, while the second selects small entropy outliers to train the first. This method enhances input asymmetry by data split. Differently, our proposed method focuses on neural network structure asymmetry.

Ablation Studies

Effectiveness of each component in ABMT. Compared with traditional clustering-based Re-ID methods, the performance improvement mainly comes from DBSCAN on reranked distance, asymmetric branches and cross-branch supervision. We use a Mean Teacher Baseline where original ResNet-50 and a K-Means++ clustering of 500 clusters are adopted. We conduct ablation studies by gradually adding one component at each time. Results are shown in Table 3. We can observe: (1) Our proposed asymmetric branches bring the most significant performance improvement during the adaptation. Moreover, as we can see from first two rows in Table 3, they can directly improve the domain generalizability of appearance signatures without target adaptation. (2) DBSCAN on re-ranked distance works better than a K-Means++ clustering of 500 clusters during the adaptation. (3) Cross-branch supervision works on asymmetric branches, which can further improve the adaptation performance. Effectiveness of asymmetric branch structure. To validate the effectiveness of our proposed asymmetric branch structure, we compare several possible structures: (1) 2 branches with different pooling methods and different depths, (2) 2 branches with same pooling methods (global average pooling) but different depths, (3) 2 branches with different methods but same depths, (4) 3 branches where the new branch is composed of 5 bottleneck blocks and global average pooling. From results given in Table 4, we can conclude that different pooling methods play a more important role in asymmetric branches.

Effectiveness of loss functions. We conduct ablation studies on loss functions used in our proposed method and report results in Table 5. The degree of influence of 3 loss functions used in our proposed method: L ce > L sce > L stri .

Can traditional decoupling methods further improve the performance? Stochastic data augmentation (teacher inputs and student inputs are under stochastic data augmentation methods) and drop out (teacher feature vectors and student feature vectors are under independent drop out operations before classifiers) are 2 widely-used methods to provide random noise, which also helps to decouple the weights between the teacher and the student. We conduct experiments with stochastic data augmentation . The results in Table 3 show that they can not further improve the UDA Re-ID performance. These methods are not designed for fine-grained Re-ID task. As UDA Re-ID performance is already very high, they can not contribute anymore.

Conclusion

In this paper, we propose a novel unsupervised crossdomain Re-ID framework. Our proposed method is mainly based on learning from noisy pseudo labels generated by clustering and Mean Teacher. A self-ensembled Mean Teacher is robust to label noise, but the coupling problem inside paired teacher-student networks leads to a performance bottleneck. To address this problem, we propose asymmetric branches and cross-branch supervision, which can effectively enhance the diversity in two aspects: appearance signature features and teacher-student weights. By enhancing the diversity in the teacher-student networks, our proposed method achieves better performance on both unsupervised domain adaptation and fully unsupervised Re-ID tasks.

Figure 1 .

 1 Figure 1. Source domain pre-training for asymmetric branched network. One ResNet bottleneck block corresponds to three convolutional layers. For UDA setting, inputs are labelled images from source training set. GAP refers Global Average Pooling, while GMP refers to Global Max Pooling. FC refers to Fully Connected layer.

Figure 2 .

 2 Figure 2. ABMT adaptation. For UDA setting, inputs are training set images from both source and target domains. For fully unsupervised setting, inputs are unlabeled images from target training set.

Stage1:

  Source domain supervised pre-training. We set λ s ce = 0.5 and λ s tri = 0.5 in Equation 1. The max epoch E pre is set to 80. For each epoch, the networks are trained R pre = 200 iterations.The initial learning rate is set to 0.00035 and is multiplied by 0.1 at the 40th and 70th epoch. For each iteration, 64 images of 16 identities are resized to 256*128 and fed into networks.

Figure 3 .

 3 Figure 3. Comparison between (a) Mean Teacher Baseline (b) Mutual Mean Teaching [8] and (c) our Mean Teacher with cross-branch supervised asymmetric branches. Teacher network is formed by exponential moving average (EMA) values of student network.

Figure 4 .

 4 Figure 4. Distance comparison between features extracted from a ResNet50 backbone on all samples in DukeMTMC-reid training set for Market → Duke task. Left: Feature distance between two teacher models in MMT and between two teacher branches in our proposed method. Right: Feature distance between teacher and student networks.

Table 2 .

 2 Comparison of unsupervised Re-ID methods (%) with a ResNet50 backbone on Market and Duke datasets. * refers to our implementation where we remove the source pre-training step. DBSCAN refers to a DBSCAN clustering based on re-ranked distance.

Figure 5 .

 5 Figure 5. Examples of retrieved most similar 5 images in Market → Duke task from MMT [8] and our proposed method. Given a query image, different identity images are highlighted by red bounding boxes, while same identity images are highlighted by green bounding boxes.

  dataset is collected in front of a supermarket in Tsinghua University from 6 cameras. It contains 12,936 images of 751 identities in the training set and 19,732 images of 750 identities in the testing set. DukeMTMC-reID [21] is a subset of the DukeMTMC dataset. It contains 16,522 images of 702 persons in the training set, 2,228 query images and 17,661 gallery images of 702 persons for testing from 8 cameras. MSMT17 [25] is a large-scale Re-ID dataset, which contains 32,621 training images of 1,041 identities and 93,820 testing images of 3,060 identities collected from 15 cameras. Both Cumulative Matching Characteristics (CMC) and mean Average Precisions (mAP) are used in our experiments.

Table 1 .

 1 Comparison of unsupervised domain adaptation (UDA) Re-ID methods (%) on medium-to-medium datasets (Market→ Duke and Duke → Market) and medium-to-large datasets (Market → MSMT and Duke → MSMT).

	UDA Methods			Market → Duke Duke → Market Market → MSMT Duke → MSMT mAP Rank1 mAP Rank1 mAP Rank1 mAP Rank1
	HHL (ECCV'18)[33]			27.2	46.9	31.4	62.2	-	-	-	-
	ECN (CVPR'19)[34]			40.4	63.3	43.0	75.1	8.5	25.3	10.2	30.2
	PCB-PAST (ICCV'19)[29]		54.3	72.4	54.6	78.4	-	-	-	-
	SSG (ICCV'19)[6]			53.4	73.0	58.3	80.0	13.2	31.6	13.3	32.2
	UDAP (PR'20)[22]			49.0	68.4	53.7	75.8	-	-	-	-
	ACT (AAAI'20)[27]			54.5	72.4	60.6	80.5	-	-	-	-
	ECN+ (PAMI'20) [35]			54.4	74.0	63.8	84.1	15.2	40.4	16.0	42.5
	MMT500 (ICLR'20)(ResNet50)[8]	63.1	76.8	71.2	87.7	16.6	37.5	17.9	41.3
	MMT700 (ICLR'20)(ResNet50)[8]	65.1	78.0	69.0	86.8	-	-	-	-
	MMT1500 (ICLR'20)(ResNet50)[8]		-	-	-	-	22.9	49.2	23.3	50.1
	ours (ResNet50)			69.1	82.0	78.3	92.5	23.2	49.2	26.5	54.3
	MMT500 (ICLR'20)(IBN-ResNet50)[8]	65.7	79.3	76.5	90.9	19.6	43.3	23.3	50.0
	MMT700 (ICLR'20)(IBN-ResNet50)[8]	68.7	81.8	74.5	91.1	-	-	-	-
	MMT1500 (ICLR'20)(IBN-ResNet50)[8]	-	-	-	-	26.6	54.4	29.3	58.2
	ours (IBN-ResNet50)			70.8	83.3	80.4	93.0	27.8	55.5	33.0	61.8
	Unsupervised methods	Market mAP Rank1 mAP Rank1 Duke
	MMT500*(ICLR'20)[8]	26.9	48.0	7.3	12.7
	BUC (AAAI'19)[16]	30.6	61.0	21.9	40.2
	SoftSim (CVPR'20)[17]	37.8	71.7	28.6	52.5
	TSSL (AAAI'20)[26]	43.3	71.2	38.5	62.2
	MMT*+DBSCAN (ICLR'20)[8] 53.5	73.1	54.5	69.5
	ours w/o Source pre-training	65.1	82.6	63.1	77.7

Table 3 .

 3 Ablation studies with ResNet50 backbone. MT-Baseline corresponds to the Mean Teacher Baseline in Figure3(a) with a ResNet-50. K-Means refers to a K-Means++ clustering whose cluster number is set to 500. AB refers to asymmetric branches. DBSCAN refers to a DBSCAN clustering[START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF].

	Source pre-training			Market → Duke Duke → Market mAP Rank1 mAP Rank1
	ResNet50				29.6	46.0	31.8	61.9
	ResNet50+AB				31.5	49.7	33.2	63.2
	Target adaptation			Market → Duke Duke → Market mAP Rank1 mAP Rank1
	MT-Baseline+K-Means			59.9	74.8	68.9	88.2
	MT-Baseline+DBSCAN			61.9	77.3	69.9	88.3
	MT-Baseline+K-Means+AB		64.7	78.1	74.8	90.5
	MT-Baseline+K-Means+AB+Cross-branch	66.4	79.9	76.8	91.7
	MT-Baseline+DBSCAN+AB		67.8	81.1	77.3	92.0
	ABMT(MT-Baseline+DBSCAN+AB+Cross-branch) 69.1	82.0	78.3	92.5
	ABMT+Stochastic data augmentation	68.8	81.2	77.6	91.7
	ABMT+Drop out			68.3	81.8	77.9	92.0
	Structure	Market → Duke Duke → Market mAP Rank1 mAP Rank1
	ABMT	69.1	82.0	78.3	92.5
	ABMT w/o different pooling 65.2	79.7	74.2	90.1
	ABMT w/o extra bottleneck	67.5	80.6	77.6	92.4
	ABMT + one more branch	68.1	80.7	76.2	90.4

Table 4 .

 4 Ablation studies on structure of asymmetric branches.

	Loss	Market → Duke Duke → Market mAP Rank1 mAP Rank1
	ABMT	69.1	82.0	78.3	92.5
	ABMT w/o L ce	52.5	69.6	57.5	79.8
	ABMT w/o L sce	66.7	79.8	77.7	92.2
	ABMT w/o L stri 64.7	78.5	75.5	91.2

Table 5 .

 5 Ablation studies on loss functions.

Code at https://github.com/chenhao2345/ABMT.


Marie Balasse
email: marie.balasse@mnhn.fr

Léo Renault-Fabregon

Henri Gandois

Denis Fiorillo

John Gorczyk

Krum Bacvarov

Maria Ivanova

Neolithic sheep birth distribution: results from Nova Nadezhda (sixth millennium BC, Bulgaria) and a reassessment of European data with a new modern reference set including upper and lower molars

Keywords: Sheep, birth seasonality, stable oxygen isotopes, modern reference, Nova Nadezhda

    Neolithic sheep birth distribution Results from Nova Nadezhda (sixth millennium BC, Bulgaria) and a reassessment of European data with a new modern reference set including upper and lower molars

Introduction

The Neolithic diffusion of domestic plants and animals across Europe from the seventh to the fourth millennia BC involved the reshaping of agropastoral systems in keeping with local constraints and resources, leading to the adaptation of practices as well as changes in animal and plant physiology [START_REF] Balasse | A key to the adaptation of Neolithic husbandry in the Orkneys: contribution of seaweed to the sheep diet at the Holm of Papa Westray, revealed through stable isotope analysis (δ 13 C and δ 18 O) of teeth[END_REF][START_REF] Bogaard | Crop manuring and intensive land management by Europe's first farmers[END_REF][START_REF] Marinova | Crop diversity and choice in prehistoric southeastern Europe: cultural and environmental factors shaping the archaeobotanical record of northern Greece and Bulgaria[END_REF][START_REF] Orton | Between the Danube and the Deep Blue Sea: Zooarchaeological Meta-Analysis Reveals Variability in the Spread and Development of Neolithic Farming across the Western Balkans[END_REF][START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF][START_REF] Salavert | Agricultural Dispersals in Mediterranean and Temperate Europe[END_REF][START_REF] Ivanova | Pioneer farming in southeast Europe during the early sixth millennium BC: Climate-related adaptations in the exploitation of plants and animals[END_REF]. In particular, the dispersal of plants and animals initially domesticated in south-eastern Anatolia towards higher latitudes led to selection in the photoperiod responsive biological cycles of imported domesticates. This led, for example, to the differential expansion across Europe of two variants of domestic barley with different flowering-time adaptations [START_REF] Jones | Population-Based Resequencing Reveals That the Flowering Time Adaptation of Cultivated Barley Originated East of the Fertile Crescent[END_REF]. The timing of animal breeding (i.e. when, in the yearly cycle, did births occur) was also impacted. This may have been especially true for sheep, which, in temperate latitudes, show strongly photoperiodic-controlled seasonal reproductive behavior, inherited from their wild ancestor [START_REF] Hafez | Studies on the breeding season and reproduction of the ewe. Part I. The breeding season in different environments. Part II. The breeding season in one locality[END_REF][START_REF] Karsch | Neuroendocrinal basis of seasonal reproduction[END_REF][START_REF] Malpaux | Control of the circannual rhythm of reproduction by melatonin in the ewe[END_REF][START_REF] Thiéry | Neuroendocrinal interactions and seasonality[END_REF]. Currently, most European sheep are short-day breeders: they enter a period of sexual activity when the length of the day begins decreasing in late summer and up to early winter, with births occurring five months later, from winter to early summer. As an additional consequence of photoperiodic control, the breeding season is further reduced in higher latitudes, with a later onset of sexual activity [START_REF] Hafez | Studies on the breeding season and reproduction of the ewe. Part I. The breeding season in different environments. Part II. The breeding season in one locality[END_REF]. The consecutive constraints with regard to the organization of agro-pastoral tasks and the availability of animal production on an annual scale, were probably just as strong in ancient agropastoral systems as in present-day husbandries [START_REF] Balasse | Environmental constraints on reproductive activity of domestic sheep and cattle: what latitude for the herder?[END_REF][START_REF] Chemineau | Seasonality of reproduction in mammals: intimate regulatory mechanisms and practical implications[END_REF]. Current research aims to clarify the pattern of sheep lambing in the Neolithic (i.e. the spread of births, the shape of the distribution, and the timing of births), on regional and European scales [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF].

Past sheep birth seasonality can be investigated through the reconstruction of the seasonal cycle record in sheep molars, using sequential analysis of the stable oxygen isotope composition (δ 18 O) of enamel [START_REF] Balasse | Determining sheep birth seasonality by analysis of tooth enamel oxygen isotope ratios: the Late Stone Age site of Kasteelberg (South Africa)[END_REF][START_REF] Balasse | Investigating seasonality and season of birth in past herds: a reference set of sheep enamel stable oxygen isotope ratios[END_REF]. A review of the data obtained from European sites dated from the sixth to the third millennia BC demonstrated a restricted breeding season for sheep, with lambing occurring over a three-to-four-month period, from late winter to early summer [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF]. Later lambing was observed at higher-latitude sites (Orkney), confirming a rapid adaptation of sheep breeding behavior to the photoperiodic conditions of North Atlantic Europe [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF]. Further work should help to refine the picture on a regional scale. Nevertheless, two methodological points were left unresolved. Firstly, some of the assemblages included in previous studies considered upper and lower third molars together [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF]. Although no major delay is expected in the timing of the formation of upper and lower molars in sheep (Milhaud and Nezit, 1991), the comparability of the results obtained from both molar rows still requires careful investigation, given the high resolution of δ 18 O sequential sampling. Secondly, the great majority of δ 18 O sequential datasets used to reconstruct sheep birth seasonality were obtained from the third molar. This choice is often dictated by the necessity to secure tooth identification in assemblages mostly composed of isolated teeth, where second and first molars may be mixed up. Yet, available modern reference sets, used to determine the season of birth, were obtained from sheep second molars [START_REF] Blaise | Seasonality and season of birth of modern and late Neolithic sheep from south-eastern France using tooth enamel δ 18 O analysis[END_REF][START_REF] Balasse | Investigating seasonality and season of birth in past herds: a reference set of sheep enamel stable oxygen isotope ratios[END_REF]. Although the lack of modern references on the third molar does not preclude an adequate description of the distribution of births (seasonality of births i.e. the spread of the lambing period and the number of peaks), the season of births (the time of the year when births occurred) cannot be directly determined [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF].

The present study addresses these two methodological points by introducing a new modern reference set; the Kemenez sheep herd (Brittany). The sample includes paired upper and lower third molars in order to test for comparability between the sequences of δ 18 O values measured in both teeth as well as to provide reference values for winter births in the third molar. This new reference set interlinks to previous ones (Blaise andBalasse, 2011 Balasse et al., 2012) through the second lower molar. The synthesis for Neolithic Europe published in [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF] is then revised through a normalization of the data obtained from upper and lower M3s. Additionally, we provide a new set of results for sheep from the archaeological assemblage from Nova Nadezhda in Bulgaria (early sixth millennium BC), mixing upper and lower molars, in order to expand data relating to the earliest stages of the introduction of sheep to Europe. Previously published data for sheep birth seasonality in the Early Neolithic of Southeast Europe include the Starčevo-Criş 1 assemblages of Măgura-Boldul lui Moş Ivănuş in the Northern Balkans (Romania; [START_REF] Balasse | Early herding at Magura-Boldul lui Mos Ivanus (early sixth millennium BC, Romania): environments and seasonality from stable isotope analysis[END_REF] and Alsónyék-Bátaszék in the Hungarian plains [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF]. In the lower Danube and the Carpathian basin where Măgura and Alsónyék-Bátaszék lie, occupying an intermediate position between the Mediterranean and the Central European bioclimatic zones, the prevailing type of husbandry is characterized by a strong predominance of ruminants among livestock, a higher relative proportion of cattle and milk exploitation as evidenced from organic residue analysis in pottery. This contrasts with the sub-Mediterranean Southern Balkans where husbandry focused primarily on caprines and where the incidence of dairy residues in pots is low [START_REF] Ethier | Earliest expansion of animal husbandry beyond the Mediterranean zone in the sixth millennium BC[END_REF]. Investigating sheep birth seasonality at Nova Nadezhda should therefore determine whether a disparity can also be discerned in the rhythms of husbandry systems, between this assemblage located in the Southern Balkans and previously acquired data in the temperate area of the Balkans.

Material

Kemenez is a small island (27 hectares) located in the Molene archipelago, off Brest and the western coast of Brittany (latitude 48°N; Figure 1). The region is subjected to an oceanic climate characterized by mild temperatures (around 7.5°C on average in winter and 15°C in summer) and high regular rainfalls. Monthly precipitation amounts exceed 100 mm from October to February and are the lowest in July [START_REF] Lambert | Human-induced river runoff overlapping natural climate variability over the last 150 years: Palynological evidence (Bay of Brest, NW France)[END_REF]. Data related to the oxygen isotopic composition of precipitation (δ 18 O P ) at the nearest Global Network of Isotopes in Precipitation (GNIP) station at Brest Plouzane is available for the years 1998-2002. Over this time period, the monthly δ 18 O P values were comprised between -5‰ and -6 ‰ in the winter, and -3 ‰ and -4‰ in the summer on average; however, interannual variability is high. The yearly amplitude of variation varied from 3.6 ‰ in 2002 to 7.5 ‰ in 1998 (IAEA/WMO, 2020). The Kemenez reference set is composed of 20 sheep from a Ouessant x Landes de Bretagne cross breed, raised on the island of Kemenez (KMZ; Molène, Finistère). In this herd, sheep are kept under free-range conditions, with no extra food supply. Uncontrolled reproduction leads to winter lambing. Over a seven-year period from 2008 to 2013 and in 2015, a total of 255 referenced births occurred, mainly from mid-January to mid-April (N= 240 or 94 %), with a marked peak over a six-week period during February until mid-March (N = 187 or 73 %). A small number of out-of-season births are also observed in summer (August) and autumn (October) (Figure 2). The specimens used for stable isotope analysis were retrieved from the semi-articulated skeletons of individuals who had died of natural causes or were slaughtered due to a pathological condition, collected between 2013 and 2017. In total, the KMZ sheep reference set is composed of 31 teeth, including four lower second molars (M2) and 19 lower third molars (M3), eight of which were paired with an upper M3 (Table 1).

The archaeological site of Nova Nadezhda, in south-eastern Bulgaria (Figure 1), is situated on a lowlying alluvial terrace in the Middle Maritsa Valley floodplain. The site covers an area of more than five hectares and consists of two low tell-like formations separated by a stream. Excavations carried out in the northwest tell revealed continuous occupation from the beginning of the Early Neolithic (c. 6000 cal BC) to the Final Chalcolithic (c. 4000 cal BC), then again in the Early Iron Age. The Early Neolithic occupation spanned approximately four centuries (6000 -5600/5500 cal BC). The archaeological features include five concentric enclosure ditches used successively from the innermost (earliest) to the outermost (latest). The latest ditch may have been combined with a palisade. One Early Neolithic building was excavated, overlapping three of the inner ditches and most probably contemporaneous with the latest ditch. The material culture is typical of the Karanovo I period [START_REF] Bacvarov | The dead and the nested pots: an early Neolithic ditch burial at Nova Nadezhda, Bulgarian Thrace[END_REF]. [START_REF] Digard | Techniques des nomades Baxtyari d'Iran[END_REF]-1982) ;Brest Plouzane (1998-2002) ;Draix (2008Draix ( -2013) ;) ;Zagreb (1975Zagreb ( -1995) ;) ;Ramnicu Valcea (2012-2016) and Edirne (2008-2016) (IAEA/WMO, 2020). The early Neolithic faunal assemblage at Nova Nadezhda is primarily made up of bones from domesticated caprines-sheep and goat. Where distinctions between the two taxa could be made, sheep are overwhelmingly more abundant than goats. It would seem that caprines were the primary focus of husbandry, although cattle herding and pig rearing were also practiced. Hunting also seems to have been an important activity, and deer (mostly red deer and roe deer) were the second most abundant taxon after caprines. The mortality profiles of sheep mandibles suggest a mixed kill-off pattern not focusing on any single product. Mandibles from the earliest age group (0-2 months) are completely absent, with slaughtering peaks at 1-2 years and again at 4-6 years. The low sample size (n= 22 mandibles) precludes any definitive assessment of caprine husbandry, but it is more than likely that the slaughter pattern was determined by the need to ensure a stable population of animals for meat, breeding, and exchange. Organic residue analysis of pottery detected animal fats in 20 % of the corpus (total n=45), among which 89 % and 11 % were identified as ruminant and nonruminant adipose fats respectively. No dairy fats were identified [START_REF] Ethier | Earliest expansion of animal husbandry beyond the Mediterranean zone in the sixth millennium BC[END_REF].

J F M A M J J A S O N D d 18 O (‰) Edirne -13 -11 -9 -7 -5 -3 J F M A M J J A S O N D d 18 O (‰) Zagreb -15 -13 -11 -9 -7 -5 -3 J F M A M J J A S O N D d 18 O (‰) Ramnicu Valcea
The sheep remains included in this study are composed of 11 lower third molars and three upper third molars from the fill of the second earliest to the latest ditches, with a majority of teeth recovered from the second latest ditch (nine teeth out of a total of 14) (Table 2).

Lower M2 Lower M3 Upper M3 KMZ Ovis 1 1** KMZ 10013 1** 1* KMZ 00018 1** 1** KMZ 20001 1 1 1 KMZ 20038 1 1 KMZ 20025 1 1 KMZ 10026 1 1 KMZ 10033 1 1 KMZ 40001 1 1 KMZ c123 1 1 KMZ c119 1 1 KMZ c122 1 KMZ 00004 1 KMZ c116 1 KMZ c147 1 KMZ 10046 1 KMZ 10051 1 KMZ c128 1 KMZ 00017 1* KMZ c127 1* TOTAL 4 19 8

Principles and methods

Stable oxygen isotope ratios (δ 18 O) from the mineral fraction (bioapatite) of vertebrate skeletons are indirectly linked to local precipitation δ 18 O [START_REF] Land | Isotopic ecology of deer bones[END_REF][START_REF] D'angela | Oxygen isotopes in living mammal's bone phosphate: further results[END_REF]. Body water tracks changes in surface water δ 18 O values through ingested water. Bioapatite precipitates in oxygen isotopic equilibrium with body water, offset by temperature-dependent fractionation held constant in homeothermic mammals [START_REF] Longinelli | Oxygen isotopes in mammal bone phospate: A new tool for paleohydrological and paleoclimatological research[END_REF]. The oxygen isotopic composition of precipitation (δ 18 O P ) varies with the origin and the trajectory of the air mass. At the location where condensation occurs it is then affected by the air temperature and precipitation amount [START_REF] Rozanski | Isotopic patterns in modern global precipitation[END_REF]. In continental Europe, the δ 18 O P values are mainly affected by the air temperature and its seasonal variations: the highest δ 18 O P values are recorded when air temperature is the highest in the summer months, the lowest when air temperature is the lowest in the winter months [START_REF] Rozanski | Isotopic patterns in modern global precipitation[END_REF]. The seasonal variations in δ 18 O P values in the present days at Nova Nadezhda and Eastern Europe are consistent with this pattern (Figure 1). However, correlation between δ 18 O P values and rainfall amount have been observed in maritime temperate regions, including in the British Isles [START_REF] Tyler | Spatial patterns in the oxygen isotope composition of daily rainfall in the British Isles[END_REF] and in particular in northern Scotland [START_REF] Fuller | Isotope hydrology of dripwaters in a Scottish cave and implications for stalagmite palaeoclimate research[END_REF]. This would apply the late Neolithic archaeological sites in Orkney. The data recorded at Brest Plouzane, the nearest GNIP station to Kemenez, similarly shows inverse correlation between the rainfall amount and δ 18 O P values. At these locations where the amount effect applies, the highest δ 18 O P values are recorded when the monthly precipitation is the lowest -in the summer months, the lowest δ 18 O P values when the monthly precipitation are the highest, in the winter months (IAEA/WMO, 2020; Figure 1). Consequently, although different controls apply to the precipitation δ 18 O values in the different study areas (Eastern Europe /North Western Europe maritime regions), the seasonal patterns are comparable at all locations.

Sheep molar crowns grow over a one to one-and-a-half-year period, during which enamel bioapatite records changes in environmental stable isotope values, through diet. A sequential sampling procedure detects cyclical variations in δ 18 O values, reflecting the seasonal cycle. As the timing of tooth growth is fixed within a species -notwithstanding some degree of variability (Milhaut & Nezit,1991) on which the influence of sex, breeding age and nutritional regime were shown to be statistically non-significant in sheep molars [START_REF] Worley | The Sheep Project (2): The effects of plane of nutrition, castration and the timing of first breeding in ewes on dental eruption and wear in unimproved Shetland sheep[END_REF] -the season of birth determines the sequence of the annual cycle recorded in a given tooth [START_REF] Bryant | Biologic and climatic signals in the oxygen isotopic composition of Eocene-Oligocene equid enamel phosphate[END_REF][START_REF] Fricke | Inter-and intra-tooth variation in the oxygen isotope composition of mammalian tooth enamel phosphate: implications for palaeoclimatological and palaeobiological research[END_REF].

Variability in the birth season may be described through the comparison of the position of the maximum value of the δ 18 O cycle in the tooth crown [START_REF] Balasse | Determining sheep birth seasonality by analysis of tooth enamel oxygen isotope ratios: the Late Stone Age site of Kasteelberg (South Africa)[END_REF]. A quantitative estimation of inter-individual variability involves modelling the seasonality in δ 18 O sequences and normalizing the distances in order to eliminate the tooth size factor [START_REF] Balasse | Investigating seasonality and season of birth in past herds: a reference set of sheep enamel stable oxygen isotope ratios[END_REF].

The δ 18 O values were measured in the carbonate fraction of enamel bioapatite. Enamel was sequentially sampled following the procedure described in [START_REF] Balasse | Determining sheep birth seasonality by analysis of tooth enamel oxygen isotope ratios: the Late Stone Age site of Kasteelberg (South Africa)[END_REF] on the posterior (M2), middle (lower M3) or anterior lobe (upper M3). Archaeological enamel powders from Nova Nadezhda were pre-treated to eliminate diagenetic carbonates (0.1 M acetic acid for 4h at room temperature, 0.1ml/mg; [START_REF] Tornero | Seasonality and season of birth in early Eneolithic sheep from Cheia (Romania): methodological advances and implications for animal economy[END_REF]. Oxidizing agents, aimed at removing organic matter, are not deemed necessary in enamel pretreatment and were shown to cause isotopic shift due to unwanted bioapatite dissolution or incorporation of exogenous carbonates [START_REF] Pellegrini | Comparing bioapatite carbonate pre-treatments for isotopicçmeasurements: Part 2 -Impact on carbon and oxygen isotope compositions[END_REF]. For this reason they were not applied. No pretreatment was applied to enamel from modern teeth. Bioapatite samples weighing 580-620 μg were reacted with 100 % phosphoric acid at 70°C in individual vessels in an automated cryogenic distillation system (Kiel IV device), interfaced with a DeltaVAdvantage isotope ratio mass spectrometer. The analytical precision within each run, estimated from six to eight analyses of our laboratory carbonate standard (Marbre LM, expected value -1,83 ‰ calibrated to the NBS19 international standard) was lower than 0.05‰ for δ 18 O values. Over the period of analysis of the Kemenez and Nova Nadezhda enamel samples, the analysis of 220 Marbre LM gave an average δ 18 O value of -1.99±0.04 ‰ (Supplementary Material 1).

The δ 18 O sequences were modelled using an equation derived from a cosine function [START_REF] Balasse | Investigating seasonality and season of birth in past herds: a reference set of sheep enamel stable oxygen isotope ratios[END_REF] and Table S3), in order to define the position (distance to enamel-root junction) of the maximum value (x 0 ). The period of the cycle (X, distance over which the isotopic record covers one annual cycle) was used to normalize the distances (x 0 /X) in order to eliminate inter-individual variability in tooth size. Each modelled δ 18 O sequence yields a x 0 /X ratio which may be directly compared for all specimens. The distribution (spread and shape) of these ratios within a population is an expression of the distribution of births. We use a circular representation of the data to reflect the cyclical nature of seasonality: January follows December, i.e. when (x 0 /X) reaches 1, it also reaches 0. Because of the small size of our archaeological samples, we applied Student statistics. Births occurring outside the 95 % confidence interval were considered as outliers. The main period of birth is defined by the 66 % confidence interval around the mean (excluding the outliers).

The season of birth in archaeological specimens is estimated by comparison with reference x 0 /X ratios obtained from modern sheep whose birth season is known; the Kemenez sheep are used as reference for x 0 /X ratios measured in the M3. This procedure assumes a similar timing of tooth development for modern and ancient sheep breeds.

Results

Results from the sequential analysis of the Kemenez sheep molars are shown in Figure 3. All δ 18 O sequences are shown individually in Supplementary Materials 2 and 3. In the lower M3s, the lowest δ 18 O value for each intra-tooth sequence varies from -4.91 ‰ to -3.33 ‰ (mean-4.20 ± 0.41 ‰) and the highest δ 18 O value varies from -2.78 ‰ to -0.11 ‰ (mean -1.47 ± 0.66 ‰). The intra-tooth amplitude of variation is 2.73 ‰ on average (from 1.21 ‰ to 4.38 ‰). The M2s delivered comparable results (Supplementary Materials 2). All sheep included in this reference set were born in 2008 or after. The oxygen isotopic composition of precipitation during these sheep's lifetime is not known. A comparison with the data available for previous years (1998)(1999)(2000)(2001)(2002) suggests that the δ 18 O values measured in the sheep enamel carbonate are higher than the δ 18 O P values, although the extant of the shift cannot be estimated due to great interannual variability in δ 18 O P values (IAEA/WMO, 2020). The amplitude of intra-tooth variation in the δ 18 O values is also muted compared to the annual amplitude in the δ 18 O P values (3.6 ‰ to 7.5 ‰; IAEA/WMO, 2020). This is explained by prolonged mineralization of enamel, leading to the attenuation of the amplitude of cyclical variations [START_REF] Balasse | Potential biases in sampling design and interpretation of intra-tooth isotope analysis[END_REF]. 3. In most cases, they exceed the analytical precision for the IRMS measurement. However, no systematic relationship may be observed between δ 18 O values measured in the upper/lower molar. In KMZ 10033 M3, the difference is great, with a mean δ 18 O value 0.84‰ higher in the upper molar (Figure 4 A and B). This shift is unexplained.

Results from the modelling of the δ 18 O sequences measured in KMZ M2s and M3s are shown in Figure 5 (and Supplementary Material 4). In three of the 19 lower M3s, the results from the modelling were considered unreliable: in KMZ10013 M3 and KMZc127 M3 because the measured δ 18 O sequence did not show a clearly identified maximum (Supplementary Material 3); in KMZ 00017M3 because it combined a rather short sequence (22 mm), a very low amplitude of variation (1.2 ‰), and also the lowest period (X= 16.11 mm) of the data set (average 22.45±3.57 ‰ in lower M3s, N= 16; Supplementary Material 4). In the second molar, (x 0 /X) ratios vary between 0.14 and 0.20. These values are compared with the previous reference data sets from Carmejane (CAR, Blaise and Balasse, 2011) and Rousay (ROU, [START_REF] Balasse | Investigating seasonality and season of birth in past herds: a reference set of sheep enamel stable oxygen isotope ratios[END_REF] on Figure 5A. In the lower M3s, (x 0 /X) ratios vary between 0.001 and 0.98. Among those, eleven (x 0 /X) ratios are included in the 66 % confidence interval with (x 0 /X) values comprised between 0.07 and 0.26, constituting the main group of births (Supplementary Material 4 and Figure 5B). When the average (x 0 /X) value is set between February and March to fit with the observed distribution of births (Figure 2), KMZ c123 and KMZ 10046, which yielded (x 0 /X) ratios of 0.001 and 0.98 (or -0.02 considering the cyclicity of this parameter) respectively, reflect early winter births (late December and January). KMZ c116 and KMZ c122, with (x 0 /X) ratios of 0.42 and 0.48 respectively, point to the occurrence of births in May and June. KMZ40001, with (x 0 /X) values of 0.77 and 0.85 in lower and upper M3s respectively, represents an October birth (Figure 5B). O values measured in each intra-tooth sequence, amplitude of intra-tooth variation (A) and mean (=(min+max)/2). Difference between upper and lower third molars (∆ = upper -lower). O value (x 0 /X) (A) in the second (M2) and (B) in the third molar (M3). Comparison with previous reference values measured in the M2 at Carmejane (CAR; squares; January/February and September births; [START_REF] Blaise | Seasonality and season of birth of modern and late Neolithic sheep from south-eastern France using tooth enamel δ 18 O analysis[END_REF] and Rousay (ROU; diamonds; April/May births, [START_REF] Balasse | Investigating seasonality and season of birth in past herds: a reference set of sheep enamel stable oxygen isotope ratios[END_REF]. Mean and 66% confidence interval. The outlier (grey) is outside of the 95% confidence interval.

Considering paired M3s, the upper M3 (x 0 /X) value is similar to or higher than in the lower M3, with an average difference of 0.073±0.04, or approximately one month (26 ±14 days) (Supplementary Material 4 and Figure 5B). In all cases, a higher (x 0 /X) ratio for the upper M3 is explained by a higher x 0 value in this tooth (+ 0.47 to +2.78) compared to the lower M3, which does not correspond to a longer period (X) in the upper M3. In some cases, the difference is actually further accentuated by a shorter period in the upper M3 (up to -2.4 mm in KMZ 10033 and KMZ 40001). This result could be indicative of a later onset of tooth growth in the upper dentition. Although acquired on a limited number of individuals (N=8) this average difference between upper and lower M3s (x 0 /X) ratio will be used as a first approximation to correct the data acquired in the upper third molar of archaeological specimens. Ovis 12, 13 and14) cannot be matched with any of the sequences measured in the lower molars (Supplementary Material 6), confirming that all teeth belong to different individuals. The position in the tooth crown of the maximum δ 18 O value normalized to the period of the annual cycle (x 0 /X) varies from 0.17 to 0.47 in lower molars; from 0.39 to 0.40 in upper molars. Considering corrected (x 0 /X) values for upper molars (-0.073), the amplitude of variation for all 14 specimens from Nova Nadezhda is 0.30, or approximately 3.6 months. The main peak of births, as defined by the 66% confidence interval, occurs over 2 months. Sheep births distribution at Nova Nadezhda is compared to previous datasets from archaeological sites dated from the sixth to the third millennia BC [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF] in Figure 8, where all previously published (x 0 /X) ratios measured in upper M3s have been normalized to the lower M3s by correcting by -0.073. The timing of sheep births in the archaeological sites is determined by comparison with the KMZ reference value for winter births in Figure 9. The timing of sheep births at Nova Nadezhda is set within the first half of the spring.

Discussion

Correspondence between the new Kemenez and previous modern sheep reference sets

Considering the lower M3s at Kemenez, when the average measured (x 0 /X) ratio is set between February and March to fit with the observed distribution of births (Figure 2), there is a good agreement between the expected shape of the distribution and the one reconstructed from the stable isotope analysis, except for one individual whose birth is set in late December (Figure 5B) while the earliest recorded births occurred in early January. Overall, this suggests minimal interindividual variability in the timing of the lower M3 growth within the sampled population. Considering the M2, when comparing the reference values obtained from the Kemenez sheep and those previously obtained from Carmejane (CAR; [START_REF] Blaise | Seasonality and season of birth of modern and late Neolithic sheep from south-eastern France using tooth enamel δ 18 O analysis[END_REF] and Rousay (ROU; [START_REF] Balasse | Investigating seasonality and season of birth in past herds: a reference set of sheep enamel stable oxygen isotope ratios[END_REF] respectively on a Préalpes du Sud breed and a Shetland cross, there is a good agreement between the three sets as long as their relative positioning in the yearly cycle is considered (Figure 5A): the (x 0 /X) ratios measured in KMZ, reflecting mainly February and March births, are correctly ranked in between the January and February births at Carmejane and the April/May births at Rousay. When setting the sequence of births within the yearly cycle using the KMZ as a reference for February/March, the Carmejane births are correctly attributed to February and September (although one belonging to the last group is wrongly attributed to August), while the Rousay births are globally set in March and April instead of April/May. This slight discrepancy may result from interbreed variability in the timing of tooth growth, although this would need to be investigated further on a higher number of individuals and breeds in the future. This creates uncertainty in the attribution of birth to a month of the year and imposes the season as a more reliable resolution scale. The main group of births at Kemenez provides values (0.07 < x 0 /X < 0.26) that will be used as reference values for winter births when defining the lambing period in past breeds.

Sheep birth distribution at Nova Nadezhda

The results obtained at Nova Nadezhda point to a main peak of births over a two-month period. This spread and shape of of births distribution are similar to those observed in wild ovine populations, as well as in domestic herds with no human control over reproduction. Observations in European mouflon populations report peaks of fertilization/births over periods of three to seven weeks [START_REF] Hadjisterkotis | Reproduction of Cyprus mouflon Ovis gmelini ophion in captivity and in the wild[END_REF][START_REF] Garel | Reproductive output of female mouflon (Ovis gmelini musimon x Ovis sp.): a comparative analysis[END_REF]. At Kemenez, where reproduction occurs without human control and males and females are kept together all year round, the main peak of births occurs over a six-week period. However, a lower number of births directly precedes and follows this main peak, and isolated births occur out-of-season in the summer and autumn (Figure 2). In domestic herds, reproduction can be controlled by separating males and females and reintroducing males at the onset of the breeding season. This practice, aimed at synchronizing births, leads to massive lambing at the very beginning of the lambing season, and the absence of out-of-season births [START_REF] Etienne | Des moutons en forêt littorale varoise[END_REF]. At Nova Nadezhda, the occurrence of births before and after the main peak does not indicate this kind of control over sheep reproduction. At this site, the non-occurrence of out-of-season births -observed elsewhere in Neolithic European assemblages [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF] and Figure 8) -may also be due to sampling bias given the small assemblage size. Birth synchronization is observed in modern industrialized as well as in traditional husbandry and may be driven by different factors. In the case of large herds, in view of the intense effort required to monitor lambing and take care of the youngest animals, synchronization allows farmers to concentrate on this task at a given time in order to reduce time conflicts with other tasks in the agropastoral calendar. For nomadism and herding systems involving mobility, synchronized births are preferable because young animals are vulnerable to long-distance moves. In modern husbandry, another objective is to be able to provide the market with large quantities of products at the most advantageous time [START_REF] Barth | The nomads of South Persia[END_REF][START_REF] Digard | Techniques des nomades Baxtyari d'Iran[END_REF][START_REF] Mace | Nomadic pastoralists adopt subsistence strategies that maximize long-term household survival[END_REF][START_REF] Etienne | Des moutons en forêt littorale varoise[END_REF]. The absence of this type of control of sheep reproduction at Nova Nadezhda may simply be due to the lack of need for it.

Although the orientation of production cannot be clearly inferred from the slaughtering profile on account of the low number of caprine dental remains, a strong focus on seasonal production (milk, tender meat) is not evidenced either.

Sheep birth season

After correction of the data obtained from upper M3s at Nova Nadezhda and all previously published European assemblages (Figures 8 and9), the main conclusions in [START_REF] Balasse | Sheep birth distribution in past herds: a review for prehistoric Europe (6 th to 3 rd millennia BC)[END_REF] are refined here. In all instances, the main peak of sheep births occurred over a period of two to three months except at Alsónyék-Bátaszék in the Hungarian plains, where the peak is spread over four months with an earlier onset of the breeding season (Figure 9). The small number of individuals included in most assemblages included in this comparison prevents from discussing further the archaeological significance of this find. The new Kemenez reference set shows that births globally occurred in the spring at sites in the Balkan area (Figure 9). Comparatively to the Balkans (Bulgaria, Hungary and Romania), the onset of the breeding season occurred later under higher latitudes (in late spring at Knap of Howar and Skara Brae in Orkney), possibly as a consequence of the photoperiodic control of the fertility cycle in sheep (Hafez, 1962).

Figures 8 and 9 compare the distribution and timing of sheep births at Nova Nadezhda in the Southern Balkans with contemporaneous sites of the Starčevo-Criş 1 culture in more northern parts of the Balkans: Măgura-Boldul lui Moş Ivănuş in Romania and Alsónyék-Bátaszék in the Hungarian plains (Figure 1). In contrast to Nova Nadezhda, where husbandry focused mainly on caprines, cattle played a more important meat-supplying role at Alsónyék-Bátaszék, where they appear in similar proportions to caprines [START_REF] Nyerges | Preliminary report on the Neolithic archaeozoological finds from Alsonyek-Bataszek, Hungary[END_REF][START_REF] Nyerges | Neolithic animal husbandry in the Tolnai-Sarköz region on the basis of the archaeozoological finds from the Alsonyek-Bataszek archaeological site[END_REF] and at Măgura, where the larger body mass of cattle compensates for the predominance of caprines (2/1) [START_REF] Bălăşescu | Omul si Mediul Animal între Mileniile VII-IV[END_REF]. The presence of dairy fat residues in ceramic vessels also attests to milk exploitation at Alsónyék-Bátaszék, whereas no dairy residues were detected at Nova Nadezhda [START_REF] Ethier | Earliest expansion of animal husbandry beyond the Mediterranean zone in the sixth millennium BC[END_REF]. In spite of different climatic settings and different herding strategies in terms of flock composition and the orientation of animal exploitation, the results do not suggest any major shift in the timing of sheep births between the sub-Mediterranean southern and temperate northern parts of the Balkans. With 14 specimens, the Nova Nadezhda sample is the largest ever published for the Early Neolithic in south-eastern Europe and for the European Neolithic as a whole. As a result, in addition to determining the timing of the main birthing period, the shape of sheep birth distribution could also be described. This level of description is probably less reliable for other sites available for comparison in Figure 8 with smaller sample sizes. It is reasonable to assume that this apparent uniformity may be challenged in the future with enlarged datasets.

In Romania, the early Chalcolithic assemblage at Cheia (early fifth millennium cal BC, Hamangia culture) and the Chalcolithic assemblages at Hârşova-tell and Borduşani-Popină (second half of the fifth millennium cal BC, Gumelniţa A2 culture; Figure 1) allow for the investigation of the stability of sheep birth seasonality over time. Over this time period in Romania, the composition of domestic stock was very different among cultural groups. At Cheia, as well as at other sites from the Hamangia culture, caprines were outnumbered by cattle but nonetheless played an important second role and pigs were rare [START_REF] Bălăşescu | Oameni şi animale. Strategii şi resurse la comunităţile preistorice Hamangia şi Boian. Bucureşti[END_REF][START_REF] Balasse | Cattle and sheep herding at Cheia, Romania, at the turn of the fifth millennium cal BC: a view from stable isotope analysis[END_REF]. In the Gumelniţa, where pigs were considerably more important but cattle and caprines were nonetheless predominant, sheep were probably exploited for meat, with a focus on early culling, and cattle mortality profiles suggest milk exploitation [START_REF] Bréhard | What's behind the tell phenomenon? An archaeozoological approach of Eneolithic sites in Romania[END_REF]. At all three Chalcolithic sites, the timing of sheep births remained similar compared to birthing at Măgura and at all Early Neolithic sites in the Balkans as a whole (Figures 8 and9). However, again, reduced sample sizes prevent more detailed comparisons of the shape of sheep birth distribution and may possibly contribute to this apparent uniformity.

Conclusions and recommendations

The Kemenez sheep provide a solid new reference set of δ 18 O sequences for winter births and the very first reference set for the sheep third molar. It also provides the first direct comparison of δ 18 O sequences measured in upper and lower sheep third molars. The results from this comparison on eight paired third molars showed a moderate shift in the positioning of the isotopic sequences in tooth crowns: the maximum δ 18 O were recorded at a greater distance from the enamel-root junction i.e. at an earlier stage of tooth formation in the upper compared to the lower third molar. The consecutive difference in the final assessment of the time of year when births occurred was estimated at one month on average (26±14 days), with a later estimation when inferred from the upper molar compared to the lower molar. This discrepancy may be considered significant or not depending on the aim pursued. In temperate Europe, given the strongly seasonal reproductive behavior of sheep, the timing of the main birthing season can be defined with a limited number of sheep teeth (winter/spring versus autumn/winter, for example, which are currently the two main patterns). When pursuing this objective, upper and lower molars may be directly combined without introducing too much of a bias in interpretation. However, beyond this goal, a more accurate definition of the shape of birth distribution (i.e., spread of the main peak of births and positioning in the yearly cycle), would lead to refined comparisons between sites and over time, but also to conclusions on the degree of control by herders of sheep reproduction. To attain this objective, it is clear that a higher number of specimens is required -the 14 sheep from Nova Nadezhda were just enough to outline birth distribution at the site. The analysis of larger datasets is recommended, and this may be attained by the combined use of upper and lower molars, and would in this case require the normalization of the results obtained from upper and lower teeth. In this regard, the correction used in the present study was obtained on a limited number of specimens (N=8) and must be refined on a larger sample size.

The KMZ reference set also provided four δ 18 O sequences measured in the M2. A comparison with other reference sets (Carmejane, [START_REF] Blaise | Seasonality and season of birth of modern and late Neolithic sheep from south-eastern France using tooth enamel δ 18 O analysis[END_REF]Rousay, Balasse et al., 2012) showed a good agreement between the three sets as long as their relative positioning in the yearly cycle was considered. Considering the absolute positioning in the yearly cycle (month of birth) and when setting the KMZ reference values at February/March, there was a slight offset in the attribution of the Rousay births (March/April instead of April/May). This is possibly due to interbreed variability in the timing of tooth growth. Alternatively, it could be related to the small sample size of all modern reference sets currently available. Either way, until this is clarified with enlarged datasets, it is recommended to use the season as the resolution scale rather than the monthly scale When compared to the Kemenez reference set, the extended dataset at Nova Nadezhda provide evidence for a peak of births in early spring at this site, and a birth distribution pattern pointing to the absence of sheep reproduction control by herders; i.e., no separation of males and females. The spring birth periods at Nova Nadezhda are similar to previous data from other Neolithic and Chalcolithic sites in the Balkans. In this area, no latitudinal and chronological trend could be observed for six sites located in the Southern Balkans, the Northern Balkans and the Hungarian plains and dated from the early sixth millennium to the second half of the fifth millennium BC. This apparent uniformity for the spread (2-4 months) and timing (early spring) of the birthing period could be challenged in the future by enlarged datasets, potentially showing different birth distributions and different timing for the main peak of births.

Figure 1 :

 1 Figure 1: Archaeological cultures mentioned in the text. Location of archaeological sites, modern reference sets and Global Network of Isotopes in Precipitation (GNIP) stations. Available modern data for average monthly oxygen isotope composition of precipitation at Altnabreac[START_REF] Digard | Techniques des nomades Baxtyari d'Iran[END_REF]-1982) ; Brest Plouzane (1998-2002) ; Draix (2008 Draix ( -2013) ;) ; Zagreb (1975 Zagreb ( - 1995) ;) ; Ramnicu Valcea (2012-2016) and Edirne (2008-2016) (IAEA/WMO, 2020).

Figure 2 :

 2 Figure 2: Distribution of births in the Kemenez sheep herd over the years 2008-2013 and 2015. Jan1: first two weeks of January; Jan2 last two weeks of January. 17 individuals are missing from this count (unknown dates of birth).

Figure 3 :

 3 Figure 3: Results from the sequential analysis of the Kemenez (KMZ) sheep second (M2) and third molars (M3). Each sample is located in the tooth crown using its distance from the enamel-root junction. The δ 18 O sequences measured in upper and lower M3s are compared in Figure 4A. The difference in the δ 18 O values measured in paired lower and upper M3s are reported in Table3. In most cases, they exceed the analytical precision for the IRMS measurement. However, no systematic relationship may be observed between δ 18 O values measured in the upper/lower molar. In KMZ 10033 M3, the difference is great, with a mean δ 18 O value 0.84‰ higher in the upper molar (Figure4A and B). This shift is unexplained.

Figure 4 :

 4 Figure 4: A-Comparison of the δ 18 O sequences measured in upper and lower third molars (M3) from Kemenez. B -Comparison of the normalized location in tooth crown of the maximum δ 18 O value (x 0 /X) in upper and lower third molars.

Figure 5 :

 5 Figure 5: Kemenez (KMZ) sheep births positioned in the yearly cycle, as reflected from the normalized location in tooth crown of the maximum δ 18

  sequential analysis of sheep M3 from Nova Nadezhda are shown in Figure 6. The δ 18 O sequences are individually shown in Supplementary Material 5. The δ 18 O values vary between -8.6 ‰ and +0.1 ‰. Intra-tooth variation follows a marked sinusoidal pattern with amplitude varying between 2.5 ‰ and 7.0 ‰ (mean 5.3 ‰). Results from the modelling of the δ 18 O sequences are shown in Figure 7 (and Supplementary Material 6). The δ 18 O sequences retrieved from the three upper molars (NND

Figure 6 :

 6 Figure 6: Results from the sequential δ 18 O analysis of third molars (M3) of 14 sheep from Nova Nadezhda. Each sample is located in the tooth crown using its distance from the enamel-root junction.

Figure 7 :

 7 Figure 7: A -Nova Nadezhda sheep births positioned in the yearly cycle, as reflected from the normalized location in tooth crown of the maximum δ 18O value (x 0 /X) in the third molar (M3). B-After correction of (x 0 /X) ratios obtained in the upper molars.

Figure 8 :

 8 Figure 8: Distribution of sheep births at Nova Nadezhda and other archaeological sites cited in the text. Data: this study, Balasse et al., 2017 and Balasse et al., 2013. All (x 0 /X) ratios measured in upper M3s have been normalized to the lower M3s by correcting by -0.073. In order to enlarge the size of the dataset, Hârşova-tell and Borduşani-Popină, contemporaneous tell sites located in close vicinity and showing great similarity in the orientation of the economy (Bréhard and Bălăşescu 2012), are considered together.

Figure 9 :

 9 Figure 9: Timing of sheep births at Nova Nadezhda and other archaeological sites cited in the text, estimated from the comparison with the modern reference set at Kemenez, defining winter births.

Table 1 :

 1 The

	code	lower/upper	piece ID	Feature	Structure
	NNd Ovis 1 M3	L	53.56	53	latest ditch
	NNd Ovis 2 M3	L	56.4	56	second earliest ditch
	NNd Ovis 3 M3	L	1005.5.8	054/2	second latest ditch
	NNd Ovis 4 M3	L	54.2.89	054/2	second latest ditch
	NNd Ovis 5 M3	L	54.1.27	054/1	second latest ditch
	NNd Ovis 6 M3	L	54.2.95	054/2	second latest ditch
	NNd Ovis 7 M3	L	54.3.55	054/3	second latest ditch
	NNd Ovis 8 M3	L	54.1.79	054/1	second latest ditch
	NNd Ovis 9 M3	L	55.1.7	055/1	third ditch
	NNd Ovis 10 M3	L	53.9	53	latest ditch
	NNd Ovis 11 M3	L	54.3.44	054/3	second latest ditch
	NNd Ovis 12 M3	U	1005.9.2	054/2	second latest ditch
	NNd Ovis 13 M3	U	1004.7.3	054/2	second latest ditch
	NNd Ovis 14 M3	U	55.2.11	55	third ditch

Kemenez sheep composing the modern reference set. (*) could not be modelled; (**) data from Balasse et al. (2017).

Table 2 :

 2 Provenance of the sheep molars from Nova Nadezhda.

Table 3 :

 3 Minimum and maximum δ 18

Acknowledgements

We thank Soizic and David Cuisnier from the Ferme insulaire de Kemenez for providing access to their farm and sheep. K. Bacvarov wishes to acknowledge partial funding from the Bulgarian Ministry of Education and Science under the National Research Program Cultural Heritage, National Memory and Social Development approved by DCM No 577 of 17 August 2018. The stable isotope analyses were conducted at the SSMIM (MNHN Paris). This work was funded by the DFG Food Cultures project (IV 101/5-1, dir. M. Ivanova-Bieg) for the study of the material from Nova Nadezhda and the "CyclOvin" MNHN ATM project (dir. M. Balasse) for the study of the Kemenez sheep. Stéphanie Bréhard helped with gathering the sheep on Kemenez Island. We thank Gaël Obein (LNE CNAM) for his help for data treatment.