A genomic atlas of human adrenal and gonad development

Background: In humans, the adrenal glands and gonads undergo distinct biological events between 6-10 weeks post conception (wpc), such as testis determination, the onset of steroidogenesis and primordial germ cell development. However, relatively little is currently known about the genetic mechanisms underlying these processes. We therefore aimed to generate a detailed genomic atlas of adrenal and gonad development across these critical stages of human embryonic and fetal development. Methods: RNA was extracted from 53 tissue samples between 6-10 wpc (adrenal, testis, ovary and control). Affymetrix array analysis was performed and differential gene expression was analysed using Bioconductor. A mathematical model was constructed to investigate time-series changes across the dataset. Pathway analysis was performed using ClueGo and cellular localisation of novel factors confirmed using immunohistochemistry. Results: Using this approach, we have identified novel components of adrenal development (e.g. ASB4, NPR3) and confirmed the role of SRY as the main human testis-determining gene. By mathematical modelling time-series data we have found new genes up-regulated with SOX9 in the testis (e.g. CITED1), which may represent components of the testis development pathway. We have shown that testicular steroidogenesis has a distinct onset at around 8 wpc and identified potential novel components in adrenal and testicular steroidogenesis (e.g. MGARP, FOXO4, MAP3K15, GRAMD1B, RMND2), as well as testis biomarkers (e.g. SCUBE1). We have also shown that the developing human ovary expresses distinct subsets of genes (e.g. OR10G9, OR4D5), but enrichment for established biological pathways is limited. Conclusion: This genomic atlas is revealing important novel aspects of human development and new candidate genes for adrenal and reproductive disorders.


Abstract
: In humans, the adrenal glands and gonads undergo distinct Background biological events between 6-10 weeks post conception (wpc), such as testis determination, the onset of steroidogenesis and primordial germ cell development. However, relatively little is currently known about the genetic mechanisms underlying these processes. We therefore aimed to generate a detailed genomic atlas of adrenal and gonad development across these critical stages of human embryonic and fetal development.
: RNA was extracted from 53 tissue samples between 6-10 wpc Methods (adrenal, testis, ovary and control). Affymetrix array analysis was performed and differential gene expression was analysed using Bioconductor. A mathematical model was constructed to investigate time-series changes across the dataset. Pathway analysis was performed using ClueGo and cellular localisation of novel factors confirmed using immunohistochemistry.
: Using this approach, we have identified novel components of adrenal Results development (e.g. , ) and confirmed the role of as the main ASB4 NPR3 SRY human testis-determining gene. By mathematical modelling time-series data we have found new genes up-regulated with in the testis (e.g. ), SOX9 CITED1 which may represent components of the testis development pathway. We have shown that testicular steroidogenesis has a distinct onset at around 8 wpc and identified potential novel components in adrenal and testicular steroidogenesis (e.g. , , , , ), as well as testis MGARP FOXO4 MAP3K15 GRAMD1B RMND2 biomarkers (e.g. ). We have also shown that the developing human SCUBE1 ovary expresses distinct subsets of genes (e.g. , ), but OR10G9 OR4D5 enrichment for established biological pathways is limited.
: This genomic atlas is revealing important novel aspects of human Conclusion development and new candidate genes for adrenal and reproductive disorders.

Introduction
The development of the adrenal gland and gonads (testes, ovaries) are two of the most important embryological events, but relatively little is known about the exact mechanisms of these processes in humans.
Development of these structures is closely linked as they arise from a shared region of intermediate mesoderm in early embryogenesis and subsequently have common biological processes, such as the ability to synthesise steroid hormones (see Figure 1 for an overview).
The primordial adrenal gland arises as a distinct structure at around 28 days post-conception (dpc) and undergoes rapid growth during late embryonic and early fetal life (Ishimoto & Jaffe, 2011). The adrenal cortex develops into definitive and fetal zones that can produce steroid hormones (cortisol, aldosterone) and large amounts of adrenal androgens, such as dehydroepiandrosterone (DHEA) (Miller & Auchus, 2011;Xing et al., 2015). Chromaffin cells derived from the sympathetic nervous system migrate into the developing adrenal gland and later coalesce to form the adrenal medulla, which secretes adrenaline and noradrenaline ( Figure 1). It has been proposed that the entire hypothalamic-pituitary-adrenal (HPA) axis is functionally intact during early fetal life (Goto et al., 2006).
In contrast, the developing gonad remains "bipotential" until around 41 dpc as it can differentiate into either a testis or ovary ( Figure 1) (Bashamboo et al., 2017;Svingen & Koopman, 2013). The bipotential gonad and adrenal gland both arise from the adrenogonadal primordium around CS13. Sex determination in the bipotential gonad occurs between CS17 and CS23 when transient SRY expression promotes the upregulation of genes, including SOX9, in the developing testis. Sex differentiation starts with the onset of steroidogenesis around CS23, which results in development of the external genitalia and phallic growth (right panel). Ovary-specific genes are thought to suppress male pathway specification. The ovary also supports germ cell expansion and entry into meiosis. The fetal adrenal cortex develops into definitive and fetal zones, producing steroid hormones and adrenal androgens. Chromaffin cells derived from the sympathetic nervous system migrate into the developing adrenal gland, merging later to form the adrenal medulla.

Amendments from Version 1
We have undertaken further qRTPCR validation of adrenal gene expression and include these data as a revised Figure 5b. We also include higher power immunohistochemistry of FOXO4 and NR5A1 in the nuclei of fetal testis Leydig cells as a revised Figure 12d. We include a dataset of differential gene expression in the testis versus ovary (Dataset 7, OSF). We now include three Supplementary Figures of further validation data to show serial upregulation of novel genes in the testis with time (qRTPCR, Supplementary Figure 1); upregulation of novel steroidogenic components with time (qRTPCR and immunohistochemistry, Supplementary Figure 2); higher powered adrenal expression data for GRAMD1B and RMDN2 in the adrenal gland as well as additional negative control data (immunohistochemistry, Supplementary Figure 3). Finally, we have extended the discussion to include a paragraph outlining several limitations of the study and possible future approaches.

