Expression analyses in Ginkgo biloba provide new insights into the evolution and development of the seed

Although the seed is a key morphological innovation, its origin remains unknown and molecular data outside angiosperms is still limited. Ginkgo biloba, with a unique place in plant evolution, being one of the first extant gymnosperms where seeds evolved, can testify to the evolution and development of the seed. Initially, to better understand the development of the ovules in Ginkgo biloba ovules, we performed spatio-temporal expression analyses in seeds at early developing stages, of six candidate gene homologues known in angiosperms: WUSCHEL, AINTEGUMENTA, BELL1, KANADI, UNICORN, and C3HDZip. Surprisingly, the expression patterns of most these ovule homologues indicate that they are not wholly conserved between angiosperms and Ginkgo biloba. Consistent with previous studies on early diverging seedless plant lineages, ferns, lycophytes, and bryophytes, many of these candidate genes are mainly expressed in mega- and micro-sporangia. Through in-depth comparative transcriptome analyses of Ginkgo biloba developing ovules, pollen cones, and megagametophytes we have been able to identify novel genes, likely involved in ovule development. Finally, our expression analyses support the synangial or neo-synangial hypotheses for the origin of the seed, where the sporangium developmental network was likely co-opted and restricted during integument evolution.

www.nature.com/scientificreports/ been found widely expressed in the vegetative tissue, in the petiole of the leaf, the young leaves and the vascular bundles (Fig. 1s,t). No signal was detected with a GibiANT sense probe ( Fig. 1u-x).

