Transcriptome analysis of pigeon milk production – role of cornification and triglyceride synthesis genes

Differentiation of the ‘lactating’ crop

Immunohistological analysis of the proliferating cells of the pigeon crop in its resting state and during nesting demonstrated the morphological changes that occur in preparation for pigeon milk production (Figure 1). As the crop changed in preparation for lactation, the number and depth of rete pegs increased and the lamina propria became progressively more extended and narrow, which increased the surface area of the crop. During lactation the crop was highly proliferative, which resulted in the accumulation and sloughing of large tracts of cornified epithelium (Figure 2). All lactating parents in this study (48 birds) fed their young every four hours over the lactation period. Histology revealed a cycle of production and turnover of cornified epithelium over the four-hour period (Figure 2). The squabs milk intake reduced gradually toward the end of the lactation period, which lasted approximately fourteen days.

Analysis of transcriptional changes over the lactation period compared to non-‘lactating’ crop revealed no differentially expressed probes at pre-hatch (time points −8 and −2), large differences at hatch (time point 0) and 2 days post-hatch (time point +2) (17.2 and 48.8% of all probes differentially expressed, respectively), and no difference above what could be expected by chance (5%) at 10 days post-hatch (time point +10) (2.7% of probes differentially expressed) (Additional file 1: Table S1). Any effect of sex was ruled out by comparing males to females at non-‘lactating’ and ‘lactating’ time points. There was no difference above what could be expected by chance (Additional file 1: Table S2).

Cornification genes are differentially expressed in the ‘lactating’ pigeon crop

Analysis of cornification-associated genes in the draft pigeon genome identified an epidermal differentiation complex (EDC) on scaffolds 1246 and 683, respectively (Figure 3). Transcriptional analysis of these EDC genes and other cornification-associated genes in the pigeon crop at time points 0 and +2 revealed differential expression of 43 genes in 0, +2 or both ‘lactating’ pigeon crops compared with non-‘lactating’ crop (Table 1). Thirteen of these genes were up-regulated and 30 were down-regulated. Notably, the majority of cornification-associated genes up-regulated in the ‘lactating’ crop were keratins, constituting eight of the thirteen up-regulated genes. Five of these eight keratins were beta-keratins and three were alpha-keratins. Conversely, eight of the 30 down-regulated cornification-associated genes were alpha-keratins, and none were beta-keratins. Phylogenetic analysis of the beta-keratins (Figure 4), which were all part of the pigeon EDC (Figure 3) separates them into several groups. Feather, claw and scale keratins share a common ancestor related to chicken beta-keratin, from which feather keratins (none differentially expressed in ‘lactating’ crop) formed their own clade, and claw (ns and 10.5 to 12-fold up-regulated) and scale keratins (ns and 3.2-fold up-regulated) formed another monophyletic clade. Putative pigeon keratins formed three more clades not containing a chicken homolog, and ORF 683_38 formed a clade of its own. GenBank IDs of keratins with the highest amino acid identity to the pigeon keratins are found in Additional file 2.

Phylogenetic analysis of the alpha-keratins separates them into type I and type II (Figure 4). Seven type I keratins and two type II keratins were down-regulated, and two type II keratins were up-regulated in ‘lactating’ crop (Table 1). Notably, all of the type I putative pigeon keratins were constrained to scaffold988, whereas the type I keratins included 15 putative genes on scaffold748 and two on scaffold988 (Figure 4). All of the chicken alpha-keratins had a closely related putative pigeon homolog.

Aside from keratins, there were several other differentially expressed cornification-associated genes in the pigeon crop. Up-regulated cornification-associated genes in the ‘lactating’ crop at 0 and +2 timepoints included transglutaminase 5 (3.1 and 2.7-fold), S100-A10 (1.3 and 1.6-fold) and the pigeon lactation-specific transcript, cp35 (5.7 and 23-fold), and its isoform cp37 (ns and 2.7-fold). Probes on exons 3, 4, 6, 12, 13, 14 and 15 of transglutaminase 4 were up-regulated (1.2 to 1.7 and 1.3 to 2.4-fold), whereas probes on exons 7, 10 and 11 were down-regulated (−1.5 and −1.4 to −1.6-fold), which suggests it is likely to be alternatively spliced in the pigeon crop. The transglutaminase genes transglutaminase 6-like and tissue transglutaminase were also down-regulated by 2-fold and 2.1 to 4.2-fold, respectively. Several of the S100 genes were down-regulated in the ‘lactating’ crop by 1.2 to 10.2-fold, including S100-A14-like, S100-A16-like, S100-A9-like, S100-A4 and S100P (Table 1). The annexin genes annexin A2, A5, A6, A7, A8, A11 and A13 were down-regulated by 1.1 to 21-fold (Table 1). The envelope precursor protein-encoding genes desmoplakin, periplakin, envoplakin, epiplakin and sciellin and cystatin A were all down-regulated in the ‘lactating’ crop by 1.3 to 2.6-fold (Table 1).