REVISED
In the developing 46,XY embryo, transient expression of the testis-determining gene, SRY, on the Y chromosome is believed to lead to upregulation of genes such as SOX9 and the process of testis determination (Bashamboo et al., 2017;Hanley et al., 2000;Hiramatsu et al., 2009;Koopman et al., 2016;Sekido & Lovell-Badge, 2008). At around 7-8 weeks post conception (wpc), Sertoli cells secrete the peptide hormone Anti-Müllerian Hormone (AMH, also known as Müllerian Inhibiting Substance, MIS) to regress Müllerian structures (uterus and upper vagina). Fetal Leydig cells develop and synthesise testosterone. Testosterone stabilises Wolffian structures and stimulates development of the external genitalia through its peripheral conversion to dihydrotestosterone ( Figure 1) (Miller & Auchus, 2011).
Ovary development in the 46,XX embryo was originally thought to be a "passive" process that occurs in the absence of SRY and other components of the testis-determining pathway. Although the morphological changes in the ovary are less striking, data from mice and from 46,XX individuals with rare forms of ovotestis development suggest that distinct pathways are up-regulated (e.g. RSPO1/WNT4/β-catenin), which may function largely to repress testis development ( Figure 1) (Bashamboo et al., 2017;Beverdam & Koopman, 2006;Bouma et al., 2007;Jameson et al., 2012;Nef et al., 2005). The ovary also has expansion of primordial germ cells starting towards the end of this early fetal period (Fowler et al., 2009).
Given the limited data currently available (Gerrard et al., 2016;Houmard et al., 2009;O'Shaughnessy et al., 2007;Ostrer et al., 2007), we aimed to generate a genomic atlas of human adrenal and gonad development between 6 and 10 wpc. This period represents a critical window for organogenesis during which many key biological events occur in tissue specification and cellular differentiation. This unique dataset is supported by our current knowledge of key genetic components in these processes and by human conditions affecting adrenal and gonadal function. By comparing datasets at different times and in different tissues, as well as by mathematical modelling of time series data, several novel insights into human adrenal and gonad development are emerging.

Tissue samples
Human embryonic and fetal tissue samples used in this study were obtained in collaboration with the MRC/Wellcome Trustfunded Human Developmental Biology Resource (HDBR). The HDBR is a tissue bank regulated by the Human Tissue Authority. Samples were collected with appropriate maternal written consent and with approval from the NRES Committee London-Fulham (REC reference 08/H0712/34+5). Material was obtained only from surgical terminations of pregnancy to avoid any potentially interfering effects of anti-progestogens on nuclear receptor gene transcription.
All materials were kept on ice and processed rapidly after collection. The adrenal glands, gonads and any control tissues needed were identified and removed by blunt dissection under a dissecting microscope. Samples were immediately immersed in RNAlater® (Ambion, Austin, TX, USA) and stored at -20 C. The age of the embryonic or fetal sample was calculated based on morphological characteristics using published guidelines (Hern, 1984;O'Rahilly & Müller, 1987) by a single experienced researcher (DG) (Supplementary Table 1). Each tissue sample was also karyotyped using standard G-banding, to determine the sex of the embryo/fetus and also to ensure there were no significant aneuploidies or chromosomal rearrangements present.
A total of 53 different organs between 42 days post conception (dpc) (Carnegie Stage 17) and 10 weeks post conception (wpc) (Fetal Stage 3) were included in the analysis. This dataset consisted of 17 adrenal glands, 20 testes and 10 ovaries, as well as 6 "control" samples (46, XY) from spine, brain, muscle, heart, kidney and liver. The adrenal gland and testis samples came from the same fetus in 11 cases across the age range. Control samples were obtained from 2 of these. An overview of all samples is provided in Supplementary Table 1.

RNA extraction
RNA was extracted using the TRIzol method. In brief, the tissue sample was removed from RNAlater® and homogenized in 1ml TRIzol reagent using a Pellet Pestle micro-grinder (Kontes/ Kimble Chase, USA). Chloroform (0.2 ml) was added to allow phase separation, and RNA was precipitated from aqueous phase by addition of 0.5 ml isopropyl alcohol. The RNA pellet was washed in 1ml ice-cold 75% ethanol and air dried before resuspension in nuclease-free water. The concentration and A260/A280 ratio was measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Witec, Littau, Switzerland). RNA was stored at -80 C prior to use.

Arrays
Expression assay was performed using the Affymetrix GeneChip Human Gene 1.0 ST array, covering approximately 36000 RefSeq transcripts. In brief, RNA concentration and quality was assessed using an Agilent Bioanalyzer. 250ng of total RNA was converted to ssDNA using the Ambion Whole Transcript Sense Target Labeling Assay and labelled and prepared for hybridisation using the Affymetrix Genechip WT Terminal Labelling and Hybridisation kits, as per the manufacturer's protocol. Hybridisation was performed for 16 hours at 45 C. Arrays were washed and stained using the Affymetrix Fluidics Station 450 and scanned on the Affymetrix Genechip 3000S scanner as per the manufacturer's instructions.
After passing quality control, CEL files were normalised using the Robust Multi-array Average (RMA) algorithm (oligo package in R, version 1.38.0). Adjustments for batch effects was performed using ComBat (sva package, version 3.22.0) (Leek et al., 2012). One adrenal sample that was an extreme outlier in several QC platforms was excluded prior to further analysis (leaving the 17 adrenal samples reported above). Statistical analysis of the microarray data was carried out in R programming environment (version 3.3.3 running under Ubuntu 16.04.2) using Bioconductor packages 'oligo' (Carvalho & Irizarry, 2010) and 'limma' (version 3.30.12) (Smyth, 2004). The hierarchical cluster dendrogram was generated using the agglomeration method ward.D2.
Differentially expressed genes were obtained using linear models and the Benjamini and Hochberg method was applied to adjust the P-values for multiple testing (Benjamini & Hochberg, 1995). The cut-off values for adjusted P-values and log 2 fold-changes were 0.05 and 1 or 2, respectively. The pheatmap (version 1.0.8) and VennDiagram packages (version 1.6.17) were used to generate the heatmaps and Venn diagrams, respectively. Pathway enrichment analysis was carried out using the Cytoscape plug-in ClueGO (version numbers 3.4.0 and 2.3.3, respectively) (Bindea et al., 2009).
Array data are available from the ArrayExpress database under the accession number E-MTAB-5525.
A model of expression changes over time A phenomenological model ("BALT model") was developed to pinpoint the time when gene expression levels change, as well as evaluating the speed and magnitude of this transition. The general form of the model is a sigmoid function relating the transcript concentration X (on a 2-log scale) to the time elapsed since conception (t, expressed in days): In the formula above, B ("basis") is the basic expression level at the start of the time course and A ("amplitude") is the amplitude of the difference between the finishing expression level and the starting one. For down-regulated genes this latter parameter takes negative values. The parameter L ("localisation") captures the inflexion point of the sigmoid, and hence can identify at which point in time a gene undergoes a change of expression level. Finally, T ("Transition") is a measure of the speed of this transition and the sub-function f was chosen so that this parameter can also be expressed in days.
The model was fitted to the data for each gene using a least square objective function (maximum likelihood) with a gradient following algorithm (Levenberg-Marquardt). With this, the values for the aforementioned parameters that best fit the data were obtained as well as an evaluation of the goodness of fit (and thus, of the validity of the model). Furthermore, robustness of the parameters were estimated using a Markov Chain Monte Carlo procedure.
It was found that most of the differentially expressed genes in our biological system could be well described by this model or simplified versions of it. Hence the procedure was seen as a way to compress high dimensional data into a smaller set of easily interpretable indicators. This modelling framework allowed genes to be grouped according to their transition times between different expression levels in a developmental process as well as the manner of this transition. This model was used predominantly to investigate time course data in the testis. The two main groups of genes that were analysed were those that followed a curvilinear upregulation in expression from the start of the dataset similar to SOX9, as well as genes that showed a sigmoidal upregulation in expression with the onset of steroidogenesis. The mid-point (in dpc) of the Carnegie or Fetal Stage was used when analysing data for the model (Supplementary Table 1).

Quantitative RT-PCR
Purified RNA from adrenal glands or testis and control tissue (liver, brain, muscle) was quantified using the NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific). First-strand cDNA was synthesised using SuperScript II Reverse Transcriptase (Invitrogen) and random primers according to the manufacturer's instructions. Taqman® probes were obtained for the following genes of  interest: MAP3K15, ASB4, TDGF1, FOXO4, NRK, CITED1

