2.1 Ethics statement
All the animal experiments in this study were conducted according to the protocol approved by the Institutional Animal Care and Use Committees of National Chung Hsing University (Taichung, Taiwan).
2.2 Eggs and animals
Pairs of adult zebra finches were purchased from a breeder in Tainan, Taiwan, and their fertilized eggs were collected for the study. The white leghorn chicken was used as the precocial bird model to avoid blocking the signal from in situ hybridization by feather pigmentation. The pathogen free fertilized chicken eggs were obtained from the farm of National Chung Hsing University. All of the eggs used were incubated at 38°C and in relative humidity 65% until the specific stages. The stages and corresponding incubation days of zebra finch embryos followed the description of Murray et al.
(Murray, et al. 2013), and the stages and corresponding incubation days of chicken embryos followed the description of Hamburger and Hamilton (Hamburger and Hamilton 1992b). For example, E8 means embryonic incubation days 8, E12 means embryonic incubation days 12, and D1 means posthatch day 1. The corresponding stages between chicken and zebra finch embryos followed the supplementary description of Abzhanov et al. (Abzhanov, et al. 2004). The chicken and zebra finch
nuclear antigen) antibody was purchased from Chemicon (CBL407), AMV-3C2 (Gag-pro antibody) antibody was from Developmental Studies Hybridoma Bank, and CDH1 (E-cadherin) antibody was from BD Biosciences (610182). DAPI
(4',6-diamidino-2-phenylindole) was used to visualize the nuclei.
2.4 Tissue total RNA isolation
In anterior dorsal (AD) skin regions, the square of four feather buds in length and three feather buds in width was dissected. In posterior dorsal (PD) skin regions, the square of five feather buds in length and two feather buds in width was dissected (Figure 4A). The dissected skin was immersed at 4°C overnight for penetration by RNALater solution (Ambion) and then stored at −20°C before isolation of total RNA.
After thawing, the samples were homogenized by MagNA Lyzer (Roche). Total RNA was extracted using the MasterPure Complete DNA and RNA Purification kit
(Epicentre). The 30 min DNase1 treatment was carried out at room temperature as described in the manual to remove the DNA thoroughly.
2.5 Quantitative PCR
To quantify the candidate gene expressions, the cDNAs were synthesized from the total RNA by QuaniTect Reverse Transcription kit (Qiagen). Each cDNA sample
containing SYBR green (KAPA SYBR FAST qPCR kit) was run on LightCycler 480 (Roche) at annealing temperature 63°C for 45 cycles. Quantification of the TATA box binding protein (TBP) RNA was used to normalize target gene expression levels. All the PCR primers are listed in Table A1.
2.6 mRNA whole mount in situ hybridization
Gene-specific fragments were amplified from RNA extracted from dorsal skins of chicken and zebra finch embryos and subsequently cloned into pGEM-T Easy vector
system (Promega, A1360). Both antisense and sense RNA probes were made by in vitro transcription according to manufacture’s instructions (Roche, Cat #11277073910).
Whole mount in situ hybridization was performed using non-radioactive in
situ hybridization according to the procedure described in Chuong et al. (Chuong, et al.
1996). PCR primers for the cDNA amplifications are listed in Table A2.
2.7 Stranded RNA sequencing
At E8, E9, and E12 zebra finch, the skin total RNAs were pooled from 7, 5, and 3 individuals, respectively. Total RNA concentrations from six libraries (E8A: E8 anterior dorsal skin; E8P: E8 posterior dorsal skin; E9A: E9 anterior dorsal skin; E9P: E9
posterior dorsal skin; E12A: E12 anterior dorsal skin; E12P: E12 posterior dorsal skin) were measured by Qubit fluorometer (Invitrogen USA), and quality was assessed by BioAnalyzer 2100 RNA Nano kit (Agilent, USA). The Illumina library construction and sequencing was conducted by High Throughput Genomics Core of the Biodiversity Research Center, Academia Sinica, Taiwan.
2.8 Data processing and reads mapping
Low-quality bases and reads were removed by using Trimomatic version 0.30 (Bolger, et al. 2014) according to the following procedure: (i) remove adaptors, (ii) remove leading and trailing bases with Phred quality score smaller than 20 (Ewing and Green 1998), (iii) scan the read with a 4-base wide sliding window, cutting when the
mapped to the genome using Tophat version 2.0.8 (Trapnell, et al. 2009), and its embedded aligner Bowtie version 2.1.0 (Langmead, et al. 2009) with the following parameters: -N 3 --read-edit-dist 3 --no-novel-juncs --library-type fr-firststrand. The normalized expression levels of genes, represented by fragments per kilobase of exon per million fragments mapped (FPKM) (Mortazavi, et al. 2008), were generated by Cufflinks version 2.1.1 (Trapnell, et al. 2013) with the following parameters:
--max-bundle-frags 1012 --multi-read-correct --library-type fr-firststrand.
2.9 Clustering analysis and identification of differentially expressed genes
A gene is said to be expressed if its FPKM value is higher than 1 in at least one of the six transcriptomes. The Pearson correlation coefficient (PCC) between the 6
expression profiles of the genes was used as the distance metric for gene expression differences (PCC > 0.8). All expressed genes were hierarchically clustered by the WPGMA (Weighted Pair-Group Method with Arithmetic mean) method using
heatmap.2 function in the “gplots” package of R (Warnes, et al. 2009). The cut-off for the cluster analysis is given in Figure A1. I identified the differentially expressed genes (DEGs) through three sets of comparisons. Gene expressions between AD and PD skin samples in E8 and E9 libraries were compared. To increase the power of detecting the DEGs with low expression, the transcriptomes of E8 and E9 AD skins were used as the AD replicate, while the transcriptomes of E8 and E9 of PD skins were used as the PD replicate. These two replicates were compared (E8A+E9A versus E8P+E9P). Here we skipped the samples in E12 because the natal down growth was stopped at E12.
The DEGs from the comparisons were computed by NOISeq (Tarazona, et al.
2011). Only the genes with q > 0.7 (odds value) were defined as DEGs (Liu, et al.
2013).
2.10 Gene set enrichment and pathway analysis
To search the possible pathways involved in natal down growth regulation, the Ensemble gene ID of the expressed genes were converted to the ID of their chicken homologs and input into g:Profiler (Reimand, et al. 2007; Reimand, et al. 2011), a web-based toolset for functional profiling of gene lists from large-scale experiments.
Biological process, cellular component, molecular function, reactome and human phenotype were used as the dataset. The p-value of the gene enrichment was corrected by Benjamini-Hochberg FDR (false discovery rate). Only the gene ontology with the corrected p-value < 0.05 was used in further analyses.
2.11 Functional studies
For the generation of proviral constructs, full-length cDNA PCR products were cloned into the pCR8/GW/TOPO Gateway entry vector (Invitrogen, Carlsbad, CA) and sequenced. The cDNAs were transferred into a Gateway compatible RCASBP-Y DV vector through an LR recombination reaction (Loftus, et al. 2001). Virus was made according to Chuong et al. (Chuong 1998) concentrated by ultra-centrifugation. For an in vivo assay, RCAS virus directing the expression of the candidate genes was injected into the leg or anterior dorsal skins in E3 chicken embryos. Samples were harvested at E12. At least three independent experiments were conducted for each candidate gene.
The primer pairs for the full-length coding sequence amplification were listed in Table