dataesr/scientific-paragraphs-categorization huggingface dataset:
Given that implantation occurs around day 9. 5-10. 5 in the macaque  52  , immediately prior to the timepoints represented here, it is possible that the trophoblast is transitioning through a similar differentiation process as observed in our rabbit atlas at GD8. The ability of the mouse and rabbit references to precisely and accurately annotate distinct cell types within the human and macaque datasets illustrates the utility of using a more comprehensive reference set when performing cell type annotation. Since the rabbit and mouse cell types were annotated consistently, both rabbit and mouse annotation models provide unique insights, which can either increase confidence in the annotation or highlight possible cross-species differences. While the rabbit-trained model more accurately identified human PGCs and epiblast cells, the mouse-trained model predicted a smoother transition of mesodermal cell types across the UMAP embedding, possibly reflective of the higher resolution of 6-hour timepoints sampled (Figure  S8B ).
Its large size also makes it easy to work with and characterize. Given the key role the yolk sac plays in early development, we elected to both characterize this developmental process and to interrogate how different cell types may interact with one another. To this end, we mined our sc-RNAseq dataset for markers genes that are key to distinguish the different layers of the yolk sac hematopoietic niche in rabbit at GD9, and spatially resolved their localisation in situ, using single-molecule resolution microscopy (RNAscope) hybridisation. Using highly-spatially resolved RNAscope experiments, we identified blood cells, which can be seen inside of vessels formed by CDH5+ endothelium (Figure  5A ). Some of the blood cells are RUNX1 -while others are RUNX1 +. AHNAK high mesothelium cells can be seen surrounding these CDH5+ endothelium cells (Figure  5A ), with the blood cells in the center.
By inputting the cell-types that make up the yolk sac blood islands, we observed that the highest score of inferred interaction pairs are between mesothelium and YS endothelium (Figure  5E ). This includes various extracellular matrix proteins like nectin, fibronectin-1 (FN1), collagen genes (Figure  S9A, B ) and important blood maturation ligands like VEGFA and VEGFC, expressed by the mesothelium. Having identified putative signaling crosstalk between rabbit mesothelium and endothelium, we next explored to what extent these interactions might be conserved in human, taking advantage of our refined annotation of the human CS7 dataset (Figure  4A ). This analysis revealed high expression of VEGFA/C signaling from human mesothelium and consequent VEGFR1(FLT1) and FLT4 expression in the YS endothelium, thus supporting the interaction discovered in our rabbit data (Figure  5D )  56.
DISCUSSION In this study, we report high-resolution morphological and molecular maps of early rabbit development, covering gestational days 7, 8 and 9. Using cutting-edge single-cell genomics approaches, we characterized the transcriptional and chromatin accessibility profiles of over 180,000 individual cells, isolated from whole rabbit embryos, across gastrulation and early organogenesis. To our knowledge, this represents the most comprehensive view into the molecular landscape of early non-rodent mammalian development to date. We leveraged these datasets, applying a novel neighborhood-based approach, to identify what developmental processes were most conserved, at the molecular level, between the rabbit and the mouse. We observed substantial differences in molecular profiles between the gut, PGCs and extra-embryonic cell types where transcriptional differences may be associated with the different embryo topologies, maternal environments, and implantation strategies of the rabbit and mouse.
Given the human gene ID for each feature, we located positions in the human genome relating to transcripts of that gene and specifically, the positions of each transcript's 3' most exon and UTR. Using the Ensembl Compara Perl API, we then obtained alignments for these sequences in the rabbit transcriptome, taking only those with maximum alignment scores. Exon and CDS annotations spanning these positions were then added to the rabbit GTF file, similar to the extended 3' ends. Additional measures were taken to ensure that the alignments extracted did not overlap existing annotations (including the added extensions) and that alignment sequences for the same transcript were proximal (i. e. not mapping to different chromosomes). As a result of these changes, the median number of genes and UMIs increased from 1875 to 2661 and from 6987 to 10,126 respectively.
61. The 300 nearest neighbors in the batch-corrected principal component analysis were considered, with a 'min. dist' parameter of 0. 9. Force-directed graph layouts were also computed on the batch corrected PCs with the 'scanpy. api. tl. draw_graph' function using the ForceAtlas2 algorithm. The whole dataset was clustered using the Leiden algorithm (scanpy. api. tl. leiden) using the same neighborhood graph constructed for the UMAP embedding. Clusters were generated using a range of resolution parameters (1 -10) to assist in making cell type annotations at. CC-BY 4. 0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprint this version posted October 6, 2022.
Cell type annotation Cell type annotations were initially predicted at single-cell resolution by training an automated annotation model on a recent single-cell atlas of mouse gastrulation and early organogenesis  26. The Rosshandler dataset consists of 430,339 cells across embryonic days 6. 5 to 9. 5, sampled at 6-hour intervals. These cells were manually curated into 87 different cell types following consultations with domain experts and a detailed examination of known marker genes. Given that the rabbit is less-well studied and has fewer known marker genes than the mouse, transferring labels across species served to provide a fast, initial annotation that could be scrutinized and validated in more detail. To transfer annotations from the mouse to the rabbit, we utilized SingleR, an automated annotation method that assigns cell type labels based on correlation in expression between reference and query cells  62.
In order for the cell type annotation model to relate the rabbit and mouse datasets, we constructed a common feature set between the atlases using one-to-one ortholog genes. For each rabbit gene, we extracted its mouse homolog from Ensembl and cross-referenced it with the Ensembl codes of the mouse atlas. Many-to-one and many-to-many genes were excluded by filtering on the orthology type also provided by Ensembl. We then trained a SingleR model on the mouse atlas (with 'trainSingleR' from the 'SingleR' package using parameters 'de. n = 50' and 'de. method = Wilcox'), providing the cell type annotations from the original study. Given the large number of cells, we opted for the pseudobulk option within SingleR, which clusters cells of the same cell type into pseudo-samples, reducing the computational work.
Rabbit-mouse integration The rabbit and mouse  26  datasets were integrated into a joint embedding using SAMap version 0. 1. 6  63. An initial mapping between rabbit and mouse features was first obtained by reciprocal BLAST aligning the rabbit and mouse transcriptomes. The transcriptomes were obtained through Ensembl and aligned using the 'map_genes. sh' script, provided with SAMap (https://github. com/atarashansky/SAMap/). This script outputs a table of sequence similarity scores between rabbit and mouse transcripts, which are used by SAMap as an initial weighting between features. Since the features of the rabbit and mouse datasets are given in terms of Ensembl gene IDs, prior to running SAMap, we linked each Ensembl transcript ID with its associated gene ID. The SAMap algorithm also utilises self-assembling manifolds  64  to align datasets and so these were computed for the rabbit and mouse atlases using Scanpy (scanpy.
tl. external. sam). We then constructed a SAMAP object, providing the rabbit and mouse SAM objects, the transcript to gene mappings and the map_genes output directory. We then ran the SAMap pipeline on this SAMAP object (SAMAP. run) to integrate the two datasets. For visualisation purposes, we re-computed a UMAP embedding on the integrated anndata object (SAMAP. adata) using a minimum distance parameter of 0. 8.
scATAC-seq and scRNA-seq integration To perform cell-type annotation of the scATAC-seq cells, integration was performed with the labeled scRNA-seq dataset for each stage separately. For each stage, dimensionality reduction was performed using addIterativeLSI, using dimsToUse=1:35 and nFeatures=55000 and clustering was performed using addClusters with resolution=1. 5 to create high-resolution per stage clusters. Imputation matrices were added using addImputeWeights. Next, scRNA-seq integration was performed using addGeneIntegrationMatrix, which compares the estimated gene activity stored in the 'GeneScoreMatrix' for scATAC-seq with the gene expression of the scRNA-seq dataset in order to perform label transfer for cell-types and add a 'GeneIntegrationMatrix' containing integrated RNA expression counts.
scATAC-seq trophoblast analysis Trophoblast and amnion cells were isolated from the total dataset by isolating cells belonging to clusters 1 to 5. To manually curate the annotation of these cells, dimensionality reduction and reclustering was performed with addIterativeLSI, using dimsToUse=1:25 and nFeatures=25000 and addClusters with resolution=1. Markers of the 'GeneScoreMatrix' and 'MotifMatrix' per cluster were used to relabel Trophoblast, Cytotrophoblast, Syncytiotrophoblast progenitors, Syncytiotrophoblast, and Amnion cells. Next, Amnion cells were excluded and dimensionality reduction was repeated on these cells with dimsToUse=1:10 and nFeatures=15000, and the UMAP was generated using nNeighbors=35 and minDist=0. 8. Trajectory analysis was performed by determining pseudotime ordering using addTrajectory over the UMAP embedding.
We next identified a set of features with which to compare neighbourhoods across the two species. Specifically, we selected the intersection of the top 2000 highly variable genes, computed independently for each dataset. The intersection of genes was chosen to avoid confounding differences in expression with technical variation resulting from the lower quality rabbit genome annotation. It ensured that genes selected for comparisons were expressed and highly variable across both datasets. We also experimented with the number of highly variable genes but found that our results changed very little above 2000 HVGs. Features were selected using the 'getScranHVGs' function from the scran R package. We excluded from this mitochondrial genes and those associated with cell cycle GO terms. Finally we selected only those which are one-to-one orthologs across the rabbit and mouse.
Our final set of features consisted of 796 genes. Using this set, we computed the mean expression profile within each neighbourhood (using the normalised and log-transformed gene counts). We also implemented a version of gene-specificity, used in previous cross-species comparisons, to account for differences in quantification and absolute values between datasets  69, 70   The gene specificity values ( ) for each gene, across all neighbourhoods, were then used to    compute the Pearson correlation between all pairs of rabbit and mouse neighborhoods forming a matrix of neighbourhood similarity. In Figures  3  and S7 , we visualised neighbourhood similarities across the rabbit and mouse UMAP embeddings. For each rabbit neighbourhood we extracted its maximum similarity score with any mouse neighbourhood (and vice versa).
We then plotted each neighbourhood according to the UMAP position of its index cell, and coloured each point by the neighbourhood's maximum correlation value. In Figures 3D, the maximum correlation values are aggregated according to the cell type annotation of each index cell, to obtain a distribution of similarity scores across each cell type. Finally, for the trajectory comparisons, subsets of neighbourhoods were extracted from both the rabbit and mouse datasets whose index cell was annotated as one of a specified set of cell types. These neighbourhoods were then plotted in the same way as Figure  3C  and repositioned or reflected to facilitate a clear visual comparison. Lines were drawn between maximally correlated neighbourhood pairs, computed in both the rabbit-mouse and mouse-rabbit directions.
Primate cell type annotation In our cell type annotation analysis we utilised the Tyser et al. 2021  25   in-vitro cultured macaque embryos, at days 10, 12 and 14. We accessed raw count data through the GEO database under accession number GSE148683 and downloaded associated metadata from the accompanying shiny app at https://www. nhp-embryo. net/. To combine the day 10, 12 and 14 samples, we performed batch correction with fastMNN with k=20. We also normalized and log-transformed the raw counts using 'computeSumFactors' and 'logNormCounts' functions in the 'scran' and 'scuttle' R package. To predict cell types in the human and macaque datasets, we applied the same approach used to annotate the rabbit atlas originally. We trained SingleR annotation models on the rabbit and mouse datasets both jointly and independently.
Datasets were filtered to include only one-to-one ortholog genes between each reference (rabbit or mouse) and query (human or macaque) dataset. As before, the training and classification steps were performed using the 'trainSingleR' and 'classifySingleR' functions from the 'SingleR' R package, with the 'aggr. ref=TRUE' option. In Figure  4 , we display the 'pruned' cell type annotations from SingleR, which have undergone a filtering step to remove low-quality assignments.
CellphoneDB analysis of Yolk Sac Haematopoiesis Cell-cell interactions were predicted using a previously curated list of ligand-receptor and receptor-receptor pairs using CellphoneDB v. 2. 0  55. Due to the exponential nature of possible interactions across the full dataset, only a subset of cell types corresponding to the rabbit yolk sac and blood islands were used, covering 'EMP', 'visceral endoderm', 'yolk sac endothelium', and 'YS mesothelium'. The visceral endoderm cell type in the CellphoneDB heat map analysis used the merged cell type of 'visceral endoderm', 'visceral ys endoderm 1', and 'visceral ys endoderm 2' to cover all possible interaction pairs between different subtypes of the visceral endoderm as well as other cell types.
Data and code availability Raw sequencing data are available on ArrayExpress with the following accessions: scRNA-seq: E-MTAB-11836; scATAC-seq: E-MTAB-11804. Links to the data and code are available at https://marionilab. github. io/RabbitGastrulation2022/. The web page also links to processed single-cell transcriptomics data in a variety of formats for loading into both R and python analysis pipelines. The data is also available to explore interactively via a web app accessible through the same link. All other data are available from the corresponding authors on reasonable request.
and Yang et al. 2021 50 single-cell datasets of early human and macaque development. The Tyser et al. 2021 SMART-seq dataset consists of 1,195 cells from a single CS7 human embryo acquired through the Human Developmental Biology Resource. We downloaded the normalised gene expression values and UMAP coordinates from http://human-gastrula. net/. Conversely, the Yang et al. 2020 dataset was generated through 10X genomics single-cell RNA-sequencing of.
Acknowledgements We would like to thank  K. Katarzyna  and the  CRUK-CI Genomics Core  for their help with the 10X Genomics sample processing;  R. Argelaguet  for bioinformatic support;  N. Wilson  for assistance with the facilities and technical support;  M. Rostovskaya  for discussions around eutherian amnion. We would also like to thank  Mark Keller  ( Gehlenborg lab ) for help in setting up the Vitessce visualization platform. M. -L. N. T. is funded by a Herchel Smith PhD Fellowship in Science. D. K. is funded by the  Wellcome Mathematical Genomics and Medicine Programme at the University of Cambridge  ( PFZH/158 RG92770 ). C. G. was funded by the  Swedish Research Council  ( 2017-06278 , administered by  Sahlgrenska Cancer Center, University of Gothenburg ). Work in the Gottgens group is supported by  Wellcome ,  Bloodwise ,  MRC ,  CRUK , by core support grants from  Wellcome  to the  Wellcome-MRC Cambridge Stem Cell Institute.
Surface leaching following N input A leaching coefficient is calculated for each fertiliser application following the adapted Burns model  (Burns, 1976; Scotter et al. , 1993). Parameters for this coefficient are estimated based on user data (see section 3. 5 in the Supplementary Material). The leachable N from the fertiliser is determined by deducting ammonia volatilisation, and then it is multiplied by the leaching coefficient to determine surface leaching.
Background data As the goal of the model is to support the construction of agricultural life cycle inventories representing contrasting agricultural situations, its implementation was supported with the inclusion of a database containing global data pertaining to pedoclimatic conditions and crop characteristics. The database provides the soil data and monthly climatic data, as well as crops and fertilisers data, necessary to run the model. Pedoclimatic data are available for each combination of country and Global Ecological Zone (FAO, 2012) worldwide, as well as for each French department. All background data are presented in the accompanying data paper  (Galland et al. , submitted) , and are directly available as a downloadable dataset. The model's list of crops aims at representing the most important global agricultural products.
Data for model testing and validation The model was implemented as a spreadsheet tool, and was tested with several datasets of agricultural itineraries, depicting contrasting agricultural situations. The agricultural datasets, described in Table  3 , were used to test the model. Observed values of leached nitrate were plotted against predicted values, and a set of non-correlated statistics were computed according to  Yang et al. (2000) , since a single criterion is insufficient for assessing model performance  (Bellocchi et al. , 2010). We thus computed the following statistics to assess the predictive quality of the model: i) the mean deviation (MD) assessing whether the model over-or underestimates the observation (Eq. 4), ii) the mean square error (RMSE) assessing the mean deviation around the prediction (Eq.
For the French dataset, the models showed acceptable results for all the statistical criteria, being close for MD and MEF to the performance of other models  (Tribouillois et al. , 2020). A value for MEF = 0. 70 was considered by these authors as "good", based on the review by  Cannavo et al. (2008). It should be noticed that for the whole dataset and the French dataset, model respectively tended to slightly over-or underestimate the observations. 3 )    3 ) For temperate conditions (n=15), IN-3 yielded the estimations closest to the measured values while for tropical (n=8) estimations I-N3 outperformed SQCB. In organic systems, I-N2 estimations were generally closer measured values). Overall I-N2 presented respectively a lower deviation (MD=-0. 2), but a higher root mean square error (RMSE=8), and a lower model efficiency (MEF=0.
It is not yet valid for crop associations, perennials, rice or cotton, all of which demand important design adaptations. In any case, the quality of parametrisation depends on the amount of knowledge and research effort made for each crop. If for main field crops such as maize, wheat, etc. the knowledge set is large, this is not the case for minor crops, and especially not for short cycle crops like market garden vegetables. For these crops, an effort would be required to improve the parameter dataset before a validation experiment should be performed. In spite of several attempts, it was not possible to find an adequate dataset for this aim. For organic fertilisers, mineralisation curves were taken from studies in temperate areas. Furthermore, curves which are more representative of tropical organic fertilisers have to be added.
Some studies, like those of  Thuriès et al. (2002)  and  Rothé et al. (2019) , may offer possibilities in the spreadsheet implementation of the model and the associated database. Lastly, a simplified approach inspired partially from the work of  Pervanchon et al. (2005)  was applied for grasslands, so that grazing excretions are considered as punctual fertilising inputs, instead of as a continuous distributed input over a specific time period. Here again, we proposed a compromise to avoid an over complexification of the model. Actually, nitrate leaching from grasslands is low when they are not over fertilised  (Vellinga et al. , 2001) , which is taken into account in the model, while the ploughing of a grass cover causes a peak of nitrogen mineralisation, source of high nitrate leaching  (Kayser et al.
To illustrate, two rotations were calculated for a range of values for some key variables, and the results show the magnitude of the model's sensitivity (Fig. 10 ). Monte Carlo or other uncertainty propagation methods where not used to test for the model sensitivity to parameters' values, because such bootstrapping-oriented methods, when used alone, are more suitable for understanding data uncertainty, not model uncertainty or model sensitivity  (Groen and Heijungs, 2017). For a true global sensitivity analysis, a bootstrapped dataset (e. g. via Monte Carlo) would have to be further processed to identify the relative contribution of each parameter to the model's outputs  (Xiao et al. , 2017) , via for instance variance decomposition methods (e. g. Puy et al. , 2021)  or regression-based relative importance analysis methods (e.
(2018)  who found that hydraulic parameters, especially during dry seasons, mostly controlled N leaching. This highlights the importance of personalised high quality data for hydraulic variables determining drainage. Even if default data exist, data collection should always be performed as much as possible, especially for soil data. The validation test showed that the presented acceptable predictive performances for temperate arable conditions and outperformed other models like the previous version of I-N  (Bockstaller et al. , 2008; Bockstaller and Girardin, 2010) , or the models used in ecoinvent (SALCA-NO3 and SQCB), but remains to be improved for tropical conditions. However, the dataset from  Benoit et al. (2016)  contained 14 points for three pedoclimatic situations with mainly wheat fields, where pedoclimatic conditions remained quite close since fields were located around Paris.
Nonetheless, this dataset was interesting because it contained conventional and organic fields, the latter having been poorly fertilised. A temperate situation with heavy organic fertilisation should be tested with the model. Likely, further efforts have to be made to test the model for more cropping systems and conditions like short cycle crops, like vegetables or grasslands. Thus, the results presented here have to be considered as a first validation attempt. Actually, as pointed out in  Keating (2020) , validation of model should be a continuous process.