Results
Generation of an atlas of gene expression in human adrenal gland, testis and ovary Human embryonic and fetal tissues of interest were obtained following ethical approval and with informed consent from the Human Developmental Biology Resource (HDBR, www.hdbr. org). RNA was extracted and gene expression profiling performed using Affymetrix ST1.0 Gene arrays. A total of 53 different organs between 42 days (6 weeks) and 70 days (10 weeks) postconception were included in the final analysis (Supplementary Table 1). This dataset consisted of 17 adrenal glands (46,XY), 10 ovaries (46,XX), and 20 testes (46,XY), as well as 6 control samples (46,XY) from spine, brain, muscle, heart, kidney and liver. Adrenal and control samples with a 46,XY karyotype were intentionally chosen so that subtle effects of Y-chromosome genes on testis determination and sex development could be analysed.
All relevant CEL files have been deposited in the ArrayExpress database (E-MTAB-5525) and the entire dataset of log 2 normalised values is included in Dataset 1.
Initial comparisons of gene expression profiles in the different samples were undertaken using principal components analysis (PCA) and by generating correlation plots (heatmaps) and cluster dendrograms (Figure 2a, Figure 3). Samples from each tissue clearly clustered together, with more similarity between testis and ovary profiles than adrenal gland (Figure 2a, Figure 3).
Global gene expression patterns reveal known diseasecausing genes and identify novel tissue-specific factors An analysis of global gene expression patterns was initially performed for the three tissues (adrenal, testis, ovary) compared to controls. Genes were identified that were up-regulated with a log 2 fold change (FC) of ≥ 2 or ≥ 1 and adjusted P-value lower than 0.05 (P≤ 0.05) (Dataset 2-Dataset 4).
An overview of the number of differentially up-regulated genes that are specific to each tissue, or overlapping, is shown in Figure 2b. As expected, there was greater overlap between testis and adrenal (due to shared steroidogenic pathways) and testis and ovary (due to their gonadal origin and presence of germ cells) than between adrenal and ovary, which are developmentally and functionally more distinct.
Analysis of factors known to be involved in steroidogenesis showed that shared components of the steroidogenic pathways overlapped in the adrenal and testis datasets (e.g. CYP11A1, HSD3B2, CYP17A1), whereas the gene responsible for conversion of delta-4 steroids (e.g. androstenedione) into testosterone, HSD17B3, was the only testis-specific steroidogenic enzyme ( Figure 2c). Adrenal-specific steroidogenic components included the enzymes needed to make aldosterone and cortisol from precursors (e.g. CYP11B1, CYP21A2) as well as adrenal-specific receptors involved in signalling pathways (e.g. MC2R [also known as the ACTH receptor], MRAP).
An analysis of known transcription factors ( Figure 2d) revealed several key genes shared between the testis and ovary involved in gonad development (e.g., WT1, EMX2, LHX9, GATA4) or in germ cell maturation (e.g., POU5F1 [also known as OCT4] and NANOG) ( Figure 2d). One key transcription factor shared between the adrenal and testis was the nuclear receptor NR0B1 (also known as DAX-1), responsible for X-linked adrenal hypoplasia congenita with male infertility (Suntharalingham et al., 2015). Another key nuclear receptor common to all three tissues was NR5A1 (also known as steroidogenic factor-1), responsible for adrenal insufficiency, testicular dysgenesis and primary ovarian insufficiency (Achermann et al., 1999;Lourenço et al., 2009;Suntharalingham et al., 2015). Therefore, the dataset seemed validated based on known factors involved in development as well as in genetic conditions in humans.
Known and novel adrenal-specific genes Initial analysis focused on genes up-regulated in the adrenal gland compared to controls ( Figure 4a, Table 1 Table 1). These genes are known core components of adrenal steroidogenesis and most are associated with monogenic disorders in humans ( Table 1). The most strongly differentially expressed transcription factors were FOXO4 (absolute FC 10), NR0B1 (DAX-1) (absolute FC 9.5) and NR5A1 (SF-1) (absolute FC 7.8) (Figure 5a). Key novel genes or those that have not been extensively studied in the adrenal gland before include MGARP (also known as OSAP), MAP3K15, ASB4, NRK, and TDGF1 ( Figure 4d). Expression of these genes was confirmed by qRTPCR and/or immunohistochemistry (Figure 4c, Figure 5b). Analysis of the heatmap for the top 50 differentially expressed genes revealed that several novel genes appear to be highly adrenal-specific (e.g., ASB4, NPR3, GRB14), whereas others are also expressed in the later testis samples, suggesting a potential novel role for these genes in steroidogenesis (see below) ( Figure 4d).
SRY is the primary driver of human testis development Studies in transgenic mice and humans with SRY deletions and translocations have suggested that SRY is the primary Y-chromosomal testis-determining gene in humans, but direct evidence is limited (Bashamboo et al., 2017;Hanley et al., 2000;Koopman et al., 2016). SRY is thought to be transiently up-regulated in the 46,XY bipotential gonad at around 42 dpc, leading to downstream expression of testis development pathways (Hanley et al., 2000). Analysis of SRY expression in early testis samples (CS18-CS19) compared to control tissues (all 46,XY) showed that SRY is the only significantly differentially-expressed protein coding Y-chromosomal gene (absolute FC 2, P-value≤0.003) ( Figure 6a and Dataset 5). Although changes are subtle, there was a clear transient expression of SRY that had occurred by the start of our dataset (around 45 dpc), with a reduction in intensity around CS23 (57 dpc) ( Figure 6a, Figure 6b).
Other genes that are differentially expressed in the early testis compared to controls are shown in the heatmap in Figure 6c, Figure 7a and Figure 7b. These include the gene encoding renin (REN) as well as mesonephric factors such as NPHS2 and WT1. Pathway-enrichment analysis using ClueGO confirmed "sex differentiation" and "renal system development" as the two major up-regulated pathways in the dataset ( Figure 7c).
Modelling dynamic changes in the testis downstream of SRY Studies in mice have shown that SOX9 is a target of SRY, and disruption of SOX9 in humans results in testicular dysgenesis (Bashamboo et al., 2017;Bishop et al., 2000;Sekido & Lovell-Badge, 2008). However, the dynamics of SOX9 expression in human testis development are still not well understood and it is not established whether SOX9 is the only potential SRY target (Hanley et al., 2000;Ostrer et al., 2007). Analysis of SOX9 across the testis series revealed an approximate 2-fold increase in gene expression levels with a curvilinear pattern ( Figure 8a).
In order to assess detailed patterns and changes in gene expression across the time series dataset, and to be able to group genes with similar dynamic expression patterns together, a phenomenological   (b) Confirmation of differential expression of several novel genes in the adrenal gland (9 wpc) by qRT-PCR. Pooled liver, brain and muscle were used as a control. GAPDH was used as a housekeeping gene.
in adult testis but also strongly in at least one other system (e.g., ANKRD18A, SLC52A3, KEL) (Human Protein Atlas). In addition, INHBB encodes the β-subunit of inhibin B, a well-established testis protein, which together with SOX9 provides validation for the dataset (Valeri et al., 2013).
Identifying known and novel components of steroidogenesis Using the "BALT" model across the time-series dataset, the expression of key known components of testis steroidogenesis showed marked upregulation around a fixed time point ("L") between 54-57 dpc (LHCGR, STAR, CYP11A1, HSD3B2, CYP17A1, HSD17B3) ( Figure 9a). Expression levels of all these enzymes showed a sigmoid-type pattern, suggesting that discrete biological events occur in this fixed window (Supplementary Table 3). The heatmap of normalised gene expression plots for these enzymes showed similar discrete changes between early testis (CS18 to CS22) and late testis (F1 to F3) for many known components involved in testosterone synthesis ( Figure 9b). Principal components analysis also demonstrated segregation of samples into two distinct groups ( Figure 10a). A subset of genes was seen to be up-regulated with age in a correlation plot (Figure 10b), and pathway-enrichment analysis strongly reflected processes involving "cholesterol", "sterol" and "steroid" metabolism or biosynthesis ( Figure 10c). Therefore, CS23 (8 weeks) was a discrete point around which the onset of testicular steroidogenesis could be modelled, and samples either side of this time point were used for "pre-steroidogenesis" ("early") versus "poststeroidogenesis" ("late") comparisons. An overview of key genes in this dataset, many of which are associated with monogenic disorders in humans, is shown in Table 3.
In order to identify potential novel core components of steroidogenesis in the adrenal gland and testis, genes were identified that were both up-regulated in late testis samples compared to early ones (post-steroidogenesis versus pre-steroidogenesis, log 2 FC≥1), and up-regulated in adrenal samples versus controls (adrenal versus controls, log 2 FC≥2) ( Figure 11a). A total of 45 overlapping genes were found, including all relevant known steroidogenic enzymes and 17 novel genes (Figure 11a-c, Table 4). Enrichment-pathway analysis showed upregulation of genes involved in "cholesterol transport and metabolism" as well as "steroid biosynthesis", further validating the analysis ( Figure 11d). The heatmap for these 45 genes confirmed highly selective expression in the fetal adrenal gland and in the fetal testis, but largely after CS23 ( Figure 11c).
An overview of the proposed function of the 17 potentially novel factors is shown in Table 5. These genes include MGARP (also known as OSAP), which may be involved in mitochondrial trafficking, MAP3K15 (ASK3) involved in cell signalling of stress, and the imprinted gene DLK1 (PREF1) implicated in growth and differentiation (Guasti et al., 2013;Naguro et al., 2012;Zhou et al., 2011). Two genes that have not been studied previously are GRAMD1B and RMDN2 (FAM82A1). These genes are clearly expressed with the onset of steroidogenesis (Figure 11b,c) mathematical model was developed that allowed expression points of any gene to be described based on the basal gene expression value ("B"), the amplitude of the change ("A"), the localisation time of the maximum transition ("L"), and the rate of transition ("T") ("BALT model") ( Figure 8b) (See methods for details).
Using this approach, several genes were found to have similar expression dynamics to SOX9 (Figure 8c, Figure 8d, Table 2 and Supplementary Table 2). Analysis of the heatmap for these factors in the entire dataset showed that a subset of these genes had predominantly testis-specific upregulation ( Figure 3a). However, currently the true biological function of these genes remains unknown. The transcription factor FOXO4 showed expression in the nuclei of Leydig cells in the fetal testis ( Figure 12d).

