Effect of sex chromosomes versus hormones in neonatal lung injury

The main mechanisms underlying sexually dimorphic outcomes in neonatal lung injury are unknown. We tested the hypothesis that hormone- or sex chromosome–mediated mechanisms interact with hyperoxia exposure to impact injury and repair in the neonatal lung. To distinguish sex differences caused by gonadal hormones versus sex chromosome complement (XX versus XY), we used the Four Core Genotypes (FCG) mice and exposed them to hyperoxia (95% FiO2, P1–P4: saccular stage) or room air. This model generates XX and XY mice that each have either testes (with Sry, XXM, or XYM) or ovaries (without Sry, XXF, or XYF). Lung alveolarization and vascular development were more severely impacted in XYM and XYF compared with XXF and XXM mice. Cell cycle–related pathways were enriched in the gonadal or chromosomal females, while muscle-related pathways were enriched in the gonadal males, and immune-response–related pathways were enriched in chromosomal males. Female gene signatures showed a negative correlation with human patients who developed bronchopulmonary dysplasia (BPD) or needed oxygen therapy at 28 days. These results demonstrate that chromosomal sex — and not gonadal sex — impacted the response to neonatal hyperoxia exposure. The female sex chromosomal complement was protective and could mediate sex-specific differences in the neonatal lung injury.

tiled together to recreate a seamless image of the entire tissue section using Aivia 9.8 (Leica Microsystems). Using Fiji (v. 1.53h), individual regions of interest of 675 µm x 675 µm were drawn to cover as much as possible of the lung parenchyma across all 5 lobes, while avoiding major vessels area and airways. The individual regions of interest were subsequently batch processed using MLI plugin for ImageJ/Fiji (4) to measure the mean linear intercept distance along each line, in each ROI, and each tissue section. The spacing between each line (both horizontal and vertical) was set to 15 pixels (or 11.25 µm) after determination that this was roughly equivalent to the radius of a single airways, hence guaranteeing an ideal sampling rate for the MLI measurements. Raw results were then analyzed using a custom-made R script (R Studio Desktop, v. 1.4.1106) and finally de-anonymized in order to assemble animal repeats and conditions prior to final plotting using GraphPad Prism (v. 8
Immunohistochemistry: Pulmonary vessel density was determined based on vWF (1:4000 dilution, abcam; cat no ab6994) immunostaining. vWF-stained vessels with external diameter less than 50 μm per high-power field were counted in 10 random non-overlapping fields (x200 magnification; n=6/group). Care was taken to avoid fields containing large airways or vessels. F4/80 antibody for macrophage (1:500 dilution, Bio Rad laboratories; catalog number:MCA497G) was used for staining the lung macrophages. Twenty random non-overlapping high-power fields were analyzed and numbers of cells were counted and the average number of cells per high power field was calculated.
RNA isolation: Total RNA from flash frozen lung samples was isolated using the ZYMO mRNA extraction kit. Sample concentration and quality was assayed using Nanodrop-8000 (Thermo Scientific, Wilmington, DE, USA) and further with the Agilent Bioanalyzer. The 260/280 and 260/230 ratios were greater than 1.8. and the RNA Integrity Number (RIN) were between 7-10 and with a range of 1 -1.5.

RNA-Sequencing Library Preparation and Sequencing.
RNA samples underwent quality control assessment using the RNA tape on Tapestation 4200 (Agilent) and were quantified with Qubit Fluorometer (Thermo Fisher). The RNA libraries were prepared and sequenced at the University of Houston Seq-N-Edit Core per standard protocols. RNA libraries were prepared with QIAseq Stranded Total RNA library Kit (Qiagen) using 1000 ng input RNA. mRNA was enriched with Oligo-dT probes attached to Pure mRNA beads (Qiagen). RNA was enzymatically fragmented for 15 minutes, reverse transcribed into cDNA, and ligated with Illumina sequencing adaptors. The size selection for libraries was performed using SPRIselect beads (Beckman Coulter) and purity of the libraries was analyzed using the DNA 1000 tape Tapestation 4200 (Agilent). The indexed libraries were pooled and sequenced Lung transcriptome data analysis: Data was trimmed using trimGalore, it was mapped using HISAT2 (5) to the UCSC mouse build mm10 and gene expression was quantified using featureCounts (6) and the GENCODE gene model (7). Data was further normalized using the RUV method, and differential gene analysis was assessed via the EdgeR method (8), using the R statistical system. We used the significance threshold of FDR-adjusted p-value<0.05, and fold change exceeding 1.5x. Enriched pathways were assessed via the overrepresentation method (as implemented by MsigDB) using the hypergeometric distribution and the Gene Ontology Biological Processes compendium, with significance achieved at FDRadjusted p-value<0.05. Signature correlation across the GTEx healthy lung samples was achieved by first considering genes conserved between mouse and human. Next, genes were converted to z-score across the human samples. Finally, the summed z-scores for a signature and human specimen was computed by adding the z-score of up-regulated genes and subtracting the z-scores of down-regulated genes. Pearson Correlation Coefficient was computed between signatures using the Python scientific library, with significance achieved at p<0.05. Pathway enrichment using Gene Set Enrichment Analysis was also assessed, with significance achieved at FDR<0.25. qRT-PCR analysis: Total RNA was extracted from the lung tissue using Direct-zol mRNA microprep kit (ZYMO Research, R2062). cDNA was prepared using RevertAid Reverse Transcriptase kit (ThermoFisher, k1691). Quantitative PCR was performed using the QuantStudio 7 Flex real-time PCR detection system (ThermoFisher) and SYBR Green (Bio-Rad, 1725275). The thermal cycling conditions used were as follows: one cycle at 95 °C for 1 min, 40 cycles at 95 °C for 15 s, and one cycle at 60 °C for 15 s. The primers used in real-time PCR test are listed in Supplemental Table 1 Relative mRNA levels were calculated using the 2 −ΔΔCT method and normalized to Tbp expression, which was used as the housekeeping gene.

Gene
Forward

Human transcriptome data mining and analysis:
When comparing and contrasting gene signatures (differentially expressed genes, DEGs), it is important to account for the fact that in real tissues there are complex underlying relationships between genes, such as genes in the same pathway or genes at different locations in a signaling cascade of distinct pathways. Also, whereas some housekeeping genes are fairly constant across individuals, other genes exhibit high variability in a population. To capture and quantify the complex underlying dependencies between genes in one individual tissue and variability between individuals, cancer bioinformatics has developed the summed z-scores methodology, applied over the cohort of whole tissue transcriptomes (10)(11)(12). Data for 578 healthy whole lung transcriptomes from human donors was downloaded from the GTEx consortium (13). Genes were next transformed to a z-score across all samples. For each gene signature and each specimen, we obtained the summed z-scores by adding the z-scores of up-regulated genes and subtracting the z-scores of down-regulated genes, via a methodology used extensively in mining of transcriptomic cohorts from cancer patients (11,14). Gene symbols were converted to mouse genome, using the Ensemble BiomaRt reference (15). Essentially, the summed z-score compares each human specimen with the average gene expression in the human cohort across all the Four Core Genotype (FCG) signatures, and then scores each human specimen based on how well it matches each FCG DEG list. With each human specimen in a cohort assigned a summed zscore, it is possible then to do a Pearson Correlation Coefficient analysis. Significant correlations between distinct FCG signatures can then indicate underlying biological relationships. Signature correlations in transcriptomics from GTEx lungs were evaluated using Pearson correlation, with significance achieved at p<0.05.
Data for blood transcriptomes from newborns was downloaded from NCBI GEO, accession number GSE32472 (17). In this cohort, babies were enrolled if gestational age was < 32 weeks and birth weight less than or equal to 1500 gms. The data was collected from 111 preterm neonates with a mean birth weight of 1029 g and mean gestational age at birth of inflammation and lung development in FCG mice exposed to hyperoxia at PND5 and PND21 ( Figure 3A). Percent of enriched GOBP pathways (upregulated and downregulated) that are immune or inflammation related in the FCG mice exposed to neonatal hyperoxia at PND5 and PND21 ( Figure 3B). The fold change of selected genes related to immune/inflammation related pathways in the RNA-Seq experiment was validated in an independent cohort of neonatal mice by qRT-PCR at PND21 (Figure 3C). Values are mean ± SD from 4 individual animals. Analysis done by 3-way ANOVA to assess the effect of treatment, chromosomal sex and gonadal sex as well as the interactions between the independent variables. Significant differences between room air and hyperoxia within genotype are indicated by *P<0.05, **P<0.01, and ***P <0.001. Significant differences between hyperoxia-exposed mice between different genotypes are indicated by # P < 0.05, ## P < 0.01, and ### P < 0.001.

Supplemental Figure 4: Distribution and overlap of gene signatures based on the
gonadal and chromosomal sex: Figure 4 shows the schematic and the number of differentially expressed genes (up-and down-regulated) in each genotype and the overlap based on gonadal or chromosomal sex. Figure 4A shows the distribution at PND5 and Pearson correlation coefficient is shown for significant correlation (p<0.05). This information is also presented in Figure 11A and reproduced here for the reader's convenience ( Figure 5A). Scatterplots indicating the correlation between activity scores of hyperoxia signatures at PND21 and birth weights patients from a BPD cohort ( Figure 5B).

Supp Fig 4
Interestingly, there was a positive correlation between the patient activity scores of chromosomal female mice and birth weight in the male and female patient cohort. Coordinates circadian rhythm and metabolic pathways. Plays a role in the suppression of the endogenous proinflammatory mechanism in the lung (23).