Félix Quinton

Loic Landrieu

Crop Rotation Modeling for Deep Learning-Based Parcel Classification from Satellite Time Series

While annual crop rotations play a crucial role for agricultural optimization, they have been largely ignored for automated crop type mapping . In this paper, we take advantage of the increasing quantity of annotated satellite data to propose the first deep learning approach modeling simultaneously the inter-and intra-annual agricultural dynamics of parcel classification. Along with simple training adjustments, our model provides an improvement of over 6.6 mIoU points over the current state-of-the-art of crop classification. Furthermore, we release the first large-scale multi-year agricultural dataset with over 300 000 annotated parcels.

Introduction

The Common Agricultural Policy (CAP) is responsible for allocating agricultural subsidies in the European Union, which nears 50 billion euros each year [START_REF]The common agricultural policy at a glance[END_REF]. Consequently, monitoring the crop types for subsidy allocation represents a significant challenge for payment agencies, which have encouraged the development of automated crop classification tools based on machine learning [START_REF] Loudjani | Artificial Intelligence and EU Agriculture[END_REF]. In particular, The Sentinels for Common Agricultural Policy (Sen4CAP) project [START_REF] Koetz | SEN4CAP Sentinels for CAP monitoring approach[END_REF] aims to provide EU member states with algorithmic solutions and best practice studies on crop monitoring based on satellite data from the Sentinel constellation [START_REF] Drusch | Sentinel-2: ESA's optical high-resolution mission for GMES operational services[END_REF]. Despite the inherent difficulty of differentiating between the complex growth patterns of plants, this task is made possible by the nearly limitless access to data and annotations. Indeed, Sentinel-2 offers multi-spectral observations at a high revisit time of five days on average, which are particularly appropriate for characterizing the complex spectral and temporal characteristics of crop phenology. Moreover, farmers declare the crop cultivated in each of their parcels every year. This represents over 10 million of annotations each year for France alone [START_REF]Registre parcellaire graphique (RPG) : contours des parcelles et îlots culturaux et leur groupe de cultures majoritaire[END_REF], all openly accessible in the Land-Parcel Identification System (LPIS). However, the sheer scale of the problem raises interesting computational challenges: Sentinel-2 gathers over 25Tb of data each year over Europe.

The state-of-the-art of yearly parcel-based crop type classification from Satellite Image Time Series (SITS) is particularly dynamic, especially since the adoption of deep learning methods [START_REF] Sainte | Satellite image time series classification with pixel-set encoders and temporal selfattention[END_REF][START_REF] Pelletier | Deep learning for the classification of Sentinel-2 image time series[END_REF][START_REF] Rußwurm | Self-attention for raw optical satellite time series classification[END_REF]. However, most methods operate on a single year worth of temporal acquisitions and ignore inter-annual crop rotations. In this paper, we propose the first deep learning framework for classifying yearly crop types from information spanning several years, as represented in Figure 1. We show that we can improve their predictions by a large margin 1 Figure 1: Multi-Year Sentinel-2 Data. Details of our area of interest for the three years studied in this article. The crop type of each parcel is represented by the color of a polygon following their contour according to the legend above. This color code is used throughout this article for all figures representing cultivated crops.

Single-Year Crop-Type Classification. Single-year crop-type classification involves the classification of the crop grown in a parcel from a single year worth of observation. Predeep learning parcel-based classification methods rely on such as support vector machines [START_REF] Zheng | A support vector machine to identify irrigated crop types using time-series Landsat NDVI data[END_REF] or random forests [START_REF] Vuolo | How much does multi-temporal Sentinel-2 data improve crop type classification?[END_REF] operating on handcrafted descriptors such as the Normalized Difference Vegetation Index. The temporal dynamics are typically handled with stacking [START_REF] Vuolo | How much does multi-temporal Sentinel-2 data improve crop type classification?[END_REF], probabilistic graphical models [START_REF] Siachalou | A hidden Markov models approach for crop classification: Linking crop phenology to time series of multi-sensor remote sensing data[END_REF], or dynamic warping method [START_REF] Belgiu | Sentinel-2 cropland mapping using pixel-based and object-based time-weighted dynamic time warping analysis[END_REF].

In conjunction with growing data availability, the adoption of deep learning-based methods has allowed for a large increase in performance for parcel-based crop classification. The spatial dimension of parcels is typically handled with convolutional neural networks [START_REF] Kussul | Deep learning classification of land cover and crop types using remote sensing data[END_REF], parcel-based statistics [START_REF] Rußwurm | Self-attention for raw optical satellite time series classification[END_REF], or set-based encoders [START_REF] Sainte | Time-Space tradeoff in deep learning models for crop classification on satellite multispectral image time series[END_REF]. The temporal dynamics are modeled with temporal convolutions [START_REF] Pelletier | Deep learning for the classification of Sentinel-2 image time series[END_REF], recurrent neural networks [START_REF] Sainte | Time-Space tradeoff in deep learning models for crop classification on satellite multispectral image time series[END_REF], hybrid convolutionalrecurrent networks [START_REF] Rußwurm | Multi-temporal land cover classification with sequential recurrent encoders[END_REF], and temporal attention [START_REF] Rußwurm | Self-attention for raw optical satellite time series classification[END_REF][START_REF] Sainte | Satellite image time series classification with pixel-set encoders and temporal selfattention[END_REF][START_REF] Yuan | Self-Supervised Pretraining of Transformers for Satellite Image Time Series Classification[END_REF].

Multiple recent studies [START_REF] Kondmann | DENETHOR: The DynamicEarthNET dataset for Harmonized, inter-Operable, analysis-Ready, daily crop monitoring from space[END_REF][START_REF] Sainte | Lightweight Temporal Self-Attention for Classifying Satellite Images Time Series[END_REF][START_REF] Sainte | Satellite image time series classification with pixel-set encoders and temporal selfattention[END_REF][START_REF] Schneider | Re] Satellite Image Time Series Classification with Pixel-Set Encoders and Temporal Self-Attention[END_REF][START_REF] Sainte | Panoptic Segmentation of Satellite Image Time Series with Convolutional Temporal Attention Networks[END_REF] have solidified the PSE+LTAE (Pixel Set Encoder & Lightweight Temporal Attention) as the state-of-the-art of crop type classification for these reasons. Furthermore, this network is particularly parsimonious in terms of computation and memory usage, which proves well suited for training on multi-year data. Finally, the code is available (https://github.com/VSainteuf/lightweight-temporalattention-pytorch, accessed on 10/10/21). For these reasons, we choose to use this network as the basis for our analysis and design modifications.

Multi-Year Agricultural Optimization. Most of the literature on multi-year crop rotation focuses on agricultural optimization, ie, improving agricultural practices aiming to improve yields. These models generate suggested rotations according to expert knowledge [START_REF] Dury | Models to support cropping plan and crop rotation decisions. A review[END_REF], handcrafted rules [START_REF] Dogliotti | ROTAT, a tool for systematically generating crop rotations[END_REF], or statistical analysis [START_REF] Voskoboynik | Modeling Crop Rotation with Discrete Mathematics[END_REF], while other models are based on a physical analysis of the soil composition [START_REF] Brankatschk | Modeling crop rotation in agricultural LCAs-challenges and potential solutions[END_REF] such as the nitrogen cycle [START_REF] Nina | Modelling optimal crop sequences using network flows[END_REF]. Aurbacher and Dabbert also take a simple economic model into account in their study [START_REF] Aurbacher | Generating crop sequences in land-use models using maximum entropy and Markov chains[END_REF]. More sophisticated models combine different sources of knowledge for better suggestions, such as ROTOR [START_REF] Bachinger | ROTOR, a tool for generating and evaluating crop rotations for organic farming systems[END_REF] or CropRota [START_REF] Schönhart | CropRota-A crop rotation model to support integrated land use assessments[END_REF]. The RPG Explorer software [START_REF] Levavasseur | RPG Explorer: A new tool to ease the analysis of agricultural landscape dynamics with the Land Parcel Identification System[END_REF] uses a second-order Markov Chain for a more advanced statistical analysis of rotations.

Given the popularity of these tools, it is clear that the careful choice of cultivated crops can substantially impact agricultural yields and is the object of meticulous attention from farmers. This is reinforced by the multi-model, multi-country meta-study of Kollas et al. [START_REF] Kollas | Crop rotation modelling-A European model intercomparison[END_REF], showing that multi-year modeling allows for a large increase in yield prediction. Consequently, we posit that a classification model with access to multi-year data will learn inter-annual patterns to improve its accuracy.

Multi-Year Crop Type Classification. Multi-year crop type classification refers to leveraging information (satellite observations, past declarations) to improve the classification of the grown crop type in agricultural parcels. Osman et al. [START_REF] Osman | Assessment of a Markov logic model of crop rotations for early crop mapping[END_REF] propose to use probabilistic Markov models to predict the most probable crop type from the sequence of past cultivated crops of the previous 3 to 5 years. Giordano et al. [START_REF] Giordano | Improved Crop Classification with Rotation Knowledge using Sentinel-1 and-2 Time Series[END_REF] and Bailly et al. [START_REF] Bailly | Crop-rotation structured classification using multi-source Sentinel images and lpis for crop type mapping[END_REF] propose to model the multi-year rotation with a second-order chain-Conditional Random Field (CRF). Finally, Yaramasu et al. [START_REF] Yaramasu | Pre-season crop type mapping using deep neural networks[END_REF] are the first to propose to analyze multi-year data with a deep convolutional-recurrent model. However, they only choose one image per year and hence do not model both inter-and intra-annual dynamics. In contrast, we propose to explicitly our model operates at both the intra-annual scale by using the sequence of yearly observation and the inter-annual scale by considering past declarations.

We list here the main contributions of this paper:

• We propose a straightforward training scheme to leverage multi-year data and show its impact on yearly agricultural parcel classification.

• We introduce a modified attention-based temporal encoder able to model both interand intra-annual dynamics of agricultural parcels, yielding a large improvement in precision.

• We present the first open-access multi-year dataset [START_REF]multi-year dataset[END_REF] for crop classification based on Sentinel-2 images, along with the full implementation of our model.

• Our code is open-source and accessible at the following repository: https://github. com/felixquinton1/deep-crop-rotation.

i 2 Materials and Methods

We present our dataset and proposed method to model multi-year SITS, along with several baseline methods to assess the performance of its components. We denote by [1, I] the set of years for which satellite observations are available to us and use the compact pixel-set format to represent the SITS. For a given parcel and a year i ∈ [1, I], we denote the corresponding SITS by a tensor x i of size C × S × T i , with C the number of spectral channels, S the number of pixels within the parcel, and T i the number of temporal observation available for the year i. Likewise, we denote by l i ∈ {0, 1} L the one-hot-encoded label at year i, denoting which kind of crop is cultivated in the considered parcel among a set L of crop types. Note that, in this article, we focus on the prediction of the main culture, ie, only one crop type per year.

Dataset

Our proposed dataset, represented in Figure 2, is based on parcels within the 31TFM Sentinel-2 tile, covering an area of 110 × 110 km 2 in the South East of France (centered around 4.31N, 46.44E in WGS84). This area is in the Auvergne-Rhône-Alpes region, a major producer of cereal with over 54 000ha of corn and 30 000ha of wheat. Extensive livestock production makes meadow the most common crop type with over 60% of declared parcels in the LPIS. The most frequent crop rotations are permanent cultures (meadows, vineyards, pasture) and alternating between corn, wheat, and rapeseed.

Our satellite time series are constituted of Sentinel-2 level 2A images. We discard the bands B01, B09, and B10 and resample the remaining 10 spectral bands to a spatial resolution of 10m per pixel with bilinear interpolation. Our data spans three years of acquisition: 2018, 2019, and 2020, with respectively 36, 27 and 29 valid entries, see Figure 3. The length of sequences varies due to the automatic discarding of cloudy tiles by the data provider THEIA [START_REF] Baghdadi | The Theia land data centre[END_REF]. We do not apply further preprocessing such as cloud removal or radiometric calibration than what is already performed by the data provider THEIA. We select stable parcels, meaning that their contours only undergo minor changes across the three studied years. We also discard very small parcels (under 800m 2 ) small or with very narrow shapes to reflect the resolution of the Sentinel-2 satellite. Each parcel has a ground truth cultivated crop type for each year corresponding to the main culture as reported by the French LPIS, whose precision is estimated at over 97% as reported by the French Payment Agency. Note that we ignore secondary cultures for parcels with multiple growth cycles. To limit class imbalance, we only keep crop types among a list of 20 of the most cultivated species in the area of interest. In sum, our dataset is composed of 103 602 parcels, each associated with three image time sequences and three crop annotations corresponding to the farmers' declarations for 2018, 2019, and 2020.

The Sentinel2Agri dataset [START_REF] Sainte | Satellite image time series classification with pixel-set encoders and temporal selfattention[END_REF], composed of parcels from the same area, is composed of 191 703 parcels. We can estimate that our selection criteria exclude approximately every other parcel. A more detailed analysis of the evolving parcel partitions across different plots could lead to retaining a higher proportion of the original parcels.

N N As represented in Table 1, the dataset is still imbalanced: more than 60% of declarations correspond to meadows. In comparison, potato is cultivated in less than 100 parcels each year in the area of interest.

1 km N N 1 km (a) (b) (c)

Pixel-Set and Temporal Attention Encoders

The Pixel Set Encoder (PSE) [START_REF] Sainte | Satellite image time series classification with pixel-set encoders and temporal selfattention[END_REF] is an efficient spatio-spectral encoder which learns expressive descriptors of the spectral distribution of the observations by randomly sampling pixels within a parcel. Its architecture is inspired by set-encoding deep architecture [START_REF] Charles R Qi | Pointnet: Deep learning on point sets for 3d classification and segmentation[END_REF][START_REF] Kussul | Deep learning classification of land cover and crop types using remote sensing data[END_REF], and dispenses us from preprocessing parcels into image patches, saving memory and computation. The Temporal Attention Encoder (TAE) [START_REF] Sainte | Satellite image time series classification with pixel-set encoders and temporal selfattention[END_REF] and its parsimonious version Lightweight-TAE (LTAE) [START_REF] Sainte | Lightweight Temporal Self-Attention for Classifying Satellite Images Time Series[END_REF] are temporal sequence encoders based on the language processing literature [START_REF] Vaswani | Attention is all you need[END_REF] and adapted for processing SITS. Both networks can be used sequentially to map the sequence of observations x i at year i to a learned yearly spatio-temporal

e i = TAE PSE x i t Ti t=1 . ( 1
)
l i observations labels x i-2 l i-2 x i-1 x i year i-2 year i-1 year i target D l i-1 l i observations labels x i-2 l i-2 x i-1 x i year i-2 year i-1 year i target D CRF l i-1 T (a) Single-year model (b) CRF model l i observations labels x i-2 l i-2 x i-1 x i year i-2 year i-1 year i target D l i-1 l i observations labels x i-2 l i-2 x i-1 x i year i-2 year i-1 year i target D l i-1 (c) Observation bypass (d) Proposed model

Multi-Year Modeling

We now present a simple modification of the PSE+LTAE network to model crop rotation. In the original PSE+LTAE approach, the descriptor e i is directly mapped to a vector of class scores z i by a Multi-Layer Perceptron (MLP). In order to make the prediction z i covariant with past cultivated crops, we augment the spatio-temporal descriptors e i by concatenating the sum of the one-hot-encoded labels l j for the previous two years j = i -1, i -2. Then, a classifier network D, typically an MLP, maps this feature to a vector z i of L class scores:

z i = D e i l i-1 + l i-2 , (2)
with [•||•] the channelwise concatenation operator. We handle the edge effects of the first two available years by defining l 0 and l -1 as vectors of zero of size L (temporal zero-padding). This model can be trained end-to-end to simultaneously learn inter-annual crop rotations along with the intra-annual evolution of the parcels' spectral statistics.. Our model makes three simplifying assumptions:

• We only consider the last two previous years because of the limited span of our available data. However, it would be straightforward to extend our approach to a longer duration.

• We consider that the history of a parcel is entirely described by its past cultivated crop types, and we do not take the past satellite observations into account. In other words, the label at year i is independent from past observations conditionally to its past labels [START_REF] Martin | Graphical models, exponential families, and variational inference[END_REF]Chap 2]. This design choice allows the model to stay tractable in terms of memory requirements.