Identification of potential secreted proteins in the testis
The analysis of genes that showed progressive upregulation in the developing testis was extended to try to identify novel biomarkers. The protein structures encoded by genes that were up-regulated in late testis samples (F1 to F3; FC≥1, Dataset 6) were reviewed for characteristics of secreted proteins, such as the presence of a cleaved signal peptide. Using this approach, the three major testis biomarkers Anti-Müllerian Hormone (AMH/MIS), inhibin α-subunit (INHA) and insulin-like 3 (INSL3) were identified (Figure 13a-c). Novel factors identified included EGFLAM, CARTPT, ADAMTS5, SCUBE1 and EPPIN ( Figure 13a). Expression of SCUBE1 was shown in seminiferous cords in the developing testis ( Figure 13d).
Ovarian development is not a passive process Although ovary development was once thought to be a less biologically active process, several studies in the mouse have shown changes in discrete sets of differentially expressed genes in the early ovary (Beverdam & Koopman, 2006;Jameson et al., 2012;Nef et al., 2005). However, it is unclear whether similar effects are also present during early human gonad development. By studying up-regulated genes in the ovary and comparing them to control samples, (Dataset 4) and by studying up-regulated genes in the testis and comparing them to controls (Dataset 3), very similar numbers of gonad-specific differentially expressed genes were found (log 2 FC≥ 1, ovary 274 versus testis 280; log 2 FC≥2, ovary 69 versus testis 57) ( Figure 14a). As expected, there was substantial overlap in genes that were co-expressed in both the ovary and the testis (log 2 FC≥1, 463; log 2 FC≥2, 68). Although the number of strongly up-regulated genes was somewhat higher in the testis compared to the ovary when observed in a volcano plot (Figure 14b and Dataset 7), in part due to Y chromosomal genes, a similar scatter of up-and down-regulated genes was seen for the ovary versus control as for the testis versus control, again reinforcing the concept of the developing ovary having discrete genetic activity (Figure 14c).
Analysis of the top 20 up-regulated genes in the ovary compared to control data is shown in the heatmap in Figure 14d. Many of these genes were also expressed to some extent in the testis, and may represent factors involved in pluripotency and germ cell development such as HIST1H2BA, NANOG and POU5F1 (also known as OCT4). Expression of several genes appeared to increase in the ovary across the time course (CS18 to F2) and to decrease in the testis, such as OR10G9, GABRG1, OR4D5, NPNT, CNT-NAP4 and NPY (Figure 14d). These genes could potentially encode novel components of an ovary-specification program. However, pathway-enrichment analysis of the 274 genes differentially upregulated in the ovary compared to controls (log 2 FC≥1, P<0.05) resulted in very few biological processes other than "germ cell development" and "type 1 interferon signalling pathway"   Figure 15a). In contrast, pathway-enrichment analysis for the 280 genes up-regulated in the testis samples versus controls (log 2 FC≥1, P<0.05) revealed many more detailed networks, with enrichment in terms such as "androgen" and "reproductive" processes ( Figure 15b). Taken together, these data suggest that current knowledge of the genetic events in the early stages of human ovary development is still very limited.