The proteases calpain-15 (ns and 6.6-fold) and calpain 9 isoform 1 (ns and 1.8-fold) were up-regulated in ‘lactating’ crop, whilst calpain-5 was down-regulated by 2.4 to 4.4-fold. Cathepsins H, C, B-like, S, Z, L and D were all down-regulated by between 1.2 and 1.8-fold.

The epithelial cell-derived antimicrobial peptide encoding gene beta defensin 5 was up-regulated by 2.3 to 4.8-fold in ‘lactating’ crop at timepoints 0 and +2, respectively (Additional file 3).

Differential expression of cornification-associated genes in the cornified epithelial cells of the pigeon crop

Triglyceride synthesis is up-regulated in the ‘lactating’ pigeon crop

Examination of lipid droplets in ‘lactating’ pigeon crop showed that lipid was present throughout the differentiated epithelium, and was perinuclear (Figure 5). To investigate whether lipids could potentially be synthesised de novo, the expression of genes linked to milk lipid synthesis in the mouse mammary gland [19] were examined in the ‘lactating’ pigeon crop. Thirty-four mouse mammary gland-linked lipid synthesis genes were differentially expressed in the pigeon crop, including 7 variants of genes investigated in the mouse study. Expression patterns of milk lipid synthesis genes were similar in pigeon crop and mouse mammary gland, although pigeon crop expressed different variants of many genes of the triglyceride synthesis and fatty acid synthesis pathways in comparison to the mouse mammary gland. In particular, the triglyceride synthesis genes Agpat1 and Dgat1 were up-regulated in the lactating mouse mammary gland compared to pregnant mouse mammary gland [19], whereas Agpat3, Agpat9 and Dgat2 were up-regulated in the ‘lactating’ pigeon crop compared to non-‘lactating’ crop. The fatty acid synthesis gene Elovl1 was up-regulated in lactating mouse [19], whereas Elovl6 was up-regulated in ‘lactating’ pigeon crop. The lactating mouse mammary gland showed up-regulation of 5 different Fabp gene variants, whereas the ‘lactating’ pigeon crop up-regulated only Fabp5. Both lactating mouse and pigeon crop showed up-regulation of the same fatty acid transporter gene, Slc27a4, the fatty acid translocase, Cd36, and down-regulation of fatty acid transporter Slc27a1.

Pigeon tissue sample collection

Thirty-two breeding pairs of King pigeons were purchased from Kooyong Squab Producers (Moama, New South Wales). They were housed in temperature-controlled cabinets (between 21–24°C) with a 12 hour light cycle (lights on 6 am), and supplied with nest bowls and nesting materials. Pigeons had ad libitum access to pigeon mix (pro-vit-min, Ivorsons, Geelong) and water. Control non-‘lactating’ pairs (ctrl, 13 birds) were culled prior to mating. Breeding pairs were culled at different lactation time points whereby squab hatch was designated as time zero. Time points pre-hatch have the prefix ‘-‘ and post-hatch have the prefix ‘+’. Specifically, breeding pairs were euthanised at 8 days pre-hatch (−8, n = 10 birds), 2 days pre-hatch (−2, n = 10 birds), at hatch (0, n = 14 birds), 2 days post-hatch (+2, n = 10 birds), and 10 days post-hatch (+10, n = 4 birds). Whole crop tissue samples were snap frozen in liquid nitrogen and separate samples of all crops were fixed in 10% neutral buffered formalin or snap frozen in optimal cutting temperature (OCT) compound for histology. Samples of pigeon crop from a time 0 pair were fixed in PaxGene (Qiagen) fixative according to the manufacturer’s instructions for laser dissection microscopy, to investigate gene expression differences between basal and proliferating cell types. Samples of other whole tissues (brain, pituitary, thymus, esophagus, trachea, proventriculus, gizzard, heart, kidney, duodenum, ileum, jejunum, pancreas, spleen, cecum, bone marrow, muscle and skin) were snap frozen in liquid nitrogen and used for the construction of a pooled tissues cDNA library. The blood and spleen of the ten day old squab were removed; the blood into vacutainers coated with EDTA dipotassium salt and the spleen into sterile media (DMEM with 10% FCS, 100 U/mL penicillin, 100 μg/mL streptomycin, 500 μg/mL fungizone).

