dataesr/scientific-paragraphs-categorization huggingface dataset:
Quality assessment and functional annotation of the Unigene dataset From the Trinity package, we obtained the contigs arranged into clusters representing unique gene sequences, called Unigenes. Using this information, we rearranged the identified transcripts into 24 866 Unigenes, which represent the full-length sequences of unique transcript (Table  1 ). According to the transcript (contigs) and Unigene numbers, we identified an average of 1. 46 transcript isoforms per Unigene across the different black-eyed pea tissues. The Unigene statistics are displayed in Table  1  and are summarized by a N50 length of 2187 bp and an average length of 1610 bp. To evaluate the quality and the effectiveness of the assembly and to analyse the representativeness of the Unigene set, we compared the transcript sequences identified from black-eyed pea with known protein sequences from Arabidopsis and soybean, and with an existing black-eyed pea EST database (available at the HarvEST webserver, www.
harvest. ucr. edu). We revealed a high ortholog proportion between black-eyed pea Unigene sequences and soybean and Arabidopsis protein sequences, respectively 83. 37 and 80. 98% alignment coverage considering the nearly full-length of transcript sequences (Table  S2 ), supporting the effectiveness and accuracy of the assembly process. Compared to the existing black-eyed pea EST database, we identified 77% of the transcript sequences identified in this study were present in the EST database (Table  S3 ), displaying a high proportion of transcript sequence similarity identified through two different approaches. In summary, the predictive transcript sequences identified using de novo assembly, displayed high homology with known black-eyed pea ESTs and higher ortholog proportion with closely related species like soybean (G.
max) than more distant species such as Arabidopsis, which assessed the quality of the de novo assembly resulting in the transcript sequence identification. To further analyse the Unigene set, we annotated the 24 866 Unigene sequences using different strategies. First, using Trinotate from the Trinity software, 19 535 Unigenes (78%) were annotated by homology search into SwissProt, PFAM and EMBL Uniprot databases. Then, sequences containing transcription factor (TF) domains were identified using the online software PlantTFcat  (Dai et al. , 2013). From this analysis, 2485 TFs were predicted to have TF conserved domains and annotated according to their TF family. A family representation of these TFs is provided in the Table  S4. According to our prediction, 10% of the Unigenes (2485 TFs/24 866 Unigenes) were predicted to belong to a TF family and as described in other species, the most abundant TF family is the zinc finger family with the C2H2 type displaying 413 members.
Others common families displayed more than 100 member such as the MYB, bHLH and AP2-EREBP families with respectively 153, 119 and 102 members  (Verdier et al. , 2008 ; see review  Qu and Zhu, 2006). Finally, to provide an overview of the biological functions and metabolic pathways, we attributed GO annotations and Mapman bin codes to the Unigene sequences using respectively Trinotate (part of the Trinity bundle) and Mercator (http://mapman. gabipd. org/web/ guest/app/mercator) tools. A Mapman bincode file is provided in Table  S5  and was formatted to be directly usable as a mapping file in PageMan or Mapman software  (Usadel et al. , 2006).
Transcript abundance estimation and identification of constitutively and differentially expressed gene sets The alignment-based quantification method RSEM was used to estimate transcript abundance of all 36 529 identified transcripts based on fragments per kilobase per transcript per million mapped reads (FPKM) (Figure  1 ). A correlation matrix and a cluster dendrogram were performed in order to validate the reproducibility of the biological replicates within all the different tissues (Figure  S2 ). This correlation analysis of global transcript abundance between plant tissues showed two main groups of organs: the seed (at 8, 10, 14 and 18 DAP) and the other plant tissues (i. e. leaf, stem, root, flower, pod) (Figure  S2 ). This observation based on the global transcript abundance could be explained by the specialized functions and differentiation during seed development that might involve different set of transcripts with differential expression in other tissues.
According to our contrasted dataset containing all major tissues and the time course of seed developmental stages, we did not identified transcripts following constitutive expressions in the major tissues and during the seed development. For this reason, we separately analysed transcripts constitutively expressed during the seed development according to the same criteria than the other tissues. We identified 74 transcripts with stable expression throughout the seed developmental stages, for seed-focused analyses (Table  S6 ). Genes displaying preferential or specific gene expression in a particular plant tissue are important in the understanding of the specialized processes occurring within these tissues. In our analysis, we identified 1197 transcripts with specific tissue expression (i. e.
In order to further validate the annotations and the expression profiles of our dataset, we looked into the specific processes occurring in different tissues by an overrepresentation analysis (ORA) of tissues specific genes in different organs (Figure  S3 ). As expected, we observed over-representation of predictable specialized processes in different plant tissues such as functional classes related to 'cell wall' and 'apoplast transport' in pod wall, 'pollen tube' in flower and 'aleurone grain' and 'storage lipid' in seed. The detailed over-representation analyses using the GO terms of all the specific genes for the different plant tissues are presented in Figure  S3.
Implementation of the VuGEA, a webserver to store, display and analyse transcriptomics data To provide open access to the data generated in this study, we generated a web accessible database, called the Vigna unguiculata Gene Expression Atlas (VuGEA, website http:// vugea. noble. org). This webserver conserves the same architecture and tools that designed for the Medicago truncatula and Lotus japonicus Gene Expression Atlas  (He et al. , 2009). All the identified sequences and corresponding gene IDs, annotations and transcripts abundance information are stored in this webserver and are accessible by a user-friendly interface. Transcripts sequences and abundance across the tissues can be obtained by a 'search' function using the black-eyed pea ID, ortholog ID (from Medicago, Arabidopsis, Lotus and soybean), gene description keywords, GO annotation or TF family name.
Transcripts can also be found based on gene sequence using a blast homology search tool from a nucleic or protein sequence as query. The current version of the webserver offers two analytical modules: a multi-transcript viewer to visualize expression profiles of multiple transcripts and a co-expression analysis tool to identify positively or negatively co-expressed transcripts in selected tissues (using Pearson or Spearman correlation methods). Finally, we linked the three legume gene expression atlases (Medicago, Lotus and black-eyed pea) to have access to orthologous gene information from any of these gene expression atlases, allowing users to directly access expression profiles of orthologous gene(s) in other legume species. This webserver also allows users to download all or selected batches of data generated in this study.
Transcriptional regulation of black-eyed pea seed development To continue to evaluate the quality of the de novo assembly and the transcript abundance quantification, together with deciphering the black-eyed pea seed development, we analysed the molecular mechanisms occurring between 8 and 18 DAP. For this purpose, we extracted from the whole transcript dataset a subset of 34 164 transcripts expressed during the seed development (33 774 transcripts displaying a CV > 30% throughout seed development). According to a figure of merit analysis (FOM), this subset was divided into five different clusters according to their expression profile patterns. Using k-means clustering analysis, we identified cluster 'S8' (15 602 transcripts) displaying a maximum of expression at 8 DAP, cluster 'S8-S10' (7505 transcripts) with a maximum expression at 8 and 10 DAP, cluster 'S14' (2324 transcripts) with a maximum expression at 14 DAP, cluster S18 (5460 transcripts) with a maximum expression at 18 DAP and finally cluster 'S8-S10-S18' (3273 transcripts) expressed throughout the whole seed development except at 14 DAP.
We performed an ORA of the functional classes in each of the different clusters to provide an overview of the molecular mechanisms occurring in different clusters. Pageman software  (Usadel et al. , 2006)  and the mapping file generated in our study (Table  S5 ) were used to analyse these clusters using a Fisher's exact test in order to identify molecular mechanisms occurring at different seed developmental phases (Figure  3 ). Summary of the different functional classes with the highest statistical significances associated with the different clusters is displayed in Figure  4  and detailed analysis of all over-representation classes is provided in Table  S8. Transcripts belonging to Cluster S8 with a peak of expression at 8 DAP mainly belong to nucleotide metabolism, DNA synthesis/chromatin structure, cell organization and cell division confirming that the seed is at the embryogenesis stage, which is marked by intense cell division.
g. RNA processing, RNA binding) and protein synthesis ribosomal protein suggesting intense RNA regulation activity at these stages. In this study, we analysed if these miRNA involved in drought stress response could serve as a molecular signal in seed desiccation mechanisms. First, we predicted putative target genes of these 44 miRNA using psRNA Target algorithm  (Dai and Zhao, 2011)  and transcript sequence dataset generated in this study. In total, 28 of the miRNAs displayed putative target genes based on our black-eyed pea sequence database with differential transcript abundance during the black-eyed pea seed development (Table  S9 ). Interestingly, two of the miRNAs, Vun295724_336_200 and Vun90329_258_100, were predicted to target auxin response factor genes (i. e. Vun_T02367, Vun_T13063 and Vun_T03046), which are putative orthologs of AtARF10 and AtARF16 (Table  S9 ).
CONCLUSION In this study, we performed de novo assembly and quantification of the black-eyed pea transcriptome in the major plant tissues. We identified, annotated and quantified the expression of 36 529 transcripts across all the major plant organs. From these dataset, we identified 24 866 unique transcript sequences, called Unigenes. All the results generated were stored into a gene expression atlas webserver, displaying a user-friendly interface and providing tools for transcriptome analysis and comparison with other legume species. Using this gene expression atlas, we provided a baseline description of the molecular mechanisms occurring during seed development. Additionally, we revealed that some of the mechanisms identified as drought-associated in the plant could be similarly regulated during the seed desiccation phase and involved the same miRNAs.
Raw read filtering and de novo assembly of contigs Quality control of short reads was performed on the raw pairedend reads to discard low quality regions using SolexaQA (Q < 17)  (Cox et al. , 2010) , and adapter sequences using cutadapt (v1. 3)  (Martin, 2011). After filtering reads, paired reads with length above 25 bp were referred to as 'clean reads'. Cleaned reads were de novo assembled using the Trinity software package (version v. 2. 0. 6) with parameter setting of '-min_kmer_cov2,' a fixed k-mer size of 25 and other default parameters  (Grabherr et al. , 2013). We used the Linux cluster server available at the bioinformatics platform of the Shanghai Center for Plant Stress Biology using 1 TB of maximum memory. The complete Trinity command line used for the assembly process was: '-seqType fq-max_memory 1000G-left reads1.
Transcript quantification The alignment-based quantification method RSEM was used to estimate transcript abundance based on FPKM  (Li and Dewley, 2011). Transcripts with abundance superior to one FPKM were considered for further analysis. To remove highly similar or redundant contigs, we merged contigs with higher sequence identity than 95% using CD-HIT-EST. Reads were mapped to transcript sequences identified from de novo assembly for each analysed tissue and stage using BWA  (Li and Durbin, 2009). Then, SAMtools flagstat was used to calculate total mapped reads, uniquely mapped reads and perfectly mapped reads  (Li and Durbin, 2009).
Annotation of Unigene sequences Full-length transcript sequences were compared against the Arabidopsis thaliana database (TAIR10) and Glycine max database (G_max_109) using BLASTX algorithm (e-value 1e-20, best hit method: max_target_seqs 1, outfmt 6). Then a Perl script was used to calculate the percentage of the Arabidopsis thaliana target and Glycine max target alignment. The EST sequences from the cowpea EST database (http://harvest-web. org/hweb/pickassy. wc) were compared to our contig assembly using BLASTN algorithm (e-value 1e-5, perc_identity 95, max_target_seqs 1, outfmt 6), followed by a Perl script to calculate the percentage of cowpea EST matching our assembled dataset. Annotation of the de novo assembly contig set was achieved using Trinotate (http://trinotate. github. io/), a software included into the Trinity suite, which performed homology search to Swis-sProt, protein domain identification (PFAM) and EMBL Uniprot databases.
Functional annotation of contigs was done by Trinotate using eggnog and GO pathway databases. Lohse et al. , 2014). Identification of orthologous genes of soybean, Lotus and Medicago was performed by sequence similarity searches using BLASTX algorithm (e-value 1e-5, bit score 50, max_target_seqs 1, outfmt 6) between our assembled contigs and the Glycine max (G_max_109), Lotus japonicus (LjDB v. 6) and Medicago truncatula (IMGA Mt4. 0) gene database. Open reading frames (ORF) prediction of the transcript sequences was performed using the Trans-Decoder software (http://transdecoder. github. io/) included into the Trinity suite.
Differential gene expression analyses Global heatmap analysis of all the black-eyed pea transcripts and correlation analyses were performed using 'pheatmap' package from R software. Distance analysis was calculated using pairwise Pearson correlation and average linkage clustering method. Tissue-specific or preferential genes were identified using two criteria: only expression values in one tissue or 20-fold more expressed in a specific tissue with respect to others. GO enrichment analyses of the tissue-specific or preferential genes were performed using 'topGO' package available on R software  (Alexa and Rahnenfuhrer, 2010). For different tissues, we visualized tissue-specific genes using a heatmap representation of z-scores. zscore = (T s ÀT a )/SD, T a is the mean transcript level for a given gene across all samples, T s is the mean transcript level in a specific samples, SD is the standard deviation across all samples.
Housekeeping genes were identified using CV, where CV (coefficient of variation) was calculated as SD/mean. Housekeeping genes were selected due to their constant expression across different tissues (i. e. CV < 10%). The seed differential expressed gene set used for seed developmental analysis was made up of transcripts displaying a CV > 30% regarding their expression throughout seed development. Then, the k-means algorithm was used to cluster seed differential expressed genes using Pearson correlation as a distance metric. Subsequently, PageMan software (http://mapman. gabipd. org/web/ guest/pageman) performed ORA of functional classes using Fisher's exact test and Bonferroni correction.
Implementation of the VuGEA web server and RNA-seq data deposition The VuGEA webserver is using the same software and hardware as the Medicago and Lotus GEA web servers  (Benedito et al. , 2008; Verdier et al. , 2013). All the data were deposited in an in-house MySQL database. The 'biologist friendly' interactive web interface was designed using PHP and Javascript programing combined with additional GNU web development packages as described by  He et al. (2009). All the features and analytical tools of LjGEA (Lotus japonicus Gene Expression Atlas) and MtGEA are available the VuGEA web server. Moreover, links between these three web servers were generated to compare expression of orthologous genes in these three legume species. Raw and normalized RNA-seq data were deposited and publicly available at the VuGEA web server.
Raw data processing of EC data was made using EdiRe Data Software (University of Edinburgh, United Kingdom) according to the methodology based on the EuroFlux protocol  (Aubinet et al. , 2000)  with a further adjustment such as the WPL correction  (Webb et al. , 1980). Quality control was performed after  Foken and Vichura (1996). Net ecosystem exchange (NEE) was obtained after taking into account the storage term which was calculated using a profile system designed to make sequential CO 2 concentration measurements (IRGA, SBA-4, PP Systems) from six heights in a 6 min cycle. Quality assessment and quality control was made according to  Foken and Vichura (1996). Gap-filling of missing data and NEE flux partitioning to gross primary production (GPP) and total ecosystem respiration (TER) was made using online EC data gap-filling and flux partitioning tools (http://www.
Modeling exercises within international projects revealed additional problems that needed a solution related, e. g. , to stomatal conductance  (Sándor et al. , 2016). Recently published findings (e. g. , simulation of temperature acclimation of respiration;  Smith and Dukes, 2012 ) also motivated our work on model development. Addressing the known issues with the Biome-BGC model required diverse directions, but this was necessary due to the complex nature of the biogeochemical cycles of the soilplant system. Our overarching aim was to provide a stateof-the-art biogeochemical model that is in the same league as currently available models, such as LPJmL, ORCHIDEE, CLM, JULES, CASA, and others. The model that we present in this study is referred to as Biome-BGCMuSo (abbreviated as BBGCMuSo, where MuSo refers to a multilayer soil module).
Soil thermodynamics Since some of the soil processes depend on the actual soil temperature (e. g. , decomposition of soil organic matter), it is necessary to calculate soil temperature for each active layer. Given the importance of soil temperature (e. g. ,  Sándor et al. , 2016) , the soil module was also reconsidered in BBGC-MuSo. The daily soil surface temperature is determined after  Zheng et al. (1993). The basic equations are detailed in  Hidy et al. (2012). An important modification in BBGCMuSo is that the average temperature of the top soil layer (with a thickness of 10 cm) is not equal to the above-mentioned daily surface temperature. Instead, temperature of the active layers is estimated based on an empirical equation with daily time steps. Two optional empirical estimation methods are implemented in BBGCMuSo.
In BBGCMuSo, the time step (TS; seconds) of the soil water module integration is dynamically changed. In contrast to the modified Biome-BGC  (Hidy et al. , 2012) , the TS is based on the theoretical maximum of soil water flux instead of the amount of precipitation. As a first step, we calculate water flux data for a short, 1 s TS for all soil layers based on their actual SWC. The most important problem here is the selection of the optimal TS: TS values that are too small result in a slow model run, while TS values that are too large will lead to overestimation of water fluxes. The TS is estimated based on the magnitude of the maximal soil moisture content change in 1 s ( SWC max ; m 3 m -3 s -1 ). Equations (  3 )-(  4 ) describes the calculation of TS: TS = 10 EXPON , (3) where EXPON = abs [LV + (DL + 3)] if DL + LV < 0 EXPON = 0 if DL + LV ≥ 0.
Senescence calculation The original Biome-BGC ignores plant wilting and associated senescence (where the latter is an irreversible process) caused by prolonged drought. In order to solve this problem, a new module was implemented to simulate the ecophysiological effect of drought stress on plant mortality. A senescence simulation has already been implemented in developed Biome-BGC  (Hidy et al. , 2012)  and it was further improved in BBGCMuSo following a different approach. In BBGC-MuSo, if the plant-available SWC decreases below a critical value, the new module starts to calculate the number of the days under drought stress. Due to low SWC during a prolonged drought period, aboveground and belowground plant material senescence is occurring (actual, transfer, and storage C and N pools) and the wilted biomass is translocated into the litter pool.
Management modules The original Biome-BGC was developed to simulate natural ecosystems with very limited options to disturbance or human intervention (fire effect is an exception). Lack of management options limited the applicability of Biome-BGC in croplands and grasslands, but also in managed forests. One major feature of BBGCMuSo is the implementation of several management options. For grasslands, grazing and mowing modules were already published as part of the developed Biome-BGC  (Hidy et al. , 2012). Since the release of the developed Biome-BGC, we improved the model's ability to simulate management. In BBGCMuSo, the user can define seven different events for each management activity, and additionally annually varying management activities can be defined. Further technical information about using annually varying management activities is detailed in Sect.
Phenology To determine the start and end of the growing season, the phenological state simulated by the model can be used  (White et al. , 1999). We have enhanced the phenology module of the original Biome-BGC, keeping the original logic and providing the new method as an alternative. We developed the so-called heat sum growing season index (HSGSI;  Hidy et al. , 2012) , which is the extension of the GSI index proposed by  Jolly et al. (2005). In the original, as well as in the modified Biome-BGC model versions, snow cover did not affect the start of the vegetation period and photosynthesis. We implemented a new, dual snow cover limitation method in BBGCMuSo. First, the growing season can only start if the snowpack is less than a critical amount (given as millimeters of water content stored in the snowpack).
New C4 photosynthesis routine In the original and modified Biome-BGC models, C 4 photosynthesis was expressed as a sub-version of C 3 photosynthesis. It was implemented in a way that only one parameter differed (photons absorbed by transmembrane protein complex per electron transported; mol mol -1 ) for C 3 and C 4 plants  (Collatz et al. , 1991). Based on the work of Di  Vittorio et al. (2010) , we implemented a new, enzyme-driven C 4 photosynthesis routine into the photosynthesis module. A new ecophysiological parameter (fraction of leaf N in PEP carboxylase) was defined in BBGCMuSo to support the C 4 routine. The ecophysiological parameter "fraction of leaf N in RuBisCO" still affects the process of C 4 photosynthesis (see Di  Vittorio et al. , 2010  for details).