• The labels of the past two years are summed and not concatenated. The order in which the crops were cultivated is then lost, but results in a more compact model.

Baseline Models

In order to meaningfully evaluate the performance of our proposed approach, we implement different baselines for classifying parcels from multi-year data . In Figure 4, we represent schematically the main idea behind these baselines and our proposed approach. Note that the choice of backbone network to handle single-year data is out of the scope of this paper.

Single-Year: M single . We do not provide the previous years' labels, and directly map the current year's observations to a vector of class scores [START_REF] Sainte | Lightweight Temporal Self-Attention for Classifying Satellite Images Time Series[END_REF].

Conditional Random Fields: M CRF . Based on the work of [START_REF] Bailly | Crop-rotation structured classification using multi-source Sentinel images and lpis for crop type mapping[END_REF] and [START_REF] Giordano | Improved Crop Classification with Rotation Knowledge using Sentinel-1 and-2 Time Series[END_REF], we implement a simple chain-CRF probabilistic model. We use the prediction of the previous PSE+LTAE, calibrated with the method of Guo et al. [START_REF] Guo | On calibration of modern neural networks[END_REF] to approximate the posterior probability p ∈ [0, 1] L of a parcel having the label k for year i : p k = P (l i = k | x i ) (see Section 3.3 for more details). We then model the second order transition probability p(l i = k | l i-1 , l i-2 ) with a three-dimensional tensor T ∈ [0, 1] L×L×L which can be approximated based on the observed transitions in the training set. As suggested by Bailly et al. , we use a Laplace regularization [START_REF] Schütze | Introduction to information retrieval[END_REF]Chap. 13] to increase robustness. The resulting probability for a given year i is given by:

z i CRF [k] = p T [l i-2 , l i-1 , :] , (3)
with the Hadamard term-wise multiplication. This method is restricted to i > 2 as edge effects are not straightforwardly fixed with padding.

Observation Bypass: M obs . Instead of concatenating the labels of previous years to the embedding e i , we concatenate the average of the descriptors of the last two years e i-1 for e i-2 :

z i obs = D obs      e i      1 2 [e i-1 + e i-2 ] if i > 1 e 0 if i = 1 0 if i = 0       . (4)
Edge effects are handled with mirror and zero temporal padding.

Label Concatenation: M dec-concat . Instead of concatenating the sum of the last two previous years, we propose to concatenate each one-hot-encoded vector l i-1 and l i-2 with the learned descriptor z i . This approach is similar to Equation 2, but leads to a larger descriptor and a higher parameter count.

Single-Year Label Bypass: M dec-one-year . In order to evaluate the impact of describing the history of parcels as the past two cultivated crops, we only concatenate the label of the previous year to the learned descriptor e i . We propose a simple training protocol to leverage the availability of observations and farmers' declarations from multiple years.

Training Protocol

Mixed-year Training:

We train a single model with parcels from all available years. Our rationale is that exposing the model to data from several years will contribute to learning richer and more resilient descriptors. Indeed, each year has different meteorological conditions influencing the growth profiles of crops. Moreover, by increasing the size of the dataset, mixed-year training mitigates the negative impact of rare classes on the performance.

We assess the impact of mixed-year training by considering I = 3 specialized models whose training set is restricted to a given year: M 2018 , M 2019 , and M 2020 . In contrast, the model M mixed is trained with all parcels across all years with no information regarding of the year of acquisition. All models share the same PSE+LTAE configuration [START_REF] Sainte | Lightweight Temporal Self-Attention for Classifying Satellite Images Time Series[END_REF]. We visualize the training protocols in Figure 5, and report the results in Table 2.

Cross-validation:

We split our data into 5 folds for cross-validation. For each fold, we train on 3 folds and use the last fold for calibration and model selection, corresponding to a train/validation/test ratio of 60%, 20% and 20% in each fold. In order to avoid data contamination and self-correlation, our folds are all spatially separated: the fold separation is done parcel-wise and not for yearly observations. A parcel cannot appear in multiple folds for different years.

Evaluation Metrics

In order to assess the performance of the different approaches evaluated, we report the Overall Accuracy (OA), corresponding to the rate of correct prediction. If we denote by N c the number of accurate predictions and N the total number of parcels, the overall accuracy writes:

OA = N c N . (5)
To address the high class imbalance, we also report the mean Intersection over Union (mIoU), defined as the unweighted class-wise average of the intersection over Union (or Jaccard distance) between the prediction and the ground truth for each class. For a given class i, IoU i is defined as the ratio between the number of elements that are both predicted and labeled by class i (the intersection, or true positives, and the number of elements that are either predicted or labeled as belonging to class i (the union). In terms of binary classification (class i vs. not class i), this translates into the following formula:

IoU i = TP i TP i + FP i + FN i , (6)
with TP i , FP i and FN i the number of true positives, false positives, and false negatives respectively. The mIoU represents the average of the IoU calculated over the K studied classes:

mIoU = 1 K K i=1 IoU i .is (7)

Results

This section presents the quantitative and qualitative impact of our design choice in terms of training protocol and architecture.

Training Protocol

Predictably, the specialized models have good performance when evaluated on a test set composed of parcels from the year they were trained, and poor results for other years, making this training procedure ill-fitted for the application at hand. In contrast, the model M mixed vastly outperformed specialized models on average over the three considered years: over 15 points of mIoU. More surprisingly, the model M mixed also outperforms all specialized models even when evaluated on the year of their training set. This implies that the increased diversity of the mixed-year training set allows the model to learn more robust and expressive representations.

In Figure 6, we illustrate the representations learned by the mixed model M mixed and the specialized model M 2020 . We remark that the parcel embeddings of the specialized model are inconsistent from one year to another, resulting in a higher overlap between classes. In contrast, the mixed year model M mixed learns year-consistent representations. This results in embedding clusters with large margins between classes, illustrating the ability of the model to learn robust and discriminative SITS embeddings.

Influence of Crop Rotation Modeling

We evaluate all models presented in Section 2.3 and Section 2.4, and provide qualitative illustration in Figure 7 . All models are trained with the mixed-year training protocol and only tested on parcels from 2020 to avoid edge effects affecting the evaluation. We give quantitative cross-validated results in Table 3. Training our model on one fold takes 4 hours, and inference on all parcels takes under 3 minutes (over 500-parcels per second). We observe that our model appreciably improved on the single-year model, with over 6 points gained in mIoU. The CRF models also increase the results to a lesser margin. We attribute this inferior performance to an oversmoothing phenomenon already pointed out by Bailly et regardless of the specificity of the observation. In contrast, our approach simultaneously models the current year's observations and the influence of past cultivated crops. M obs barely improves the quality of the single-year model. While this model has indeed access to more information than M single , the same model is used to extract SITS descriptors for all three years. This means that the model's ambiguities and errors will be the same for all three representations, which prevent M obs from largely improving its prediction. Our approach injects new information into the model by concatenating the labels of previous years, which is independent of the model's limitations. Our method is more susceptible to the propagation of mistakes in the farmers' declarations but provides the most significant increase in performance in practice. Lastly, we concatenate both past label vectors to keep information about the order in which past crops were cultivated and observe a slight decrease in performance. The increase in model size can explain this. We conclude that this order is not crucial information for our model conditionally to the observation of the target year. Lastly, the model's performance with only the declaration of the last year performs almost as well as our model with two years worth of crop declarations. This suggests that yearly transition rules are sufficient to capture most inter-year dynamics, such as permanent culture. Alternatively, our two-year scheme may suffer from sharp edge effects with only three years' worth of data. Only a quantitative analysis over a longer period may resolve this ambiguity. On average, our M dec model obtains an mIoU of 84.7% and overall accuracy of 98.1% on the training set.

We report the confusion matrix of M dec in Figure 8, and its performance for each crop in Table 4. We also compute ∆ = IoU(M dec ) -IoU(M single ) the gain compared to the singleyear model IoU(M single ), as well as the ratio of improvement ρ = ∆/(1 -mIoU(M single )). This last number indicates the proportion of IoU that have been gained by modeling crop rotations. We observe that our model provides a large performance increase across all classes but four. The improvement is particularly stark for classes with strong temporal stability such as vineyards. To further this analysis, we arrange the crop types into three groups according to the crop grown in 2018 and the number of observed class successions over the 2018-2019-2020 period:

• Permanent Culture. Classes within this group are such that at least 90% of the • Structured Culture. A crop is said to be structured if, when grown in 2018, over 75% of the observed three-year successions fall into 10 different rotations or less and is not permanent. It contains Rapeseed, Sunflower, Soybean, Alfalfa, Leguminous, Flowers/Fruits/vegetables, and Potato.

• Other. All other classes.

We report the unweighted class average for these three groups in Table 5. Predictably, our approach considerably improves the results for permanent cultures. Our model is also able to learn non-trivial rotations as the improvement for structured classes is also noticeable. On average, our method also improves the performance for other nonstructured classes, albeit to a lesser degree. This indicates that our model can learn multi-year patterns not easily captured by simple rotation statistics.

Model Calibration

Crop mapping can be used for various downstream applications, such as environmental monitoring, subsidy allocation, and price prediction. These applications carry crucial economic and ecological stakes and hence benefit from properly calibrated prediction. A prediction is said to be calibrated when the confidence (ie, the probability associated to a given class) of the prediction corresponds to the empirical rate of correct prediction: we want 90% of the prediction with a 90% confidence to be accurate. This allows for more precise risk estimation and improves control on the rate of false positives / negatives. Deep learning methods such as ours are notoriously badly calibrated. However, this can be corrected with the simple technique proposed by Guo et al. [START_REF] Guo | On calibration of modern neural networks[END_REF]. This method, called temperature scaling, consists in minimizing the discrepancy between the predicted confidence and the rate of errors binned into quantiles (we chose here 15 bins) on the validation set by adjusting the temperature parameters in the last softmax layer [6, Chap. 2.4]. As represented in Figure 9, we can improve the calibration and observe a 43% decrease of the Expected Calibration Error (ECE) at a small computation cost.

Discussion

In this paper, we set out to develop a deep learning method to leverage both the inter-and intra-annual dynamics of crop growth for crop mapping. We propose to enrich the learned spatio-temporal features with the last two declared cultures. Our experiments show that this simple method leads to an appreciable increase in performance compared to models operating on data drawn from a single year. Our method outperforms other approaches We quantize the predicted confidence into 100 bins for visualization purposes. For a perfectly calibrated prediction, the blue histogram would exactly follow the orange line. We observe that a simple post-processing step can considerably improve calibration.

such as CRF smoothing or observation stacking. This improvement can be observed for most crop types, including those with rotation patterns beyond permanent culture. We now discuss the limitations of our method and of our analysis.

Choice of Backbone Network.

Our method can be adapted to any network with a distinct classifier module mapping a spatio-temporal learned feature vector to a predicted vector of class scores. However, the choice of spatio-temporal encoder (backbone) is out of the scope of this article. While it may be relevant to explore the effect of our modification on other architectures, we limited our investigation to the PSE+LTAE as it is the current state-of-the-art network for crop type mapping by a large margin.

Operational Setting. We showed that training our model with samples from all available years leads to considerably improved results. However, this scenario is not compatible with the operational setting of crop monitoring, in which payment agencies may want to detect erroneous declarations before all farmers' declarations have been received. Instead, we use the same setting as the vast majority of work in parcel classification and whose task is to classify parcels after the year is over [START_REF] Zheng | A support vector machine to identify irrigated crop types using time-series Landsat NDVI data[END_REF][START_REF] Vuolo | How much does multi-temporal Sentinel-2 data improve crop type classification?[END_REF][START_REF] Siachalou | A hidden Markov models approach for crop classification: Linking crop phenology to time series of multi-sensor remote sensing data[END_REF][START_REF] Belgiu | Sentinel-2 cropland mapping using pixel-based and object-based time-weighted dynamic time warping analysis[END_REF][START_REF] Kussul | Deep learning classification of land cover and crop types using remote sensing data[END_REF][START_REF] Rußwurm | Self-attention for raw optical satellite time series classification[END_REF][START_REF] Pelletier | Deep learning for the classification of Sentinel-2 image time series[END_REF][START_REF] Sainte | Time-Space tradeoff in deep learning models for crop classification on satellite multispectral image time series[END_REF][START_REF] Rußwurm | Multi-temporal land cover classification with sequential recurrent encoders[END_REF][START_REF] Yuan | Self-Supervised Pretraining of Transformers for Satellite Image Time Series Classification[END_REF][START_REF] Kondmann | DENETHOR: The DynamicEarthNET dataset for Harmonized, inter-Operable, analysis-Ready, daily crop monitoring from space[END_REF][START_REF] Sainte | Lightweight Temporal Self-Attention for Classifying Satellite Images Time Series[END_REF][START_REF] Sainte | Satellite image time series classification with pixel-set encoders and temporal selfattention[END_REF][START_REF] Schneider | Re] Satellite Image Time Series Classification with Pixel-Set Encoders and Temporal Self-Attention[END_REF][START_REF] Sainte | Panoptic Segmentation of Satellite Image Time Series with Convolutional Temporal Attention Networks[END_REF].

As the Sentinel-2 mission was only operational starting in 2017, we only have access to full-year coverage since 2018. This means that we only have three years' worth of data at the time of writing this paper. In our opinion, this prevents us from a realistic setting in which the last year is withheld from the training set. Indeed, the inter-year meteorological variations between two years are typically too great to test for a third year and reasonably expect good results, as corroborated with preliminary experiments not shown in this paper. As more Sentinel-2 data become available, we will be able to evaluate our approach in a more realistic setting.

Scope of the Study. Given the large amount of data involved and the complexity of data collection, we have limited our analysis and our proposed open-access dataset to a single area of the French Metropolitan territory. While nothing in our method is specific to this area, some of our analysis may be biased by the preponderance of stable cultures such as vineyards in this area. To confirm the generality of our conclusions, we would require a dataset with parcels taken from regions across the world with various meteorological conditions and agricultural practices. This task is complicated by the lack of harmonization between LPIS regarding open-access policy and even nomenclature. We hope that our results will encourage mapping agencies worldwide to release multi-year LPIS in open-source to help constitute a truly global dataset, allowing the community to assess the spatial generalizability of state-of-the-art methods. We also limit ourselves to predicting the main culture in each parcel while ignoring cases with multiple growth cycles within one year. This may be particularly detrimental to its application in subtropical regions.

Applicability of our Model. By requiring the last two grown crops to classify a parcel, our method cannot be applied to areas where the LPIS is not easily accessible. Furthermore, our training setting requires only selecting stable parcels. This can be easily obtained from the LPIS if it also contains information about the extent and position of each parcel, as is the case for the French LPIS. As a consequence, the effect of parcels with changing contours is out of the scope of our investigation.

Conclusions

We explored the impact of using multi-year data to improve the quality of the automatic classification of parcels from satellite image time series. We showed that training a deep learning model from multi-year observations improved its ability to generalize and better precision across the board. We proposed a simple modification to a state-of-the-art network to model both inter-and intra-year dynamics. This resulted in an increase of +6.3% of mIoU when compared to models operating on single-year data. The effect is most substantial for classes with strong temporal structures, but also impacts other crop types. We also showed how simple post-processing could improve the calibration of the models considered.

Finally, we release both our code and our data. We hope that our promising results will encourage the SITS community to develop methods modeling multiple time scales simultaneously and release more datasets spanning several years.

Figure 2 :

 2 Figure 2: Area of Interest. The studied parcels are taken from the 31TFM Sentinel-2 tile, covering 110 × 110 km and containing over 103 602 parcels meeting our size, shape, and stability criteria. (a) Large view of the tile. (b) Detail of the area.

Figure 3 :

 3 Figure 3: Intra-Year Dynamics. Evolution of two areas across three seasons of the year 2020. The top parcels contain mainly meadow parcels, while the bottom one comprises more diverse crops. The aspect of most parcel drastically changes across one year's worth of acquisition, corresponding to different phases in the growth cycle. (a) Winter. (b) Spring. (c) Summer

Figure 4 :

 4 Figure 4: Multi-Year Modeling. Different approaches to model crop rotation dynamics: (a) the model only has access to the current year's observation; (b) a chain-CRF is used to model the influence of past cultivated crop; (c) the model has access to the observation of the past two years; (d) proposed approach: the model has access to the last two declared crops.

Figure 5 :

 5 Figure 5: Training Protocol. A single model is trained with parcels taken from all three years (a), and three specialized models whose training set only comprises observation for a given year (b).

  al. : CRFs tend to resolve ambiguities with the most frequent transition

Figure 6 :

 6 Figure 6: Learned Representations. Illustrations of the learned SITS representations of the mixed-year model M mixed (a) and the specialized M 2020 (b).T-SNE algorithm is used to plot in 2D the representation for 100 parcels over 10 classes and 3 years. We observe that M mixed produced cluster of embeddings that are consistent from one year to another, and with clearer demarcation between classes.

Figure 7 :

 7 Figure 7: Qualitative Illustration. Detail of the area of interest with the ground truth in (a) and the qualification of the prediction in (b) with correct prediction in blue and errors in red.

  l e y W i n t e r S u n f l o w e r V i n e y a r d S o y b e a n S o r g h u m A l f a l f a O a t W i n t e r L e g u m i n o u s M i x c e r e a l s F l o . F r u . V e g . O a t S u m m e r P o t a t o B a r l e y S u m m e r W

Figure 8 :

 8 Figure 8: Confusion Matrix. Confusion matrix of the prediction of M dec for the year 2020. The area of each entry corresponds to the square root of the number of predictions.

Figure 9 :

 9 Figure 9: Model calibration. Empirical rate of correct prediction by predicted confidence.We quantize the predicted confidence into 100 bins for visualization purposes. For a perfectly calibrated prediction, the blue histogram would exactly follow the orange line. We observe that a simple post-processing step can considerably improve calibration.

Table 1 :

 1 Crop distribution. We indicate the number of parcels declarations in the LPIS for each class across all 103 602 parcels and all 3 years.

	Class	Count	Class	Count
	Meadow	184 489 Triticale	5114
	Maize	42 006	Rye	569
	Wheat	27 921	Rapeseed	7624
	Barley Winter	10 516	Sunflower	1886
	Vineyard	15 461	Soybean	6072
	Sorghum	820	Alfalfa	2682
	Oat Winter	529	Leguminous	1454
	Mixed cereal	1061	Flo./fru./veg. 1079
	Oat Summer	330	Potato	230
	Barley Summer 538	Wood pasture 425
	descriptor e i :

Table 2 :

 2 Quantitative evaluation. Performance (mIoU and OA) of the different specialized models M 2018 , M 2019 , M 2020 and of the mixed-years model M mixed evaluated on each year individually and all available years simultaneously with 5-fold cross-validation. The best performances are shown in bold. Boxed values correspond to evaluations where the training and evaluation sets are drawn from the same year. The mixed-year model performs better for all years, even compared to specialized models.In the rest of the paper, we use mixed year training for all models.

	Model	2018 OA mIoU	2019 OA mIoU	2020 OA mIoU	3 years OA mIoU
	M 2018	97.0	64.7	90.3	45.5	90.8	43.4	92.7	49.1
	M 2019	88.9	39.5	97.2	70.1	88.7	40.1	91.6	48.0
	M 2020	91.4	44.2	93.7	51.8	96.7	67.3	93.9	54.0
	M mixed	97.3 69.2	97.4 72.2	96.8 68.7	97.2 70.4

Table 3 :

 3 Performance by model. Performances (mIoU and OA) of the models M single , M obs , M CRF , and M dec tested for the year 2020. Our proposed model M dec achieve higher performance than M single with a 6.3% mIoU gap.

	Model	Description	OA	mIoU
	M single	single-year observation	96.8 68.7
	M obs	bypassing 2 years of observation	96.8 69.3
	M CRF	using past 2 declarations in a CRF 96.8 72.3
	M dec-one-year concatenating last declaration only 97.5 74.3
	M dec-concat	concatenating past 2 declarations	97.5 74.4
	M dec	proposed method	97.5 75.0

Table 4 :

 4 Performance by class. IoU per class of our model M dec for the year 2020, as well as the improvement ∆ compared to the single-year model M single , and the ratio of improvement ρ. All values are given in %, and we sort the classes according to decreasing ratios ρ.

	Class	IoU	∆	ρ	Class	IoU	∆	ρ
	Wood Pasture	92.4 +48.2 86.3 Oat Summer	52.8 +3.6 7.0
	Vineyard	99.3 +1.4	68.7 Rapeseed	98.3 +0.1 6.6
	Alfalfa	68.7 +23.9 49.9 Maize	95.7 +0.2 6.3
	Flo./Fru./Veg. 83.4 +14.5 46.5 Wheat	91.9 +0.3 3.9
	Meadow	98.4 +0.9	36.9 Barley Summer 64.3 +1.1 3.1
	Leguminous	45.2 +14.6 21.1 Potato	57.1 +0.5 1.2
	Rye	54.7 +6.4	12.4 Sunflower	92.2 -0.1	-0.3
	Oat Winter	57.7 +4.5	9.7	Sorghum	56.6 -0.2	-0.4
	Triticale	68.7 +2.6	7.8	Soybean	91.8 -0.2	-3.1
	Mix. Cereals	31.0 +5.1	6.8	Barley Winter	92.8 -0.6	-8.5

Table 5 :

 5 Improvement Relative to Structure. Classwise IoU and mean improvement of our model compared to the single-year model according to the rotation structure of the cultivated crops.

	Category	mIoU mean ∆
	Permanent 97.3	16.9
	Structured 77.7	7.6
	Other	66.6	2.3


Charlotte Vajhøj

Benjamin Schmid

Ania Alik

Ronald Melki

Karina Fog

Bjørn Holst

Tina Charlotte Stummann

Establishment of a human induced pluripotent stem cell neuronal model for identification of modulators of A53T αsynuclein levels and aggregation

Inhibiting formation or promoting degradation of α-synuclein aggregates are among the therapeutical approaches under investigation as disease-modifying treatment strategies for Parkinson's disease. To support these developments, several in vitro models based on seeded α-synuclein aggregation have been established in immortalized cell lines and murine primary neurons. Here, we report on a humanized model with a reproducibility and throughput that enables its use in supporting target identification and validation in pharmacological research. A human induced pluripotent stem cell (iPSC) line was genetically modified to express HA-tagged α-synuclein with the point mutation in position 53 from Alanine to Threonine (A53T) under an inducible system and differentiated into cortical neurons expressing neuronal markers and exhibiting spontaneous activity. Intracellular α-synuclein aggregation was triggered by exposure to exogenous added fibrillated recombinant wild-type human αsynuclein fibrils91 and demonstrated by several endpoints; the formation of Triton-insoluble SDS-soluble α-synuclein, biochemically in a fluorescence resonance energy transfer based aggregation assay and by immunocytochemistry of phosphorylated α-synuclein positive puncta. We demonstrate the feasibility of upscaling the iPSC neuron production for drug discovery and that the model has a suitable dynamic range allowing for both detection of increased and decreased α-synuclein aggregation. Moreover, gene modulation is feasible using siRNAs, making the model suitable for genetic screening for modulators of α-synuclein aggregation. Data on effects of USP8, USP13 and USP9X knockdown on α-synuclein expression and aggregation contradicts published data from immortalized cell lines and murine systems. This highlight the importance of including humanized neuronal models in the confirmation of biological mechanisms in specific variations of Parkinson's disease.

Introduction

Targeting α-synuclein is currently one of the most pursued disease-modifying treatment strategies for Parkinson's disease. Toxic α-synuclein species have been shown to disrupt important cellular functions, and inhibiting formation or promoting degradation of α-synuclein aggregates are among the therapeutical approaches under investigation [START_REF] Fields | Targeting Alpha-Synuclein as a Therapy for Parkinson's Disease[END_REF]. The development of cellular models, where addition of exogenous pathological α-synuclein fibrils to immortalized cell lines or murine primary neurons leads to misfolding, aggregation and phosphorylation of endogenous α-synuclein, has provided valuable model systems capturing key pathological phenotypes of Parkinson's disease [START_REF] Luk | PNAS-2009-Luk-20051-6.pdf[END_REF][START_REF] Yamasaki | Parkinson's disease and multiple system atrophy have distinct -synuclein seed characteristics[END_REF][START_REF] Volpicelli-Daley | Exogenous α-synuclein fibrils induce Lewy body pathology leading to synaptic dysfunction and neuron death[END_REF].

Models based on dopaminergic neurons derived from Parkinson's patient iPSCs have demonstrated higher α-synuclein protein levels and increased phosphorylation of α-synuclein at serine residue 129 (pS129) compared to controls [START_REF] Devine | Parkinson's disease induced pluripotent stem cells with triplication of the α-synuclein locus[END_REF][START_REF] Ryan | Isogenic Human iPSC Parkinson's Model Shows Nitrosative Stress-Induced Dysfunction in MEF2-PGC1α Transcription[END_REF]. The limited reports on spontaneous aggregation in these models show data on colocalization of Thioflavin T with α-synuclein and the presence of Triton insoluble α-synuclein [START_REF] Ryan | Isogenic Human iPSC Parkinson's Model Shows Nitrosative Stress-Induced Dysfunction in MEF2-PGC1α Transcription[END_REF][START_REF] Mazzulli | α-Synuclein-induced lysosomal dysfunction occurs through disruptions in protein trafficking in human midbrain synucleinopathy models[END_REF][START_REF] Brazdis | Demonstration of brain region-specific neuronal vulnerability in human iPSC-based model of familial Parkinson's disease[END_REF][START_REF] Prots | α-Synuclein oligomers induce early axonal dysfunction in human iPSC-based models of synucleinopathies[END_REF]. The robustness of the findings remains to be confirmed.

Recently, recombinant α-synuclein fibrils have been demonstrated to increase endogenous pS129 α-synuclein and aggregation in human iPSC cortical neurons from healthy individuals [START_REF] Gribaudo | Propagation of α-Synuclein Strains within Human Reconstructed Neuronal Network[END_REF][START_REF] Bieri | LRRK2 modifies α-syn pathology and spread in mouse models and human neurons[END_REF]. Aggregation of α-synuclein in the iPSC models was shown by assessment of Triton insoluble α-synuclein by fractionation, an assay with a very low throughput, by the use of Thioflavin T, exhibiting enhanced fluorescence upon binding to amyloid fibrils and hence is not specific for α-synuclein, and by immunostaining for pS129 α-synuclein as proxy maker for aggregation. Assays quantifying the α-synuclein aggregation directly would be preferable in the search for therapeutical approaches modulating this feature of Parkinson's disease. In this paper, we take advantage of gene-editing to establish an iPSC lines with doxycycline inducible HA-tagged α-synuclein A53T expression, exogenous fibrils to seed intracellular α-synuclein aggregation and a fluorescence resonance energy transfer (FRET)-based α-synuclein aggregation assay. Together, this enabled us to develop a robust and reliable human neuronal model suitable for drug discovery genetic screens for modulators of α-synuclein aggregation.

Materials and methods

iPSC generation and culturing

The derivation of the iPSC line BIONi010-C from an apparently healthy 18 years old male has previously been described [START_REF] Rasmussen | Transient p53 suppression increases reprogramming of human fibroblasts without affecting apoptosis and DNA damage[END_REF]. Gene-editing was carried out using TALEN technology to generate the BIONi010-C-24 iPSC line used in the present study. The cDNA sequence of the SNCA construct with A53T mutation and a single HA-tag (ATGGATGTATTCATGAAAGGA CTTTCAAAGGCCAAGGAGGGAGTTGTGGCTGCTGCTGAGAAAACCAAACAGGGTGTGGCAG AAGCAGCAGGAAAGACAAAAGAGGGTGTTCTCTATGTAGGCTCCAAAACCAAGGAGGGAGT GGTGCATGGTGTGACAACAGTGGCTGAGAAGACCAAAGAGCAAGTGACAAATGTTGGAGGA GCAGTGGTGACGGGTGTGACAGCAGTAGCCCAGAAGACAGTGGAGGGAGCAGGGAGCATTG CAGCAGCCACTGGCTTTGTCAAAAAGGACCAGTTGGGCAAGAATGAAGAAGGAGCCCCACA GGAAGGAATTCTGGAAGATATGCCTGTGGATCCTGACAATGAGGCTTATGAAATGCCTTCT GAGGAAGGGTATCAAGACTACGAACCTGAAGCCTACCCATACGATGTTCCAGATTACGCT TAA, A53T mutation in bold and underlined, HA-tag in italic and underlined) was designed with an N-terminal NcoI site and a C-terminal KpnI site. A plasmid, containing the two homology arms, a neomycin (Neo) resistance gene and a tetracycline/doxycycline responsive element, was used as backbone vector for the SNCA-HA construct (TRE-TIGHT-EGFP-backward donor, a gift from Rudolf Jaenisch corresponding to Addgene plasmid # 22077; http:// n2t.net/addgene:22077; RRID:Addgene_22077). The SNCA-HA cDNA was cloned into this vector using the two mentioned restriction enzymes. To make the line inducible for doxycycline, a plasmid with the reverse TET transactivator (m2rtTA) was used (AAVS1-Neo-M2rtTA, a gift from Rudolf Jaenisch corresponding to Addgene plasmid # 60843; http://n2t. net/addgene:60843; RRID:Addgene_60843). The m2rtTA plasmid contained a puromycine (Puro) resistance gene enabling double selection for Neo and Puro. To target the AAVS1 locus, two TALENs were used (AAVS1-TALEN-L, a gift from Danwei Huangfu corresponding to Addgene plasmid # 59025; http://n2t.net/addgene:59025; RRID:Addgene_59025 and AAV-S1-TALEN-R, a gift from Danwei Huangfu corresponding to Addgene plasmid # 59026; http://n2t.net/addgene:59026; RRID:Addgene_59026).

The iPSCs were cultured in E8 medium on Matrigel (Corning Bioscience) coated 6w-dishes and detached using Accutase. When they reached a density of 70-90%, a total of 1.5 million cells were co-nucleofected with 2 μg of both TALEN plasmids, 2 μg of the plasmid containing the m2rtTA and 2 μg of the SNCA-HA plasmid. For the nucleofection, the P3 Primary Cell Kit (Lonza) was used and run with the program CA167 following to the manufacturer's instructions. The nuclofected iPSCs were transferred back to a Matrigel-coated 100 mm dish in E8 medium supplemented with 1:200 diluted Revita cell supplement (Gibco). 24 h post-nucleofection, cells were subjected to double selection with 1 μg/mL puromycin (Invitrogen) for 2 days and 250 μg/mL neomycin (Gibco) for 4 days. After selection, cells were allowed to recover for one week. Resistant colonies were then picked and expanded. DNA from the clones was extracted and a PCR was run to check for correct integration of the plasmids into the AAVS1 locus using primer pairs specific for the M2rtTA-plasmid (ccagaggcggcgcagaagccag and atcttgttcaatggccgatcccat) and the SNCA-HA-plasmid (ccagaggcggcgcag aagccag and ccgtgggcttgtactcggtcat). The BIONi010-C-24 lines is commercially available from the EBiSC collection (European Collection of Authenticated Cell Cultures).

The BIONi010-C-24 iPSCs were cultured in Essential 8 (Gibco A15169-01) media on Matrigel (BD 354277) coat 6w-dishes and split with EDTA (Ambion AM9206G). We generated mother stock vials frozen in nitrogen tanks. From each of these, we generate a large bank of iPSC vials frozen at 140˚C to be used for starting new differentiation. In practice, this enables differentiating to always start from the same passage number.

α-synuclein fibrillar assemblies

Recombinant wild-type human α-synuclein was purified and assembled into the fibrillar polymorphs fibrils91 as previously described [START_REF] Ghee | PA700, the regulatory complex of the 26S proteasome, interferes with α-synuclein assembly[END_REF][START_REF] Makky | Nanomechanical properties of distinct fibrillar polymorphs of the protein α-synuclein[END_REF][START_REF] Rey | α-Synuclein conformational strains spread, seed and target neuronal cells differentially after injection into the olfactory bulb[END_REF]. Assembly was monitored using Thioflavin T as previously described [START_REF] Bousset | Structural and functional characterization of two alpha-synuclein strains[END_REF]. Monomeric alpha-synuclein concentration was determined spectrophotometrically using the extinction coefficients 5960 M -1 �cm -1 at 280 nm. Fibrillar assembly concentration was determined by subtracting the concentration of α-synuclein in the supernatant fraction from the initial protein concentration after centrifugation at 100,000 g for 30 min. The resuspended fibrils91 (in PBS) were fragmented prior to addition to cell cultures by sonication for 20 min in 2 mL Eppendorf tubes using a Vial Tweeter powered by ultrasonic processor UIS250v (250W, 2.4kHz; Hielscher Ultrasonic, Teltow, Germany), aliquoted, flash frozen in liquid nitrogen and stored until use at -80˚C. Quality control was performed before and after: Fibrils91 were imaged by TEM after negative staining with Uranyl acetate using a Jeol 1400 transmission electron microscope after adsorption on 200 mesh, carbon coated, electron micrograph grids. The images were recorded using a Gatan Orius CCD camera (Gatan Inc. Pleasanton, CA, USA). Fibrils91 (1.4 mg/ml equivalent monomer concentration) were also digested at 37˚C in PBS with Proteinase K (3.8 ug/ml) (Roche). Aliquots were removed at different time intervals following addition of the protease (0, 1, 5, 15, 30 and 60 minutes) and transferred into Eppendorf tubes with proteinase K inhibitors (phenylmethylsulfonyl fluoride). The samples where dried using speed vacuum and further solubilized by addition of pure Hexafluoroisopropanol (HFIP). After evaporation of HFIP, the samples were resuspended in laemmli buffer, heated 10 minutes at 70˚C and analyzed by Tris-Glycine SDS-PAGE (15%) after Coomassie blue staining.

iPSC neuronal differentiation, α-synuclein seeding and siRNA exposure A modified version of a previously published neuronal differentiation protocol was used [START_REF] Shi | Human cerebral cortex development from pluripotent stem cells to functional excitatory synapses[END_REF]. Differentiation was initiated day 0 by dual SMAD inhibition (100 nM LDN193189 (Sigma-Aldrich SML0559) and 10 μM SB431542 (Sigma-Aldrich S4317)) in N3 media (50% DMEM-F12 (Thermo Fisher Scientific 31331-028), 50% Neurobasal medium (Thermo Fisher Scientific 21103-049), 0.5% N2 (Thermo Fisher Scientific, 17502-048), 0.5% B27 (Thermo Fisher Scientific 17504-044), 0.5 mM GlutaMAX supplement (Thermo Fisher Scientific 35050-061), 0.5% NEAA supplement (Thermo Fisher Scientific 11140-050), 50 μM 2-mercaptoethanol (Thermo Fisher Scientific 21985023), 2.5 μg/mL insulin (Sigma-Aldrich 19278)) of 500,000 iPSCs per cm 2 plated on Matrigel coated dishes. The cells were split at day 7, 12 and 17 and re-plated as on day 0, except that from day 12 and onwards poly-L-ornithine/laminin coat was used and that LDN193189 and SB431542 were withdrawn from day 13. The neuronal progenitor cells were frozen day 21 to generate a neuronal progenitor bank. Upon thawing the neuronal progenitors were allowed to proliferate for 4 days, plated day 25 and subjecting to neuronal maturation medium (N3 with 20 ng/mL BDNF (R&D Systems 248-BD), 10 ng/mL GDNF (R&D Systems 212-GD), 200 μM L-Ascorbic Acid (Sigma-Aldrich A5960), and 500 μM db-cAMP (Sigma-Aldrich D0627) from day 26 and onwards. Medium containing maturation factors was used within a week and changed every 2-3 days. Day 32, the cells were re-plated into assay plates at a density of 167,000 cells/cm 2 , except that 333,000 cells/cm 2 were used for measurement of spontaneous calcium oscillation. 96 well formats were used for all assays except that 6 wells were used for fractionation and western blot samples. At all timepoints of differentiation, cells were split using accutase (Thermo Fisher Scientific A11105) and 10 μM ROCK inhibitor (Sigma-Aldrich Y0503) was added overnight upon thawing and splitting. On day 33, doxycycline (1 μg/mL unless otherwise indicated, Stemgent #04-0016) was added to the cultures to induce expression of HA-tagged A53T α-synuclein. Horizon Discovery Accell siRNAs (1 μM) were added day 36 and 43 and left on for 72h before the media was refreshed. Pools of four different siRNAs were used as non-targeting control (D-001910-10) and to target USP8 (E-005203-00), USP13 (E-006064-00) and USP9X (E-006099-00). Custom made siRNAs were used for SCNA (sense sequences for SCNA siRNA_1, 2 and 3 were GAUGUAUUCAUGAAA GGACUU, GGGUGUUCUCUAUGUAGGCUU and GAGCAAGUGACAAAUGUUGUU). On day 46, fragmented recombinant full-length wild-type human α-synuclein fibrils91 were thawed and immediately added to the iPSC neurons. The fibrils were left on for 96h before changing the media. Cells were harvest for analysis on day 60 unless indicated otherwise.

Real-time quantitative polymerase chain reaction (RT-qPCR)

SYBR Green Cell-to-Ct kit (Ambion 4402957) was used for RNA isolation and reverse transcription to cDNA. RT-qPCR was performed with TaqMan 1 Fast Advanced Master Mix and TaqMan primers purchased from Thermo Scientific. The 2 -ΔΔCt method was applied to expression data and genes of interests (POU5F1 Hs04260367_gH, TUBB3 Hs00801390_s1, MAPT Hs00902190_m1, MAP2 Hs00258900_m1, RBFOX3 Hs01370653_m1, CAMK2A Hs00947041_m1, SLC17A7 Hs00220404_m1, SLC17A6 Hs00220439_m1, GAD1 Hs01065892_m1, SLC6A1 Hs01104475_m1, SLC6A11 Hs01117194_m1, SLC32A1 Hs00369773_m1, SNCA Hs01103383_m1) were normalized to geometric mean of housekeeping genes (GUSB Hs00939627_m1, RPL13A Hs04194366_g1, HSP90AB1 Hs04194340_g1) and mean of undifferentiated iPSC cDNA samples.

Immunocytochemistry

Cells were fixed in 4% paraformaldehyde for 15 min at room temperature, blocked 20 min in KPBS (10 mM phosphate buffer, 150 mM NaCl, and 3.4 mM KCl) supplemented with 0.5% BSA, 0.1% Triton X-100, and 5% normal swine serum (Jackson ImmunoResearch, 014-000-121) and incubated overnight at 4˚C in blocking buffer containing primary antibody (Oct-4 1:500 Stemgent 09-0023, Nestin 1:5000 Millipore MAB5326, Pax6 1:1000 BioLegend 901301, Ki67 1:200 Millipore MAB4190, βIII-Tubulin 1:200 GeneTex GTX631836, total tau 1:1000 Dako A0024, synapsin-I 1:200 Abcam ab64581, GFAP 1:2000 Dako z0334, GABA 1:1000 Sigma A2052, Doublecortin 1:200 Abcam ab18723, pS129 α-synuclein 1:500 Abcam ab51253, HA-tag 1:1000 Santa Cruz sc-7392, MAP2 1:5000 Abcam ab5392). Subsequently, they were incubated with secondary antibodies (Alexa Fluor Secondary Antibodies from Invitrogen) and Hoechst 34580 (BDbiosciences 565877) for 1 hour at room temperature, mounted in Prolong Gold Antifade Mountant (Thermo Scientific, P36934), and visualized on Leica DM5500 B fluorescent microscope (Leica, DE) or using a Nikon Eclipse Ti 4.10 microscope combined with a Yokogawa CSU-X1 spinning disk confocal scanner and a Hamamatsu UltraVIEW VoX C9100-50 EMCCD camera.

Intracellular calcium measurements

Neurons were incubated for 1h with the fluorescent calcium indicator FLIPR Calcium 5 (Molecular Devices R8186) diluted 3x more than providers recommendation in Hanks balanced salt solution (ThermoFisher Scientific, 14175) with 20 mM Hepes and 1.26 mM CaCl2. Image based measurements of cytosolic calcium kinetics were conducted at 37˚C and with a 1 Hz sampling frequency (excitation and emission wavelengths of 480, and 540 nm, respectively) using the FDSS7000 system (Hamamatsu, Japan). The instrument readout is the total fluorescent intensity per well and data are reported as ratio to the first sampled data point.

Fractionation and western blotting

Whole cell lysates were generated by lysis for 30 min in FNN0011 buffer (ThermoFisher Scientific) supplemented with protease and phosphatase inhibitor (Roche, 11697498001, 04906837001). Lysate was centrifuged (18.000xg, 20 min, 4˚C), and the supernatant was collected and stored at -80˚C until western blotting. Fractionation samples were generated by lysing the cells from a 6-well in 100μl 1% Triton X-100 in 50 mM Tris, 150 mM NaCl (pH 7.6) with protease and phosphatase inhibitor (Roche, 11697498001, 04906837001) for 15 min on ice followed by sonication using a Branson Sonifier 250. Ultracentrifuge at 100.000g for 30 min at 4 ˚C resulted in the "Triton fraction" supernatant. The pellet was washed once by adding and removing buffer again without dissolving pellet, resuspend in 100μl SDS buffer (1% SDS in 50 mM Tris, 150 mM NaCl (pH 7.6) with protease and phosphatase inhibitor), sonicated and ultra-centrifuged at 100.000g for 30 min at RT. This supernatant was the "SDS fraction". Both supernatants were frozen at -80˚C until further processing. 15 μl of the fractions were loaded per well of the polyacrylamide gels. All western blots were performed with standard methods. Recombinant tau protein ladder was included in some blots (Sigma T7951). Nitrocellulose membranes were placed in 98˚C PBS for 5 min before blotting with α-synuclein antibodies for full epitope retrieval. Primary antibodies were βIII-Tubulin 1:400 GeneTex GTX631836, total tau 1:1000 Dako A0024, synapsin-I 1:500 Abcam ab64581, β-actin 1:50.000 Sigma A5441, GAPDH 1:1000 Abcam ab9485, pS129 α-synuclein 1:2000 Abcam ab51253, HA-tag Abcam ab9110, total α-synuclein 1:5000 Thermo Fisher MA1-90346-4B12, USP13 1:500 Proteintech 16840-1-AP, USP8 1:250 Cell Signaling 11832, USP9X 1:500 Cell Signaling 365353. Secondary antibodies were Alexa Flour 680 from Invitrogen or 800CW IRDye from LI-COR. Blots were scanned on a LI-COR Odyssey CLx Instrument and data extracted using LI-COR Imaging Studio software.

pS129-synuclein and pS129-pS129 aggregation α-synuclein FRET assays

Neurons were washed twice in PBS, before sucking dry and freezing at -80˚C until further analysis of the samples. Analysis was performed with commercially available time-resolved fluorescence resonance energy transfer (FRET) kits from Cisbio and a PHERAstar Microplate Reader (BMG LABTECH). One kit was the pS129 α-synuclein assay that is based on a total αsynuclein antibody and a pS129 α-synuclein antibody (Cisbio #6FSYNPEG), detecting only phosphorylated α-synuclein. The other was the pS129-pS129 aggregation α-synuclein assay that is based on one pS129-synuclein antibody labelled either with FRET donor or with acceptor (Cisbio #6FCUS000, please contact Delphine.Jaga@PERKINELMER.COM to get access to this kit), detecting only phosphorylated oligomers or larger α-synuclein assemblies.

Cell viability assays and total α-synuclein quantification

Cell viability was determined according to the manufacture's instruction using CellTiter-Glo Luminescent Cell viability assay quantifying the cellular ATP content (Promega G7570). Alternative, resazurin metabolism was used to determine viability by incubating the neurons for 1h at 37˚C in neuronal maturation medium containing 10μg/ml resazurin sodium salt (Sigma R7017) and subsequently measuring the fluorescence with excitation wavelength at 530-560nm and emission wavelength at 590nm using an Envision 2103 Multilabel Reader. Total αsynuclein protein quantification was determined according to the manufacture's instruction using human α-synuclein electrochemiluminescence immunoassay kit (Meso Scale Diagnostics K151TGD).

Curve fitting and statistical analysis

Arithmetic was conducted in Excel (Microsoft) and curve fitting and statistical analysis in Prism GraphPad. Graphs are shown as mean±SEM ("n" = independent neuronal differentiations, "m" = technical replicates).

Results

Establishment of a doxycycline inducible A53T α-synuclein expressing iPSC line and its differentiation into spontaneously active neurons

We generated the iPSC line BIONi010-C-24 with doxycycline inducible expression of A53T αsynuclein carrying the disease-causing mutations A53T and a HA-tag. We decided for A53T α-synuclein, which is more aggregation prone than wildtype α-synuclein [START_REF] Conway | Accelerated in vitro fibril formation by a mutant α-synuclein linked to early-onset Parkinson disease[END_REF], to support the establishment of an α-synuclein seeding with a robust assay window. The HA-tag was added to be able to differentiate endogenously formed fibrils from the seed material.

De novo generated fibrils91 seed α-synuclein aggregation in iPSC neurons expressing A53T α-synuclein

De novo generated recombinant full-length wild-type human α-synuclein fibrils91 were quality controlled and fragmented (S7 Fig) and added to human iPSC neurons expressing either endogenous wild-type α-synuclein (i.e. neurons not treated with doxycycline) or endogenous wild-type α-synuclein plus overexpressed HA-tagged human A53T α-synuclein (i.e. neurons RBFOX3 (NeuN), glutamatergic markers CAMK2A, SLC17A7, SLC17A6 and GABAergic markers GAD1, SLC6A1, SLC6A11 and SLC32A1 in iPSC neurons (day 60 of differentiation) versus undifferentiated iPSCs (day 0). Data are mean±SEM (n = 3, m = 4) and analyzed by Students unpaired t-test, ns. non-significant, � p<0.05, �� p<0.01, ��� p<0.001 and ���� p<0.0001. (F) Changes in intracellular calcium in iPSC neurons (day 52/53 of differentiation) following application of the neurotransmitters, glutamate (300 μM glutamate +10 μM glycine), NMDA (40 μM) and GABA (100 μM), or depolarization of the membrane potential by addition of extracellular potassium chloride (25 mM). Representative traces are shown. (G) Spontaneous calcium oscillations in iPSC neurons (day 59/60 of differentiation) were eliminated by addition (timepoint marked with triangle) of 1 μM tetrodotoxin (grey traces) but not buffer (black traces). Representative traces are shown. Exposing the doxycycline treated neurons to fibrils91 and harvesting them 1 week (day 53) and 2 weeks (day 60) after seed addition resulted in a fibrils91-concentration and time dependent increase in the FRET signal in the pS129 α-synuclein assay (Fig 2D). Analysing the samples in the pS129 α-synuclein assay as well as the pS129-pS129 α-synuclein aggregation assay resulted in similar curves indicating that the two assays are interchangeable (Fig 2E). Seeding was only observed in doxycycline treated cells demonstrating that the exogenous fibrils91 did not contributing to the signal in the assay, and that the endogenous wild-type α-synuclein alone was not enough to drive the templated aggregation in the cells.

To confirm the presence of seeded α-synuclein aggregation with an alternative method, we performed fractionation by ultra-centrifugation separating the cell lysates of neurons exposed to PBS, 0.1 μg/mL fibrils91 or 1 μg/mL fibrils91 into Triton-soluble and Triton-insoluble SDSsoluble fractions. Western blots for total α-synuclein and pS129 α-synuclein in the Triton-soluble fraction from doxycycline treated neurons revealed a 14 kDa and a 16 kDa band corresponding to full length wild-type endogenous and HA-tagged A53T α-synuclein, respectively (Fig 2F and S1 Raw images). Western blotting for pS129 α-synuclein in doxycycline treated neurons showed a fibrils91-concentration dependent increase in the Triton-insoluble SDS-soluble fractions (Fig 2F and S1 Raw images), which was absent in non-doxycycline treated neurons (S1 Raw images). This showed that the increase was not due to phosphorylation of the fibrils91 or formation of phosphorylated endogenous wild-type α-synuclein aggregates. Aggregation of HA-tagged A53T α-synuclein was confirmed by HA-tag blotting showing a fibrils91 concentration dependent increase in HA-band intensities in the Triton-insoluble SDS-soluble fraction in doxycycline treated neurons (Fig 2F and S1 Raw images). Western blots for total αsynuclein in the Triton-insoluble SDS-soluble fraction from doxycycline treated neurons further confirmed a fibrils91 concentration dependent increase in α-synuclein intensities (S1 Raw images), although some of signal comes from the fibrils91 as evident in the Triton-insoluble SDS-soluble fraction from fibrils91 exposed non-doxycycline treated neurons (S1 Raw images).

Based on the data from the fractionation and FRET-based α-synuclein assays, pS129 αsynuclein can be used as a proxy marker for α-synuclein aggregation. Immunocytochemistry for this epitope was therefore used to confirm that the templated aggregation was in neurons (Fig 2G). In line with the data from the fractionation and FRET-based α-synuclein assays, no pS129 α-synuclein signal was detectable in non-doxycycline treated neurons or neurons not exposed to fibrils91.

A53T α-synuclein seeding can be modulated

We aimed at generating an iPSC model of α-synuclein aggregation suitable for target validation and screening for modulators of α-synuclein aggregation. It was therefore crucial to demonstrate that the model has a suitable dynamic range and that gene modulation is feasible. . This demonstrates that our standard experimental paradigm, using 1 μg/mL doxycycline, allows for both detection of increased and decreased aggregation.

To address the feasibility of gene modulation, we showed that pS129 α-synuclein aggregation could be reduced by addition of three different SCNA siRNAs but not by non-targeting siRNAs (Fig 3D). RT-qPCR gene expression analysis confirmed knockdown of SCNA mRNA (Fig 3E) and western blotting confirmed reduction of α-synuclein expression by SCNA siR-NAs and no effect of non-targeting siRNAs (Fig 3F and 3G and S1 Raw images). There was correlation between α-synuclein mRNA knockdown, α-synuclein protein knockdown and αsynuclein aggregation (Fig 3D , 3E and3G).

Enzymes of the ubiquitination system affect A53T α-synuclein aggregation

USP8, USP13 and USP9X are all enzymes of the ubiquitination system that are reported to modulate α-synuclein expression or aggregation [START_REF] Alexopoulou | Deubiquitinase Usp8 regulates α-synuclein clearance and modifies its toxicity in Lewy body disease[END_REF][START_REF] Liu | Ubiquitin specific protease-13 independently regulates parkin ubiquitination and alpha-synuclein clearance in alpha-synucleinopathies[END_REF][START_REF] Rott | α-Synuclein fate is determined by USP9X-regulated monoubiquitination[END_REF]. These genes were therefore knocked down using pools of four siRNAs for each gene to investigate if our model could capture previous findings. Quantifications of western blots for the enzymes showed statistically significant reduction of USP13, USP8 and USP9X protein levels upon siRNA treatment (Fig 4A

Discussion

We aimed at generating an iPSC model to investigate α-synuclein aggregation in a robust and reproducible system with a good throughput. The first step to reach that goal was to generate the gene-editing iPSC BIONi010-C-24 line with doxycycline inducible HA-tagged α-synuclein A53T expression. Differentiation of 5 million iPSCs resulted in 170±70 million neuronal progenitors, which could be cryopreserved. Thawing and further differentiation provided enough cells for seeding 10,200 x 96-wells (50.000 cells/well day 32), demonstrating the feasibility of upscaling the iPSC neuron production for drug discovery.

Three neuronal progenitor batches were generated, which could be matured to neurons with high reproducibility. Upon maturation, the majority of the cells were of neuronal identity exhibiting clear neuronal morphology and expression at protein and mRNA levels of several pan-neuronal markers. Although the expressed tau isoform indicates embryonic identity, these postmitotic neurons responded to depolarization of the plasma membrane, expressed functional neurotransmitter receptors and showed robust spontaneous calcium oscillations indicative of their functional maturity. Overall, the human identity and the reproducible neuronal differentiation qualify the model for use in pharmaceutical neuroscience.

Overexpression of HA-tagged human A53T α-synuclein without fibrils91 treatment did not result in detectable aggregation in our model although reported by others [START_REF] Ryan | Isogenic Human iPSC Parkinson's Model Shows Nitrosative Stress-Induced Dysfunction in MEF2-PGC1α Transcription[END_REF][START_REF] Mazzulli | α-Synuclein-induced lysosomal dysfunction occurs through disruptions in protein trafficking in human midbrain synucleinopathy models[END_REF][START_REF] Brazdis | Demonstration of brain region-specific neuronal vulnerability in human iPSC-based model of familial Parkinson's disease[END_REF][START_REF] Prots | α-Synuclein oligomers induce early axonal dysfunction in human iPSC-based models of synucleinopathies[END_REF]. The fibrils91. The blots revealed a fibrils91-concentration dependent increase in α-synuclein in the Triton-insoluble SDS-soluble fractions. Representative blots of three experiments. (G) Immunocytochemistry confirmed that seeding was in neurons. Representative images of neurons stained for pS129 α-synuclein, HA-tag, MAP2 and Hoechst. Scale bars represent 100 μm. Neurons treated with fibrils91, but not doxycycline, were negative for pS129 α-synuclein demonstrating that the pS129 antibody does not stain the seed material (upper left). Neurons treated with doxycycline, but not fibrils91, were negative for HA-tag as the antibody was titrated down to show areas of high protein density (lower left). Neurons treated with fibrils91 and doxycycline revealed pS129 αsynuclein and HA-tag positive puncta (right). A rare non-neuronal cell positive for HA-tag stain is visible in the image. No seeding was seen in non-neuronal cells. Representative images from one out of three experiments. https://doi.org/10.1371/journal.pone.0261536.g002 reasons for the difference could be many; the use of a patient genetic background, the detection assays used, antibody specificity, the seed material, the culture conditions, etc.

Exposing the doxycycline treated neurons to recombinant human wild-type fibrils91 resulted in increased FRET signals in pS129 α-synuclein and pS129-pS129 α-synuclein aggregation assays, increased numbers of pS129 α-synuclein positive puncta in immunocytochemistry and the appearance of Triton-insoluble SDS-soluble HA-tagged α-synuclein. The αsynuclein aggregation FRET-based assay provides a high throughput endpoint addressing aggregation directly, which is preferable in the search for therapeutical approaches modulating this feature of Parkinson's Disease. Seeding of endogenous wild-type α-synuclein was not detectable in our system, although reported previously in studies using iPSC derived neurons https://doi.org/10.1371/journal.pone.0261536.g004 [START_REF] Gribaudo | Propagation of α-Synuclein Strains within Human Reconstructed Neuronal Network[END_REF][START_REF] Bieri | LRRK2 modifies α-syn pathology and spread in mouse models and human neurons[END_REF]. Absence of endogenous seeding can be considered an advantage in a drug discovery setting as it provides a control, demonstrating that the fibrils91 themselves are not detected by the assays.

Our standard experimental paradigm, using 1 μg/mL doxycycline, allowed for both detection of increased and decrease aggregation as demonstrated by the dependency of α-synuclein aggregation on intracellular A53T α-synuclein levels. Furthermore, the fibrils91 induced seeding could be opposed using SNCA siRNAs. This confirms that gene modulation is achievable in the model, and that it has a suitable dynamic range to detect both increased and decreased aggregation, making it suitable to be used for genetic screening for aggregation modulators.

The USP8, USP13 and USP9X enzymes of the ubiquitination system have been reported by others to modulate α-synuclein expression or aggregation [START_REF] Alexopoulou | Deubiquitinase Usp8 regulates α-synuclein clearance and modifies its toxicity in Lewy body disease[END_REF][START_REF] Liu | Ubiquitin specific protease-13 independently regulates parkin ubiquitination and alpha-synuclein clearance in alpha-synucleinopathies[END_REF][START_REF] Rott | α-Synuclein fate is determined by USP9X-regulated monoubiquitination[END_REF]. Liu et al. describe increased α-synuclein ubiquitination and clearance upon USP13 shRNA injection in substantia nigra of mice [START_REF] Liu | Ubiquitin specific protease-13 independently regulates parkin ubiquitination and alpha-synuclein clearance in alpha-synucleinopathies[END_REF]. In contrast, we found no effect of USP13 knockdown on monomeric α-synuclein protein levels. Liu et al. also report that USP13 knockdown protects from cell death induced by lentiviral α-synuclein overexpression in vivo. In contrast, we find increased formation of aggregation upon seeding pointing to a non-protective effect of USP13 knockdown. USP9X knockdown in SH-SY5Y cells is described to promote the formation of monoubiquitinated αsynuclein species, thus increasing the steady state levels of total α-synuclein and enhancing the formation of toxic α-synuclein inclusions upon proteolytic inhibition [START_REF] Rott | α-Synuclein fate is determined by USP9X-regulated monoubiquitination[END_REF]. Knockdown of USP9X in our model based on iPSC neurons resulted in decreased monomeric α-synuclein levels and reduced seeded aggregation. USP8 knockdown in SH-SY5Y cells reduced α-synuclein levels and protected against α-synuclein expression induced toxicity (rough eye phenotype) in Drosophila flies [START_REF] Alexopoulou | Deubiquitinase Usp8 regulates α-synuclein clearance and modifies its toxicity in Lewy body disease[END_REF]. In our iPSC neuron model, USP8 knockdown did not affect monomeric α-synuclein levels while increasing seeded aggregation. The discrepancies between these findings may be due to factors such as the animal species, α-synuclein species, in vivo versus in vitro system, the cell type, the protein expression levels, off-target effects of siRNAs/ shRNAs, etc. This strongly highlights the importance of confirming biological findings in independent model systems. In short, our data indicate that USP8 and USP13 facilitate degradation of aggregated α-synuclein in seeded neurons without affecting levels of monomeric αsynuclein (data from neurons not exposed to fibrils).

We succeed in generating a reproducible assay with a robust seeding window suitable for genetic screens. The cortical differentiation protocol was selected due to its high reproducibility, HA-tagged α-synuclein was used to be able to differentiate endogenously formed fibrils and A53T α-synuclein was used to increase the assay window as it is more aggregation prone than wildtype α-synuclein [START_REF] Conway | Accelerated in vitro fibril formation by a mutant α-synuclein linked to early-onset Parkinson disease[END_REF]. It remains to be shown whether the same assay build on dopaminergic neurons and non-tagged wildtype α-synuclein would give identical results to understand better translatability into humans and the broader Parkinson's patient cohort.

In overall conclusion, our study shows that doxycycline induced HA-tagged A53T α-synuclein expression can be used to generate a robust humanized seeding model. The model exhibits two fundamental characteristics for genetic screens and target validation in pharmacological research: a high reproducibility and a remarkable versatility for throughput studies.

sequence: The Neo integration PCR (left) detects the correct insertion of the construct containing the reverse tetracycline transactivator (M2rtTA). The Puro integration PCR (middle) detects the insertion of the construct with the SNCA cDNA. The absence of the wildtype sequence of the AAVS1 locus (right) confirms that both alleles of the AAVS1 locus contain an integration and that no contaminating wild-type cells are present. X indicate samples that are not part of this study. (B) Normal karyotype was shown by G-banding: The iPSCs were incubated with Colcemid (Gibco) for 1.5 h and split with Accutase. Next, cells were incubated with 0.075M KCl for 30 minutes at 37˚C and afterwards fixed with a mixture of 25% acidic acid and 75% methanol. Fixed cells were placed at -20˚C overnight and then shipped for G-band karyotyping at the Institute of Medical Genetics and Applied Genomics, University of Tu ¨bingen. A normal male karyotype was found, 46,XY (C) Immunocytochemistry demonstrated expression of the pluripotency markers Oct-4, Sox2 and Tra-1-60: The cells were fixed with ice cold methanol, blocked and permeabilized with PBS containing 2% BSA and 0.1% Triton-X-100, and incubated with the primary antibodies over night at 4˚C (anti-Oct4, polyclonal goat, Abcam ab27985, RRID:AB_776898, 1:100; anti-Sox2, polyclonal rabbit, Abcam ab97959, RRID:AB_2341193, 1:100; anti-Tra-1-60, monoclonal mouse, Merck MAB4360, RRID:AB_2119183, 1:100). The cells were washed and incubated with the secondary antibody and Hoechst for 1 hour (donkey anti-goat Alexa fluor 488, Invitrogen A11055 RRID: AB_2534102, 1:200; donkey anti-mouse Alexa fluor 488, Life technologies A21202 RRID: AB_141607, 1:200; donkey anti-rabbit Alexa fluor 488, Life technologies A21206 RRID: AB_2535792, 1:200). The cells were loaded with mounting solution containing DAPI from Invitrogen and investigated by fluorescence microscopy. (D) Pluripotency was supported by the capability of trilineage differentiation in vitro, as confirmed by flow cytometry showing expression of the ectodermal marker PAX6 (5662388 BD Bioscience, Cy5, 1:50) and Sox1 (561592 BD Bioscience, PE, 1:50), the mesodermal marker CD56 (555518 BD Bioscience, APC, 1:25) and CD34 (555822 BD Bioscience, PE, 1:25) and the endodermal marker Sox17 (562594 BD Bioscience, AF, 1:50) and CD184 (555974 BD Bioscience, PE, 1:25): The iPSCs were split with accutase into a well of a 12-well plate with stem cell medium on matrigel in different densities: 200,000 cells/cm2 for ecto-and endoderm and 50,000 cells/cm2 for mesoderm. For ectodermal differentiation, the medium was changed to neural induction medium (50% DMEM F12 and 50% Neurobasal medium, 1X B27 without retinoic acid, 1X N2 supplement, 1X glutamax, 1X Pen/Strep (all Gibco), 10 μM SB431542, 0.1 μM LDN193189 (both from Selleckchem)) on day one. The medium was changed every day until day six. For endodermal differentiation, the medium was changed to MCDB131-1 medium (MCDB131 medium, 1.5 g/L NaHCO3, 1X glutamax, 1X Pen/Strep (all Gibco), 10 mM glucose (Sigma), 0.5% BSA) on day one including 3 μM CHIR99021 (Selleckchem) and 100ng/mL Activin A (Cell Guidance Systems). On day two, CHIR99021 was withdrawn and MCDB131-1 medium with activin A was changed every day until day six. For mesodermal differentiation, the medium was replaced by mesodermal induction medium (APEL medium (Gibco), 25 μg/mL Activin A (Cell Guidance Systems), 30 ng/mL BMP4 (Peprotech), 50 ng/mL VEGF (peprotech), 1.5 μM CHIR99021 (Selleckchem)), which was left on the cells for two days. On day three, the medium was changed to vascular specification medium (APEL medium, 50 ng/mL VEGF, 10 μM SB431542 (Selleckchem)), which was then changed every day until day six. The Staining Buffer Set from Invitrogen was used for flow cytometry analysis. iPSCs were detached with Accutase for 5-10 minutes and 200,000 cells were fixed in 0.5 mL of the fixation/permeabilization buffer according to the guidelines. After 30 minutes incubation at RT, the cells were washed three times with permeabilization buffer and then incubated for 45 minutes with 100 μL of antibody. Cells were washed additionally three times and

  The line was made by inserting two constructs into the AAVS1 safe locus site: One allele contained the SNCA-A53T-HA with a tetracycline responsive element and the other allele contained the reverse TET transactivator (Fig 1A). Addition of doxycycline to the pluripotent iPSCs confirmed inducible expression of HA-tagged A53T α-synuclein (Fig 1B and S1 Fig). Quality control experiments demonstrated purity, normal karyotype and pluripotency of the line (S1 Fig).

Fig 1 .

 1 Fig 1. Establishment and differentiation into neurons of iPSC line BIONi010-C-24 with doxycycline inducible A53T αsynuclein. (A) Illustration of the doxycycline (DOX) inducible SNCA expression system: One allele contained the m2 version of the reverse TET transactivator (M2rtTA) under regulation of the CAG promoter. When doxycycline binds the tetracycline reverse transactivator (rtTA), its 3D structure changes facilitating binding to the tetracycline responsive element (TRE) leading to the expression of SNCA-A53T-HA. Neomycin (Neo) and puromycine (Puro) resistance genes enable double selection. SA = splice acceptor, 2A = self-splicing protein. (B) Western Blot analysis for total α-synuclein in pluripotent iPSCs exposed to doxycycline for 2 days showed that the protein was expressed in an inducible manner. (C) Immunocytochemistry of iPSC neurons (day 60 of differentiation) showing expression of the pan-neuronal filament markers tau and βIII-tubulin and the synaptic protein synapsin-I in most cells. A subpopulation of GABA positive cells was present as well as a few cells positive for the neuroprogenitor markers nestin. Scale bars. 200 μm. Representative images from three independent neuronal differentiations. (D) Western Blot analysis for tau, βIII-tubulin and synapsin-I expression in undifferentiated iPSCs (day 0) and iPSC neurons (day 59 of differentiation). Representative blots of three independent neuronal differentiations. (E) RT-qPCR for gene expression of pluripotency marker POU5F1 (Oct-4), pan-neuronal markers TUBB3 (βIII-tubulin), MAPT (tau), MAP2 (microtubule associated protein-2) and

  https://doi.org/10.1371/journal.pone.0261536.g001 treated with doxycycline) as illustrated in Fig 2A. The experimental paradigm was as outlined in Fig 2B unless otherwise indicated. The expression of endogenous wild-type α-synuclein increased from day 46 to day 60. Doxycycline treatment of the cells resulted in stable overexpression of HA-tagged human A53T α-synuclein reaching levels of 16-, 9-and 7 folds above endogenous wild type protein on day 46, 53 and 60, respectively (Fig 2C). Overexpression of HA-tagged human A53T α-synuclein without fibrils91 treatment did not result in aggregation detectable by the pS129 α-synuclein FRET assay, the pS129-pS129 α-synuclein aggregation FRET assay, immunocytochemistry or fractionation (Fig 2D, 2E and 2G; S1 Raw images).

  First, α-synuclein aggregation was shown to increase with increasing concentration of doxycycline (Fig 3A) e.g. increasing A53T α-synuclein expression (Fig 3B), without measurable effect

Fig 2 .

 2 Fig 2. Exogenous fibrils91 seed α-synuclein aggregation in iPSC neurons expressing A53T α-synuclein. All graphs are shown as mean±SEM and ���� p<0.0001. (A) Overview of the in vitro seeding model. (B) Outline of seeding experiments. Day 21 neuronal progenitors were thawed, expanded and split before changing from neuronal expansion medium to neuronal maturation medium day 26. A final split was performed day 32 followed by addition of doxycycline (dox) day 33. Fibrils91 seeds were added day 46 and left on for 96h before the media was refreshed. Cells were harvested for analysis day 60 unless otherwise indicated. (C) Electrochemiluminescence immunoassays showed upregulation of endogenous α-synuclein (αSyn) levels during neuronal differentiation. Doxycycline treatment resulted in stable A53T α-synuclein overexpression. Data were analyzed by one-way ANOVA followed by Sidaks's multiple comparison test comparing doxycycline treated with non-treated for each timepoint (n = 3, m = 4). (D) FRET assays for pS129 α-synuclein showed a fibrils91-concentration and time dependent increase in the signal in doxycycline treated cells. Data were analyzed by one-way ANOVA followed by Dunnett's multiple comparison test comparing to fibrils91 untreated samples for doxycycline treated and non-treated, respectively (n = 4, m = 6). (E) FRET assays for pS129 α-synuclein and aggregated pS129-pS129 α-synuclein on identical sample sets showed strong correlation. Data were analyzed by one-way ANOVA followed by Dunnett's multiple comparison test comparing to fibrils91 untreated samples for doxycycline treated and non-treated, respectively (n = 4, m = 6). (F) Western blotting for pS129 α-synuclein, HA-tag and total α-synuclein in the Triton-soluble (Triton bar) and Triton-insoluble SDS-soluble (SDS bar) fractions from doxycycline treated neurons exposed to PBS, 0.1 μg/mL or 1 μg/mL

  and S1 Raw images). None of the siRNAs induced cytotoxicity as determined by total protein content determinations (Fig 4B). USP13 and USP8 knockdown led to increased α-synuclein aggregation (Fig 4C) without affecting α-synuclein protein levels (Fig 4D). USP9X knockdown reduced α-synuclein aggregation and protein levels (Fig 4C and 4D).

Fig 3 .

 3 Fig 3. A53T α-synuclein seeding can be modulated. All graphs are shown as mean±SEM, ns = non-significant, � p<0.05, �� p<0.01 and ���� p<0.0001. (A) All samples were seeded with 0.1 μg/mL fibrils91. FRET assays for pS129-pS129 α-synuclein aggregation showed a doxycycline-concentration dependent increase in the signal. (B) Electrochemiluminescence immunoassay showed a doxycycline-concentration dependent increase in A53T αsynuclein expression. (C) Measuring the cellular ATP content revealed no cytotoxicity from addition of 0.25-4 μg/mL doxycycline. (D) All samples were seeded with 0.1 μg/mL fibrils91. FRET assays for pS129-pS129 α-synuclein aggregation showed a doxycycline dependent increase in the signal, which was reduced by addition of SCNA siRNAs, but not a pool of four non-targeting (NT) negative control siRNAs. (E) Relative gene expression analysis confirmed knockdown of SCNA mRNA by the SCNA siRNAs but not the NT siRNAs. No MAP2 off-targets effects were seen. (F) Western blotting for α-synuclein and HA-tag. Representative blots of two experiments. (G) Quantification of α-synuclein blotting in F) confirmed a doxycycline dependent increase α-synuclein expression. Both wild-type endogenous α-synuclein (14 kDa) and HA-tagged A53T α-synuclein (16 kDa) were reduced by addition of SCNA siRNAs but not a pool of four non-targeting (NT) negative control siRNAs. Data in A), B) and C) were analyzed by one-way ANOVA followed by Dunnett's multiple comparison test comparing to doxycycline untreated samples (n = 3, m = 6). Data in D), E) and G) were analyzed by one-way ANOVA followed by Sidaks's multiple comparison test comparing the groups indicated with horizontal black bars (D. n = 3, m = 6; E. n = 3, m = 4; G. n = 2, m = 2). https://doi.org/10.1371/journal.pone.0261536.g003

Fig 4 .

 4 Fig 4. Enzymes of the ubiquitin system affect A53T α-synuclein aggregation. (A) Quantification of western blots (shown in S1 Raw images) showed reduction of USP13, USP8 and USP9X upon addition of specific pools of four siRNAs. The non-targeting (NT) negative control siRNAs were also pools of four siRNAs. (B) Quantification of total protein content in western blot samples showed no significant differences between treatments. (C) All samples were seeded with 0.1 μg/mL fibrils91. FRET assays for pS129-pS129 α-synuclein aggregation showed increased aggregation in response to USP8 and USP13 siRNAs and decreased aggregation in response to USP9X siRNAs. SCNA siRNAs was included as positive control and non-targeting (NT) siRNAs as negative control. (D) Electrochemiluminescence immunoassay showed down-regulation of α-synuclein protein levels in doxycycline treated neurons in response to SCNA and USP9X siRNAs. Data are shown as mean±SEM. Data in (A) and (B) were analyzed by one-way ANOVA followed by Dunnett's multiple comparison test (n = 2, m = 2). Data in (C) and (D) were analyzed by one-way ANOVA followed by Sidaks's multiple comparison test (n = 3, m = 6). ns = non-significant, � p<0.05 and �� p<0.01 and ���� p<0.0001.

PLOS ONE | https://doi.org/10.1371/journal.pone.0261536 December 21, 2021

The BIONi010-C-24 line was made as part of the EBiSC project (grant agreement No 115582) that received funding from the Innovative Medicines Initiative Joint Undertaking (IMI2 JU), resources of which are composed of financial contribution from the European Union and EFPIA companies' in-kind contribution.

All relevant data are within the paper and its Supporting information files.

Competing interests: NO authors have competing interests.

Supporting information