The animal model used in this study had been successfully constructed in our previous study, the flowchart of the experimental procedure is displayed in Figure 1 (Zhang et al., 2019). Briefly, all animal procedures were performed in accordance with the guidelines of the Animal Ethics Committee of Zhejiang University. Timed mating was performed in virgin Sprague–Dawley females weighing 250 ± 10 g (mean ± SD) obtained from Slack Experimental Animal Center of Chinese Academy of Sciences (Shanghai, China). After confirmation of pregnancy by vaginal smear plug, the females were randomly divided into two groups with an ad libitum diet of standard laboratory chow or a 50% food-restricted diet determined by the quantification of normal intake in ad libitum-fed rats from day 1 of pregnancy until parturition. The offspring delivered by rats with normal food intake were defined as the control group, while the offspring delivered by food-restricted rats with a birth weight below the 10th percentile of controls were defined as IUGR. To avoid hormonal cycle disturbances, only male pups were further investigated. All pups were then weaned onto standard chow at 3 weeks. At 9 weeks of age, the IUGR and control offspring were randomly divided into two groups. The normoxia group was placed in a normal room air chamber, while the hypoxia group was placed in the hypoxia chamber with the fraction of inspired oxygen of 11% for 2 weeks. These subdivisions yielded two subgroups per group: control-normoxia (n = 8), control-hypoxia (n = 8), IUGR-normoxia (n = 8), and IUGR-hypoxia (n = 8). Rats in all the groups were fed with standard chow ad libitum and were housed in separate cages maintained with constant temperature and humidity. The cages and soda lime in the chambers were changed daily to prevent the accumulation of excess ammonia and CO₂, respectively. Samples from the liver were harvested and snap-frozen in liquid nitrogen, and five samples were selected randomly each from control-normoxia, control-hypoxia, IUGR-normoxia, and IUGR-hypoxia for subsequent processing.

Total RNA was isolated from the liver using TRIzol reagent (Invitrogen Life Technologies), followed by the determination of the quality and integrity using a NanoDrop spectrophotometer. First-strand cDNA was synthesized using random hexamers and SuperScript II, while second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The second-strand cDNA was then purified and treated using terminal repair and ligation primers. Products were purified (AMPure XP beads) and quantified using the Agilent high-sensitivity DNA assay on a Bioanalyzer 2100 system (Agilent). After selection using AMPure XP beads, the second-strand cDNA containing U was degraded by the USER enzyme. PCR amplification was then performed, followed by sequencing of the library using Illumina HiSeq sequencing by Whbioacme Co. Ltd.

First, we excluded the samples that emerged as clear outliers using the hierarchical clustering method. Second, we applied WGCNA with the expression profile using the R package of “WGCNA” (Langfelder and Horvath, 2008). Pearson’s correlations were performed for all pair-wise genes, followed by the construction of a correlation matrix. The correlation matrix was then transformed into an adjacency matrix (also known as the scale-free network) using an appropriated soft thresholding powers β. The β value was set from 1 to 20, followed by the calculation of the scale-free fit index and mean connectivity for each β value. The β value was determined when the scale-free fit index was up to 0.85, and the β value that had the highest mean connectivity was considered to be the most appropriated one. The topological overlap matrix (TOM) was then obtained on the basis of the adjacency matrix. Finally, the average linkage hierarchical clustering was conducted according to the TOM-based dissimilarity measure in order to divide all genes into several co-expression modules having gene module sizes above 30. We then calculated Pearson’s correlations of the eigengenes after defining the first principal component of a given module to further analyze the module. The modules whose eigengenes were highly correlated were merged into one module, with Pearson’s correlation being higher than 0.75.

The clinical traits of our samples included IUGR and hypoxia. The correlation between modules and traits was calculated, and the modules having a positive correlation with IUGR and hypoxia were considered as playing roles in the pathogenesis of the disease. On the other hand, genes in modules having a positive correlation with a normal trait were indispensable for maintaining normal biological functions. Thus, we extracted gene modules having the highest correlation with IUGR and hypoxia for subsequent studies.

Gene significance (GS) and module membership (MM) were introduced to the correlation of a given gene with a clinical trait and module eigengene, respectively. Genes in clinical-related modules should have high values and preferable correlations of GS and MM.

We used clusterProfiler R package (version 3.15.3) to perform Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis in order to explore the biological functions of genes in the clinical related modules (Shannon et al., 2003). The p-value < 0.01 and Benjamini–Hochberg adjusted p-value (Benjamini and Hochberg 1995) <0.01 were considered to be statistically significant.

Generally, the hub genes in the co-expression network should have high connectivity with the whole module and clinical trait. For identifying the hub gene related with IUGR-hypoxia, we extracted gene clusters that were enriched in certain GO terms from the WGCNA network to construct subnetworks after GO enrichment analysis of the IUGR-hypoxia-related module. The CytoHubba Cytoscape plugin was then used to identify hub genes (Chin et al., 2014). This study regarded the genes with high Maximal Clique Centrality (MCC) values as being hub genes because MCC has a better performance in predicting hub genes.

We validated the expression of hub genes using qRT-PCR from liver tissue in order to obtain further evidence for the significance of hub genes in the module of interest.