Expression analyses of BELL1 Ginkgo homologues: GibiBEL1 and GibiBEL1-2. Expression analy-
ses of the two BELL1 homologues, GibiBEL1 and GibiBEL1-2, in developing ovules show restricted expression patterns for each copy. At S1 and S2, GibiBEL1 is expressed at the base of the ovule (Fig. 2a,b). At S5, GibiBEL1 is expressed in the abscission zone ( Fig. 2c). At S7, when the ovule is competent for fertilization, GibiBEL1 is detected in the pollen chamber (Fig. 2c) and at S8, in the megaspore mother cell and jacket cells once they have formed (Fig. 2d,e). At S10, GibiBEL1 is strongly expressed in the abscission zone, and detected in the nucellus and endotesta cells of the integument (Fig. 2f). No GibiBEL1 expression was found in the young developing pollen cone or in the blade of the young leaf (Fig. 2g,h). No signal was detected with a GibiBEL1 sense probes ( Fig. 2i-l). Unlike GibiBEL1, GibiBEL1-2 is expressed in the nucellus from the early stages of ovule development (S4; Fig. 2m) with this expression restricted to the megaspore mother cell, once it develops (S7-8; Fig. 2n,o). Gibi-BEL1-2 is also expressed in the jacket cells (S10, Fig. 2p) and in the abscission zone (Fig. 2q,r). No expression was detected in the pollen cones or the leaves (Fig. 2s.

t).
Expression analyses of Ginkgo homologues GibiKAN, GibiUCN, GibiUCN2 and GbC3HDZs. Ini-tially, in S2, GibiKAN is expressed throughout the ovule primordia and the funiculus (Fig. 3a). Later, at S3, GibiKAN is expressed in the region that will become the nucellus (Fig. 3b) and it is maintained as the nucellus develops, 5 (Fig. 3c). At this stage GibiKAN is also expressed in the apical region of the integument (Fig. 3c). In S7 and 8, GibiKAN is also expressed in the integument when the integument begins to close the micropyle (Fig. 3d,e). These expression patterns in the integument and nucellus are maintained, and additional expression is detected in the jacket cells at S10 (Fig. 3f). GibiKAN is also expressed in microspores and pollen grains (Fig. 3g). In vegetative tissues, GibiKAN is highly expressed throughout leaf development and its expression does not appear polar (Fig. 3h). Sense probes show no expression (Fig. S10).
The two UNICORN homologues, GibiUCN and GibiUCN2, show low levels of expression throughout ovule development ( Fig. 3i-t). As the integument becomes distinct from the nucellus, GibiUCN is specifically expressed in the apical region of the integument forming the micropyle (Fig. 3j). Both paralogues are strongly expressed in the tapetum and in the nearly mature pollen grains ( Fig. 3o-u). We did not detect expression of either homologue in the blade of young leaves (Fig. 3n,v). Sense probes show no expression (Fig. S10).
From the five paralogues identified for the C3HDZip genes in Ginkgo, GbC3HDZ1 to 5 49 , we were able to assess the expression of four paralogs GbC3HDZ1 to 4 (Fig. 4, Supplementary Fig. S11). At S2, GbC3HDZ1 is expressed in the ovule primordia (Fig. 5a); and in S4 and S6, in the young developing nucellus (Fig. 4b,c) specifically in the region where the megaspore will develop (Fig. 4c). GbC3HDZ1 is expressed in the adaxial side of the integument, in the region that delimits the integument and nucellus (Fig. 4c). This expression is maintained in the adaxial side of the integument and in the jacket cells, S9 (Fig. 4d). At S10, in the fleshy integument, we did not detect expression of GbC3HDZ1 (Fig. 4e), but it is expressed in the tapetum of the microsporangium and in the nearly mature pollen grains (Fig. 4f). GbC3HDZ1 is detected in the leaf and petiole vasculature and appears adaxial in young developing leaves (Fig. 4g). GbC3HDZ1 is not detected in the blade of well-developed leaves (Fig. 4h). GbC3HDZ1 sense probes show no expression (Fig. S10).
In ovules at S4, no expression of GbC3HDZ2 was detected ( Fig. 4i) but in S8, as soon as the megaspore and the jacket cells start to develop, expression is detected (Fig. 4j). This expression is maintained as the ovule matures, S9 (Fig. 4k). Later, at S10 after pollination, GbC3DZ2 is found expressed in the jacket cells (Fig. 4l,m) and throughout the fleshy seed coat (Fig. 4n). GbC3DZ2 expression is detected in the tapetum and the nearly mature pollen grains ( Fig. 4o) and in young developing leaves (Fig. 4p) becoming restricted to the vascular bundles and the adaxial side of the well-developed leaves (Fig. 4q).
GbC3HDZ3 is expressed similarly to GbC3HDZ1 in the young developing nucellus (Fig. 4r), the jacket cells, the tissue that will form the megaspore, the adaxial side of the integument in the region in close contact with the nucellus (Fig. 4t-v), in the tapetum and pollen grains (Fig. 4x), and throughout the vegetative tissue including vascular bundles (Fig. 4y,z). GbC3HDZ4 is expressed at S11 in well-developed ovules in the inner region of the integument, the endotesta (Fig. S11a-e) and the pollen grains (Fig. S11f). No expression of GbC3HDZ4 is detected in the leaf (Fig. S11g).
Transcriptome assembly. A de novo reference transcriptome of Ginkgo was generated from RNAs isolated from young ovule (S4), integument, megagametophyte, collar (dissected from ovules in S9), pollen cone and leaf. Using Trinity software 86,050 transcripts were obtained, with an average GC content of 41.52% and a maximum assembled contig length of 18,726 bases. To improve the quality of the assembly, the contigs were mapped to the initial assembly with ABySS. This gives a final total of 53,970 transcripts (Table 1). Based on read coverage, the E90N50 statistic was ~ 1.8 Kb (Fig. S6), the reference transcriptome contained 86.9% of the conserved Embryophyte genes using BUSCO annotation (Fig. S7).
Samples were compared with a Principal Component Analysis (PCA), which shows that the integument and the megagametophyte are the most dissimilar samples in the data set in terms of gene expression (Fig. S8a). A hierarchical cluster analysis was performed to better understand the similarities within samples. The resulting dendrogram shows that the integument is the most distinct sample with longer branch distance (Y-axis) but it is more similar to the megagametophyte (Fig. S9). Collar, leaf, pollen cone, and young ovules form another cluster (Fig. S9).    Fig. 5a, Table S1). Differentially expressed genes in the integument were filtered by statistical significance (FDR p ≤ 0.05) and a comparison of all tissues against integument was performed. We found that most of the DE genes, that belong to the ovule genetic network seem to be similarly upregulated in all tissues except for GibiANT (Fig. 5b). Subsequently, to focus on genes with a larger change (log2FC ≤ − 2 and ≥ 2), we added a Fold Change threshold which detected 2137 DE genes (Fig. 5c). None of the genes in the ovule regulatory network passed this filter.