All work using animals was conducted in accordance with the Australian Code of Practice for the Care and Use of Animals for Scientific Purposes (7th edition), and in accordance with institutional animal ethics guidelines (CSIRO AAHL Animal Ethics Committee).

Pigeon splenocyte stimulation

The squab spleen was minced through a 70 μm filter using a syringe plunger into 15 mL phosphate buffered saline (PBS). The blood was diluted in PBS, and the cell suspensions were layered slowly over the same volume of Lymphoprep (Axis-Shield, Oslo, Norway). After centrifugation the cells were removed from the interface of the gradient and washed twice with 50 mL PBS + 10% FCS. The cells were seeded on a 24-well plate (Nunc) at 5 × 10⁵ cells/mL. To each well 10 μg/mL concanavalin A (Astral Scientific) was added and the cells were incubated at 37°C in 5% CO₂. After 24 hours the cells were pelleted and re-suspended in 1 mL TRIreagant RT (Molecular Research Center) for RNA extraction and synthesis of the immune library.

cDNA library synthesis

Double-stranded cDNA (dscDNA) libraries were syn-thesised from poly(A)⁺ RNA extracted from the female crop tissue samples and the pooled tissues and immune cells. Total RNA was extracted with TRIreagant RT (Molecular Research Center) according to the manufacturer’s instructions, with an additional high salt precipitation buffer step to remove glycoproteins. Poly(A)⁺ RNA was twice purified from total RNA using Dynal beads (Life Technologies), according to the manufacturer’s instructions. 200 ng poly(A)⁺ RNA was fragmented at 70°C for 30 s in a 20 μl reaction containing 18 μl fragmentation solution (100 mM Tris–HCl pH 7.0, 100 mM ZnCl₂). The reaction was stopped on ice with 2 μl of 0.5 M EDTA pH 8.0 and 28 μl of 10 mM Tris–HCl pH 7.5 and cleaned up with an RNeasy column (Qiagen). The sequencing libraries were synthesised from the fragmented poly(A)⁺ RNA according to the Roche cDNA Rapid Library Preparation Method Manual, Rev. Jan 2010 beginning section 3.2. Briefly, fragmented RNA was reverse transcribed into cDNA and made double-stranded using a mixture of DNA polymerase I, ligase and RNase H and blunt ended with T4 DNA Polymerase (Roche cDNA synthesis system).

High throughput sequencing, assembly, microarray design and annotation

Sequencing libraries were prepared from the dscDNA libraries using a Rapid Library Preparation Kit (Roche). Sequencing beads were generated with a SV-emPCR Kit (Roche) and sequenced on a 454 GS FLX using the titanium chemistry (Roche). Each sample was sequenced in a separate region. The raw reads from all regions were combined (430654) and assembled with Newbler v2.3 shotgun assembler (Roche). The resulting contigs (10463) and remaining Singletons (71997) were used to design unique microarray probes with OligoArray 2.1 [36]. The microarray probes were annotated via the source contigs or read sequences by a series of BLAST searches using an E-value of 10^-3 as cut-off for all searches. The first search used BLASTX [37] with all sequences against a local copy of the non-redundant protein database (dated 11 April 2012). All non-matched sequences were then used in a BLASTN query [38] against the non-redundant nucleotide database (dated 13 April 2012). Finally, the remaining unmatched sequences were used as queries in a TBLASTX search [37] against the nucleotide database.

Identification of putative pigeon cornification-associated full-length genes

Cornification-associated genes were identified by literature search, and Raw 454 reads or assembled contigs were used as local megaBLAST [38] queries against the Columba livia draft genome [39] to identify in which scaffold each gene of interest was present. The scaffolds of interest were then submitted to a Hidden Markov Model gene prediction program (FGENESH, Softberry; http://linux1.softberry.com/berry.phtml) using parameters for chicken (aves) to identify predicted full-length gene sequences. Where scaffolds were too large to be processed by FGENESH, the region of the scaffold with the BLAST match was submitted. Microarray probes were mapped to the predicted gene sequence by local BLAST. Non-redundant, non-overlapping microarray probes matching predicted gene coding sequences were identified by megaBLAST against the predicted gene sequences and the Columba livia draft genome sequence.

Phylogenetic analysis of pigeon alpha and beta keratins

Phylogenetic trees were constructed separately for alpha- and beta-keratins. The evolutionary relatedness was inferred using the Minimum Evolution method [40]. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) was calculated [41]. The tree was drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the JTT matrix-based method [42] and are in the units of the number of amino acid substitutions per site. The ME tree was searched using the Close-Neighbor-Interchange (CNI) algorithm [43] at a search level of 1. The Neighbor-joining algorithm [44] was used to generate the initial tree. All positions containing alignment gaps and missing data were eliminated in pairwise sequence comparisons (Pairwise deletion option). Phylogenetic analyses were conducted in MEGA4 [45] after alignment in ClustalX [46].

