Tissue collection
Tissues were collected from fetuses at 80 (80 d, n = 4), 96–100 (100 d, n = 4), 120 (120 d, n = 4), 130 (130 d, n = 4), and 142–144 (145 d, n = 4) days of gestation and on the first (1 d, n = 4) day after delivery. Each group included one set of twin fetuses. None of the ewes showed any signs of impending labor. For collection of fetal tissues, ewes were killed with 20 ml of Euthasol solution (7.8 g pentobarbital and 1 g phenytoin sodium; Virbac AH, Fort Worth, TX) administered intravenously, the fetus was quickly removed, and the fetal brain was removed. Brain dissection was performed immediately, using sterile instruments and gloves. First, the head was separated from the backbone cutting cranial to the axis of the vertebrae. Second, the posterior cranial bones were separated using forceps. With a scalpel we cut the pituitary stalk and optic nerves to release the brain, which was then deposited onto a sterile container or surface. Holding the brain upside down, the hypothalamus was dissected, followed by both hippocampi and the brainstem. Finally, a section of the frontal cortex was dissected. Fetal tissues were rapidly frozen in liquid nitrogen and stored at -80°C.
The use of animals in this project was approved by the University of Florida Institutional Animal Care and Use Committee.
RNA extraction and preparation
RNA was extracted from the cortex, hippocampus, hypothalamus and brainstem using Trizol (Invitrogen, Carlsbad, CA) and following the manufacturer’s directions. The RNA was resuspended in RNAsecure, and stored at -80°C in aliquots until use. For microarray analysis, 20 ug of RNA was DNase treated using the Turbo RNase-free DNase kit (Ambion, Foster City, CA), the concentration determined with a Nanodrop spectrophotometer (ND-1000, ThermoFisher, Wilmington DE) and the integrity of the RNA was measured using an Agilent Bioanalyzer, 2100 model. One μg of the DNase-treated RNA was labeled with Cyanine 3 (Cy3) CTP with the Agilent Quick Amp kit (5190–0442, New Castle, DE) according to their methodology, purified with the Qiagen RNeasy kit (Valencia, CA) as according to Agilent’s revision of the Qiagen protocol as shown in the Quick Amp kit protocol except that the microcentrifugation spins were performed at room temperature instead of 4°C. The resulting labeled cRNA was analyzed with the Nanodrop spectrophotometer, and the specific activities and the yields of the cRNAs were calculated. The labeled cRNA was stored at -80°C until use.
Microarray hybridization
This was performed following protocols from Agilent. Briefly, 600 ng of each labeled cRNA was fragmented and then mixed with hybridization buffer using the Agilent gene expression hybridization kit. These were applied to sheep 8 × 15 K array slides (Agilent 019921), containing 8 arrays with 15,208 oligomers with a length of 60 bases and hybridized at 65°C for 17 h at 10 rpm. A total of 3 slides (24 arrays) were employed per each brain region (six gestational ages times four replicates per gestational age). The arrays were washed, dried, stabilized, and scanned with an Agilent G2505B 2 dye scanner at the Interdisciplinary Center for Biotechnology Research at the University of Florida. Features were extracted with Agilent Feature extraction 9.1 software.
Microarray data
The limma package was employed to import the raw data into R (http://www.r-project.org), perform background correction and normalize the data using the quantile normalization method [31]. Control probes and low expressed probes were filtered out, retaining for further analysis the probes that were at least 10% brighter than the negative controls on at least four arrays.
Statistical analysis
The Bayesian Estimation of Temporal Regulation (BETR) algorithm was used to identify the differentially expressed genes (DEG) for each brain region at a False Discovery Rate (FDR) < 0.05. The first gestational age (80 days) was considered as baseline measurement and compared to the subsequent time points, to correlate the differential expression between time points. The algorithm was applied using the BETR package for R software. A detailed explanation of the mathematical model can be found in Aryee et al. [32]. This method returns the probabilities of differential expression for each gene in the data set. Genes with a probability higher than 99.99% were considered as DEG.
Supervised weighted gene co-expression network analysis (WGCNA)
The DEG for each brain region was subjected to signed-WGCNA. Rows of each data set were collapsed forming an average expression of the genes with the same official symbol, in order to obtain unique identifiers for each gene in the working data set. The automatic method was employed for block-wise network construction and module detection. The co-expression similarity was raised to a soft thresholding power (β) of 12 to calculate adjacency. The adjacency for the signed network is defined as a
ij
= |(1 + cor(x
i
,x
j
))/2|β[33]. The resulting modules for each network were related with the gestational age to identify modules, or clusters, of co-expressed genes with increasing or decreasing expression pattern. Gene significance (GS) was defined as the correlation of i-th gene with a temporal pattern. Module membership (MM) was defined as the correlation of the i-th gene respects its corresponding module (the higher is the MM the more connected is the i-th gene with the other genes of the corresponding modules). The correlation coefficient of MM and GS was measure for each module, plotting MM versus GS. Higher correlation between MM and GS means that genes that are highly associated with the temporal pattern are also the central elements of the given module [33] (Additional file 3: Table S1). Modules that showed the highest correlation between MM and GS, either with increasing pattern (positive GS) or decreasing pattern (negative GS), were selected for Enrichment analysis. The top positively correlated module for each brain region was: the turquoise module for cortex and brainstem, the brown module for hippocampus and the blue module for hypothalamus. The top negatively correlated module for each brain region was: the blue module for cortex and brainstem, the yellow module for hippocampus and the turquoise module for hypothalamus (Additional file 4: Figure S3).
All the analyses were performed with the WGCNA package for R software. More details about the methodology for WGNA for network construction can be found at: http://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/.
Enrichment analysis
The bioinformatics tool applied for this purpose was WebGestalt (WEB-based GEne SeT AnaLysis Toolkit). This program is designed for functional genomic studies that generate large number of gene list. Then, it organizes the large set of genes based on common functional features, like GO categories or biochemical pathways [8]. The list of human official symbols for the genes composing the top modules for each network was submitted for enrichment biological processes and KEGG pathway analysis, selecting H. Sapiens as the organism of interest. The statistical method employed was the hypergeometric test, adjusting the p-values for Benjamini & Hochberg method, selecting the minimun of 5 genes per category and using the human genome as reference set. Significantly enriched Kegg pathways (p < 0.01) were compared to identify common enriched pathways between all brain regions.
Network analysis
A gene network was constructed for each brain region using CytoScape version 3.0 [13] through the GeneMania plugin [14], which was used to infer network data. The set of functional association data between genes was downloaded from the Homo sapiens database. The list of human official symbols for the genes of interest was input into the GeneMania plugin to retrieve the corresponding association network. The association data employed was protein-protein and protein-DNA interaction. Then, the resulting networks were merged together using the Merge Networks tool of Cytoscape 3.0.
Quantitative Real-Time (qRT)-PCR Validation
The mRNA samples extracted from the four brain regions at the six different gestational ages (4 fetuses/gestational age; 96 samples in total) were converted to cDNA with a High Capacity cDNA Archive kit using the methodology recommended by the kit manufacturer (Applied Biosystems, Foster City, Calif., USA). The newly synthesized cDNA was stored at -20°C until qRT-PCR was performed.
A total of 18 genes were selected for further validation by qRT-PCR. The cDNA from corresponding brain region was used to measure the mRNA of the following genes: Brainstem: CD34, CD109, CD44, CD5 and CD9; Hippocampus: CD3 gamma (CD3G), CD3 delta (C3D) and CD3 epsilon (CD3E); Cortex: CSF-1 CSF1R, IL34, CD11B, CD81, FCGRIIB, IL10, TGFB, CD24 and MBP.
Relative expression of selected genes were determined using primers (Sigma-Aldrich, St Louis, MO) and Sybr Green PCR Master Mix (Applied Biosystems, Foster City, CA). Primers were designed with Primer Express software (Applied Biosystems. Primers for CD109, CD44, CSF1, CSF1R, IL34, CD24 and MBP were designed from the corresponding bovine mRNA. Primers for CD34, CD5, CD9, CD3G, CD3D, CD3E, CD81, CD11B, FCGR2B, IL10 and TGFB were designed from the ovine mRNA. Primers sequences and accession numbers are reported in Additional file 11: Table S3. All primer pairs had efficiencies greater than 95%. The abundance of B-actin mRNA was determined in each brain region sample, using primers and VIC Taqman probes designed from the ovine B-actin sequence and Taqman qRT-PCR master mix (Applied Biosystems, Foster City, CA). All samples were run in triplicate for each gene and for B-actin. Relative mRNA expression of each gene was calculated by determining change in threshold cycle (ΔCt) between the mean Ct for each gene and the mean Ct for B-actin mRNA from the same sample. The effect of gestational age on each gene was analyzed by ANOVA using the ΔCt values. Pairwise comparisons among gestational age means were done by Duncan’s procedure. Data were graphed as the mean fold change in mRNA relative to the 80d group; fold change in each sample was calculated as 2-ΔΔCt, where ΔΔCt is the difference between ΔCt in each sample and the mean ΔCt in the control group. For all statistical analyses, the criterion for achieving statistical significance was P < 0.05.