The dissected liver was homogenized, followed by the isolation of RNA using TRIzol reagent (Invitrogen) in accordance with the manufacturer’s instructions. The obtained total RNA was then converted to cDNA using a reverse transcription reagent kit (TaKaRa). Real-time quantitative PCR was performed using a 7500 Real-Time PCR Detection System (Applied Biosystems) using a PCR kit containing SYBR Green (TaKaRa). The obtained results were expressed in relative expression using the comparative 2^{−∆∆ Ct} method normalized by the housekeeping gene GAPDH. Statistical significance was calculated using Student’s t-test by SPSS 18.0 (SPSS, Inc., Chicago, IL, United States). p < 0.05 was considered to indicate a statistical difference (Table 1).

The mean pulmonary arterial pressure, right ventricular hypertrophy, the media wall area thickness, and pulmonary vascular endothelial function assessment had been shown in our previous study (Zhang et al., 2019), which provided evidence that our model was built successfully.

Since non-varying genes are usually regarded as background noise, in this study, we filtered the genes by variance, and the top 50% (5,000) genes with larger variance were chosen for subsequent analyses. Sample clustering excluded the outlier samples, thereby leaving a total of five control-normoxia samples, five IUGR-normoxia samples, four control-hypoxia, and five IUGR-hypoxia samples. The final results after sample clustering revealed intra-group consistency and intergroup variability, as shown in Figure 2A.

In this study, the soft thresholding was determined as eight on the base of the scale-free fit index and mean connectivity (Figures 2B,C). Accordingly, the correlation matrix was transformed to the adjacency matrix and then converted to a topological overlap matrix to improve the sensitivity of co-expression identification. A total of eight co-expression modules with different genes were generated and displayed with different colors based on the average linkage hierarchical clustering and module merging (Figure 3A). The modules included 85 genes in a midnight blue module, 1,641 genes in a blue module, 447 genes in a tan module, 219 genes in a green–yellow module, 179 genes in a purple module, 591 genes in a magenta module, and 1,761 genes in a turquoise module. The gray module contained genes that could not be divided into any co-expression modules.

Figure 3B shows the topological overlap matrix among all genes analyzed by WGCNA; light color represents low overlap, while the progressively darker red color represents a higher overlap.

The associations between modules and trait were investigated using Pearson’s correlation analysis in the WGCNA, and the obtained results are shown in Figure 4A. The results indicated that the blue module displayed highest correlation with IUGR-hypoxia (r = 0.84, p = 4e-06). The GS vs. MM values of all member genes of the blue module are shown in the scatterplot (Figure 4B). The GS and MM values were highly correlated in the blue module (cor = 0.81, p < 1e-200), indicating that genes in the blue module were highly significantly associated with clinical traits, and thus, they were suitable for further analyses.

The GO and KEGG enrichment analyses were performed to confirm the biological themes of genes in the blue module and obtain further insights into the underlying biological pathways behind IUGR-hypoxia. The top 15 significant terms of GO and KEGG are exhibited in Figures 5A,B, respectively.

Enriched GO-biological process (BP) terms for the top IUGR-hypoxia module were mainly about metabolic process (Figure 5A), such as “fatty acid metabolic process” (gene count = 87, p = 1.24E-20), “lipid modification” (gene count = 52, p = 1.44E-12), and “fatty acid catabolic process” (gene count = 32, p = 6.57E-13).

The KEGG enrichment results, mainly on lipid metabolism, were similar to those obtained from GO-BP (Figure 5B). The top KEGG terms included “fatty acid degradation” (gene count = 23, p = 7.20E-13), “fatty acid metabolism” (gene count = 22, p = 6.61E-09), “PPAR signaling pathway” (gene count = 22, p = 2.46E-06), and “biosynthesis of unsaturated fatty acids” (gene count = 13, p = 2.62E-06).

We extracted the genes enriched in two of the top 15 significant GO-BP terms in the top IUGR-hypoxia module, namely, “fatty acid metabolic process” and “lipid modification” and used them to construct two subnetworks of the weighted co-expression network, respectively. The obtained MCC value of genes predicted by the CytoHubba plugin was used to search for hub genes, and the most central genes of the subnetwork were screened out, as shown in Figure 6. The top 10 genes obtained using the MCC values in the pathway of the fatty acid metabolic process were Hadha, Acsm3, Cyp2j4, Mapk14, Cyp2j10, Acads, Lipc, Hacl1, Cyp2f4, and Acat1, while the top 10 genes in the pathway of lipid modification were Stard4, Acadvl, Mapk14, Ecil, Acadl, Etfdh, Fam126b, Acsl5, Agtrla, and Hacl1.

The expression of hub genes between IUGR-hypoxia and IUGR-normoxia are summarized in Figure 7A. The downregulated hub genes in IUGR-hypoxia included Cyp2f4, Agtr1a, Eci1, Acsl5, Acat1, Acadl, Lipc, Hadha, Acads, and Acadvl, while the upregulated hub genes in IUGR-hypoxia included Acsm3, Cyp2j4, Mapk14, Cyp2j10, Stard4, Fam126b, and Hacl1. The qRT-PCR analysis was conducted in IUGR-normoxia and IUGR-hypoxia liver tissues in order to determine the aberrant expression of IUGR-hypoxia-related hub genes. The obtained results indicated that the mRNA expression of Cyp2f4, Lipc, and Acadl was significantly decreased in IUGR-hypoxia liver tissues when compared with IUGR-normoxia liver tissues, whereas Hacl1 mRNA expression was increased in IUGR-hypoxia liver tissues (Figures 7B,C). These results were in concordance with the WGCNA data.