Laser dissection microscopy and RNA amplification

PaxGene fixed time 0 pigeon crop of a female and male breeding pair were dehydrated through fresh ethanol and xylene using an automated processor (Leica), and embedded in paraffin according to the PaxGene manufacturer’s instructions. Sections of 4 μm were cut by microtome and floated on to laser dissection slides (Leica #11505158 membrane slides PEN-membrane 2 um). The cornified crop epithelial cells and the basal cells of 5 serial sections of each crop were laser dissected using a Leica LMD6000 machine and collected by gravity into 500 μl PCR tubes. The dissected cells were dissolved in QIAzol by pipetting up and down, and RNA was extracted using the RNeasy Lipid Tissue kit according to the manufacturer’s instructions, and eluted in 30 μl water. RNA was quantified using a Bioanalyzer RNA Pico chip, and an equal amount of RNA of each of the four samples was used for two rounds of RNA amplification using an Ambion MessageAmp II aRNA Amplification Kit, according to the manufacturer’s instructions.

Microarray hybridisation, scanning and data pre-processing

RNA was extracted from whole frozen pigeon crop tissue according to the manufacturer’s instructions (Qiagen RNeasy Lipid Tissue kit). RNA quality and quantity was measured using a Bioanalyser RNA Pico chip and 5 μg of this RNA was used to synthesise first-strand cDNA with oligo_dt primer according to the manufacturer’s instructions (Invitrogen SuperScript SuperMix) which was then purified using a PCR purification kit (Qiagen). cDNA was synthesised and purified from whole crop RNA and from amplified laser dissected sample RNA. All cDNA samples were labelled with Cy3 using a Roche One-Color DNA Labelling Kit according to the manufacturer’s instructions. The labelled microarray probes were re-suspended with a sample tracking control and hybridisation buffer and loaded on 12-plex 135 k custom pigeon microarrays (ArrayExpress ID A-MEXP-2257). These were hybridised for 20 hours in a NimbleGen Hybridisation Station (Roche) at 42°C and then washed using the NimbleGen wash buffer kit (Roche) according to the manufacturer’s instructions. Each subarray was scanned at 2 μm on autogain with a NimbleGen MS200 microarray scanner (Roche). Sample tracking controls and control spots were used to autoalign a grid over each subarray using NimbleGen MS200 Software (Roche).

Microarray normalisation and statistical analysis

Robust Multichip Average (RMA) analysis [47] was used to background correct and normalise spot signal intensity. To compare datasets hybridised to different slides, the data were subjected to the non-parametric CombatR algorithm to remove batch effects [48]. The datasets were exported into GeneSpring (Agilent) and differentially expressed genes were identified using an unpaired Welch t-test assuming unequal variances with a Benjamini and Hochberg post-hoc test, with a false discovery rate of p = 0.05. The comparison of cell layers from laser dissected RNA omitted the post-hoc test as there were only two samples per group. All microarray data has been deposited into ArrayExpress (accession numbers E-MTAB-1317 for whole crop and E-MTAB-1318 for laser dissected cell layers).

Proliferating cell nuclear antigen immunohistochemistry

Formalin fixed pigeon crop was dehydrated through ethanol and xylene and embedded in paraffin. Sections of 4 μm were dewaxed in xylene and rehydrated through ethanol. Antigen retrieval was performed using the Dako PT Link (97°C for 30 min) while immersed in Target Retrieval Solution High pH (Dako). Following retrieval the sections were quenched with hydrogen peroxide. Sections were then incubated for 1 h with primary antibody (PCNA, 1:8000; Abcam ab29). This was followed by a Mouse linker (EnVision Flex+, Dako) for 15 minutes to enhance the staining. Horseradish peroxidase conjugated secondary antibody (Envision Flex+, Dako) was then applied for 20 minutes. Sections were stained with 3-amino- 9-ethylcarbazole (AEC) substrate chromogen (Dako) for 10 min, and counterstained with Lillie-Mayer’s haematoxylin.

Oil Red O staining and confocal microscopy

70 μm sections of formalin-fixed ‘lactating’ crop tissue were sectioned by vibrating microtome (Leica VT1200S). Sections were stained with Oil Red O according to the method of Lillie and Ashburn [49] and nuclei were labelled with DAPI for 15 min. Following a water wash, sections were mounted with Vectashield (Abacus, Australia). Samples were imaged sequentially for each dye with a Leica SP5 confocal microscope (Leica Microsystems, Sydney).