Transcription factors differentially expressed in integument.
We focused our transcriptome analyses on transcription factors (TF) which are known to control transcription levels and act as major developmental switches. 134 DE genes were detected as TF and the differential expression of each of these TF within tissues was also compared (Fig. 5d). Of these TFs, compared to other tissues, 21 are found to be largely upregulated in the integument (Table S2) and there are 97 down regulated transcription factors (Table S3). By comparing the results of the samples of young ovules (Fig. S10) and integument, we detected genes that are expressed throughout integument development (from early stages of the ovule to the mature integument) suggesting that there are 13 throughout integument development (Fig. 5e, Table S4).

Differentially expressed FANTASTIC FOUR homologues. Among the 21 transcription factors
upregulated in the integument, the FANTASTIC FOUR (FAF) gene family stood out as they are known to repress WUSCHEL genes in Arabidopsis 45 . To understand the relationships among these transcripts, a detailed phylogenetic analysis of this family of transcription factors was performed. We were able to identify one sequence as a FAF homologue, referred herein as GibiFAF (Fig. 6a). We identified a duplication event occurred before the diversification of angiosperms giving rise to clades FAF1/2 and FAF3/4. In addition, two Brassicaceae-specific duplication events were detected in each clade, resulting in the clades FAF1, FAF2, FAF3 and FAF4 respectively (Fig. 6a). With expression studies in Ginkgo, we found that GibiFAF expression is restricted to the integument throughout S4 to S9 of ovule development (Fig. 6b,c). GibiFAF does not appear to be expressed in the pollen cones or leaves (Fig. 6d,e).