Discussion
Adrenal development and gonad development are two of the most fundamental biological processes. However, current knowledge about the genetic mechanisms underlying these events is derived largely from studies in mice (Beverdam & Koopman, 2006;Inoue et al., 2016;Jameson et al., 2012;McClelland et al., 2015;McClelland & Yao, 2017;McDowell et al., 2012;Nef et al., 2005;Xing et al., 2015) and few data are currently available from humans (Fowler et al., 2009;Gerrard et al., 2016;Houmard et al., 2009;O'Shaughnessy et al., 2007), especially during the first trimester or across all three tissues. Of the published data, the recent study from Gerrard et al. provides the most detailed insight into early human adrenal and testis expressed genes currently available, but this work included just two samples from pooled RNA between CS18 and CS22 (Gerrard et al., 2016). Given the potential biological differences between species, as well as a lack of detailed time course data for all three tissues, our aim was to generate a detailed "Atlas" of genomic events across a critical period in human embryonic and fetal development (6 to 10 wpc).
Initial analysis focused on global gene expression patterns in the adrenal, testis and ovary compared to control tissues. Control      (a) Immunohistochemistry for GRAMD1B in human fetal testis at 9 wpc. NR5A1 (SF-1) was used to highlight Leydig cells (green). DAPI was used to counterstain nuclei (blue) and to highlight the outer capsule. Scale bars, 50 µm (top panels) and 20 µm (bottom panels). (b) Immunohistochemistry of RMDN2 performed as above. Scale bars, 50 µm (top panels) and 20 µm (bottom panels). (c) Immunohistochemistry of GRAMD1B and RMDN2 in the fetal adrenal gland at 9 wpc. NR5A1 (SF-1) was used to highlight the definitive zone and fetal zone cells (green). DAPI was used to counterstain nuclei (blue) and to highlight the outer capsule. Scale bars, 100 µm. (d) Immunohistochemistry for FOXO4 in human fetal testis at 11 wpc. NR5A1 (SF-1) was used to highlight Leydig cells (green). DAP-I was used to counterstain nuclei (blue). Scale bar, 100 µm (top panels) and 20 µm (bottom panels). (b) Immunohistochemistry for AMH (MIS) in human fetal testis at 9 wpc. NR5A1 (SF-1) was used to highlight Leydig cells (green). DAPI was used to counterstain nuclei (blue) and to highlight the outer capsule. Scale bars, 100 µm (top panels) and 20 µm (bottom panels) (c) Immunohistochemistry for INSL3 in human fetal testis at 9 wpc. NR5A1 (SF-1) was used to highlight Leydig cells (green). DAPI was used to counterstain nuclei (blue). Scale bar, 50 µm (d) Immunohistochemistry for SCUBE1 in human fetal testis at 9 wpc. NR5A1 (SF-1) was used to highlight Leydig cells (green). DAPI was used to counterstain nuclei (blue). Scale bars, 100 µm (top panels) and 20 µm (bottom panels).  samples were chosen to represent a range of developmental tissues (neuroectoderm, mesoderm, endoderm) across the study period. Approximately 0.5-1% of all genes were differentially expressed in each tissue at log 2 FC≥2 (absolute fold change 4) and 3-5% of all genes at log 2 FC≥1 (absolute fold change 2). The strong overlap between the adrenal and testis partly reflected shared components of steroidogenesis, whereas the overlap between testis and ovary reflected the shared origins of these tissues as well as the presence of germ cells.
Initial validation of our dataset came from known genes involved in distinct biological processes or monogenic disorders in humans. For example, two of the most strongly expressed transcription factors were the nuclear receptors NR0B1 (DAX-1), which was found in the adrenal and testis, and NR5A1 (SF-1), which was found in all three tissues. Pathogenic variants in NR0B1 cause X-linked adrenal hypoplasia and impaired spermatogenesis in males (Suntharalingham et al., 2015), whereas pathogenic variants in NR5A1 are associated with adrenal insufficiency, testicular dysfunction and primary ovarian insufficiency (Achermann et al., 1999;Lourenço et al., 2009;McElreavey & Achermann, 2017;Suntharalingham et al., 2015). These findings, together with the identification and distribution of steroidogenic components and the detection of known gonadal regulators (e.g. EMX2, LHX9, TCF21, WT1, GATA4) and germ cell pluripotency factors (e.g POU5F1, NANOG) provided strong validation for our dataset (Svingen & Koopman, 2013).
Analysis of genes that were differentially expressed in the adrenal gland provided further support for our experimental approach. The six genes most strongly expressed all represent major core components of adrenal steroidogenesis (CYP17A1, CYP11A1, SULT2A1, STAR, MC2R, CYP11B1) (Miller & Auchus, 2011). All have associations with human adrenal disorders except for sulfotransferase 2A1, which is involved in sulphation of dehydroepiandrosterone (DHEA) to dehydroepiandrosterone sulphate (DHEAS), the main steroid synthesised by the fetal adrenal zone. Novel factors identified in the adrenal gland include adrenal-specific factors such as ASB4 and NPR3, and two kinases (MAP3K15, NRK) that are also expressed in the testis. FOXO4 and TBX3 emerged as novel transcription factors that have not previously been implicated in adrenal development. Of note, several genes associated with adrenal insufficiency in humans were not differentially expressed during early development; namely NNT, TXNRD2, AAAS, MCM4 and SGPL1. Most of these factors have ubiquitous expression and may have a pathogenic effect through oxidative stress pathways or through the accumulation of toxic metabolites. These conditions may be precipitated by postnatal stress and affected children very rarely present with adrenal insufficiency in early infancy, unlike other conditions associated with differentially expressed genes where children often present in the neonatal period ( (Hanley et al., 2000). SRY is already upregulated at the start of our study (45 dpc) and falls with the onset of steroidogenesis around 56dpc. By comparing early testis samples (CS18-CS22) with 46,XY controls, SRY was the only protein coding Y-chromosome gene detected. Although it is widely assumed that SRY is the only Y-chromosomal testis-determining factor, this provides the first direct evidence in humans that this is the case.
Other genes up-regulated in the early testis show an overlap with mesonephric development and factors involved in nephrotic syndrome (e.g., NPHS2, WT1). Of note, the highest differentially expressed gene encodes renin (REN), which supports data from previous studies of mouse testis development (Beverdam & Koopman, 2006;Bouma et al., 2007;Inoue et al., 2016;Nef et al., 2005). Renin is a hormone typically secreted from the juxtaglomerular apparatus in the kidney that regulates angiotensin and aldosterone, thereby modulating renal sodium reabsorption and arterial blood pressure. Its role in the testis is unclear, as is its expression in extra-adrenal tissue. Furthermore, it has been proposed that SRY-dependent expression of renin-angiotensin system genes, including REN itself, is a major factor responsible for differences in blood pressure between men and women (Araujo et al., 2015;Prokop et al., 2012).
It is well-established that SOX9 is a target of SRY in testis determination, and SOX9 expression rises in the early human testis following expression of SRY (Hanley et al., 2000;Ostrer et al., 2007;Sekido & Lovell-Badge, 2008). Although data from mice and humans suggest that SOX9 is in itself sufficient to trigger downstream testis pathways (Bishop et al., 2000;Kim et al., 2015) it is unclear if other SRY targets exist, or what direct targets of SOX9 might be. Attempts have been made to investigate this using ChIP-Chip in mice, but this approach is challenging given the amount of fetal tissue needed at the correct stage of development (Bhandari et al., 2012;Li et al., 2014).
An alternative approach is to study the pattern of SOX9 expression across a time-series, and to identify other genes with similar expression profiles. Rather than using a correlative approach, we developed a mathematical model that could be used to capture subtle changes in our dataset more quantitatively. Using this approach, coupled with subsequent review of gene expression in different fetal and adult tissues, we have identified a small group of genes that are up-regulated with SOX9 in the testis (e.g., CITED1, ZNF280B, PRPS2).
Of these, CITED1 (encoding Cbp/p300 interacting transactivator with Glu/Asp rich carboxy-terminal domain 1) is of greatest interest, as it is a transcriptional co-regulator, and strongly expressed in the adult testis, epididymis and pituitary gland.  Plisov, 2005;Rodriguez et al., 2004). Although Cited1 was identified as a target of Sry in one mouse ChIP-CHIP study (Li et al., 2014), little is known about Cited1 in mouse testis. In contrast, the related protein Cited2 influences adrenal and gonad development in the mouse, potentially through the regulation of Sf1 as well as of Sry itself (Buaas et al., 2009;Combes et al., 2010). CITED2 is ubiquitously expressed in humans and does not show differential expression in our datasets, and pathogenic variants in CITED2 have not been found in patients with adrenal and adrenogonadal disorders, suggesting that CITED1 might be a more important coregulator in humans (Ferraz-de-Souza et al., 2009). Amongst the other key genes identified, ZNF280B has been proposed to encode a negative regulator of p53 in prostate cancer cells (Gao et al., 2013), whereas PRPS2 encodes a testis specific phosphoribosylpyrophosphate synthetase that may function as an anti-apoptotic factor (Lei et al., 2015;Taira et al., 1990).
The advantage of having multiple serial samples for the testis is that the BALT model could be used to accurately pinpoint the onset of steroidogenesis. Using this approach a distinct sigmoid-shaped upregulation of all the known genes involved in testicular steroidogenesis was seen between 54 and 57 dpc, with marked increases in fold-change. This discrete change in activity allowed us to divide the datasets into "pre-steroidogenic" and "post-steroidogenic" stages.
In order to identify novel components involved in steroidogenesis, genes were identified that were up-regulated in the testis poststeroidogenesis compared to pre-steroidogenesis, and also upregulated in the adrenal gland compared to control samples. A total of 45 genes were identified, including all known core components of steroidogenesis and many factors involved in cholesterol biosynthesis and oxidative stress.
The 17 "novel" genes identified included MGARP (also known as OSAP), DLK1 (PREF1), and MAP3K15 (ASK3). MGARP has been implicated in mitochondrial trafficking and steroidogenesis (Jia et al., 2014;Jin et al., 2012;Zhou et al., 2011), but may also be a key target of OCT4 during reprogramming to pluripotency (Tiemann et al., 2014). DLK1 is an imprinted gene that is involved in inhibiting adipocyte differentiation, which has previously been shown to be expressed in the outer undifferentiated zone of the male rat adrenal gland and signals to the adrenal capsule to regulate zonation (Guasti et al., 2013). MAP3K15 seems to be a novel kinase linked to steroidogenesis, which has very recently been shown to be differentially expressed in the adult adrenal gland (Bergman et al., 2017). Of note, the Map3k15/Ask3 knockout mouse is hypertensive and studies have linked this kinase to the response to osmotic stress and blood pressure regulation in the mouse kidney, through interactions with WNK1 or WNK4 (Maruyama et al., 2016;Naguro et al., 2012). In our dataset, MAP3K15 was predominantly expressed in steroidogenic tissue, and not in the kidney. Similar findings are seen in RNA expression data for MAP3K15 in the Human Protein Atlas. Taken together, these findings suggest that adrenal effects of MAP3K15 may be an important area to focus future research on.
Other potential novel components of steroidogenesis with unknown function include FOXO4, GRAMD1B and RMDN2 (FAM82A1). FOXO transcription factors play important roles in development, cancer and metabolic homeostasis (Coomans de Brachène & Demoulin, 2016). FOXO4 has been implicated as a negative regulator of cell proliferation in several tumours, but a link to steroidogenesis has not previously been established (Li et al., 2016). GRAMD1B has been shown to be enriched in mouse fetal Leydig cells following CRHR1 agonist stimulation (McDowell et al., 2012), and has been identified in other expression studies of mouse fetal Leydig cells (Inoue et al., 2016;McClelland & Yao, 2017). We have shown by immunohistochemistry that there is clear upregulation of GRAMD1B in the cytoplasm of NR5A1positive steroidogenic cells with the onset of steroidogenesis in humans. Similar upregulation of expression of RMDN2, a putative microtubule regulator, was also seen in these cells. However, the exact functional role of both GRAMD1B and RMDN2 in steroidogenesis is currently unknown.
In addition to steroidogenesis, the fetal testis has important biological roles through the secretion of proteins or hormones that regulate Müllerian regression (e.g., AMH), endocrine feedback (e.g., inhibin B) and testis descent (e.g. INSL3) (Valeri et al., 2013). These circulating proteins are also proving to be important biomarkers of testis integrity and function postnatally. Using a strategy to identify genes encoding secreted proteins, all three major factors were identified, namely AMH, inhibin (α-subunit) and INSL3. Potential novel secreted proteins include EPPIN, which has been extensively studied as a novel approach to male contraception (O'Rand et al., 2016), and SCUBE1, a secreted cell surface glycoprotein belonging to the EGF superfamily. SCUBE1 was first identified from vascular endothelial cells and is released by activated platelets. Plasma levels of SCUBE1 have been measured as a potential biomarker in several conditions, such as gastric cancer, hypertension, ischaemic heart disease, and stroke (Özkan et al., 2013). Higher levels of SCUBE1 have been reported in adult males compared to females in some, but not all studies, although few control data exist (Ulusoy et al., 2012). An acute elevation in SCUBE1 has been demonstrated in experimental models of rat testis torsion, but it is unclear whether testicular SCUBE1 contributes significantly to circulating levels of SCUBE1 in humans (Turedi et al., 2015).
Although ovary development was originally viewed as a relatively passive process compared to testis development, studies in mice have shown that a similar number of differentially expressed genes are up-regulated in both organs during early development (Beverdam & Koopman, 2006;Bouma et al., 2007;Nef et al., 2005). Our data show that this is true in human gonad development too, with 274 and 280 genes differentially expressed in the ovary and testis, respectively, at log 2 FC>1, and 69 and 57 genes at log 2 FC>2. However, pathway-enrichment analysis of the tissue-specific genes revealed virtually no annotated pathways in the ovary, whereas tissue-specific pathways were more abundant in the testis. As expected, amongst the differentially-expressed ovary factors are several key genes involved in germ cell development, e.g. NANOG, POU5F1, TFAP2C, ZFP42 (REX1). In addition, several genes involved in neurotransmission were higher in the developing ovary but were attenuated in the developing testis (OR10G9, GABRG1, OR4D5, NPNT, CNTNAP4 and NPY). These include two olfactory receptors of unknown function, a GABA receptor isoform subunit, neuropeptide Y and contactin associated protein like 4 (CNTNAP4), which has recently been implicated in GABAergic synaptic transmission (Karayannis et al., 2014). Clearly, our understanding of very early human ovary development is still limited.
Although we feel this is an important study in human development, there are several limitations. Firstly, whole tissue samples were used for RNA extraction and RNA analysis. Using this approach transcriptomic signals from small populations of cells might be relatively attenuated. For example, the interstitial/Leydig cell population in the testis is only a relatively small population of the total number of cells. We have addressed this by studying relative changes in gene expression over time, which is independent of the total signal strength, and confirming key findings by immunohistochemistry, although the proportion of cell populations many also change with time. Future studies of single cell transcriptomics will be able to address this in more detail. Secondly, these studies use microarray analysis rather than RNA-Seq. Microarray analyses are restricted by the specificity and number of probes and the arrays used in this study did not address differently spliced isoforms. Future studies using RNA-Seq might address both these points.

Conclusion
We have developed an important dataset describing the genomic landscape across three different tissues and spanning critical stages in human embryonic and fetal development. We have identified several important elements previously determined from studies in mice (e.g. Renin in testis), but also several novel aspects that may be specific to humans (e.g. CITED1 rather than CITED2 in early testis). By studying time-series changes and differential expression patterns, we have discovered novel components of these important biological processes and interesting new pathways for further investigation. The dataset is highly validated based on known human disorders of adrenal and gonad function and will provide important insight to assign likely biological significance to new genes identified in whole exome and whole genome studies of patients with undiagnosed conditions. Finally, from the limited number of examples discussed above, novel insights from our genomic atlas could fuel many aspects of translational research, including biological sex-differences in humans, vulnerability to early life insults and programming, blood pressure regulation and stem cell/germ cell development.

Data availability
Array data are available from the ArrayExpress database under the accession number E-MTAB-5525.

Consent
Human embryonic and fetal tissue samples used in this study were obtained in collaboration with the MRC/Wellcome Trust-funded Human Developmental Biology Resource (HDBR). The HDBR is a tissue bank regulated by the Human Tissue Authority. Samples were collected with appropriate maternal written consent and with approval from the NRES Committee London-Fulham (REC reference 08/H0712/34+5).

Author contributions
IdV and FB were involved in bioinformatic analysis, study design and manuscript preparation. AJD and RP undertook immunohistochemistry and protein characterisation. LL undertook sample preparation, RNA extraction and RTPCR. MB and SS were involved in data modelling and bioinformatic analysis. MH facilitated microarray studies and study design. DG was involved in sample coordination, dissection and staging, and study design. JCA was involved in study design, interpreting clinical and biological aspects of the work, data analysis and manuscript preparation. All authors were involved in the revision of the draft manuscript and have agreed to the final content.

Competing interests
No competing interests were disclosed. The authors have addressed all the comments satisfactorily.

Grant information
No competing interests were disclosed.

Competing Interests:
We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard. We are happy to approve publication of the revised version of this manuscript.
No competing interests were disclosed.

Competing Interests:
We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard. Research manuscript entitled by del Valle "A genomic atlas of human adrenal and gonad development" et is an ambitious research effort and it addresses many puzzling questions concerning embryonic and al. early fetal gonad and adrenal gland development in human. The authors need to be thanked for taking up such a challenge and congratulated for being able to create a nicely flowing and very well written article from an enormous amount of data. The article is indeed a pleasure to read and the microarray data is presented in a clear and easily digestible manner. The article makes a substantial contribution to our canon of knowledge concerning early events in embryonic and fetal testicular, ovarian, and adrenal gland development. While it is able to corroborate earlier findings by others, it also provides an extensive amount of new data and sheds light on species-dependent differences in organogenesis. Despite these achievements, there are concerns that should be addressed in order to improve the quality of the data and make the presented arguments more scientifically sound.
There are two major concerns that should be addressed to a satisfactory extent. Firstly, the main concern of the manuscript is data validation. Novel findings based on microarray data should always be confirmed using another method, such as qPCR and/or immunohistochemical staining (as the authors have done for the adrenal data: Fig4 and Fig 5). At many points in the manuscript (p8, p11, p25, p27) the authors propose that co-enrichment with known markers associated with distinct biological processes or monogenic disorders validates the data. This is, however, at best merely correlative data and provides only weak validation. It would greatly improve the quality of the research if the authors would validate the data at least for the most important candidate genes using qPCR and/or immunohistochemistry. The authors have included a good selection of immunohistochemical stainings but they are partially a cause for concern in their own right (see below) plus fail to provide sufficient validation for the data. The latter statement might sound harsh but in order to validate microarray data the authors need to provide at least two stainings: one that shows low/no expression at a specific time point (say CITED1 in the early testis) and high expression in another time point (CITED1 in late testis). Including this kind of data would be sufficient to address this concern and would provide the minimum acceptable validation for the microarray findings. Strictly speaking IHC does not validate microarray data as it visualizes proteins, not RNA, but at the end it is the protein expression that we are interested in and in most cases microarray/RNAseq/qPCR are just tools to identify the differentially expressed mRNAs that are translated into proteins.
The second major concern is the quality of immunohistochemical stainings. Particularly lack of negative control stainings and difficulty to interpret the data due to inappropriate magnification and/or absence of high magnification insets. Please include negative control stainings at least for the novel markers. Provide high magnification inset for the following figures: Fig.4c top and middle panel, Fig.12c and Fig.12d. Due to poor quality, please consider providing an alternative image of higher magnification or after cropping (+inset) for Fig.12d and Fig.13d top panel.
Minor points: P8, left column, row 12 from top. For Figure 2c FC≥2 was used, whereas for 2d FC≥1. Therefore it is a bit misleading to use term "…similar analysis…" in this case.
P8, right column, row 1 from top: check spelling. Related to referee question #5: Is the data presented in Figure 14b and 14c available somewhere to have a look at the differentially expressed genes (cf. Figure 4b/ Table 1)? P26, left column, rows 5-7 from the bottom. The provided data only indicates that GRAMD1B and RMDN2 are expressed in 9 wpc testis. Because only one time point was included in the analysis, it does not demonstrate of these genes/enhanced expression of the protein products, upregulation as stated in the text.

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate? I cannot comment. A qualified statistician is required.

Are the conclusions drawn adequately supported by the results? Partly
No competing interests were disclosed.

Competing Interests:
Referee Expertise: JT: developmental endocrinology, spermatogenesis; J-AM: stem cells, spermatogenesis We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 17 Oct 2017 , UCL GOS Institute of Child Health, UK John Achermann We thank the reviewers for their constructive comments. We have addressed these as shown in bold below and modified the revised manuscript accordingly.
There are two major concerns that should be addressed to a satisfactory extent. Firstly, the main concern of the manuscript is data validation. Novel findings based on microarray data should always be confirmed using another method, such as qPCR and/or immunohistochemical staining (as the authors have done for the adrenal data: Fig4 and Fig 5). At many points in the manuscript (p8, p11, p25, p27) the authors propose that co-enrichment with known markers associated with distinct biological processes or monogenic disorders validates the data. This is, however, at best merely correlative data and provides only weak validation. It would greatly improve the quality of distinct biological processes or monogenic disorders validates the data. This is, however, at best merely correlative data and provides only weak validation. It would greatly improve the quality of the research if the authors would validate the data at least for the most important candidate genes using qPCRand/or immunohistochemistry. The authors have included a good selection of immunohistochemical stainings but they are partially a cause for concern in their own right (see below) plus fail to provide sufficient validation for the data. The latter statement might sound harsh but in order to validate microarray data the authors need to provide at least two stainings: one that shows low/no expression at a specific time point (say CITED1 in the early testis) and high expression in another time point (CITED1 in late testis). Including this kind of data would be sufficient to address this concern and would provide the minimum acceptable validation for the microarray findings. Strictly speaking IHC does not validate microarray data as it visualizes proteins, not RNA, but at the end it is the protein expression that we are interested in and in most cases microarray/RNAseq/qPCR are just tools to identify the differentially expressed mRNAs that are translated into proteins.
We agreed that further validation would be ideal at both the RNA level (qRTPCR) and protein level (immunohistochemistry), but we tried to focus on the key new genes described and were limited by the specificity of antibodies for IHC for some genes/protein of interest. We did focus largely on IHC data as we felt this provided validation at the protein level as well as RNA, but agree we could have included more validation across the time series (some of which we omitted because of space).
To address these issues, we have: Antonarakis SE, Dermitzakis ET, Nef S. http://www.biorxiv.org/content/early/2017/09/18/190264). We feel that these are additional novel, independent data to show the potential link between SOX9 and CITED1 in testis development.

4)
We now provide sequential IHC data for GRAMD1B and RMDN2 before and after the onset of steroidogenesis to show an increase in protein expression in the cytoplasm of Leydig cells (Supplementary Figure 2b).
Please note that in the original manuscript we have included links to the Human Protein Atlas for many key novel factors. This resource provides IHC data for adult tissues as well as graphical data for qPCR datasets (e.g. FANTOM5, GTEx).
The second major concern is the quality of immunohistochemical stainings. Particularly lack of negative control stainings and difficulty to interpret the data due to inappropriate magnification and/or absence of high magnification insets. Please include negative control stainings at least for the novel markers. Provide high magnification inset for the following figures: Fig.4c top and middle panel, Fig.12c and Fig.12d. Due to poor quality, please consider providing an alternative image of higher magnification or after cropping (+inset) for Fig.12d and Fig.13d top panel.
Negative control data are particularly important to assess binding of secondary antibodies and non-specific effects. We include representative negative data in Supplementary Figure 3b Minor points: P8, left column, row 12 from top. For Figure 2c FC≥2 was used, whereas for 2d FC≥1. Therefore it is a bit misleading to use term "…similar analysis…" in this case.
We have changed this to "analysis".
We are unclear what this refers to. We agree that these data should have been presented better. We have repeated the experiment again three times in triplicate and present data with error bars. We have replaced this part of the Figure 5b.
Related to referee question #5: Is the data presented in Figure 14b and 14c available somewhere to have a look at the differentially expressed genes (cf. Figure 4b/ Table 1)?
We have uploaded this dataset into OSF as Dataset 7. We apologise that this was not included.
P26, left column, rows 5-7 from the bottom. The provided data only indicates that GRAMD1B and RMDN2 are expressed in 9 wpc testis. Because only one time point was included in the analysis, it does not demonstrate upregulation of these genes/enhanced expression of the protein products, as stated in the text.
We have now included qRTPCR and IHC to show up-regulation over time in Supplementary Figure 2a,  The comprehension of the mechanisms underlying adrenal and gonadal development in the first trimester of intrauterine life is essential for the understanding of human congenital disorders. Most of the information available on normal embryonic and fetal development of the adrenals and gonads derives from studies in experimental models. This study provides crucial insight into the comprehension of these developmental processes in humans. In fact, by performing microarray studies on tissues obtained from accurately staged embryos and fetuses to assess time-series changes and differential expression patterns, del Valle and colleagues have constructed an atlas of RNA expression in the adrenals, testes and ovaries from 6 to 10 weeks post-conception, a critical period in the differentiation of these organs. The authors validate their findings using quantitative RT-PCR and immunohistochemistry of genes known to be important in adrenal and gonadal development. The expression profiles of known genes were corroborated whereas unknown profiles of genes in the gonads and adrenals have been described. These findings set the bases for a wide variety of future studies addressing the physiology and pathophysiology of human adrenal and gonadal development.
We have minor comments: We have minor comments: Title and text: We have concerns about the use of the phrase "genomic atlas" in this study. In fact, genomic atlases, e.g. cancer genome atlas, refer to the comprehensive description of variants of the genome, whereas this work makes a comprehensive description of the transcripts present in the three organs at different developmental stages. The phrases "transcriptomic atlas" or "atlas of gene expression" (as used on page 5) seem more adequate.
Methods and Results: Carnegie and fetal stages are based on the morphological development of the embryo, and are not directly dependent on either age or size. This concept should be more clearly expressed (e.g. Supplementary Table 1 and text in Methods) in order to avoid that the non-specialised reader interpret that there are precise cut-offs (e.g. CS18=44-47.5 days, etc.).
Methods and Results: Tissue transcriptome analysis by microarray is a powerful technique. Nonetheless, like other techniques based on RNA extraction from whole tissue, one limitation is the confounding effect due to modifications in RNA levels in single components whose relative proportion in the tissue varies or differs. For instance, STAR was 73-fold higher in the adrenal tissue probably because most of the tissue sample was composed of steroidogenic cells, whereas a similar amount of testis sample would only have a few Leydig cells. This does not mean that one cortical adrenal cell expresses more STAR than one Leydig cell. Or, if Sertoli cells represent 80% of testicular tissue on week 6 and 50% on week 10, microarray results from whole testicular tissue will show a lower expression level of a given Sertoli cell RNA even if each Sertoli cell expresses the same amount of the given RNA. This limitation should be considered in the interpretation of the results: up/down-regulation give the idea of regulated increase/decrease of a given RNA expression in one cell rather than increase/decrease of RNA levels in whole tissue, which could be due to the relative proportion of cells expressing that given RNA. Therefore, the words up-regulation and down-regulation may be misleading, and may need to be replaced (e.g. "more/less abundant in whole testicular/ovarian/adrenal tissue").
Methods: We do not have expertise to assess the BALT model used to pinpoint the time when gene expression levels change. A mathematician/statistician should be consulted.
Results: We would have expected SRY and SOX9 to appear as a testis differentially expressed genes on Figure 2. Could the authors discuss why they did not select them for this figure?
Figure 6: A differential expression of SRY on the adrenal is unexpected. Could the authors comment on this?
Discussion: page 25, "The role of SRY as the primary testis-determining gene has been known for 25 years, following transgenic mouse studies as well as the discovery of duplication of SRY in individuals with 46,XX testicular disorders of sex development (DSD)". Do the authors mean translocation of SRY, or is it really duplication?

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate? I cannot comment. A qualified statistician is required.
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Yes
No competing interests were disclosed.

Competing Interests:
We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 17 Oct 2017 , UCL GOS Institute of Child Health, UK John Achermann We thank the reviewers for their positive and constructive comments. We have addressed these as shown in bold below and modified the revised manuscript accordingly.

Minor comments:
Title and text: We have concerns about the use of the phrase "genomic atlas" in this study. In fact, genomic atlases, e.g. cancer genome atlas, refer to the comprehensive description of variants of the genome, whereas this work makes a comprehensive description of the transcripts present in the three organs at different developmental stages. The phrases "transcriptomic atlas" or "atlas of gene expression" (as used on page 5) seem more adequate.
We agree that in the age of genome sequencing there might be confusion and that other terms might be more accurate, but this term has been used to describe developmental datasets before so we would prefer to leave the title as it is.
Methods and Results: Carnegie and fetal stages are based on the morphological development of the embryo, and are not directly dependent on either age or size. This concept should be more clearly expressed (e.g. Supplementary Table 1 and text in Methods) in order to avoid that the non-specialised reader interpret that there are precise cut-offs (e.g. CS18=44-47.5 days, etc.).
We have added a note in the methods to state that staging is based on morphological characteristics.

characteristics.
Methods and Results: Tissue transcriptome analysis by microarray is a powerful technique. Nonetheless, like other techniques based on RNA extraction from whole tissue, one limitation is the confounding effect due to modifications in RNA levels in single components whose relative proportion in the tissue varies or differs. For instance, STAR was 73-fold higher in the adrenal tissue probably because most of the tissue sample was composed of steroidogenic cells, whereas a similar amount of testis sample would only have a few Leydig cells. This does not mean that one cortical adrenal cell expresses more STAR than one Leydig cell. Or, if Sertoli cells represent 80% of testicular tissue on week 6 and 50% on week 10, microarray results from whole testicular tissue will show a lower expression level of a given Sertoli cell RNA even if each Sertoli cell expresses the same amount of the given RNA. This limitation should be considered in the interpretation of the results: up/down-regulation give the idea of regulated increase/decrease of a given RNA expression in one cell rather than increase/decrease of RNA levels in whole tissue, which could be due to the relative proportion of cells expressing that given RNA. Therefore, the words up-regulation and down-regulation may be misleading, and may need to be replaced (e.g. "more/less abundant in whole testicular/ovarian/adrenal tissue").
We have approached this study on the "whole tissue" level and talk about up-regulation or down-regulation of genes in this regard. However, it is an important point to consider, so we have added a paragraph in the discussion to address some of the limitations of this work and included this point.
Results: We would have expected SRY and SOX9 to appear as a testis differentially expressed genes on Figure 2. Could the authors discuss why they did not select them for this figure?
We used a cut-off value of log2 FC >1 (i.e., two-fold) for this analysis and both factors were just below. As can be seen in Fig. 6b and 8a even changes within the tissue are subtle. SRY may have some extra-gonadal expression, for example in the adrenal and brain. We mentioned the potential role of SRY in the regulation of renin in the discussion. This is an interesting area for future research.
The inv(9)(p11q13) balanced translocation is quite a common finding in the general population (1.2-2.3%). Although data originally suggested this could be linked to infertility, more recent publications have shown that differences in 1800 infertile people were not statistically significantly higher ( . Several studies are also biased as this change can increase with age and individuals being investigated for fertility were often older. There were no specific outlying effects of these samples in our analysis, so we felt comfortable including them in the dataset. including them in the dataset. Discussion: page 25, "The role of SRY as the primary testis-determining gene has been known for 25 years, following transgenic mouse studies as well as the discovery of duplication of SRY in individuals with 46,XX testicular disorders of sex development (DSD)". Do the authors mean translocation of SRY, or is it really duplication?
. The development of the adrenal glands and gonads represent two significant embryological events. The development of both tissues is closely linked as they both derive from the same region of intermediate mesoderm in early embryogenesis. Given this shared history it is not surprising that they share many processes such as steroidogenesis in common. However, until the current study most of our understanding of spatial and temporal gene expression profiles in these tissues has come from mouse studies which are not always applicable to humans. This study by Achermann and colleagues finally addresses this fundamental gap in our knowledge. The authors extracted RNA from adrenal, testis, ovary and control tissues between 6-10 weeks (wpc) in human embryos. The RNA was then post coitum subjected to an Affymetrix array and subsequent differential gene expression analysis and mathematical modelling to identify temporal changes across the dataset. In addition, they used pathway analysis and confirmed the cellular localization of novel genes using immunohistochemistry.
The study confirmed the expression profiles of many known genes such as and identified many NR5A1 potential novel genes. Interestingly, they showed that the only protein coding Y-chromosome gene detected in the developing human testis (CS18-CS22) was the gene. Twenty-five years since the SRY discovery of as the sole primary testis determining factor we now have further proof in human tissue SRY that this is the case. Novel components of the adrenal pathway have been identified as have genes that are potentially up-regulated by such as , and . In addition, potential novel SOX9 CITED1 ZNF280B PRPS2 genes were associated with adrenal and testicular steroidogenesis (eg. , and ).

FOXO4 MAP3K15
RMND2 Of interest was the finding that while a similar number of differentially expressed genes were found in the testis and ovary, pathway-enrichment analysis failed to reveal any annotated pathways in the later. The study did reveal some interesting candidate ovarian genes such as and , some of which OR10G9 OR4D5 appear to be involved in neurotransmission. This clearly highlights our lack of understanding of early ovarian development, and suggests that additional methods may be required to find novel genetic pathways important in this tissue.
The public availability of this data will be a great future resource for any researchers interested in human adrenal or gonadal development. In addition, this database may allow researchers to isolate novel genes involved in gonadal and adrenal disorders, many of which currently lack a defined genetic aetiology. Finally, while the authors used a highly comprehensive set of tissues, expression analysis using microarray platforms can have several drawbacks compared to technologies such as RNA-Seq. microarray platforms can have several drawbacks compared to technologies such as RNA-Seq. Microarrays have been shown to have poorer sensitivity for lowly-abundant RNA compared to RNA-Seq and may also miss expression data for poorly annotated genes and alternative transcripts. We therefore look forward to similar studies on human fetal adrenal and gonads using RNA-Seq and ultimately single cell fate mapping approaches as have been achieved in mouse.

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and is the work technically sound? Yes

If applicable, is the statistical analysis and its interpretation appropriate? Yes
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Yes
No competing interests were disclosed.

Competing Interests:
We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 17 Oct 2017 , UCL GOS Institute of Child Health, UK John Achermann We thank the reviewers for their positive comments and outlining some limitations. We have addressed these as shown in bold below and modified the revised manuscript accordingly.
Finally, while the authors used a highly comprehensive set of tissues, expression analysis using microarray platforms can have several drawbacks compared to technologies such as RNA-Seq. Microarrays have been shown to have poorer sensitivity for lowly-abundant RNA compared to RNA-Seq and may also miss expression data for poorly annotated genes and alternative transcripts. We therefore look forward to similar studies on human fetal adrenal and gonads using RNA-Seq and ultimately single cell fate mapping approaches as have been achieved in mouse.
These comments provide a very useful summary of potential drawbacks. We have included a paragraph in the discussion outlining some limitations of the study and future directions.
No competing interests were disclosed. Competing Interests:

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 17 Oct 2017 , UCL GOS Institute of Child Health, UK John Achermann We thank the reviewers for their constructive comments. We have addressed these as shown in bold below and modified the revised manuscript accordingly.
This reviewer has only a few comments to be considered.
The mathematical model and part of the statistics used are complex and not easy to follow for a non-expert. It should therefore be reviewed by an expert statistician.
Most of the statistics used for comparing expression between samples and heatmaps use a standard R package. The model developed is also a fairly standard approach and similar profiles were obtained using soft and hard clustering analysis (where genes can be clustered to multiple or single clusters), but the model provided more detailed analysis over time.
It is not entirely clear if data from XX adrenals were compared with data from XY adrenals and if there was any difference.
The adrenal samples were intentionally chosen to be from 46,XY embryos/fetuses so that they could be compared directly with testis without the complicating effects of any sex-chromosome gene dosage. We agree that sex differences in extra gonadal tissues are an interesting area to study but are beyond the scope of this manuscript.
In Fig. 13 D, lower part, SCUBE1 expression seems to be localized to flat cells in the periphery of the chords. Have the authors considered that these cells may be peritubular cells (PTC) rather than Sertoli cells or gonocytes? PTC may be identified by their expression of aSma and PDGF-Ra.
Thank you. At this stage of development the cells at the boundary of the cords are not particularly flattened like classic peritubular myoid cells. Certainly the cells containing SCUBE1 are fairly round and in some cases staining is seen throughout the cytoplasm and in some more centrally positioned cells. As SCUBE1 is potentially a secreted protein we wondered if the peripheral localisation was to do with this. However, we take the reviewers point that these could be developing PTM cells. We have not been able to address this further at present but will do so if we explore the role of SCUBE1 in the future.
No competing interests were disclosed. Competing Interests: