Influence of Aromatase Absence on the Gene Expression and Histology of the Mouse Meibomian Gland

Abstract

Purpose.: We hypothesize that aromatase, an enzyme that controls estrogen biosynthesis, plays a major role in the sex-related differences of the meibomian gland. To begin to test this hypothesis, we examined the influence of aromatase absence, which completely eliminates estrogen production, on glandular gene expression and histology in male and female mice.

Methods.: Meibomian glands were obtained from adult, age-matched wild-type (WT) and aromatase knockout (ArKO) mice. Tissues were processed for histology or the isolation of total RNA, which was analyzed for differentially expressed mRNAs by using microarrays.

Results.: Our results show that aromatase significantly influences the expression of more than a thousand genes in the meibomian gland. The nature of this effect is primarily sex-dependent. In addition, the influence of aromatase on sex-related differences in gene expression is predominantly genotype-specific. However, many of the sex-related variations in biological process, molecular function, and cellular component ontologies, as well as in KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, are remarkably similar between WT and ArKO mice. The loss of aromatase activity has no obvious effect on the histology of meibomian glands in male or female mice.

Conclusions.: Our findings demonstrate that aromatase has a significant impact on gene expression in the meibomian gland. The nature of this influence is sex-dependent and genotype-specific; however, many of the sex-related variations in gene ontologies and KEGG pathways are similar between WT and ArKO mice. Consequently, it appears that aromatase, and by extension estrogen, do not play a major role in the sex-related differences of the mouse meibomian gland.

Introduction

Sex steroids exert a significant impact on the meibomian gland. ^1–3 These hormones regulate glandular function, modulate the quality and/or quantity of lipids produced by this tissue, and influence the formation of the tear film's lipid layer. ^1–3 Sex steroids may also prevent (e.g., androgens) or promote (e.g., estrogens) meibomian gland dysfunction (MGD) and evaporative dry eye disease. ^1–3  

Of particular interest, the differential action of sex steroids, as with other sebaceous glands, ^4,5 may contribute to the sex-related differences known to exist in the morphology, gene expression, neutral and polar lipid patterns, and secretory output of the meibomian gland. ^2,6–17 Consistent with this proposal is the observation that androgens appear to mediate almost 30% of the sex-associated variations in gene expression of the mouse meibomian gland. ⁷  

It is quite possible that sex-specific aromatase activities may also play a role in the sex-related differences of the meibomian gland. Aromatase is a cytochrome P450 enzyme that catalyzes the formation of estrogens from androgens, ^18–20 and contributes to a number of sex-specific differences throughout the body. ^21–25 This steroidogenic enzyme is present in numerous tissues, ^18-20 including the meibomian gland. ²⁶ In the absence of aromatase, the synthesis of estrogens is completely eliminated. ^27–29  

We hypothesize that aromatase, through its control of estrogen biosynthesis, is a critical factor in the sex-related differences of the meibomian gland. To begin to test this hypothesis, we examined the effect of aromatase absence on gene expression in male and female meibomian glands. We also explored the impact of aromatase loss on tissue morphology. These studies were made possible by our use of aromatase knockout (ArKO) mice, which were generated by the targeted disruption of exon IX in the cyp19 gene and produce no estrogen. ³⁰ In addition, we used both Affymetrix and CodeLink microarrays in these studies, because these platforms have different gene identification capabilities. ^31–35 The use of different microarrays to analyze the same samples has been proposed to provide a more meaningful transcriptome assessment. ³³  

Materials and Methods

Animals and Tissue Collections

Breeding pairs of C57BL/6J ArKO heterozygous mice were obtained from Dr Orhan Oz (University of Texas Southwestern Medical Center, Dallas, TX). Animals were shipped to Charles River Breeding Laboratories (Wilmington, MA) for initial quarantine, health monitoring, and serology, and then forwarded to the Animal Facilities of the Schepens Eye Research Institute (Boston, MA). Mice were housed and bred in constant-temperature rooms with fixed light/dark intervals of 12 hours' duration. Offspring were genotyped according to a published protocol. ³⁰ In brief, genomic DNA was isolated from tails by using a GenElute Mammalian Genomic DNA Miniprep Kit (Sigma-Aldrich, St. Louis, MO). The DNA was amplified by PCR with a Hybaid OMN-E thermocycler (Thermo Electron Corp, Burlington, Ontario, Canada) by using exon IX gene primers (forward: GTGACAGAGACATAAAGATCG; reverse: GTAAATTCATTGGGCTTAGGG) and neo gene primers (forward: ATCAGGATGATCTGGACGAAGA; reverse: CCACAGTCGATGAATCCAGAA). The PCR conditions were 1 cycle (3 minutes at 94°C), 35 cycles (40 seconds at 94°C, 30 seconds at 55°C, 45 seconds at 72°C) and 1 cycle (5 minutes at 72°C) and amplicons were evaluated on 2.5% agarose gels. Band sizes were 220 bp from wild-type (WT) mice, 170 bp from ArKO mice, and both fragments from heterozygotes. 

When indicated, mice were killed by CO₂ inhalation and either meibomian glands were removed for molecular biological procedures, or lids and globes were obtained for histology. For molecular studies, the meibomian glands were excised from the upper and lower lids under direct visualization with a biomicroscope. This surgical method involved making a small incision near the inner corner of the eyelid, separating skin and subcutaneous tissue from the inner to outer aspect of the lid, and then removing skin from the meibomian glands by cutting at the mucocutaneous junction. After these steps, the palpebral conjunctiva was removed from the meibomian glands, and the glands were dissected from the remaining tissue by starting at the outer lid corner and carefully avoiding an adjacent vein. Glandular samples were prepared by combining tissues from 5 mice/sex/group. Three such pooled sample preparations were made for each tissue/sex/group and then processed for RNA analysis. 

For histological analyses, tissues were fixed in 10% buffered formalin, dehydrated, embedded in methacrylate, cut into 3-μm sections, and stained with hematoxylin and eosin or with periodic acid-Schiff, according to previously described techniques. ³⁶ Sections were obtained from different tissue areas, all separated by minimal distances of at least 30 μm. Slides were observed with an Olympus BH-2 light microscope (Olympus, Center Valley, PA) and micrographs were taken with a Nikon Eclipse E800 (Nikon, Tokyo, Japan) and SPOT camera (MicroVideo Instruments, Inc., Avon, MA). 

All studies with mice were approved by the Institutional Animal Care and Use Committee of The Schepens Eye Research Institute and adhered to the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research. 

Molecular Biological Procedures

Total RNA was extracted from tissues by using TRIzol reagent (Invitrogen Corp., Carlsbad, CA) and further purified with RNAqueous spin columns (Ambion, Austin, TX). Glandular RNA samples were treated with RNase-free DNase (Invitrogen), analyzed spectrophotometrically at 260 nm to determine concentration and examined with a RNA 6000 Nano LabChip (Agilent Technologies, Palo Alto, CA) and an Agilent 2100 Bioanalyzer (Agilent Technologies) to confirm RNA integrity. The RNA samples were then stored at −80°C until processed by using two different procedures. 

The first procedure to evaluate gene expression involved the use of CodeLink UniSet Mouse 20K I Bioarrays (n ∼ 20,000 genes/array; Amersham Biosciences/GE Healthcare, Piscataway, NJ). The RNA samples were processed for CodeLink Bioarray hybridization, according to published procedures. ³⁷ Briefly, cDNA was synthesized from RNA (2 μg) with a CodeLink Expression Assay Reagent Kit (Amersham) and then purified with a QIAquick purification kit (Qiagen, Valencia, CA). Samples were dried, and cRNA was produced with a CodeLink Expression Assay Reagent Kit (Amersham), recovered with an RNeasy kit (Qiagen) and quantified with an UV spectrophotometer. Fragmented and biotin-labeled cRNA was incubated and shaken (300 rpm shaker) on a CodeLink Bioarray at 37°C for 18 hours. Following this time period, the Bioarray was washed, exposed to streptavidin-Alexa 647, and scanned by using ScanArray Express software and a ScanArray Express HT scanner (Packard BioScience, Meriden, CT) with the laser set at 635 nm, laser power at 100%, and photomultiplier tube voltage at 60%. Scanned image files were evaluated by utilizing CodeLink image and data analysis software (Amersham), which yielded both raw and normalized hybridization signal intensities for each array spot. The intensities of the ∼20,000 spots on the Bioarray image were standardized to a median of 1. Normalized data, with signal intensities exceeding 0.50, were analyzed with sophisticated software (Geospiza, Seattle, WA). This comprehensive software also generated gene ontology and z score reports. These ontologies covered biological processes, molecular functions, and cellular components and were organized according to the guidelines of the Gene Ontology Consortium (http://www.geneontology.org/GO.doc.html). ³⁸  

The second procedure to evaluate, and in part to verify, the differential expression of various mRNAs, involved the hybridization of each cRNA (20 μg) sample to a GeneChip Mouse Genome 430A 2.0 Array (Affymetrix, Santa Clara, CA) according to the manufacturer's protocol. Hybridized GeneChips were scanned with an Affymetrix Model 700 Scanner and expression data files were produced from array images using Affymetrix Microarray Suite 4.0 software. GeneChip data were standardized by selecting the default scaling in Affymetrix GeneChip Operating Software, which generates a trimmed mean intensity of 500 for each GeneChip microarray. Normalized data with a quality value of 1.0 were then analyzed with GeneSifter software. 

Gene expression data were assessed without log transformation, and statistical analysis of these data was performed with Student's t-test (two-tailed, unpaired) by using GeneSifter software. Our statistical analyses were not tailored for multiple comparisons. Genes that were expressed in analogous directions in different groups were identified by using “Gene ID” accession numbers and an intersector program (Geospiza; www.public.genesifter.net). Data from the individual CodeLink and Affymetrix arrays are freely accessible for download through the National Center for Biotechnology Information's Gene Expression Omnibus (http://www.genesifter.net/web/dataCenter.html) via series accession number GSE5878. 

Results

Effect of Aromatase Absence on Gene Expression in Male and Female Meibomian Glands

To determine the effect of aromatase absence on gene expression in the meibomian gland, tissues were obtained from adult, age-matched male and female WT (male = 159 ± 6 days; female = 150 ± 7 days) and ArKO (male = 151 ± 7 days; female = 152 ± 8 days) mice (n = 5/sex/group; n = 3 separate experiments). Glands were pooled according to sex and group, processed for the isolation of total RNA, and analyzed for differentially expressed mRNAs by using CodeLink Bioarrays and Affymetrix GeneChips. Microarray data were evaluated with Geospiza bioinformatics software. 

Our results demonstrate that aromatase significantly influences gene expression in the meibomian gland. As shown in Tables 1 and 2, both aromatase presence (i.e., WT) and absence (i.e., ArKO) are associated with a significant up- and downregulation of more than a thousand genes. Genes with the greatest expression differences in terms of ratios included those increased in female WT (Affymetrix [A]: phosphopantothenoylcysteine decarboxylase, 6.5↑; CodeLink [C]: secreted frizzled-related protein 2, 12.2↑) and ArKO (A: hydroxymethylbilane synthase, 4.6↑; C: palladin, cytoskeletal-associated protein, 2.3↑), as well as in male WT (A: phosphopantothenoylcysteine decarboxylase, 5.9↑; C: Crystallin, α B, 3.2↑) and ArKO (A: β defensin 8, 14.6↑; C: tryptophan 2,3-dioxygenase, 1.97↑) meibomian glands. 

Table 1. 

 

Influence of Aromatase Absence on Gene Expression in the Female and Male Mouse Meibomian Gland

Group	Microarray	Genes
		WT > ArKO	ArKO > WT	Total
Female	Affymetrix	783	950	1733
	CodeLink	719	254	973
Male	Affymetrix	856	388	1244
	CodeLink	483	321	804
Male + Female	Affymetrix	837	872	1709
	CodeLink	435	220	655

 The expression of listed genes was significantly (P < 0.05) upregulated between the comparative groups.

Table 2. 

 

Influence of Aromatase on Gene Expression Ratios in the Mouse Meibomian Gland

Accession No.	Gene	Ratio	P value	Ontology
F, WT > ArKO
NM_031188	Major urinary protein 1	3.89	0.0073	Pheromone binding
NM_022984	Resistin	2.89	0.0145	Response to insulin stimulus
NM_134066	Aldo-keto reductase family 1, member C18	2.68	0.0046	Progesterone catabolic process
BM230253	Centromere protein N	2.48	0.0499	Chromosome, centromeric region
NM_054077	Proline arginine-rich end leucine-rich repeat	1.96	0.0453	Cell aging
NM_013459	Complement factor D (adipsin)	1.96	0.0301	Proteolysis
F, ArKO > WT
BB000512	Hydroxymethylbilane synthase	4.62	0.0021	Porphyrin biosynthetic process
BB490338	Calponin 3, acidic	4.59	0.0002	Actomyosin structure organization
AI019214	Transcription elongation factor B, polypeptide 1	4.44	0.0240	Transcription, DNA-dependent
BC002148	Fatty acid binding protein 4, adipocyte	4.33	0.0221	Cytokine production
AK018666	Cysteine-rich repeat-containing protein CRIM1	3.78	0.0406	Regulation of cell growth
AK014590	Kinesin family member 5B	3.53	0.0163	Mitochondrial transport
M, WT > ArKO
BB329391	Phosphopantothenoylcysteine decarboxylase	5.91	0.0002	Coenzyme A biosynthetic process
C77434	Lanosterol synthase	5.16	0.0283	Steroid biosynthetic process
BM115860	Mortality factor 4 like 1	3.87	0.0026	Double-strand break repair via homologous recombination
C85064	V-ATPase E2 subunit	3.74	0.0090	ATP catabolic process
AJ002522	Myosin heavy chain 2X, partial	3.09	0.0110	Muscle myosin complex
BG141806	Ribosomal protein L27a	2.91	0.0267	Translation
M, ArKO > WT
AJ300674	βdefensin 8	14.57	0.0107	Defense response
AF425084	Serine (or cysteine) peptidase inhibitor, clade B, member 6c	3.14	0.0030	Negative regulation of endopeptidase activity
AV006000	Rho GTPase activating protein 5	2.31	0.0456	Positive regulation of mesenchymal cell proliferation
BC003267	cDNA sequence BC003267	2.16	0.0235	Metal ion binding
NM_011726	X-linked lymphocyte-regulated 3A	2.02	0.0029	Biological process
Z67747	Zinc finger protein 62, transcript variant 1	1.88	0.0363	Transcription, DNA-dependent

 Relative ratios were calculated from Affymetrix data by comparing the degree of gene expression in meibomian glands from male (M), female (F), WT, and ArKO groups. Genes listed had a signal intensity of ≥ 100 in at least one group, a comparative P value (between treatments) of < 0.05 and a known identity.

As anticipated, many genes identified as differentially expressed by the Affymetrix and CodeLink arrays were unique to each platform (Table 3). For example, Affymetrix (ratio ≥ 1.80, signal intensity ≥ 500 in at least one group) and CodeLink (ratio ≥ 1.80, signal intensity ≥ 5 in at least one group) platforms identified 72 and 8 genes, respectively, that showed significant differences between female WT and ArKO meibomian glands. Of these genes, only one, small proline-rich protein 2A (ArKO, 2.1-fold↑), was the same in both platforms. Similarly, if all male and female platform-specific data were combined, then only 2 genes, with ratios and intensity thresholds as above, were similar in Affymetrix and CodeLink platforms. These genes were adipsin (WT, 2.0- and 1.9-fold↑) and resistin (WT, 2.7- and 1.9-fold↑). Our findings are consistent with other investigations, which report little agreement between Affymetrix and CodeLink arrays in the detection of differential gene expression. ^31–34,39 However, although these arrays appear to measure different things, ³³ most gene expression changes revealed by each of the platforms are believed to be biologically correct. ^32,33 Given that all our data from the Affymetrix and CodeLink arrays are available in GEO (series number), we will focus predominantly on the Affymetrix data in this Results section. 

Table 3. 

 

Differences between Affymetrix and CodeLink Platforms in the Identification of Differentially Expressed Genes in WT and ArKO Meibomian Glands

Platform	Database Genes	Genes with Ratio ≥ 1.8	Same Response in Other Array
Affymetrix
F, WT vs ArKO	6112	197	10
M, WT vs ArKO	6273	64	1
CodeLink
F, WT vs ArKO	4969	18	1
M, WT vs ArKO	4946	2	1

 Database genes are the number of genes in Affymetrix and CodeLink platforms with signal intensities ≥ 500 and ≥ 5 units, respectively. These threshold intensities, which occur in regions of minimal variation on scatter plots (data not shown), include approximately 40% of the total genes (range: 12,564 to 14,536) detected with the Affymetrix (intensity range: 2.8 to 25,650 units) and CodeLink (intensity range: 0.01 to 1,127 units) arrays. Genes with Ratio ≥ 1.8: total number of genes significantly up- and downregulated by ≥ 1.8-fold on the designated platform; Same Response as in Other Array: number of genes from the “Ratio” column that responded in the same way in the database of the other platform.

The nature of the aromatase impact on meibomian gland gene expression is sex-dependent. As demonstrated by analyses of both Affymetrix and CodeLink data, fewer than 5% of the genes significantly influenced by aromatase were the same in males and females (Table 4). In other words, more than 95% of the aromatase-linked genes are sex-specific. 

Table 4. 

 

Comparative Influence of Aromatase and Sex on Gene Expression in the Mouse Meibomian Gland

Microarray	Comparison 1	Comparison 2	Total Genes	Common Genes	Common/Total (%)
Affy	F, WT > ArKO	M, WT > ArKO	1540	75	4.9
CL	F, WT > ArKO	M, WT > ArKO	1167	35	3.0
Affy	F, ArKO > WT	M, ArKO > WT	1269	51	4.0
CL	F, ArKO > WT	M, ArKO > WT	562	11	2.0
Affy	WT, M > F	ArKO, M > F	2526	508	20.1
CL	WT, M > F	ArKO, M > F	1818	408	22.4
Affy	WT, F > M	ArKO, F > M	1624	174	10.7
CL	WT, F > M	ArKO, F > M	2415	320	13.3

 Following analyses of Affymetrix (Affy) and CodeLink (CL) data, common genes were identified by using the Geospiza Intersector program. The expressions of genes in the “Total” column were significantly (P < 0.05) upregulated in the designated comparison groups.

The absence of aromatase exerted a significant effect on numerous biological processes, molecular functions, and cellular components. As shown by z score analyses, this enzyme deficiency stimulated the expression of male and female genes related to actin binding and non–membrane-bounded organelles, and suppressed genes associated with catalytic activity and protein transport (Table 5). Aromatase absence also influenced many KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways (e.g., ArKO, pentose phosphate pathway, 2.2–3.8↑ in male and female). Of particular interest, most of the affected ontologies and KEGG pathways were sex-dependent (Table 6). For example, analysis of Affymetrix data identified 185 and 85 biological process ontologies (≥30 genes/ontology) that were significantly altered in meibomian glands of females and males, respectively. However, only 22 of these ontologies were the same in each sex. Similarly, evaluation of CodeLink data demonstrated that almost all ontologies (≥20 or 30 genes/ontology) influenced by aromatase absence were different in meibomian glands of females and males (data not shown). 

Table 5. 

 

Similar Gene Ontologies Influenced by Aromatase Absence in Female and Male Meibomian Glands

Ontology	WT Genes ↑	ArKO Genes ↑	WT z Score	ArKO z Score
Cytoplasm	431, 436	445, 188	7.54, 4.82	2.91, 1.60
Cellular metabolic process	372, 380	380, 168	4.54, 3.03	0.66, 1.18
Catalytic activity	258, 259	233, 110	4.43, 2.40	−0.9, 0.58
Ligase activity	33, 28	21, 13	4.25, 2.37	0.37, 1.65
ATP binding	80, 93	77, 33	3.1, 3.83	1.1, 0.68
Intracellular transport	53, 49	56, 25	3.82, 2.5	3.07, 2.28
Intracellular organelle lumen	101, 114	104, 41	3.32, 3.76	1.62, 0.43
Protein transport	63, 59	60, 30	3.54, 2.25	1.71, 2.13
DNA metabolic process	41, 38	34, 15	3.18, 2.05	0.73, 0.63
Neurological system process	19, 27	32, 15	−3.34, −2.42	−2.12, −1.02
Actin binding	15, 22	30, 12	0.45, 1.91	3.71, 2.03
Non–membrane-bounded organelle	115, 121	153, 64	1.99, 1.42	3.7, 2.15
Cellular macromolecular complex subunit organization	20, 28	33, 19	0.29, 1.79	2.43, 3.23
Cellular macromolecular complex assembly	17, 25	29, 18	0.1, 1.74	2.24, 3.47
Cation transmembrane transporter activity	13, 23	13, 4	−1.68, 0.11	−2.34, −2.02

 Specific biological process, molecular function, and cellular component ontologies, with some of the highest and lowest z scores, were chosen following the analysis of nontransformed data. A z score is a statistical rating of the relative expression of genes, and shows how much they are over- or underrepresented in a given gene list. ⁷⁰ Positive z scores represent a greater number of genes meeting the criterion than is expected by chance, whereas negative z scores represent fewer genes meeting the criterion than expected by chance. ⁷⁰ The z scores with values >2.0 or less than −2.0 are significant and are highlighted in bold print. Criteria for inclusion in this table were an ontology containing ≥27 genes, a z score >2.0 or less than −2.0, and similar responses in both female and male meibomian glands. The first number in each column was from female tissues, and the second from male tissues.

 WT Genes ↑, number of genes upregulated in WT tissues; ArKO Genes ↑, number of genes upregulated in ArKO tissues; z Score, specific score for the upregulated genes in the WT and ArKO tissues.

Table 6. 

 

Different Gene Ontologies Influenced by Aromatase Absence in Female and Male Meibomian Glands

Ontology	WT Genes ↑	ArKO Genes ↑	WT z Score	ArKO z Score
Female
Coenzyme metabolic process	23	13	5.05	0.92
Catabolic process	113	103	4.97	1.9
mRNA binding	13	7	4.42	0.95
Posttranscriptional regulation of gene expression	26	18	4.37	1.22
Lyase activity	18	10	4.13	0.51
Cellular ketone metabolic process	52	34	4.01	−0.4
Organelle membrane	80	48	3.98	−2.12
Mitochondrial part	44	30	3.92	−0.07
Integral to membrane	134	174	−2.89	−2.11
Molecular transducer activity	36	76	−4.19	−0.49
Receptor activity	25	61	−4.46	−0.69
Extracellular region	35	79	−4.7	−0.84
Cellular component organization or biogenesis	153	205	1.67	4.12
Cytoskeletal protein binding	26	45	0.75	3.74
Glycosaminoglycan binding	5	17	−0.57	3.57
Immune response	14	51	−2.65	3.56
Protein binding	264	337	1.29	3.35
Protein targeting	17	26	1.33	3.08
Basolateral plasma membrane	10	23	−0.15	3.05
Response to chemical stimulus	62	116	−1.86	2.95
Cytoskeleton	63	91	0.78	2.83
Male
Cellular component organization or biogenesis	188	72	3.91	0.56
Protein modification process	129	41	3.81	−0.5
Regulation of protein kinase activity	34	5	3.62	−1.3
Protein kinase activity	44	7	3.54	−1.53
Intracellular signal transduction	95	32	3.5	0.08
Enzyme binding	63	24	3.38	1.1
Cell junction	46	11	3.31	−0.71
Nuclear lumen	96	38	3.25	0.94
Protein complex binding	30	8	2.6	−0.26
Extracellular region	55	24	−2.97	−2.27
RNA binding	33	28	−0.13	3.32
Cellular component assembly at cellular level	42	26	1.57	2.8
Cation binding	149	86	0.28	2.66
Golgi apparatus	53	31	1.25	2.29
Endomembrane system	42	27	0.56	2.21
Hydrolase activity, acting on acid anhydrides, in phosphorus-containing anhydrides	34	22	0.57	2.15
Ion transport	40	11	−0.35	−2.01
Plasma membrane	168	54	1.14	−2.48
Signal transducer activity	70	20	−0.67	−2.64
Transporter activity	45	10	−0.51	−2.7
Cellular response to stimulus	206	65	1.38	−2.89

 Specific biological process, molecular function, and cellular component ontologies, with some of the highest and lowest z scores, were selected following the analysis of nontransformed data. Criteria for inclusion in this table were an ontology containing ≥ 20 genes, a z score >2.0 or less than −2.0 (in bold print), and different responses in female and male meibomian glands.

For those ontologies that were the same in males and females, most of the genes were also different. For instance, in the protein transport ontology, the expressions of 123 and 89 genes were significantly influenced in females and males, but only 5 of the genes were the same between sexes. Another example was the cellular macromolecular complex assembly ontology, in which 46 genes were significantly affected in the female, and 43 in the male, but only 1 gene was the same in both sexes. 

The ontologies (≥20–30 genes/ontology) identified by the Affymetrix and CodeLink arrays, were, for the most part, unique to the specific platform (data not shown). In those instances in which ontologies overlapped, the inherent gene populations were commonly different between the arrays (Table 7). 

Table 7. 

 

Differences in Ontology Gene Profiles Identified with Affymetrix and CodeLink arrays

Comparison	Ontology	Genes – Affy	Genes – CL	Same – Affy + CL
WT > ArKO	Posttranscriptional regulation of gene expression	34	11	3
WT > ArKO	Transferase activity, transferring acyl groups	34	10	1
ArKO > WT	Mitotic cell cycle	63	20	3
ArKO > WT	Proteasome	13	4	1
M > F	RNA splicing	150	91	43
M > F	Helicase activity	33	16	8
F > M	Chemotaxis	56	73	22
F > M	Polysaccharide binding	39	41	16

 The comparison groups include genes from WT (male plus female), ArKO (male plus female), male (M, WT plus ArKO), and female (F, WT plus ArKO) mice. The first group in each comparison had a significantly increased z score for the given ontology. The ontologies are from biological processes, molecular functions, and cellular components categories. The “Genes – Affy” and “Genes – CL” columns reflect the number of genes, identified by the Affymetrix and CodeLink microarrays, that were significantly up- and downregulated in each ontology. The “Same – Affy + CL” column represents the number of genes that show the same pattern of response in both Affymetrix and CodeLink platforms.

Overall, our data show that aromatase significantly influences gene expression in the meibomian gland. The nature of this effect, though, is predominantly sex-dependent. 

Influence of Aromatase Absence on Sex-Related Differences in Meibomian Gland Gene Expression

To determine whether aromatase absence influences sex-related differences in meibomian gland gene expression, we identified the female and male gene profiles of WT and ArKO meibomian glands and then compared the sex-associated patterns of WT and ArKO mice. 

Our findings show that sex significantly influences the expression of thousands of genes in WT and ArKO meibomian glands (Table 8). Genes linked to X and Y chromosomes demonstrated the greatest differences in terms of ratios (Table 9), but more than 96.6% of the affected genes were linked to other chromosomes (data not shown). In addition, as we discovered during our WT and ArKO comparisons (above), most genes found to be differentially expressed between males and females were platform-specific (Table 10). Thus, Affymetrix and CodeLink arrays both revealed significant sex-related differences in the expression of meibomian gland gene populations, but the gene identities were most often unique to the platform (Table 10). 

Table 8. 

 

Impact of Aromatase Absence on Sex-Related Gene Expression in the Mouse Meibomian Gland

Group	Microarray	Genes
		M > F	F > M	Total
WT	Affymetrix	1844	861	2705
	CodeLink	1014	1399	2413
ArKO	Affymetrix	1271	968	2239
	CodeLink	1216	1458	2674
WT + ArKO	Affymetrix	2371	1261	3632
	CodeLink	1766	2398	4164

 The expression of listed genes was significantly (P < 0.05) upregulated between the comparative groups.

Table 9. 

 

Effect of Aromatase on Sex-Related Differences in Gene Expression in Meibomian Glands

Accession No.	Gene	Ratio	P Value	Ontology
WT, M > F
BB667072	DEAD (Asp-Glu-Ala-Asp) box polypeptide 3, Y-linked	335.89	0.0009	ATP catabolic process
NM_012011	Eukaryotic translation initiation factor 2, subunit 3, structural gene Y-linked	185.55	0.0020	GTP catabolic process
AF127244	Smcy	23.65	0.0003	Chromatin modification
BB742957	Ubiquitously transcribed tetratricopeptide repeat gene, Y chromosome	20.08	0.0067	Chromatin modification
BF537798	Metastasis-associated lung adenocarcinoma transcript 1	5.16	0.0261	Biological process
BE197934	Keratin 14	5.07	0.0011	Epithelial cell differentiation
WT, F > M
L04961	Xist (X inactive specific transcript)	202.1	0.0073	Dosage compensation, by inactivation of X chromosome
BG806300	Domesticus antisense RNA from the Xist locus, complete sequence	113.26	0.0094	Dosage compensation, by inactivation of X chromosome
NM_031188	Major urinary protein 1	4.2	0.0064	Pheromone binding
NM_007753	Carboxypeptidase A3, mast cell	2.64	0.0004	Proteolysis
NM_134072	Aldo-keto reductase family 1, member C14	2.63	0.0006	Metabolic process
BC009165	Thyroid hormone responsive SPOT14 homolog	2.32	0.0020	Transcription, DNA-dependent
ArKO, M > F
BB667072	DEAD (Asp-Glu-Ala-Asp) box polypeptide 3, Y-linked	384.4	0.0004	ATP catabolic process
NM_012011	Eukaryotic translation initiation factor 2, subunit 3, structural gene Y-linked	225.87	0.0002	GTP catabolic process
BB742957	Ubiquitously transcribed tetratricopeptide repeat gene, Y chromosome	118.68	0.0011	Chromatin modification
AJ300674	β-defensin 8	68.91	0.0086	Defense response
AF127244	Smcy	14.25	0.0000	Chromatin modification
BE197934	Keratin 14	3.83	0.0092	Epithelial cell differentiation
ArKO, F > M
L04961	Xist (X inactive specific transcript) mRNA for open reading frame	273.5	0.0001	Dosage compensation, by inactivation of X chromosome
BG806300	Domesticus antisense RNA from the Xist locus, complete sequence	108.04	0.0000	Dosage compensation, by inactivation of X chromosome
NM_007894	Eosinophil-associated, ribonuclease A family, member 1	3.69	0.0075	Nucleic acid binding
BC026054	Neuronal pentraxin 2	2.82	0.0053	Sugar binding
BB031119	Src-like adaptor	2.71	0.0108	Protein binding
BE994609	Myelin basic protein	2.68	0.0171	Myelination

 Ratios were calculated from Affymetrix data by comparing the degree of gene expression in meibomian glands from M, F, WT, and ArKO groups. The genes listed had a signal intensity of ≥100 in at least one group, a comparative P value (between treatments) of <0.05, and a known identity.

Table 10. 

 

Differences between Affymetrix and CodeLink Platforms in the Identification of Sex-Related Differences in Gene Expression of WT and ArKO Meibomian Glands

Platform	Database Genes	Genes with Ratio ≥ 1.8	Same Response as in Other Array
Affymetrix
WT, F vs M	6263	366	20
ArKO, F vs M	6221	290	26
CodeLink
WT, F vs M	5093	331	22
ArKO, F vs M	5022	233	23

 Database genes are the number of genes in Affymetrix and CodeLink arrays with signal intensities ≥500 and ≥5 units, respectively. These threshold intensities include about 40% of the total genes (range: 12,717 to 14,486) detected with the Affymetrix (intensity range: 1.03 to 20,266 units) and CodeLink (intensity range: 0.02 to 980.3 units) platforms. Genes with ratio ≥1.8 are the total number of genes significantly up- and downregulated by ≥1.8-fold on the designated platform. Same Response as in Other Array, number of genes from the “Ratio” column that responded in the same way in the database of the other platform.

The effect of aromatase absence on sex-related differences in gene expression was genotype-specific. As detailed in Table 4, between 78% and 80% of the genes expressed more highly in the male, as compared with female, meibomian glands were different between WT and ArKO mice. Similarly, between 87% and 89% of those genes upregulated in females, relative to males, were different between WT and ArKO mice. 

However, despite these specific gene differences, many, although not all, sex-related variations in ontology categories were remarkably similar between WT and ArKO mice (Table 11). The 17 highest z scores of biological process ontologies (Affymetrix, ≥30 genes/ontology) in the WT male were also all increased in the ArKO male. These included an upregulation of RNA splicing and processing, DNA replication, and cell-cycle processes in both genotypes. In an analogous manner, a number of the increased z scores in biological process ontologies of female glands were elevated in both WT and ArKO mice. These included antiapoptosis, tissue development, cell proliferation, and protein complex biogenesis (Table 11). These similarities in many WT and ArKO ontology responses were also observed in the molecular function and cellular component categories (Table 11; ≥20–46 genes/ontology), as well as in KEGG pathways (e.g., female, focal adhesion, z scores = 2.74 and 6.29; male, cell cycle, z scores = 2.66 and 3.09). Of interest, female ArKO mice also had significant increases in a variety of immune function ontologies, whereas male ArKO mice had increases in mitotic ontologies (Table 12). There were a number of sex-related parallels in the CodeLink data that were similar to those described above for Affymetrix ontologies (e.g., male, glycerophospholipid metabolism, z scores = 3.50 and 2.09). 

Table 11. 

 

Sex-Related Differences in Ontologies in Meibomian Glands of WT and ArKO Mice

Ontology	M Genes ↑	F Genes ↑	M z Score	F z Score
RNA splicing	73, 53	21, 12	8.48, 7.66	2.23, −0.76
RNA processing	106, 73	29, 20	7.99, 6.66	1.26, −1.24
Membrane-bounded organelle	916, 621	408, 402	7.37, 6.08	4.13, 0.07
Nucleotide binding	291, 205	109, 88	7.1, 6.47	2.06, −1.37
Cytoplasm	884, 594	410, 425	6.41, 4.96	4.8, 2.33
ATP binding	196, 134	52, 49	6.1, 5.06	−1.09, −2.24
DNA replication	37, 38	18, 9	2.96, 5.79	2.44, −0.82
Nuclear lumen	205, 136	72, 81	5.51, 4.16	0.76, 0.78
Catalytic activity	551, 274	232, 217	4.59, 3.75	1.63, −1.42
Cell cycle process	89, 77	32, 25	2.14, 4.29	0.02, −1.91
Protein kinase activity	81, 56	21, 27	3.83, 3.28	−0.77, −0.03
Extracellular matrix	19, 7	27, 49	−2.56, −3.37	3.3, 8.25
Extracellular region	79, 36	93, 147	−7.44, −7.78	2.38, 7.57
Cell proliferation	91, 68	68, 98	−3.23, −1.83	2, 5.12
Glycosaminoglycan binding	6, 5	14, 20	−2.42, −1.64	2.99, 4.92
GTPase activity	14, 11	2, 15	−0.86, −0.22	4.9, 2.28
Pattern binding	7, 7	16, 21	−2.48, −1.28	3.3, 4.71
Growth factor binding	11, 5	14, 17	−0.47, −1.26	3.63, 4.48
Antiapoptosis	21, 14	18, 25	−0.56, −0.5	2.46, 4.08
Tissue development	97, 43	61, 75	−0.9, −3.63	2.28, 3.39
Protein complex biogenesis	62, 44	41, 42	0.42, 0.68	3.08, 2.49
Endoplasmic reticulum	111, 73	65, 68	1.21, 0.71	3.08, 2.49

 Biological process, molecular function, and cellular component ontologies, with some of the highest z scores, were selected following the analysis of nontransformed data. Criteria for inclusion in this table were an ontology containing ≥20 genes, a z score >2.0 (in bold print), and similar responses between WT and ArKO mice. The first number in each column was from WT tissues, and the second from ArKO. M Genes ↑, number of genes upregulated in M tissues; F Genes ↑, number of genes upregulated in F tissues; z score, specific score for the upregulated genes in the M and F tissues.

Table 12. 

 

Expression of Various Immune Function and Mitotic Ontologies in Meibomian Glands of Male and Female ArKO mice

Ontology	M Genes ↑	F Genes ↑	M z Score	F z Score
Immune Function
Regulation of immune system process	34	67	−1.37	6.82
Positive regulation of immune system process	20	49	−1.89	6.12
Immune system process	56	97	−2.93	5.4
Regulation of immune response	17	39	−1.54	5.3
Immune response	29	58	−2.13	5.2
Activation of immune response	11	24	−0.78	4.56
Positive regulation of immune response	12	30	−1.67	4.54
T-cell activation	14	28	−0.99	4.27
Defense response	24	55	−3.38	4.01
Lymphocyte activation	18	36	−1.75	3.94
Immune effector process	12	30	−2.04	3.94
Leukocyte activation	21	40	−1.9	3.82
Regulation of leukocyte activation	14	24	−0.82	3.34
Leukocyte differentiation	15	25	−0.96	3.12
Regulation of lymphocyte activation	12	21	−0.91	3.01
Lymphocyte differentiation	12	19	−0.48	2.96
Innate immune response	17	24	−0.47	2.84
Regulation of cytokine production	9	23	−2.24	2.81
Cytokine production	12	25	−2.01	2.65
Response to cytokine stimulus	22	19	−0.49	2.64
Regulation of defense response	15	22	−1.06	2.15
Mitosis
Cell cycle	103	45	4.1	−1.03
Mitotic cell cycle	56	22	3.48	−0.92
Regulation of cell cycle process	29	9	3.27	−0.9
Negative regulation of cell cycle	30	11	3.24	−0.47
Cell cycle phase	54	19	3.15	−1.52
M phase	41	11	2.93	−2.02
Regulation of cell cycle	47	25	2.56	0.17
M phase of mitotic cell cycle	29	10	2.51	−1.03
Mitosis	27	10	2.19	−0.93
Nuclear division	27	10	2.19	−0.93
Negative regulation of cell proliferation	27	33	−0.14	2.97

 Listed biological process ontologies were selected following the analysis of nontransformed Affymetrix data. The z scores with values >2.0 are highlighted in bold print. Criteria for inclusion in this table were an immune function or mitotic ontology containing ≥30 genes in the meibomian glands of M and F ArKO mice. One immune function ontology (i.e., positive regulation of leukocyte activation) was also upregulated in F WT mice. Meibomian glands of M WT mice expressed genes associated with a downregulation of cell cycle ontologies (i.e., ↑ cell cycle arrest and negative regulation of cell cycle; ↓ epithelial cell proliferation and regulation of epithelial cell proliferation). Similar findings related to immune function and mitotic ontologies in males and females were observed with CodeLink data (not shown).

Although a number of sex-associated differences in ontologies were identical between WT and ArKO mice, the specific genes within these ontologies were often not the same. For example, analysis of multiple biological process ontologies in Affymetrix and CodeLink databases showed that about 70% of the genes in a given ontology were different between WT and ArKO mice. Similarly, if comparing various Affymetrix and CodeLink ontologies found to be upregulated in WT and ArKO mice in both platforms, fewer than 9% of the genes were the same within each ontology. 

Impact of Aromatase Absence on the Histological Appearance of Male and Female Meibomian Glands

To assess whether aromatase absence alters the structure of meibomian glands, tissues (n = 4/sex/strain) from adult, age-matched male and female WT (male = 127 ± 0 days; female = 125 ± 2 days) and ArKO (male = 122 ± 5 days; female = 123 ± 2 days) mice were processed for histological evaluation. As shown in the Figure, our analyses revealed no effect of aromatase deficiency on the glandular histology of either male or female mice. In addition, no abnormalities or inflammation were observed in any of the gland sections. 

Figure. 

 

Representative light micrographs of male and female WT and ArKO meibomian glands. There are no obvious histological differences between WT and ArKO glands. Sections were stained with hematoxylin and eosin. All panels are ×40 magnification.

Discussion

The present study demonstrates that aromatase has a significant impact on gene expression in the mouse meibomian gland. The nature of this effect is primarily sex-dependent and genotype-specific. However, many of the sex-related variations in ontologies and KEGG pathways are similar between WT and ArKO mice. Given this finding, it appears that aromatase, and by extension estrogen, do not play a major role in the sex-related differences of the meibomian gland. Of particular interest, the loss of aromatase activity has no obvious effect on meibomian gland histology. 

Our results show that the presence or absence of aromatase significantly influences the expression of more than a thousand genes in meibomian glands of male and female mice. The loss of aromatase, and consequently estrogen, led to significant changes in the expression of numerous ontologies related to biological processes, molecular function, and cellular components, as well as to KEGG pathways, in both male and female glands. These included an upregulation of genes related to actin binding, the pentose phosphate pathway and non–membrane-bounded organelles in ArKO tissues, and a corresponding downregulation of genes linked to such ontologies as catalytic activity and protein transport. The ability of aromatase to significantly affect gene transcripts has also been found in the ovary. ⁴⁰ Further, our finding that aromatase affects gene expression in the male meibomian gland is consistent with investigations showing that estrogen action is extremely important in other male tissues. ^41,42  

We discovered that the nature of aromatase's effect on meibomian gland gene expression is sex-dependent. As demonstrated by both Affymetrix and CodeLink data analyses, fewer than 5% of the genes significantly influenced by aromatase were the same in both males and females. Put another way, more than 95% of the aromatase-linked genes were sex-specific. Similarly, the majority of affected ontologies and KEGG pathways were sex-dependent: most differences between female WT and ArKO glands were not the same as those between male WT and ArKO tissues. This sexual dimorphism in response to aromatase is not uncommon, and has also been reported for aromatase/estrogen regulation of the liver, adipose tissue, brain, and serum leptin levels. ^22,43–47 How this sexual difference in responsiveness may affect meibomian gland physiology in males and females remains to be elucidated. 

In these studies, we used both Affymetrix GeneChips and CodeLink Bioarrays, because we and others have identified significant differences between these platforms in their ability to detect differential gene expression. ^31–34 Indeed, as found again in this investigation, there is little agreement between these platforms concerning the lists of the differentially expressed genes. ^31–35 This low concordance in gene identification seems to be due to intrinsic variations in platform design, including differences in probe length and content, deposition technology, labeling approaches, hybridizing protocols, image segmentation, signal detection, background correction, data normalization, and data mining, ^32–34,39 as well as to the innate instability of lists of significantly changed genes based on P value cutoffs. ⁴⁸ The consequence is that CodeLink and Affymetrix arrays, both of which have proven reproducibility and accuracy, appear to measure different things. ³⁴ Most gene expression changes revealed by each of the platforms are believed to be biologically correct, and these differences cannot be attributed to technological variations. ^32,33 For these reasons, researchers have suggested that a more meaningful transcriptome assessment may be obtained by analyzing the same sample with different microarray platforms. ³³ The genes contained in the intersection of Affymetrix and CodeLink lists can be used as reliable biomarkers, whereas the genes in the union can be used to identify biological pathways. 

Our analyses demonstrate that the effect of aromatase absence on sex-related differences in gene expression was genotype-specific. Thus, up to 90% of the sex-associated genes in meibomian glands of WT mice were different from those in tissues of ArKO mice. However, despite these specific gene differences, many sex-related variations in gene ontologies and KEGG pathways were remarkably similar between WT and ArKO mice. For instance, ontologies, such as RNA splicing and ATP binding, were significantly increased in male WT and ArKO glands, whereas those linked to growth factor binding and endoplasmic reticulum were upregulated in female WT and ArKO tissues. Given these similarities, it appears that aromatase does not play a predominant role in the sex-related differences of the mouse meibomian gland. This conclusion is consistent with a recent finding that estrogen action does not seem to make a major contribution to the sex-associated differences in gene expression of the meibomian gland. ⁷  

Our observation that female and male ArKO mice have significant differences in the expression of a variety of immune function and mitotic ontologies was intriguing, but the underlying bases have yet to be determined. These variations may represent sexually dimorphic responses to the known increases in the serum concentrations of testosterone, androstenedione, prolactin, follicle-stimulating hormone, luteinizing hormone, or leptin in ArKO mice, ^30,45,49,50 or to hormone-independent processes, ^51,52 or to the effects of Y-linked genes or X inactivation in the ArKO strain. ^53–56  

Some limitations in this study should be considered. First, it is possible that tiny pieces of adjacent lid tissue adhered to the meibomian glands during isolation under the dissecting microscope; however, such possible contamination would be similar for all tissues, given that the same individual obtained glands from every lid segment in the same manner. Second, because lid glands were pooled, we do not know whether differences in gene expression exist between the meibomian glands of upper and lower lids. Such differences could theoretically occur, given that glandular anatomy varies between the upper and lower lids of humans. ⁴ Third, our data represent mRNA levels. We have yet to determine whether these transcripts are translated into active proteins. 

Our tissue evaluations revealed no effect of aromatase absence on the meibomian gland histology of male or female mice. In addition, no abnormalities or inflammation were observed in any of the gland sections. This finding is in stark contrast to a report that describes a significant B-cell accumulation in salivary glands of male and female ArKO mice, and concludes that estrogen deficiency results in a lymphoproliferative autoimmune disease resembling Sjögren's syndrome. ⁵⁷ MGD is known to occur in this autoimmune disease, ³ but there is no evidence of MGD in ArKO mice. 

The lack of influence of complete aromatase and estrogen absence on the meibomian gland is also noteworthy, because this tissue is a large sebaceous gland. Estrogens are known to cause a significant reduction in the size, activity, and lipid synthesis of sebaceous glands in multiple species. ^{2,3,4,58–64} In fact, estrogen was once called the prototype agent for the suppression of sebum production, ⁶³ and for years was used to decrease sebaceous gland function and secretion in humans. ^{5,61,62,65–67} A mechanism proposed for this hormone action is that estrogen stimulates the intra-sebocyte release of lysosomal enzymes, leading to premature cellular destruction and reduced sebum output. ^63,68 Another suggested mechanism is that estrogens attenuate testosterone uptake, interfere with testosterone's conversion to dihydrotestosterone, and antagonize androgen activity in the sebaceous gland. ^60,63,69 Consequently, one might anticipate that loss of estrogen might affect the histology of the meibomian gland. Given that such an effect was not found, it may simply reflect the heterogeneity known to exist in the sebaceous gland control mechanisms in different parts of the body. ⁵⁸  

Acknowledgments

The authors thank Orhan Oz (Dallas, TX), Evan Simpson (Clayton, Victoria, Australia) and Britt Bromberg (New Orleans, LA), Alexia Thomas (Dallas, TX), Michael J. Lombardi, Patricia Rowley and Marie Ortega (Boston, MA), and Adam M. Papallo (Cambridge, MA) for their assistance. 

Footnotes

 Supported by National Institutes of Health Grant R01EY05612, the Margaret S. Sinon Scholar in Ocular Surface Research fund, and postdoctoral fellowship awards from Allergan, Japan, Bausch & Lomb, and the Japan Eye Bank Association.

Footnotes

 Disclosure: R. Rahimi Darabad, None; T. Suzuki, None; S.M. Richards, None; R.V. Jensen, None; F.A. Jakobiec, None; F.R. Zakka, None; S. Liu, None; D.A. Sullivan, None