Discussion
Unlike angiosperms, in Ginkgo, we found that the expression patterns of the WUS homologue is not only in the nucellus but also in the integument, pollen cone, and leaf ( Fig. 1a-h). In gymnosperms, the Gnetum homologue (GgWUS/WOX5), exhibits expression patterns like those of monocots, in lateral organ primordia, as well as in the nucellus 48 . In the fern Ceratopteris richardii, CrWOXB a WUS-RELATED homeobox promotes cell divisions in the gametophytic generation and organ development in the sporophytic generation 53 . In land plants, all members of WOX gene lineage are mainly known for their function in meristem identity [20][21][22]54 . However, GbWUS expression is found in the basal region of the integuments. In Ginkgo ovules, the expression patterns we detected could be associated with the meristematic activity of the pachychalaza region of the ovule (Fig. 1a-f). Shifts in the expression patterns of this gene lineage in seed plants may be linked to major morphological differences or to changes in the cis-regulatory regions, as the protein sequence seems highly conserved in seed plants 55 . We did not find GibiANT expression in young developing integuments. However, the expression pattern of ANT varies in ovules of different gymnosperms lineages. In gymnosperms Pinus thumbergii, and Gnetum parvifolium, expression analyses in young developing ovules show expression in the nucellus and integument 47 . In Gnetum gnemon and Ginkgo (Fig. 1m-t), expression is detected only in the micropyle at a pre-pollination stage (Fig. 1n) 13 . ANT in the fern Ceratopteris richardii, CerANT, is expressed in the sperm, in the archegonial neck canal just before fertilization (gametophyte structure), and in the fertilized egg, (i.e., the zygote), and in the fiddlehead (sporophyte) 56 . The expression detected in the pollen grains, suggests that ANT homologues were retained in gymnosperms as key factors in the development of the mega and the microspores, gametophyte development, similar to what is found in ferns (Fig. 1i,n). Overall, the ancestral function of ANT seems to be in cell division as it is present in active cell division regions and in young developing tissue throughout land plants 47,56 .
In Ginkgo, we found GibiBEL1 and GibiBEL1-2 expressed in megaspore and pollen grains (Fig. 2) similar to expression of the only Gnetum gnemon homologue, Melbel1, detected in the nucellus 13 . Loss of function of PpBELL1 in Physcomitrella patens moss generates bigger egg cells unable to form embryos, suggesting that BELL1 has been key in facilitating the diversification of land plants (embryophytes 57 ). This suggests that BEL function in the proper formation of the spores, may be conserved in mosses and gymnosperms (Fig. 2) 57 . Notably,   (Figs. 1d, 2e). We did not find any polar (abaxial) expression of GibiKAN in Ginkgo ovules in particular (Fig. 3a-h). KAN genes are expressed in the micropylar region of the integument in gymnosperms, suggesting differences in the proximal-distal development of these ovules compared to angiosperms (Fig. 3a-h) 13 . In the lycophyte Selaginella moellendorffii, three KAN specific homologues are expressed throughout sporangium development 60 . The expression patterns in the megaspore are conserved between S. moellendorffii and gymnosperms. KAN genes are generally known for their function in establishing abaxial organ polarity in land plants 37,[40][41][42] . This function is likely conserved in ferns 60 and in monocot homologues 61,62 . This allow us to hypothesize that the ancestral function of KAN genes is in the development of the sporangium and that this function is conserved in lycophytes and gymnosperms.
It seems that the abaxial-adaxial polarity function is not conserved in the integument of gymnosperms as UCN homologues are expressed only in the nucellus and apical portion of the integument. Intriguingly, both UCN and KAN homologues in Gnetum gnemon and in Ginkgo, are expressed in the tips of the integuments which suggesting: (1) the interaction between UCN and KAN may be conserved in this region; (2) their function in gymnosperms may be more in establishing the proximal-distal axis; and (3) this indicates major developmental differences between gymnosperm and angiosperms ovules (Fig. 3i-u) 13,45 .
Interestingly, GbC3HDZ1 and 3 are also expressed in the adaxial side of the integument, likely involved in the separation of the nucellus and integument in the pachychalazal region (Fig. 4, Supplementary Fig. S9). Notably, GbC3HDZ1 and 3 expression is not only adaxial in the integument but also at the base of the ovule. In Ginkgo, previous studies revealed expression in the leaf primordia 63 (Fig. 4g,p,y). C3HD-Zip homologues are expressed in the sporangia of the lycophyte Selaginella moellendorffii and the fern Psilotum nudum 49 . In vascular plants C3HD-Zips are involved in vasculature development, also observed here in the stalk 49 . However, the ovules of Ginkgo are not vascularized (Fig. 4, Supplmentary Fig. S9). The data available so far suggests that sporangia development could be the ancestral function of this gene lineage 49 .
The main sources of diversity and changes underlying evolution are alterations in the expression of genes encoding transcriptional regulators 64 . We focused on differentially expressed (DE) genes annotated as transcription factors (Fig. 5c,d, Tables S3, S4, Fig. S8) 65 .
We have identified a gene upregulated in the integument transcriptome related to the FANTASTIC FOUR (FAF; Fig. 6a), a plant-specific gene family with four paralogues in Arabidopsis: FAF1 to 4 51 (Table S2). FAF1 and 2 proteins, are known for their ability to regulate the size of the shoot apical meristem and expression in the embryo (Fig. S12) 51 ; this function in the meristem is linked to its ability to repress WUS 51 . Our Maximum Likelihood analysis shows that there are three duplication events. One before the diversification of all angiosperms giving rise to two clades: FAF1/2 and FAF3/4 corresponding to a whole genome duplication (WGD) event ε 66 . In addition, there is a Brassicaceae-specific duplication event in each of these clades that corresponds to the α and β WGD Table 1. Statistics for Ginkgo reference transcriptome. The initial assembly was improved with a re-assembly method using AbySS.  www.nature.com/scientificreports/ events specific to Brassicales 66 . Gymnosperms are pre-duplication homologues (Fig. 6a). Our expression analyses in Ginkgo indicate that GibiFAF is expressed at higher levels in the integument (Fig. 6b,c) and neither in the pollen cone nor in the leaf (Fig. 6d,e) corroborating the analysis of DE genes (Fig. 6d). It is not yet clear whether FAF directly represses WUS in Arabidopsis as their expression overlaps, but, GibiFAF and GibiWUS expression only overlap in the integument of Ginkgo (Figs. 1, 6) suggesting that, GibiFAF is likely a novel regulator of integument development in Ginkgo. To determine if this function is conserved in other species, further studies are needed. Beyond understanding morphological and developmental patterns of Ginkgo ovule, our results also provides molecular evidence on the origin of the seed.
(1) Expression patterns do not appear to be wholly conserved between angiosperms and gymnosperms (Fig. 7), but the main function of the gene(s) may be conserved. The expanded expression of GbWUS, at base of the ovule and in the basal portion of the integument, indicates that this region has persistent meristematic function. GbWUS expression, additionally, provides molecular support for the interpretation of Ginkgo integument as pachychalazal, where the chalaza domain extends upward from base of the ovule. The expression of GbC3HDZ1 in the adaxial basal region of the integument, indicates that its role in repressing the meristematic activity of GbWUS 12 may have occurred early during the evolution of the seed.  (Table S5) 67 . BEL1 and KAN expression in Ginkgo and Gnetum ovules, are expressed comparably in the nucellus at the sporogenous stage (Fig. 7, Table S5), however, in Gnetum, it occurs prior to pollination, whereas in Ginkgo it occurs during pollination (5) Molecular analyses available in land plants show that integument genes are expressed during sporangia development (in lycophytes, ferns, Ginkgo and Gnetum) suggesting that the integument developmental network was co-opted from a sporangia development network. (6) The outcomes of these studies, together with recent molecular studies, provide additional molecular evidence supporting the synagial/neo-synangial hypotheses, by showing the expression patterns of integument genes in both micro-megasporangium and in the apical region of the integument of Gnetum and Ginkgo 12,13 . Indeed, the data available to date, suggests, that the sporangia development genes were co-opted for the development of the integument and that the integuments have evolved according to the synangial/ neo-synangial hypothesis.
It is enticing to speculate that apical-basal expression patterns reflect the integumentary lobes envisioned in the neo-synangial hypothesis. With WUS in the base of the integument and the nucellus, it is not clear what mechanism accounts for the sterile integument. Recent reports suggest this could be due to BEL1 repression of SPL/NZZ 12 . Future studies of SPL/NZZ homologues in gymnosperms could provide further molecular support for the synangial/neo-synangial origin of the seed.

Methods
Expression analyses by in situ hybridization. The WUSCHEL homologue was previously identified with phylogenetic analysis 42 (GenBank accession number: FM882128). Other homologues were identified with a BLAST amino acid search using Arabidopsis sequences as query (Table S6). Ginkgo sequences were identified from the OneKP database (Table S6; https:// db. cngb. org/ onekp). A BLAST search was performed in the genome as well, but no hits were retrieved (PLAZA database: https:// bioin forma tics. psb. ugent. be/ plaza/ versi ons/ gymnoplaza/). The relationships of these sequences were previously shown with maximum likelihood analyses 13,50 . There are five homologues of Class III HD-Zip (C3HDZ) genes in Ginkgo, which have been previously reported 49 . However, the synthesis of the probe for one of the paralogues, GbC3HDZ5 was not effective; thus, we will present results for GbC3HDZ1 to 4.
Plant material was collected from the NYBG grounds (Accession number: 1353/97) and immediately fixed in FAA (FAA; 3.7% formaldehyde: 5% glacial acetic acid: 50% ethanol). Our characterization of the expression patterns begins around S4 of ovule development for most of these genes (i.e., GbWUS, GibiANT, GibiBEL1-2, GibiUCN, GibiUCN2, and GbC3HDZ2 and 3). This is because collection of ovules at early stages is highly variable as they are covered by the bracts of the short shoots. Only GibiBEL1, GibiKAN and GbC3HDZ1 were assessed starting at S2. After a 4-h incubation in FAA, samples were dehydrated in a standard ethanol series, then transferred to fresh Paraplast. The samples were sectioned on a Microm HM3555 rotary microtome. DNA templates for RNA probe synthesis were obtained by PCR amplification of 280-480 bp fragments. To ensure specificity, the probe templates were designed outside of conserved domains (Fig. S2, Table S7). Sense probes were used as negative controls. The fragments were cleaned using QIAquick PCR purification Kit (Qiagen, Valencia, CA, USA). Digoxigenin labeled RNA probes were prepared using T7 polymerase (Roche, Switzerland), murine RNAse inhibitor (New England Biolabs, Ipswich, MA, USA) and RNA labeling mix (Roche, Switzerland) according to the protocol of each manufacturer. The RNA in situ hybridization was performed according to Ambrose et al. 68 (Fig. S3). Raw data quality was assessed using FastQC (v 0.11.5; Andrews, 2010). Sequence adapters and low-quality reads (Phred score < 5) were removed using Trimmomatic (v 0.36) with all default parameters 69 . Transcripts were assembled using AbySS (v 2.0.2) 70 and the Trinity (v 2.8.4) software pipeline 71 for comparison (Fig. S4). Because of better statistics, we continued to work with the Trinity assemblies (Fig. S5). An initial reference transcriptome was assembled de novo from all RNA samples and all contigs ≥ 200 nucleotides length. The quality of the transcriptome assembly was assessed based on the calculated E90N50 contig length (E90N50 ~ 1.8 Kb; Fig. S6). The initial reference transcriptome was annotated using DIAMOND (v 0.9.13) 72 . To identify possible contaminants, Ginkgo contigs were searched against bacterial and fungal databases mainly associated with soil and plants, sequence databases compiled from UniProt (www. unipo rt. org). Sequences with an identity ≥ 50% were removed from the reference transcriptome (N = 2656). This initial transcriptome was re-assembled to improve the assembly stats using AbySS 70 , the quality of the transcriptome was assessed with contig length and BUSCO annotation (Fig. S7) 73 , the resulting assembly was used for the following steps. Long open reading frames (ORF) were predicted using TransDecoder (v 3.0.0) 71 . For gene annotation, the contigs of Ginkgo were searched in several databases of sequence coding land plant proteins (Amborella trichopoda: AMTR1.0_13333, Arabidopsis thaliana: TAIR10_3702, Capsicum annuum: ASM51225v2, Ginkgo biloba: NCBI:txid3311, Gnetum montanum: NCBI:txid3381, Oryza sativa: IRGSP-1.0, Picea abies: NCBI:txid3329, Selaginella moellendorfii: v1.0_88036, Vitis vinifera: 12X_29760; available through Ensembl and Plaza for gymnosperms; Table 1).
To interpret the overall structure of these samples in terms of the gene expression, a Principal Component Analysis (PCA) was performed using the normalized TPM values, as it allows to better interpret the variation of high-dimensional interrelated dataset (with high number of variables) and to detect major differences between samples. PCA was performed using the Python packages: sklearn, seaborn, and bioinfokit (v 2.0.2) Thus, to better understand the similarities within samples a dendrogram was obtained by performing a hierarchical clustering of the samples using a 'complete linkage' method (Fig. S8). Dendrogram was obtained using the SciPy package on Python (v 1.5.0).

Transcriptome abundance (RSEM) and expression level (EBSeq) analyses.
These analyses were carried out following the pipeline previously proposed 74 . Sequence reads from the different plant tissues were aligned to the reference transcriptome using Bowtie2 (v 2.4.2) 75 and RSEM (RNA-Seq by Expectation Maximization; v 1.3.0) was used to obtain estimates of transcripts abundance for all transcripts 76 . The resulting expression levels are calculated in terms of Transcripts Per Million (TPM). Transcripts were considered to be differentially expressed between integuments and the other tissues, when TPM was ≥ 0.95 for at least a single tissue and the fold change (log2FC) was ≤ − 2 and ≥ 2 with an FDR p ≤ 0.05 (Fold Discovery Rate). To identify the corresponding Gene Ontology (GO) terms, the differentially expressed genes were further analyzed with Blast2GO (v 5.2.5; Fig. S9). Data analyses and results were plotted using Matplotlib v 3.4.2 and Seaborn v 0.8.1 Python libraries (Fig. S4).
Identification of Ginkgo homologues and maximum likelihood analyses for gene lineages of interest. One of the genes found in the transcriptome analyses to be putatively involved in integument development in Ginkgo is similar to the Arabidopsis gene FANTASTIC FOUR 3 (FAF3; AT5G19260). To reconstruct the evolution of the FANTASTIC FOUR gene family, we used the four Arabidopsis paralogues (AT4G02810, AT1G03170, AT5G19260, AT3G06020) as a query to perform an amino acid BLAST search in seed plants, using the Phytozome and OneKP databases. A total of 88 sequences were compiled and aligned using the online version of MAFFT (v 7) 77 . Three Selaginella sequences were used as outgroups to root the tree (LGDQ_scaffold_2012011; JKAA_scaffold_2181098; ZFGK_scaffold_2040141).
Phylogenetic analyses using the nucleotide sequences were performed with RaxML-HPC2 BlackBox 78 . The newly isolated sequence was deposited in GenBank (accession OK255713).

Data availability
The data underlying this article are available in the GenBank Nucleotide Database with accessions provided in the methods and supplemental material. Additional data underlying this article is available upon request to the corresponding author.