text:
Small RNA and degradome profiling reveals miRNA regulation in the seed germination of ancient eudicot Nelumbo nucifera

Background MicroRNAs (miRNAs) play important roles in plant growth and development. MiRNAs and their targets have been widely studied in model plants, but limited knowledge is available concerning this small RNA population and their targets in sacred lotus (Nelumbo nucifera Gaertn.). Results In this study, a total of 145 known miRNAs belonging to 47 families and 78 novel miRNAs were identified during seed germination using high-throughput small RNA sequencing. Furthermore, some miRNA families which have not yet been reported in monocot or eudicot species were detected in N. nucifera, indicating that these miRNAs was divergence from monocots and core eudicots during evolution. Using degradome sequencing, 2580 targets were detected for all the miRNAs. GO (Gene Ontology) and KEGG pathway analyses showed that many target genes enriched in “regulation of transcription” and involved in “carbohydrate”, “amino acid and energy metabolism”. Nine miRNAs and three corresponding targets of them were further validated by quantitative RT-PCR. Conclusions The results present here suggested that many miRNAs were involved in the regulation of seed germination of sacred lotus, providing a foundation for future studies of sacred lotus seed longevity. Comparative analysis of miRNAs from different plants also provided insight into the evolutionary gains and losses of miRNAs in plants. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-3032-4) contains supplementary material, which is available to authorized users.

Background
MicroRNAs (miRNAs) are approximately 21 nucleotides endogenous non-coding RNAs that play important roles in the regulation of gene expression at the post-transcriptional level, resulting in the cleavage or translation repression of target mRNAs [1,2]. In plants, miRNA-guided gene regulation is involved in multiple developmental processes, such as floral development [3], male and female sterility [4,5], leaf growth [6], and RNA metabolism [7]. Most plant miR-NAs have perfect or near-perfect sequences complementary with their target mRNA to regulate gene expression at the post-transcriptional level [8]. Thus, identifying miRNAs and their targets, and elucidating their regulatory mechanisms are critical to understand the development processes.
Recently, degradome sequencing has been successfully used to analyse the miRNA targets. Combined with small RNA sequencing, the pace of miRNA and their targets identification have been greatly improved in plants [9]. Other studies have also suggested that wholescale miRNA identification and targets analysis might reveal regulatory networks in plant growth and development [10]. Lately, degradome sequencing has been used to find target genes on Arabidopsis, rice, cotton and so on [11][12][13]. However, the miRNA targets in sacred lotus (Nelumbo nucifera Gaertn.) have only been predicted by bioinformatics method, and little information has been experimentally validated [14,15].
Sacred lotus is a perennial aquatic plant with ornamental, edible, medicinal value and phylogenetic importance [16]. In Asia, sacred lotus has been cultivated for more than 7000 years due to its beautiful flowers, edible rhizomes and seeds [17,18]. Recently, the genome sequences of the sacred lotus variety 'China Antique' have been sequenced, which is available for further identifying miRNAs in N. nucifera [19,20].
As an ancient basal eudicots, sacred lotus seed is one of the oldest directly dated seeds for~1300 year-old. Previous studies reported that 16 old fruit of 'China Antique' were collected from the same lakebed and 84 % of them can be germinated in~3 days [21,22]. The proteins which were heat stable and ROS removal were reported to be important for seed longevity [21][22][23]. ROS adversely affects cellular proteins and enzymes and renders them inactive. Some thermal proteins, including superoxidae dismutase (SOD) and HSP were demonstrated to reduce deleterious ROS accumulation for improving the seed longevity [22,23]. Because of the sacred lotus seeds remain viable after hundreds of years, resistance and damage repair in seeds must be extraordinarily effective [24]. Therefore, sacred lotus seeds could be an excellent model for studying the adaptive mechanisms of seeds [25].
Seed germination is the first step for plant growth and plays a pivotal role in seedling establishment. It is a complex physiological and biochemical process and previous studies suggested a role for miRNAs in seed germination [26][27][28][29]. Using small RNA sequencing, the miRNA expression patterns in rice embryos at 0, 12 and 24 h after imbibition showed that osa-miR159f, osa-miR166l and osa-319b were predominantly expressed in rice seeds [27]. In soybean, gma-miR1530 and gma-miR1536 have been implicated in the conversion between photosynthesis and lipid accumulation by regulating their target genes of transketolase and carboxylase, respcetively [28]. In maize, 115 know miRNAs were identified in the seeds at 24 h after imbibition, and the regulation of the genes targeted by these miRNAs was involved in the early stage of seed germination [29]. However, little miRNAs were known about the seed germination of sacred lotus. Therefore, identification of novel miRNAs and elucidation of their functions in the seed germination will help us to understand the long living of sacred lotus seed. Furthermore, as a basal eudicot, identification of the sacred lotus miRNA and comparison with other angiosperms will aid the study of conservation and evolution of miRNAs in plants.
To comprehensively investigate small RNAs and their targets and provide some information for further understanding the miRNA-mediated regulation network in the seed germination of sacred lotus, five small RNA libraries and a degradome library were constructed in this study. The profiling of the miRNAs and their target genes provides some insights into the regulatory pathways governing gene expression during the seed germination, and the data presented here will lay a foundation for future studies of sacred lotus seed longevity.

Overview of small RNA library sequencing
To elucidate the regulatory roles of miRNAs during the seed germination of sacred lotus, five small RNA libraries were constructed from germinating seeds at 0 h, 12 h, 24 h, 36 h and 72 h after imbibition and then sequenced by a Hiseq2,500 platform. After filtering the low-quality reads, including reads <18 nt in length, 5′ adaptor contaminates, reads with polyA, and reads without 3′adaptors, ranging from 10,245,283 clean reads (at 72 h) to 14,412,672 clean reads (at 0 h) were obtained. More than 60 % of the clean reads were successfully mapped to the sacred lotus reference genome (Table 1). Moreover, 3,347,363, 4,084,098, 3,575,988, 4,002,600 and 3,497,201 clean reads were unique mapped to the genome. Other unannotated unique sequences were used for novel miRNAs prediction and further analysis (Table 1).
In all of the five libraries, the small RNA (sRNA) sequences were ranged from 18 to 28 nt in length, with the majority were 19 to 24 nt. And the 24 nt sRNAs were the most abundant within the total sRNA, which approximately accounted from 20 % to 35 % in all of the five sRNA libraries, followed by 21 nt sRNAs, the typical length of canonical miRNAs (Additional file 1: Figure S1).
The common and specific of total and unique small RNA sequences were compared between the adjacent libraries (0 vs12h, 12 h vs 24 h, 24 h vs 36 h and 36 h vs 72 h), showing that more than 75 % of the total sRNAs were shared (Additional file 1: Figure S2). However, only a small fraction (approximately 12.5 %~14 %) of unique sRNA shared in adjacent libraries, suggesting that these small RNAs were low in abundance but highly in diverse.

Identification of miRNAs and compared with other plants
Conserved miRNA families have been identified in many plant species and have important functions in plant development [29]. In the present study, 145 known miR-NAs (belonging to 47 families) were identified (Table 2 and Additional file 2: Table S1). The unannotated small RNA tags were used to predict novel miRNAs and their secondary structures of the precursors were further detected by MFOLD or RNAfold to remove the false positive (Additional file 1: Figure S3). Finally, 78 novel miRNAs were identified in sacred lotus (Additional file 1: Figure S3, Additional file 2: Table S2). The largest miRNA family size identified was Nnu-miR396 comprising 15 members, following by Nnu-miR169 and Nnu-miR393 possessed 12 and 8 members, respectively. Other miRNA families such as Nnu-miR397, Nnu-miR403, Nnu-miR529 and Nnu-miR827 possessed only one member detected in sacred lotus during seed germination of sacred lotus (Table 2 and Additional file 2: Table S1).
A comparison of the identified known miRNA families from sacred lotus with other three different plants (Arabidopsis, rice and Amborella) showed that the number of miRNA families widely varied (Table 2) [30,31]. In sacred lotus, most of the families were represented by from 1 to 8 variants, except for miR396 (15) and miR169 (12). More than 10 variants of miR169 family were also detcted in Arabidopsis and rice ( Table 2). When comparing 25 miRNA families of sacred lotus miRNAs to those in 10 other angiosperm species, five families (miR165, miR403, miR828, miR837 and miR2111) were only present in eudicot species. And other two families (miR1432 and miR2275) were found in basal angiosperm and monocots but not detected in core eudicot species (Fig. 1). MiR529 was identified in basal eudicots (Eschscholzia californica and N. nucifera) and most of the monocots, but no homologues were found in core eudicots (Fig. 1).
To validate the predicted novel miRNAs, stem-loop RT-PCR was performed for six novel miRNAs (novel_ miR_2, novel_miR_13, novel_miR_29, novel_miR_32, novel_miR_36 and novel_miR_66) with relatively high abundance at all the five stages of seed germination (Fig. 2, Additional file 1: Figure S3 and Additional file 2: Table S2). As a result, all the selected six novel miRNAs were found to be expressed in the seed germination (Fig. 2b).

Differentially expressed miRNA at different stages of seed germination
To compare the different miRNA expression profiles during the seed germination of sacred lotus, the differential expression analysis of the miRNAs was performed between the four stages with 0 h (CK) of germination, based on the normalized read count for each identified miRNA (Additional file 2: Tables S3 and S4). And a total of 107 differentially expressed (DE) miRNAs were found at the four stages during sacred lotus seed germination. Then cluster analysis was conducted to further elucidate the expression patterns of all known miRNAs, based on the criteria (at least one comparison has a fold change value ≥ 2.0 or ≤ −2.0 with a p-value < 0.05) (Fig. 3a). Compared with the 0 h of seed germination, many DE miRNAs were detected in the other four stages (Additional file 2: Table S4). Furthermore, eight DE miR-NAs were differentially expressed in all the four RNA libraries, while 26 DE miRNAs were expressed in at least two of the four small RNA libraries, and 17 miRNAs were identified only once in the four samples (Fig. 3b). These results indicated that a larger number of miRNAs were differentially expressed during all five stages of seed germination. Compared with germinating seed at 0 h, a total of 35 DE miRNAs were up-regulated, while 16 miRNAs were down-regulated (Table 3). For instance, Nnu-miR156c-5p and Nnu-miR159a were significantly expressed at 36 h (Table 3). Nnu-miR165a-5p, Nnu-miR166c, Nnu-miR166d and Nnu-miR2111b showed the highest expression levels at 24 h, while Nnu-miR396b-3p showed maximum expression level at 0 h ( Table 3). Most of the novel miRNAs showed low expression levels and some of them were only found at one stage during the seed germination. However, novel_mir_13, novel_-mir_33 and novel_mir_39 were highly expressed during all the five stages (Additional file 2: Table S4) and stemloop RT-PCR also confirmed the expression of novel_-mir_13 during sacred lotus seed germination (Fig. 2b).

Degradome sequencing analysis
In this study, a degradome library was constructed from the germinating seed of sacred lotus, followed by degradome sequencing. After discarding the low quality sequences, a total of 13,852,945 clean reads with 4,362,200 unique clean reads was remained. And 99.53 % (13,787,820) of the total sequences, corresponding to representing 99.35 % (4,333,854) of the unique sequences, were mapped to the reference genome to identify the fragments of degraded mRNAs. Cleavage of target mRNAs of the miRNAs were identified by degradome sequencing (Fig. 4 and Additional file 1: Figure S4). Generally, the target genes of kown miR-NAs were mostly transcription factors [32], and similar results were also obtained for the novel miRNAs (Additional file 2: Tables S5 and S6). For example, both Nnu-miR159b and Nnu-miR319c targeted mRNA encoding the MYB transcription factor. Furthermore, Nnu-miR156a, Nnu-miR160a, Nnu-miR160a-5p and Nnu-miR169a targeted mRNA encoding squamosa promoter-binding-like protein 17, auxin response factor 18, NAC transcription factor and nuclear transcription factor Y subunit A-10, respectively (Fig. 4, Additional file 1: Figure S4 and Additional file 2: Table S5). The targets of novel miRNAs were also identified by degradome sequencing and many of them were transcription factors, including GAMYB (novel_mir_13), MYB (novel_mir_32), NAC (novel_mir_36) and TCP (novel_mir_27) (Additional file 2: Table S6).

Functional annotation of the target genes
To further perceive the putative function of the predicted target genes, Gene Ontology (GO) analysis was performed. And twenty-five different biological processes, fifteen different cellular components and ten different molecular functions were predicted (Additional file 1: Figure S5 and Additional file 2: Table S7). The most significantly enriched GO terms were involved in the three main categories, "protein phosphorylation", "nucleus", and "ATP binding", followed by "regulation of transcription", "membrane" and "DNA binding" (Fig. 5a and Additional file 1: Figure S5). These results suggested that regulation of transcription played a vital role in the process of the seed germination of sacred lotus. Moreover, all the targets could be used for KOG classifications which were classified into 23 groups (Additional file 1: Figure S6). The most dominative group was the "general function prediction only", followed by the "transcription" (Additional file 1: Figure S6).
In addition, KEGG analysis was further conducted to elucidate the biological pathways of the DE miRNA target genes, which were classified into five groups including "organismal system", "metabolism", "genetic information processing", "environment information processes" and "cellular processes" (Fig. 5b and Additional file 2: Table  S8). The categories of "carbohydrate metabolism", "amino acid metabolism", "energy metabolism" and "biosynthesis of other secondary metabolism" were the most enriched pathways, which were all in the group of "metabolism", suggesting that the miRNA regulating the expression of metabolic genes played crucial roles in the seed germination.

Quantitative RT-PCR validation of miRNAs and their target genes
In this study, real time quantitative RT-PCR (qRT-PCR) was performed to confirm the expression levels of nine significantly DE miRNAs and three respective predicted target genes, obtained from the high-throughput sequencing ( Fig. 6 and Additional file 2: Table S5). Three miRNAs (Nnu-miR393a-5p, Nnu-miR396b-5p and Nnu-miR397) showed the similar expression pattern in terms of read abundance during the five stages of seed germination ( Fig. 6 and Additional file 2: Table S3). For instance, Nnu-miR156a showed the highest expression levels at 0 h, and then gradually declined from 12 h to 72 h (Fig. 6g). Moreover, the expression levels of Nnu-miR393a-5p, Nnu-miR394a, Nnu-miR396b-5p, Nnu-miR397 and Nnu-miR2111a were up-regulated and maximized at 24 h and then dramatically decreased to relatively low levels at 36 h and 72 h (Fig. 6). Some of the miRNAs showed the minor discrepancies between sequencing data and the qPCR results, which partially due to the differences in the sensitivity, specificity and algorithm between the two techniques. The expression levels of some of the predicted target genes (SPL16, TCP2 and FBX) were selected for validation using qRT-PCR, and multiple of target genes showed contrasting expression patterns with their corresponding miRNAs ( Fig. 6 and Additional file 2: Table  S5). For example, TCP2 (NNU_00026-RA) was the predicted target gene of Nnu-miR319c which expression level was sharply declined to the lowest level at 24 h, then gradually increased and remained at a relatively high level at 72 h. On the contrary, the target (TCP2) of the Nnu-miR319c was increased at the 24 h and then gradually declined at 72 h ( Fig. 6h and k). Moreover, Nnu-miR156a also showed a negative relationship with its target SPL16 ( Fig. 6g and j). The results suggested that some miRNAs might be invovled in the regulation of seed germination of sacred lotus [33]. However, the expression of FBX showed partially negative relationship with Nnu-miR394a and both of them showed high expression levels at 24 h ( Fig. 6i and l). Thus, it is likely that the inverse relationship between miRNAs and their targets might be limited during specific stages of germination, because miRNAs are not the only regulatory factors affecting their targets [34]. Since miRNAs have been shown to also act by translation inhibition in plants, the fact that the RNA expression level of a target  is not negatively correlated with miRNA expression level does not mean that the target is not real.

Putative microRNA editing in sacred lotus
The miRNA editing patterns showed a significant discrepancy among the germinating seeds at 0 h, 12 h, 24 h, 36 h and 72 h (Fig. 7a,b). For example, the most dominant nucleotide substitution type were at 12 h and 24 h with T to C (more than 35 % in both libraries); while the most common nucleotide substitutions were at 0 h and 36 h with A to G and C to T (Fig. 7b). At the nucleotide positions of 4 and 11, the miRNA editing events showed the specifically lowest number (less than 5 number of editing) (Fig. 7a, b and Additional file 2: Table S9). These results were consistent with previous reported miRNA editing in rice grain filling [35].
To validate the occurrence of miRNA editing events in sacred lotus, two editing types (T to A, and U to C) were examined in Nnu-miR160a and Nnu-miR319c, respectively (Fig. 7c). Using PCR amplification, precursor miRNA sequences from genomic DNA and mature miRNA sequences from cDNA were cloned and sequenced, respectively. Comparison analysis of the sequences further confirmed that the miRNA editing events really existed in sacred lotus (Fig. 7c).

Discussion
As non-coding RNAs, miRNAs play key roles in plant growth and development. Many studies of miRNAs have been reported in model plants such as rice, Arabidopsis and maize and so on. Sacred lotus is a basal eudicot with ornamental, agricultural, medicinal importance. Previous studies on miRNAs in sacred lotus were only in leaves or flowers and for bioinformatics searching without experimentally validation [14,15]. Phased small interfering RNAs (phasiRNAs) which regulated by miRNAs were also identified in Chinese scared lotus by bioinformatics analysis [36]. However, the number of miRNAs identified in sacred lotus is very few and little is known about the miRNAs involved in the seed germination of N. nucifera. In this study, 148 known miRNAs (belonging to 47 families) and 78 novel miR-NAs were identified in sacred lotus via deep sequencing and the results were validated by experimental approaches.
Most of the sacred lotus miRNAs identified in the study were present in monocots and at least one of the core eudicot species, indicating their existence in the ancestor of the two large angiosperm clades. Eleven miRNA families were conserved among monocots, basal angiosperm, basal eudicots and core eudicots, suggesting that these miRNAs were involved in important gene regulation pathways in plants (Fig. 1). However, some miRNAs might become lost or evolve within some angiosperm lingeages [30]. MiR165 and miR828 were found in Arbidopsis and two basal eudicot species, but not reported in monocots, while miR529 was present in monocots but not found in the core eudicot species (Fig. 1). That may be because the last common ancestor of monocots and eudicots likely had miR529, with the gene appearing near the speciation of that last common ancestor, but the gene was lot in the ancestor of core eudicots [30]. In addition, most of the target mRNAs of miR165 and miR166 as well as miR156 and miR529 were the same, respectively [33,34]. The loss of miR165 or miR529 in monocot or core eudicot species indicated the distinct evolution of these species. Moreover, miR1432 and miR2275 were only present in monocots and N. nucifera, suggesting that they might be ancient miRNAs in angiosperm. Further investigation of these miRNAs and their targets would give some insights into the functional divergence and evolution in angiosperm.
Using small RNA sequencing, many novel miRNAs with low abundance are discovered from various developmental stages. In this study, 78 novel miRNAs were identified and some of them were validated in the sacred lotus during seed germination (Fig. 2). The expression profiles of these novel potential miRNAs during the different stages might provide valuable information for seed germination and seed longevity in sacred lotus. In addition, a large number of miRNA editing events were found in the seed germination of sacred lotus. Analysis of the nucleotide positions of editing sites showed that the 4 and 11 sites had less editing frequency and positions of 6 and 7 sites had more editing frequency (Fig. 7a). The results were in accordance with rice grains, while in rice young panicles, the two sites had more editing events [5,35]. These results indicated that miRNA editing was diversity in different tissues and species.
The computational analysis using the cleaveland pipeline revealed a total of 2580 targets for these known and novel miRNAs by degradome sequencing in the present study. Our results showed that most miRNA families had multiple target genes, while several miRNAs from different miRNA families shared some target genes. GO analysis showed that these targets mainly involved in "protein phosphorylation" and "regulation of transcription" (Additional file 1: Figure S5). For KEGG pathway, "carbohydrate metabolism" was enriched in the germinating seeds of sacred lotus (Fig. 5b), indicating that many carbohydrate metabolism related genes were regulated by miRNAs during the seed germination in sacred lotus. "std" represents normalized expression level (Normalised expression, RPM) of miRNAs "*" and "**" represnent significant differentially expression between 0 h and other hours, *p < 0.05, **p < 0.01 In many plant species, plant hormones including Gibberellins (GAs) and ABA, can promote and inhibit seed germination. And some studies also indicated that cross talk between auxin and GA or ethylene might affect seed germination [29,37]. In Arabidopsis thaliana, overexpression of MIR160 (35S:MIR160) resulted in reduced ABA sensitivity during germination. Transcriptome analysis of germinating ARF10 and the transgenic plants expressing ARF10 with silent mutations (mARF10) revealed that many ABAresponsive genes were overexpressed in germinating mARF10 seeds [38]. Furthermore, miR160 has been shown to regulate the expression of ARF17 which is another auxin response transcription factor [39]. In this study, the Nnu-miR160 also targeted the ARF17 (NNU_16091-RA) and the expression of the three miRNAs of the miR160 families was up-regulated at 24 h during the seed germination (Additional file 2: Tables S3 and S5). These results indicated that overexpression of miR160 might promote seed germinating by decreasing the ABA levels in the seed germination of sacred lotus. Auxin signal pathways are also regulated by miR393 which targets the auxin receptor TIR1 in Arabidopsis [40]. Two TIR1 gene homologues (NNU_01291-RA and NNU_06904-RA) were also detected for the targets of Nnu-miR393 in our degradome dataset (Additional file 1: Figure S4 and Additional file 2: Table S5). The expression level of Nnu-miR393a-5p was also validated by qRT-PCR Fig. 6 Validation of differentially expressed miRNAs and some of their targets during the seed germination in sacred lotus by qRT-PCR a Nnu-miR157a, b Nnu-miR168a-5p, c Nnu-miR393A-5p, d Nnu-miR396b-5p, e Nnu-miR397, f Nnu-miR2111a, g Nnu-miR156a, h Nnu-miR319c, i Nnu-miR394a, j Target of Nnu-miR156a (NNU_02036-RA, SPL16), k Target of Nnu-miR319c (NNU_00026-RA, TCP2), l Target of Nnu-miR394a (NNU_22183-RA, FBX) (Fig. 6c). The down-regulation of some components in auxin signal transduction by miRNAs might be a regulatory step to decrease ABA sensitivity in mature seeds and to switch to the germination mode [29]. MiR159 has been reported to target GAMYB transcription factors which interact with GA-response elements [4]. In Arabidopsis seeds, miR159 regulates the abundance of MYB33 and MYB101 mRNAs during germination in response to ABA. These two MYB transcription factors were positive regulators of ABA responses and overexpressed miR159 would render plants hyposensitive to ABA [41,42]. In addition, the two target genes of Nnu-miR159 also participate in aleurone vacuolation which was in GA-induced pathways during seed germination. In this study, Nnu-miR159 targeted three GAMYB genes and their expression were highly up-regulated at 36 h during seed germination (Additional file 2: Table S3).
The Nnu-miR156 family members targeted 7 Squamosa Promoter Binding proteins like (SPL) plant specific transcription factors in the present study (Additional file 1: Figure S4 and Additional file 2: Table S5). The expression levels of Nnu-miR156a and its target gene SPL16 (NNU_02036-RA) were also validated ( Fig. 6g and j). In plants, SPLs were reported to be involved in diverse developmental processes including leaf development, phase change and flowering in plants [43]. Overexpression of miR156 in Arabidopsis and rice suppressed the expression of SPL genes and reduced apical dominance, delayed flowering time and affected the juvenile to adult transition [43][44][45]. MiR156 was also reported to be involved in seed  [46,47]. However, unlike the maize [43], the expression level of miR156 in sacred lotus was down-regulated after 72 h of the seed germination (Additional file 2: Table S2). That may be a factor for the long living of sacred lotus seed and viable for nearly thousands years.
Nuclear factor Y (NF-Y) was documented to affect the seed development and was targeted by miR169 in many plants. In the degradome dataset in the present study, NFYA2 (NNU_25577-RA), NFYA3 (NNU_11871-RA) and NFYA10 (NNU_09768-RA) were the target genes of miR169 in sacred lotus (Additional file 1: Figure S4 and Additional file 2: Table S5). In Arabidopsis, overexpression of the target of miR169 (NF-YA5) resulted in ABA hypersensitivity during seed germination [48]. Thus, upregulation of most of the miR169 genes in sacred lotus might reduce the ABA responsiveness during seed germination (Additional file 2: Table S2).
In plants, seed germination is promoted by appropriate environmental conditions including light, temperature, and nutrient availability [49]. Reactive oxygen species (ROS) are largetly generated during the seed germination and reducing the ROS and minimizing cell damage might affect the seed germination. MiR408 which targets L-ascorbate oxidase (NNU_10795-RA) might scavenge the oxidative species produced during stress. In the present study, all the miRNAs of Nnu-miR408 family were significant differentially expressed, indicating their roles in seed germination (Table 3 and Additional file 2: Table S3). PPR proteins that regulate gene expression in mitochondria and chloroplasts were also regulated by some miRNAs, including Nnu-miR168, Nun-miR2673, novel_mir_7 and novel_mir_40. These miRNAs might affect the pathways of photosynthesis and lipid accumulation during the seed germination [50].
Taken together, we proposed potential roles for miR-NAs playing during the early stages of seed germination in the sacred lotus (Fig. 8). These miRNA mediated the regulation of gene expression in many processes, including morphological changes, developmental process, metabolism and responsive to stress process (Fig. 8). These results indicated that diverse and complex miRNAs were involved in the seed germination. Further functional research and target analyses of the conserved and novel miRNAs could provide additional clues to the different regulation of gene expression during seed germination. The results in this study will lay a foundation for further investigating of long living of sacred lotus.

Conclusions
In this study, combined with small RNA sequencing and degradome sequencing, 145 known miRNAs and 78 novel miRNAs with 2580 targets were identified during seed germination of sacred lotus. GO and KEGG pathway analyses showed that many target genes enriched in regulation and metabolism. Compared with other angiosperm species, some miRNA families were present or absent in monocots and core eudicots, indicating the functional divergence and evolution in angiosperm. These results suggested that many miRNAs were involved in the regulation of seed germination of sacred lotus and provided insights into the evolutionary gains and losses of miRNAs in plants.

Plant materials
Sacred lotus (N. nucifera Gaertn.) seeds used here are maintained in our laboratory. The mature seeds of 3 years old were germinated in distilled water under 16-h light photoperiod at 28°C in a phytotron. The seeds at 0, 12, 24, 36 and 72 h after imbibition was collected and immediately frozen in liquid nitrogen, then stored in -80°C until further use.

RNA isolation and small RNA library construction
Total RNA was extracted from germinating seeds at 0, 12, 24, 36 and 72 h after imbibition using TRIZOL reagent (TianGen, China) according to the manufacturer's instructions. At least five germinating seeds were mixed for each sample. RNA-free Dnase I (Promega) was then used to remove DNA contamination from total RNA for 30 min at 37°C. RNA quantity was determined by a Qubit Fluometer. The RNA purity and integrity was detected by Agilent 2100 Bioanalyzer.
Small RNAs of 18~30 nt were purified from the 5 μg of total RNA by 15 % polyacrylamide gel electrophoresis (PAGE) as described previously [51]. A pair of Solexa adaptors was ligated to the 5′ and 3′ termini of the small RNAs. And then reverse transcription was performed by Super-Script II Reverse Transcriptase (Invitrogen) and amplified using 15 PCR cycles to produce sequencing library. The final quality of the cDNA libraries were examined with the Agilent 2100 Bioanalyzer. The small RNA fragments were finally sequenced by Hiseq2, 500 platforms.

Bioinformatics analysis of the miRNAs
The raw reads from the Illumina sequencing were first filtered to remove low-quality reads which including reads in length <18 nt, 5′ adaptor contaminates, reads with polyA, and reads without 3′ adaptor. And then adaptor sequences were trimmed to get clean reads. All unique sequences were aligned against the sacred lotus genome (http://lotus-db.wbgcas.cn/) to map the sequences using SOAP [52,53].
Subsequently, the clean sequences were used to search the RepeatMasker and Rfam database to annotate repeats, rRNA, tRNA, scRNA, snoRNA and snRNA. After removing the sequences belonging to rRNA, tRNA, scRNA, snoRNA, snRNA and repeat sequence tags, the remaining sequences were used to BLAST against miRbase 20.0 (http://www.mirbase.org) to identify conserved miRNA [54]. The remaining unannotated small RNA tags were used to predict novel miRNAs using the miRNA prediction algorithm MIREAP (http://sourceforge.net/projects/mireap). This program was utilized to identify novel miRNA candidates by exploring the secondary structures, Dicer cleavage sites and the minimum free energy (MFE) of the unannotated small RNA tags that could be mapped to the genome. The parameters for identifying potential miRNAs or pre-miRNAs as described previously [5,33]. And the secondary structures of potential precursor miR-NAs were constructed by MFOLD 3.2 [55].
To identify the expression patterns of miRNAs among the five different stages, the read count of each identified miRNA was normalized to reads per million (RPM) using the following formula: normalized expression = count of miRNA/total count of clean reads × 10 -6 . If one miRNA have no reads in a sample, the expression was arbitrarily set to 0.01 for further calculation [56]. Fold change of log2 scale value and P value were calculated using the formula as described [57]. The expression of |fold change| >2 and P value <0.01 was as significant differential expression. All the significantly differential expression profiles were clustered into different groups by Cluster 3.0 with default parameters.

Degradome library construction and data analysis
The degradome library was constructed as described previously [13]. Poly (A) RNA was extracted from total RNA of mixture of the five samples (1 μg of total RNA for each) using the Oligotex kit (Qiagen). With 5′monophosphates, polyadenylated transcripts were ligated to RNA adapters which consisted of a Mme I recognition site at its 3′ end. Then, first-strand cDNA was reverse-transcribed using oligo d (T) and amplified for five PCR cycles. The PCR product was purified and digested with Mme I and was then ligated to a 3′ double DNA adapter. After amplifying 18 PCR cycles, the PCR product was gel-purified for deep sequencing. The sliced miRNA targets were identified according to the Cleave-Land pipeline [13]. The potential targets were as followed: alignments with scores less than four (G: U pairs scored 0.5), and no mismatches were found at the site between the 10th and 11th nucleotides of the miRNAs.
The target genes were annotated using the information from Lotus-DB (http://lotus-db.wbgcas.cn/) [47]. The functional classification and pathways based on the three database: KOG (Clusters of Orthologous Groups of proteins), KO (KEGG Ortholog database) and GO (Gene Ontology), using BLAST with a cutoff E-value of 10 −5 . GO enrichment analysis was implemented by the TopGO R package and GO terms with corrected FDR < 0.05 were considered significantly enrichment [58]. The miRNA editing sites located on the known miRNA sequences were identified according the method as described previously [35]. The RNA editing level was calculated by the ratio of the reads support the mismatch in the site to the total reads detected on this site. The common and specific miRNAs with RNA editing during the different time points were also analyzed. Genomic DNA and total RNA were isolated from germinating seeds of sacred lotus to validate miRNA editing. Precursor miRNA sequences were amplified from DNA, and the corresponding mature miRNA sequences were also amplified from cDNA which was reverse-transcribed using stem-loop RT-PCR. Then the PCR products were purified, cloned and sequenced with a minimum of six clones for DNA products as well as at least twelve clones for miRNA PCR products. All primers used in this study are listed in Additional file 2: Table S10.
Detection of potential miRNA using stem-loop RT-PCR and qRT-PCR Total RNA was extracted from sacred lotus germinating seeds at different stages. Then, RNase-free DNase I (Promega) was used to remove the DNA contamination by incubating at 37°C for 30 min. For each miRNA, approximately 2 μg of total RNA was reverse-transcribed using miRNA-specific stem-loop primers in a 20 μl of reaction volume using a Fermentas RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA) [59]. The reactions were incubated for 30 min at 16°C, followed by pulsed RT of 60 cycles at 30°C for 30s, 42°C for 30s and 50°C for 1 s and finally the reactions were terminated at 70°C for 5 min [59]. The cDNA for the miRNA target was generated using 2 μg of total RNA and Oli-godT18 primer with the RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA). The novel miRNAs were cloned using specific primers and the PCR products were detected by gel electrophoresis. All the primers were listed in Additional file 2: Table S10.
Real time qRT-PCR analysis of the miRNA and their targets was performed using the FastStart Universal SYBR Green Master Mix (Roche) on the Ste-pOne plus PCR platform (Applied Biosystems). The qRT-PCR reactions were conducted with the following protocol: 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 56°C for 30s and 72°C for 15 s. Scared lotus NnEF1a (GI: 226897264) was used as an endogenous control. To avoid non-specific amplification, melting curve was carried out for each PCR product. The expression level of the miRNAs and their targets in different samples were calculated by comparative 2 -△△ CT method [60].

Availability of supporting data
All the sequencing data were deposited in the NCBI Short Read Archive (SRA) database under the accession number SRP070743.

Additional files
Additional file 1: Figure S1. Size distribution of small RNAs in all of the five libraries from germinating seeds at 0 h, 12 h, 24 h, 36 h and 72 h. Figure S2. Venn chart for total sRNA (A) and unique sRNAs (B) between the adjacent pairwise libraries. Figure S3. The secondary structures of novel miRNA precursors in sacred lotus. The mature miRNAs are in red and miRNA* in blue. ("." represent baase mismatches, "("represent base matches). Figure S4. Target plats (T-plots) of identified known and novel miRNA targets using degradome sequencing. The red lines indicate signatures consistent with miRNA-directed cleavage. Figure S5. GO classification of target transcripts for all identified miRNAs in sacred lotus. Figure S6. KOG function classification of the miRNA targets in sacred lotus. (PDF 3161 kb) Additional file 2: Table S1. Detailed information of the known miRNAs identified in sacred lotus during seed germination. Table S2. Detailed information of the novel miRNAs identified from sacred lotus during seed germination.  Table S5. The targets of known miRNA identified in sacred lotus during the seed germination. Table S6. The targets of novel miRNA predicted in sacred lotus during the seed germination. Table S7. GO enrichment analysis for all the target genes. Table S8. KEGG pathway enrichment analysis for the targets of differentially expressed miRNAs. Table S9. The common and specific miRNA editing in sacred lotus during the seed germination. Table S10. Primers of miRNAs and targets in sacred lotus for qRT-PCR and miRNA editing in this study. (XLSX 203 kb) Funding


extfieldsofstudy: Biology Medicine

text:
Successive requirement of Glass and Hazy for photoreceptor specification and maintenance in Drosophila

ABSTRACT Development of the insect compound eye requires a highly controlled interplay between transcription factors. However, the genetic mechanisms that link early eye field specification to photoreceptor terminal differentiation and fate maintenance remain largely unknown. Here, we decipher the function of 2 transcription factors, Glass and Hazy, which play a central role during photoreceptor development. The regulatory interactions between Glass and Hazy suggest that they function together in a coherent feed-forward loop in all types of Drosophila photoreceptors. While the glass mutant eye lacks the expression of virtually all photoreceptor genes, young hazy mutants correctly express most phototransduction genes. Interestingly, the expression of these genes is drastically reduced in old hazy mutants. This age-dependent loss of the phototransduction cascade correlates with a loss of phototaxis in old hazy mutant flies. We conclude that Glass can either directly or indirectly initiate the expression of most phototransduction proteins in a Hazy-independent manner, and that Hazy is mainly required for the maintenance of functional photoreceptors in adult flies.

Introduction
Cell differentiation is typically controlled by networks of transcription factors, which gradually shift during development from specifying organ and cellular identity to activating terminal gene expression in mature cells. Such a network acts during the formation of photoreceptors (PRs) in Drosophila. [1][2][3] During early eye development, an evolutionarily conserved set of transcription factors, called the retinal determination network (RDN), specifies eye field identity in the eye imaginal disc. RDN genes are both necessary and sufficient for eye formation, and thus can induce formation of ectopic eyes when misexpressed in other imaginal discs. [4][5][6][7][8] Recently, we have shown that PR differentiation is critically regulated by the zinc finger transcription factor Glass. Glass provides a genetic link between the early-acting RDN, terminal differentiation transcription factors, and genes functioning in mature PRs, such as those involved in the phototransduction cascade (Bernardo-Garcia et al. 9 ). In the absence of Glass, PR precursors retain a neuronal identity but fail to express any PR markers. Therefore, PR precursors require Glass for differentiating into functional lightsensing cells.
A direct target of Glass is the homeobox transcription factor Hazy. Ectopic expression of Glass is sufficient to induce expression of Hazy and some phototransduction proteins. Similarly, ectopic expression of Hazy can only induce a subset of phototransduction components. However, when co-expressed, Glass and Hazy can ectopically induce most of the phototransduction cascade, suggesting that both Glass and Hazy act synergistically during PR development. 9 Thus, the combinatorial action of these 2 transcription factors appears to play a central role in directing PR precursors toward a terminal differentiation program.
Here, we investigate the regulatory interaction between Glass and Hazy. We found that, while Glass is able to activate its own promoter, Hazy does not appear to auto-activate its own expression nor that of Glass. Also, by analyzing hazy mutants we disentangle the individual roles of Glass and Hazy in regulating the expression of phototransduction genes. Interestingly, we found that Hazy is particularly relevant to ensure the continued expression of phototransduction proteins in adult PRs. Young hazy mutants correctly express most of the phototransduction components, and show a similar attraction to white light as wildtype flies. By contrast, the expression of most phototransduction proteins is reduced in old hazy mutants, and they fail to show phototaxis. Together with previous results, our data suggest that Glass and Hazy are required for different tasks and at different steps in PR development. During early eye development, Glass contributes to specifying PR identity. [9][10][11][12] Later, during terminal differentiation, both Glass and Hazy activate genes that are required for the maturation of functional PRs. 9,12,13 Finally, Hazy is required for maintaining expression of phototransduction genes, and thus ensures the continued functionality of adult PRs.

Results and discussion
Hazy is a direct target of Glass in all visual organs in the fly PR development in the eye imaginal disc starts after the passage of the so-called morphogenetic furrow (MF), which sweeps across the disc from the posterior edge toward anterior, initiating the formation of ommatidia. RDN genes are present in the eye disc prior to the MF, while the proneural gene atonal is transiently expressed at the MF. 1,2 Glass expression is initiated at the MF, and maintained in differentiating cells in the retina. 14,15 In contrast, the expression of Hazy starts later during pupation, after all PRs have been specified. 3,13 We and others have shown that the expression of a hazy(wt)-GFP reporter in compound eye PRs depends on Glass binding to 2 sites in the hazy promoter. 9,12 In addition to the compound eye, flies also have PRs in the ocelli (3 separate eyes located at the top of the adult head) and in the larval eye, also termed Bolwig's organ. Since both Glass and Hazy also play a role in the development of these organs, 9,11,14,[16][17][18][19] we hypothesized that Glass might similarly activate hazy in PRs outside the compound eye.
The hazy(wt)-GFP reporter, which we have previously used to study the expression of hazy in the compound eye, was also expressed in PRs of the ocelli and Bolwig's organ (Figs. 1A, A 0 , D-D 000 ), reflecting the expression pattern of the endogenous Hazy protein. 13,17 When the hazy(wt)-GFP transgene was placed in glass mutant background, GFP expression was lost in the ocelli and the Bolwig's organ (Figs. 1B, E-E 000 ). The wild-type embryonic Bolwig's organ consists of 12 PRs, and can be identified in stage 14 embryos because of the co-expression of Kruppel and Fasciclin 2 (Fas2). 17,20 Interestingly, in glass mutants we only found 4 Kruppel-positive cells in the Bolwig's organ (Fig. 1E 000 ), indicating an early defect in larval eye formation. In addition, a hazy(gl1,2mut)-GFP reporter in which the Glass binding sites were mutated was not expressed in the ocelli nor in the Bolwig's organ (Figs. 1C, F-F 000 ).
Together, these results corroborate our findings for the compound eye. 9 We conclude that the hazy promoter is directly bound and activated by Glass in all PRs in the 3 visual organs of Drosophila.

Glass can auto-activate its own expression
We have previously shown that ectopically expressing Glass or Hazy induces the expression of some phototransduction proteins in the central nervous system (CNS). Co-expressing Glass and Hazy displays a synergistic effect on the induction of phototransduction components. Not only the genes that are activated by either Glass or Hazy alone become ectopically expressed, but also additional phototransduction proteins are induced, suggesting that Glass and Hazy function together in a coherent feed-forward loop. Glass activates the expression of Hazy and together they are able to activate the expression of more target genes than either Glass or Hazy could activate on their own. 9 Here we tested additional regulatory interactions between Glass and Hazy.
For this, we used a glass-DsRed and a hazy(wt)-GFP reporter. 9,21 glass-DsRed larvae expressed nuclear DsRed in the Bolwig's organ and in Glass-expressing cells in the brain (Figs. 2A-A 000 ), while hazy(wt)-GFP was expressed exclusively in the Bolwig's organ PRs (Figs. 1D, 2B), whose axons could be seen projecting into the brain (Fig. 2B 0 ). Thus, both reporters faithfully reflect the expression patterns of Glass and Hazy. 13,14 We performed flip-out experiments in which we ectopically induced either Glass or Hazy in the CNS of third instar larvae, and found that Glass was able to activate the glass-DsRed reporter in a subset of cells (Figs. 2C-C 000 ). Therefore, Glass may be able to maintain its own expression by auto-regulation. In contrast, ectopic expression of Hazy did not induce the expression of the hazy(wt)-GFP reporter (Figs. 2D, D 0 ), suggesting that Hazy cannot activate its own expression. Similarly, we did not observe ectopic expression of Glass in the CNS after ectopically expressing Hazy (Figs. 2E, E 0 ), suggesting that Glass expression is not activated by Hazy.

Glass can initiate the expression of most phototransduction proteins independently of Hazy
Both Glass and Hazy are required for the expression of PR genes. 9,13 Since we have shown that Glass directly activates hazy, and that inducing the expression of Hazy partly rescues the glass mutant phenotype, 9 it would be possible that Glass mainly relies on Hazy for activating the expression of phototransduction proteins. To test this, we examined the individual role of Hazy.
We found that young hazy mutant flies-less than one day old-failed to express Rhodopsin 6 (Rh6) and Transient receptor potential-like (Trpl), but did express correctly most of the phototransduction proteins that we tested: Rhodopsin 1 (Rh1), G protein a q subunit (Gaq), No receptor potential A (NorpA), Transient receptor potential (Trp), Inactivation no afterpotential D (InaD) and Arrestin 1 (Arr1) (Figs. 3A-P). [22][23][24][25][26][27][28][29][30] Of these, after 10 d the levels of Rh1, Gaq, NorpA, Trp and InaD were strongly reduced, and only Arr1 expression appeared unchanged (Figs. 3a-p). Thus, our results indicate that Hazy is . Samples were stained with antibodies against GFP (green), and Elav (used as a neuronal marker, magenta). The hazy(wt)-GFP reporter was expressed in the ocelli in wild-type (A, A'), but not glass mutant background (B). A hazy(gl1,2mut)-GFP reporter in which the 2 Glass binding sites were mutated was not expressed in the ocelli (C). (D-F) In the case of the Bolwig's organ, this is a larval eye that develops from the optic placode during embryogenesis. Stage 14 embryos were stained with antibodies against GFP (green), Fas2 (red) and Kruppel (blue). At this stage, the developing Bolwig's organ is located dorsally, still in contact with the surface, and can be identified both because of its position and the co-expression of Fas2 and Kruppel., 17,20,36 Similar to the ocelli, the hazy(wt)-GFP reporter was expressed in the Bolwig's organ in wild-type (D), but not glass mutant background (E). Also, hazy(gl1,2mut)-GFP was not expressed in the Bolwig's organ (F). For each image, the 3 channels from a close-up of the Bolwig's organ were separated and are shown below in grayscale (D 0 -F 000 ). Scale bars represent 20 mm in D 0 -F 000 ; 30 mm in A 0 , B-F; and 100 mm in A.
not required for initiating the expression of most phototransduction proteins, but that it plays an important role in maintaining the differentiated state of PRs. These results contrast with an earlier description of the hazy mutant phenotype. 13 The differences between our findings and this report may be explained because we analyzed young and old flies separately.
Since most phototransduction proteins are expressed in the retina of young hazy mutants, we infer that Glass does not mainly act via Hazy for initiating the expression of most phototransduction proteins. In the case of Rh6, Trpl, Rh1, Gaq, NorpA, Trp, InaD and Arr1, we found that all these genes contain putative Glass binding sites in their regulatory sequences (the GAARCC motif, which is present in either their promoter or their introns). 31 Therefore, it would be possible that Glass either directly or indirectly activates the expression of these phototransduction components.

Hazy is not required for white light detection in young flies
To further assess the impact of their age-dependent loss of phototransduction genes we analyzed the phototactic behavior of young versus old hazy mutant flies. It has been previously reported that hazy mutants fail to detect light due to the absence of many phototransduction proteins. 13,16 However, our finding suggests that young hazy mutant flies express a set of genes sufficient for the phototransduction machinery to detect white light.
Adult wild-type flies display a positive phototactic behavior. In a 2-choice assay, they move toward the light. This preference for light is very high in newly eclosed flies and shows some reduction when the flies get older. 32 10 day old flies are still able to distinguish between light and darkness (Fig. 4). In contrast, glass mutant flies do not show phototaxis from the day they eclose (Fig. 4). 33 This is in agreement with our previous finding that glass mutants do not express any of the proteins in the phototransduction cascade. 9 In the case of hazy mutants, we observed normal phototactic behavior in young flies, comparable with that of wildtype. However, at the age of 10 d they differed from wild-type and displayed the same disability to distinguish between light and darkness as glass mutants (Fig. 4). These experiments are in agreement with our antibody analysis above, and show that young hazy mutants are able to detect white light, but lose this ability over time. Thus, Hazy is required for the maintenance of PR functionality.  Hazy. (A, B) We used Hoechst (blue), which labels cell nuclei, as a counterstain to analyze the expression pattern of the glass-DsRed and hazy(wt)-GFP reporters in the CNS of third instar larvae. The glass-DsRed reporter was expressed the nuclei of some cells in the brain (green, A). A close-up to the right shows that those neurons endogenously expressing Glass (red/magenta) also co-express the reporter (green, A 0 ). These two channels are shown separately to the right in grayscale (A, 00 A 000 ). The hazy(wt)-GFP reporter is exclusively expressed in PRs (green, B). A grayscale image to the right shows GFP labeling the axonal projections of the PRs in the brain (arrows, B 0 ). (C-E) In flip-out experiments we ectopically induced either Glass or Hazy expression in clones labeled with nuclear b-galactosidase (bGal). We stained the CNS of third instar larvae with antibodies against bGal (red/magenta); either DsRed, GFP or Glass (green) and with Hoechst (blue). We found that Glass ectopically induced the expression of the glass-DsRed reporter in the ventral nerve cord (C, C 0 ; channels are also shown separately in grayscale in C, 00 C 000 ). By contrast, Hazy did not ectopically induce the hazy(wt)-GFP reporter (D, D 0 ) nor Glass (E, E 0 ). Scale bars represent 20 mm in A 0 , C 0 -E 0 ; and 80 mm in A-E.

Conclusions
Glass and Hazy play important roles in PR specification and maintenance. 9,11,13,34 Here we have shown that Glass activates the expression of hazy in all PRs in the 3 visual organs of Drosophila-the compound eye, the ocelli and the Bolwig's organ-by binding to the same sites in the hazy promoter. Also, in agreement with a previous report, Glass is able to auto-activate its own expression, 14 but we found no evidence that Hazy either activates glass or auto-activates its own expression. Together, our data favor a feed-forward mechanism in which Hazy acts downstream of Glass. 9,35 The notion that Glass and Hazy function through a coherent feed-forward loop to activate PR genes is further supported by our previous findings that inducing the expression of Hazy partly rescues the glass mutant phenotype, and also that co-expressing Glass and Hazy together ectopically induces more phototransduction proteins than either Glass or Hazy alone. 9 However, these findings were based on ectopic expression of both transcription factors in CNS cells in which they are normally not expressed and where we were not able to control their expression levels. Our analysis of the hazy mutant phenotype rather suggests that, by the end of terminal PR differentiation, Glass can induce the expression of most phototransduction components even in the absence of Hazy. Hazy itself is one of the targets of Glass and is required to maintain the expression of most phototransduction components throughout adult life. Thus, for the initial induction, Glass either directly activates PR-specific genes, or it interacts with other transcription factors to induce them. Hazy might be required for the initial induction of only a few genes, such as Rh6 or Trpl, but it is essential for the maintenance of most phototransduction components. Therefore, we anticipate that further research on the direct targets of Glass will reveal novel mechanisms for activating the expression of PR genes. Also, it will be interesting to explore further how Hazy functions in PR fate maintenance. Particularly, it remains unresolved whether Hazy is only required for maintaining PR gene expression or if it also functions to suppress PR dedifferentiation or degeneration. Whisker lines indicate the maximum and minimum data point that are closer than 1.5 interquartile range of its nearest quartile. Circles indicate outliers. We used Welch's t-test for comparing the PIs between groups (n D 7 per age and genotype) and to zero. Significance levels represent p > 0.05 (not significant, n.s.), p 0.05 ( Ã ), p 0.01 ( ÃÃ ), and p 0.001 ( ÃÃÃ ). In a 2-choice assay, groups of wild-type flies of every age showed positive phototaxis, which decreases with age (indicated by positive PI values, which were significantly different from zero). glass mutants were photoneutral at all ages (their PIs were not significantly different from zero). Newly eclosed hazy mutants showed positive phototaxis, not different from that of wildtype flies (p D 0.67, median wild-type PI D 0.83). Five day old hazy mutants and wild-type flies show a decreased positive phototaxis, but their PIs are not different from each other (p D 0.30, median wild-type PI D 0.42). Ten day old hazy mutants were photoneutral, with their PIs comparable to zero (p D 0.08) or to glass mutants (p D 0.56), and significantly different from wild-type (median wild-type PI D 0.22).

Fly stocks
In the present work w 1118 (courtesy of R. Stocker) was used as wild-type control, and w 1118 ; hazy hazy was used to study the hazy mutant phenotype. The hazy hazy mutant allele was provided by C. Desplan, 13 and was isogenized by crossing it to Canton-S flies for 6 generations. Other stocks that we used are: hazy(wt)-GFP, 9 hazy(gl1,2mut)-GFP, 9 gl 60j (Bloomington Drosophila Stock Center, No. 509), 11 UASglass 9 and UAS-hazy (courtesy of J. Bischof). 37 As a reporter for glass we used glass-DsRed flies (courtesy of S. Kim), which are also called glass5.2-RHS. 21 Flipout misexpression experiments were performed as described previously 9 by using hsFLP; tub(FRTcassette)Gal4, UAS-lacZ.nls (courtesy of E. Piddini).
Flies were raised at 25 C in a 12:12 hours light-dark cycle on cornmeal medium supplemented with molasses, fructose and yeast.

Antibody stainings
Immunohistochemistry was performed as previously described. 9,17,38 Antibodies against proteins in the phototransduction cascade were kindly provided by N. Colley and S. Britt. To compare the expression of phototransduction proteins in control and hazy mutants, head sections from both genotypes were taken simultaneously, stained together on the same slide, and imaged with identical settings on a Leica SP5 confocal microscope.

Phototaxis assay
Our phototaxis analysis was based on a previous protocol. 32 Briefly, newly eclosed flies were transferred each day to vials containing fresh food, which we used to stage them. For each experiment, we tested an average of 25.2 flies (ranging from 20 to 33). These flies were kept for 10 min in darkness, and then placed without anesthesia into a T-maze under red light conditions. Our set-up consisted of 2 tubes connected to each other, where a single LED illuminated from the end of one of the tubes with white light (SOLAROX, Germany, No. 50008300001). Light intensity was moderate: we measured 418.0739 mW/cm 2 with a spectrometer (the emitted light spectrum possessed 2 intensity peaks: the first peak was at 444 nm with an intensity of 1.494 mW/cm 2 /nm and half-widths of 16 nm, and the second peak was at 585 nm with an intensity of 2.768 mW/cm 2 /nm and half-widths of 61.5 nm). We allowed the flies to move freely between both tubes for 2 minutes. Then, we counted the flies in the illuminated tube (L), in the dark tube (D), and those in the intersection between the 2 tubes (M). The preference index (PI) was calculated from the formula PI D (L ¡ D)/(D C L C M). Data were analyzed in R with Welch's t-test.

Disclosure of potential conflicts of interest
No potential conflicts of interest were disclosed.


extfieldsofstudy: Chemistry Medicine

text:
Deficiency of ATP6V1H Causes Bone Loss by Inhibiting Bone Resorption and Bone Formation through the TGF-β1 Pathway

Vacuolar-type H +-ATPase (V-ATPase) is a highly conserved, ancient enzyme that couples the energy of ATP hydrolysis to proton transport across vesicular and plasma membranes of eukaryotic cells. Previously reported mutations of various V-ATPase subunits are associated with increased bone density. We now show that haploinsufficiency for the H subunit of the V1 domain (ATP6V1H) is associated with osteoporosis in humans and mice. A genome-wide SNP array analysis of 1625 Han Chinese found that 4 of 15 tag SNPs (26.7%) within ATP6V1H were significantly associated with low spine bone mineral density. Atp6v1h+/- knockout mice generated by the CRISPR/Cas9 technique had decreased bone remodeling and a net bone matrix loss. Atp6v1h+/- osteoclasts showed impaired bone formation and increased bone resorption. The increased intracellular pH of Atp6v1h+/- osteoclasts downregulated TGF-β1 activation, thereby reducing induction of osteoblast formation but the bone mineralization was not altered. However, bone formation was reduced more than bone resorption. Our data provide evidence that partial loss of ATP6V1H function results in osteoporosis/osteopenia. We propose that defective osteoclast formation triggers impaired bone formation by altering bone remodeling. In the future, ATP6V1H might, therefore, serve as a target for the therapy of osteoporosis.

Introduction
Vacuolar-type H + -ATPase (V-ATPase), a highly conserved enzyme, mediates acidification in organelles by pumping protons and coupling with the energy of ATP hydrolysis. V-ATPase-dependent organelle acidification is necessary for intracellular processes such as protein sorting, zymogen activation, receptor-mediated endocytosis, and synaptic vesicle proton gradient generation. V-ATPase consists of a peripheral catalytic V1 domain and a membrane proton channel V 0 domain [1]. Subunits of V-ATPase play important roles in the biological and physiological functions of osteoclasts [2] and the dysfunction of these subunits is usually associated with osteopetrosis and other bone abnormalities.

Ivyspring
International Publisher T-cell immune regulator 1 (TCIRG1), also known as V-ATPase V 0 subunit a3, resides in the ruffled border of osteoclasts and controls osteoclast-mediated extracellular acidification. Mutations in TCIRG1 gene cause autosomal recessive osteopetrosis in humans [3][4][5][6], and mice deficient in Tcirg1 (or Atp6i) have hypocalcemia and osteopetrorickets [7]. The V0D2 subunit of V-ATPase is required for efficient preosteoclast fusion. Atp6v0d2-deficient mice have increased bone mass and enhanced bone formation [8][9][10]. The C subunit of the V 1 domain of V-ATPase regulates F-actin ring formation during osteoclast activation [11][12][13] and the reversible dissociation of V-ATPases [14]. Consequently, Atp6v1c1 deficiency severely impairs osteoclast acidification and bone resorption, whereas osteoclast differentiation is unaffected [11]. The V1B2 subunit of V-ATPase is associated with F-actin and facilitates the recruitment of V-ATPase complexes to the osteoclast ruffled border during polarization and bone resorption [15,16]. Recent findings have shown that mutations in the ATP6V1B1 gene cause Zimmermann-Laband syndrome and dominant deafness-onychodystrophy syndrome, which might be related to impaired assembly of the V1 subcomplex of ATP6V 1B1 [17,18].
Osteoporosis is a common metabolic bone disease characterized by reduced bone mineral density (BMD) and increased risk of osteoporotic fractures. Genetic factors have been closely associated with the risk of osteoporosis inciting research on the genes involved in osteoclasts functions in osteoporosis [19,20]. One of our recent bivariate genome-wide association studies (GWAS) implicated ATP6V1G1 as a novel pleiotropic gene affecting BMD in osteoporosis patients [21]. This study suggested that other subunits of V-ATPase could also be associated with osteoclast function and bone disorders. ATP6V1H is a small subunit of V-ATPase that connects the V1 and V 0 domains; its role in bone development and bone disorders remains unknown. Further GWAS screening of the osteoporosis population helped us identify ATP6V1H as an additional subunit of V-ATPase involved in osteoporosis. To define the function of ATP6V1H and its relationship with osteoporosis, we generated Atp6v1h knockout mice using the CRISPR/Cas9 technique and defined the bone phenotypes in this mouse model. We observed that the deficiency of ATP6V1H caused a lower bone turn over and bone loss by inhibiting osteoclasts formation and bone formation at the same time; bone formation was reduced more than bone resorption, resulting in a net bone loss. This imbalance in bone homeostasis occurred from altered interactions between osteoclasts and osteoblasts through the TGF-β1 pathway. 1625 Han Chinese were divided into three groups with different Z-scores (Table 1). We chose -1SD to -2SD of Z-score as a threshold to separate controls from osteopenia group with low bone mass and osteoporosis group. A Z-score below -2 was considered as a significant risk for osteoporosis. A GWAS analysis involving 281,533 SNPs [22,23] was performed and fifteen SNPs targeted to the intronic regions of the ATP6V1H gene were selected. These variants were compared among the three groups. Multiple testing identified 4 of them (26.7%) having a nominal significant association with spine BMD, P<0.05 ( Table 2).

ATP6V1H variants and low bone mass in a human population.
Bone loss and decreased bone formation in Atp6v1h +/mice. We used the CRISPR/Cas9 technique, designed guide RNA (gRNA) to target exon 2 of the Atp6v1h gene, and created a 5bp (CGAGG) deletion and one base replacement (T>G). The mutations caused an early stop codon (Fig.1A). The gRNA has less chance of off-target events according to the results from MIT CRISPR design and analysis website (http://crispr.mit.edu/). We tested 7 genes (Aldh16a1, Tbrg1, Wdr4, Atp2b3, Gcsh, Zscan, and Kcnt11) with 4%~5% predicted occurrence ratio of off-target events, but none of them showed off-target effects. All homozygous Atp6v1h -/embryos died before embryonic day 10.
Atp6v1h +/heterozygous mice were born at the expected Mendelian ratio and exhibited a normal growth rate. There were no obvious differences in weight or size between wild-type and heterozygous mice. IHC staining and tartrate-resistant acidic phosphatase (TRAP) staining indicated that ATP6V1H was localized in osteoclasts. The positive staining of ATP6V1H in osteoclasts became weak in heterozygous mice (Fig. 1B). TRAP staining for osteoclasts showed an equal number of osteoclasts and same size of osteoclasts in the distal femurs of Atp6v1h +/mice compared with control littermates (Fig.1 C). High-resolution micro-CT analysis in calvaria (Fig. 1D), femurs ( Fig.1E and Table 3) and lumbar vertebrae (Table 3), however, demonstrated obvious decreased bone mass in Atp6v1h +/mice. Ponceau red staining of calcified bone tissues further showed decreased bone volume and trabecular thickness ( Fig.1 F). Transmission electron microscope (TEM) analysis found that Atp6v1h +/osteoclasts had normal ruffled borders (Fig. 1G). The serum CTX-I and TRAP, which reflected osteoclastic resorptive activity, were not altered in Atp6v1h +/mice (Fig. 1H,  I). Serum calcium concentrations were slightly higher ( Fig. 1J), and serum phosphate levels were lower in Atp6v1h +/mice compared to the wild-type (Fig. 1K).
The changes of bone loss and decreased bone formation co-existed in Atp6v1h +/mice. This suggests that the osteoporotic changes of Atp6v1h +/mice might result from an imbalance of osteoclasts and osteoblasts activities.
Impaired osteoclast formation and function in Atp6v1h +/mice. To characterize the osteoclasts of Atp6v1h +/mice in vitro, we stimulated primary bone marrow cells (BMs) with receptor activator of nuclear factor kappa-B (NF-κ B) ligand (RANKL) and macrophage colony-stimulating factor (M-CSF). The number of TRAP-positive Atp6v1h +/osteoclasts was reduced and the size of Atp6v1h +/osteoclasts was not significantly changed when the cells were cultured in a general 96-well plate (Fig.3A). When BMs were cultured on an osteoassay microplate with the same RANKL/M-CSF stimulating conditions, we observed fewer and smaller osteoclasts in the Atp6v1h +/group, and the total resorptive area, as well as the osteoclast number normalized to the resorptive area, were also decreased (Fig. 3B). Consistent with the above results, the number of pits per bone slice (4mm×4mm) was decreased as was the average pit area of Atp6v1h +/osteoclasts compared with those of osteoclasts from wild type mice (Fig. 3C); this indicates that the Atp6v1h defect influenced bone resorptive activity. Comparing the formation of actin ring between wild-type and Atp6v1h +/osteoclasts detected no difference (Fig. 3D).    There was no significant difference between wild-type and Atp6v1h +/-mice. n=9. J and K. The serum level of calcium was higher and that of phosphate was lower in Atp6v1h +/-mice than in wild-type. n=8. All data represent mean ±s.d. ** P<0.01, *P<0.05. NS: not significant difference.
The protein level of ATP6V1H in wild-type osteoclasts gradually increased during days 0~4 of RANKL stimulation, whereas the level of ATP6V1H in Atp6v1h +/osteoclasts was lower than that of wild-type osteoclasts throughout (Fig.  3E). Knockdown of Atp6v1h by siRNA reduced the formation of TRAP+ osteoclasts obtained from BMs stimulated with M-CSF and RANKL (Fig. 3F).
These findings suggest that Atp6v1h deficiency affects the formation of osteoclasts and the resorptive ability of osteoclasts by impairing the V-ATPase activity and increasing the intracellular pH.
Decreased interaction between osteoblast and osteoclast in Atp6v1h +/mice. The lower rate of bone formation prompted us to study whether the Atp6v1h deficiency also affected the osteoblasts directly. The immunohistochemical staining did not detect ATP6V1H protein in osteoblasts from wild-type and Atp6v1h +/mice, although comparable mRNA levels of Atp6v1h could be detected in the cranial bone of new born mice and primary cultured osteoblasts of wild-type and Atp6v1h +/mice (Fig. 4A). The screening of the mRNA levels of osteogenic markers in cranial bones and osteoblasts detected lower Alp, Col1a1 and Bsp expression in osteogenic-induced osteoblasts (Fig.  4A). When the mRNA level of Atp6v1h was downregulated in primary cultured osteoblasts or the osteoblast MC3T3-E1 cell line by Atp6v1h siRNA, only the levels of Ocn and Opn were altered (Fig. 4B), whereas the Atp6v1h +/osteoblasts had the normal ability of calcified node formation (Fig. 4C). Thus we believe that a cell-autonomous effect of ATP6V1H haploinsufficiency on the osteoblast function is unlikely, and the osteoblastic effect is probably secondary to the reduced osteoclast formation and function.
Since TGF-β1 activation is pH sensitive [24,25] and TGF-β1 could regulate the functions of osteoclasts and osteoblasts [26,27], we hypothesized that the changes in intracellular pH in Atp6v1h +/osteoclasts alter the level and activity of TGF-β1. Indeed, the protein level of TGF-β1 in Atp6v1h +/osteoclasts was reduced more than half, whereas in Atp6v1h +/osteoblasts, it was comparable to that of wild-type cells (Fig. 4D). The active TGF-β1 level in sera was also reduced in Atp6v1h +/mice (Fig. 4E). Using ELISA to evaluate the secretion of active TGF-β1 from osteoclasts and osteoblasts, we observed that Atp6v1h +/osteoclasts cultured on bone slices secreted less TGF-β1 (Fig. 4F), and that the osteogenic-induced Atp6v1h +/osteoblasts produced more TGF-β1 (Fig.  4G). The lower levels of TGF-β1 have been reported to downregulate the differentiation of mesenchymal stem cells into osteoblasts. Thus the reduced number and size of osteoblasts in Atp6v1h +/mice might result from a lower rate of osteoblasts formation.   Alizarin red S and ALP staining showed that calvarial osteoblasts from wild-type and Atp6v1h +/-mice formed mineralized nodules but no difference existed between them. n≥3. D. Western results showed the decreased TGF-β1 in Atp6v1h +/-osteoclasts and unchanged TGF-β1 in Atp6v1h +/-osteoblasts, respectively. n=3. E. ELISA detection revealed the decreased active TGF-β1 in sera of Atp6v1h +/-mice. n≥8. F. The active TGF-β1 level in the supernatant of Atp6v1h +/-osteoclasts cultured on the bone slices was also decreased. Control: control medium. n≥3. G. Atp6v1h +/-osteoblasts with osteogenic induction secreted more active TGF-β1 by ELISA detection. No difference existed between wild-type and Atp6v1h +/-osteoblasts under a non-osteogenic induction condition. n≥3. H. Western blot results showed the decreased RANKL and unchanged OPG in Atp6v1h +/-osteoblasts. n≥3. I. ELISA results revealed the less soluble RANKL (sRANKL) in sera and the supernatant of Atp6v1h +/-osteoblasts with osteogenic induction. n≥3. J. OPG level was increased in sera of Atp6v1h +/-mice and the supernatant of Atp6v1h +/-osteoblasts with osteogenic induction. n≥3. K. Comparison of the number of osteoclasts in the co-culture system of osteoblasts and osteoclasts. When co-cultured with wild-type osteoblasts, wild-type and Atp6v1h +/-monocytes formed the same amount TRAP+ osteoclasts. When co-cultured with Atp6v1h +/-osteoblasts, fewer TRAP+ osteoclasts were formed in Atp6v1h +/-group. n≥3. All data represent mean ±s.d. ** P<0.01, *P<0.05. NS: not significant difference. RANKL/OPG pathway was also analyzed to determine whether a lower bone formation rate downregulated the osteoclasts function. The protein level of RANKL was decreased in Atp6v1h +/osteoblasts (Fig. 4H). The serum levels of RANKL and OPG were respectively decreased and increased in Atp6v1h +/mice (Fig. 4I, J). The secretion of RANKL and OPG from Atp6v1h +/osteoblasts following osteogenic induction was also respectively decreased and increased (Fig. 4I, J). The levels of RANKL and OPG from primary cultured osteoblasts were unchanged under non-osteogenic inducing conditions. In a co-culture system using different combinations of monocytes and osteoblasts from wild-type and heterozygous mice, co-culture of wild-type osteoblasts with wild-type and Atp6v1h +/-monocytes induced the same amount of TRAP-positive osteoclasts, whereas Atp6v1h +/osteoblasts induced fewer osteoclasts from Atp6v1h +/monocytes ( Fig.4 K).

Discussion
Several genes are known to influence bone mass, including the vitamin D receptor, estrogen receptor, and collagen type I alpha I gene [19,20]. In the present study, we report that defects of ATP6V1H are associated with bone matrix loss in mice and humans.
V-ATPase plays important roles in various cellular functions, such as acidification, endocytosis, and zymogen activation. The subunits of V-ATPase previously known to be critical for osteoclast function include those encoded by TCIRG1, ATP6V0D2, ATP6V1C1, and ATP6V1B2. Dysfunctions of these subunits manifest in osteopetrosis or osteopetrosis-related changes. The mechanisms involve changes in extracellular acidification, a fusion of preosteoclasts, and cytoskeletal F-actin assembly of osteoclasts [3,8,11,17,28,29]. In this context, we expected that defects in ATP6V1H would impair osteoclast function; indeed, we discovered that ATP6V1H was highly expressed in osteoclasts and the expression level of ATP6V1H was decreased in Atp6v1h +/osteoclasts. The reduced level of ATP6V1H in osteoclasts did not affect the formation of F-actin rings or ruffled borders of osteoclasts, but inhibited the number and size of TRAP+ osteoclasts, impaired V-ATPase activity, increased intracellular pH, and impaired osteoclastic bone resorption. Thus it appears that the key role of ATP6V1H is to promote osteoclast formation and maintain local acidification for bone resorption, rather than to affect osteoclast structure or movement.
We also anticipated that mutations of ATP6V1H homologs would impair osteoclast function and result in osteopetrosis, but instead detected that the in vivo effect of Atp6v1h +/mice was bone loss. Atp6v1h +/mice showed decreased bone density and serum phosphate level and increased serum calcium level. Bone formation and bone calcification were decreased with fewer and smaller osteoblasts and reduced ALP activity. One possibility to explain the bone loss in Atp6v1h +/mice involves an imbalance of bone resorption and bone formation, altering bone homeostasis through coordinated cellular and molecular events. When the decrease in bone formation exceeds the decrease in bone resorption, the net effect is osteoporosis or osteopenia. Indeed, the impaired bone remodeling and a negative bone balance, causing profound osteoporosis, have been reported in other gene knockout mice [30][31][32], emphasizing the complex interactions and the delicate balance required for bone remodeling.
Although Atp6v1h mRNA could be detected in cranial bones, in primarily cultured osteoblasts, and in the osteoblastic cell line MC3T3-E1, we did not detect differences in the expression of osteogenic genes between wild-type and heterozygous mice, or when Atp6v1h was knocked down by siRNA. Thus, the decreased bone formation is probably secondary to osteoclasts instead of a cell-autonomous effect of ATP6V1H haploinsufficiency in osteoblasts. Research on ATP6V0D2 also revealed the importance of interactions between osteoclasts and osteoblasts [8,9]. Thus, ATP6V1H is another subunit of V-ATPase that affects the coordination of osteoclast and osteoblastic functions as well as bone formation.
There are several pathways by which osteoclasts interact with osteoblasts, which arise from mesenchymal stem cells. TGF-β1 may stimulate the proliferation of mesenchymal stem cells and induce their differentiation into osteoblasts [33,34]; the low concentration of TGF-β1 promotes the maturation of osteoclasts [35]. These process couple bone resorption with bone formation [26,36,37]. Our data showed decreased levels of active TGF-β1 in the sera of Atp6v1h +/mice and the supernatant of Atp6v1h +/osteoclasts cultured with bone slices. The reduced level of TGF-β1 may explain the differences between impaired Atp6v1h +/osteoclasts in vitro and unexpected bone loss in Atp6v1h +/mice in vivo. Reduced levels of TGF-β1 negatively regulate osteoblast differentiation, resulting in the fewer and smaller osteoblasts and finally causing the decreased bone formation in Atp6v1h +/mice. Since the decrease in bone formation was greater than the decrease in osteoclast function, the final in vivo net effect was the bone loss (Fig. 5). pH is one of the important factors that can activate TGF-β [24,25]. The increased intracellular pH resulted in decreased secretion of TGF-β1 from Atp6v1h +/osteoclasts. We propose that TGF-β1 is the key factor that regulates the ATP6V1H-dependent osteoporosis changes.
The RANKL/OPG signaling pathway is important for the osteoclast differentiation. TGF-β1 level may trigger the expression of RANKL in osteoblasts. High concentrations of TGF-β1 upregulate the expression of OPG, which may inhibit RANKL-induced osteoclast differentiation [37]. In Atp6v1h +/mice, the amount of RANKL and OPG secreted by osteoblasts was not changed when the cells were cultured in a non-osteogenic medium, but Atp6v1h +/osteoblasts generated reduced RANKL and enhanced OPG signal under osteogenic conditions. This suggested that Atp6v1h +/osteoblasts sent decreased signals to induce osteoclast formation only under calcifying conditions. Thus the bone remodeling of heterozygous mice was impaired.
Our population data identified alleles of ATP6V1H gene as novel risk factors for osteopenia/ osteoporosis implying ATP6V1H alleles might also influence the risk for osteopenia. The relationship between ATP6V1H and bone density suggests that Atp6v1h serves as a regulator of osteoclast formation and function. In summary, the ATP6V1H gene could be an important target for therapeutic interventions aimed at preventing bone resorption and treating osteoporosis [38]. ; TGF-β1 also promotes the maturation of osteoclasts (OC). The increased intracellular pH resulted in less secretion of TGF-β1 from Atp6v1h +/osteoclasts. Decreased TGF-β1 negatively regulates osteoblasts differentiation, results in the fewer and smaller osteoblasts and finally causes the decreased bone formation in Atp6v1h +/-mice. The low level of TGF-β1 is associated with reduced RANKL and increased OPG signal, sending decreased signals to induce the osteoclasts formation. The decreased bone formation is more than decreased osteoclast function. Thus, the final in vivo net effect is bone loss.

Methods
Population study. A total of 1625 Chinese volunteers were recruited from Midwestern Chinese Han adults living in Xi'an and Changsha cities. The study was approved by Xi'an Jiaotong University. Signed informed consent documents were obtained from all study participants before entering the study. The following subjects were excluded from the study: Subjects with chronic diseases and conditions that might potentially affect bone mass, structure, or metabolism; subjects taking anti-bone resorptive or bone anabolic agents/drugs, such as bisphosphonates. BMD values at spine were measured using Hologic 4500W machines (Hologic Inc., Bedford, MA, USA). The coefficient of variation (CV) values of the dual-energy X-ray absorptiometry (DXA) measurements for spine was approximately 1.01% [22,23]. Three groups were generated based on Z-score. Group 1 represented those with a Z-score above -1, Group 2 between -1 and -2 and Group 3 below -2.
Genomic DNA was extracted from peripheral blood leukocytes and SNP genotyping was performed using Genome-Wide Human SNP Array 6.0 (Affymetrix, Santa Clara, CA, USA) according to the manufacture's protocol. Data management and analyses were conducted using the Genotyping Command Console. The minimum call rate was above 95%. A total of 15 tag SNPs within the intronic region of ATP6V1H gene were included in subsequent association analyses. PLINK was used to test for association and meta-analysis and calculations were done using the METAL software package. A raw P value of < 0.05 in our study was considered nominally significant.
Atp6v1h knockout mouse generation using CRISPR/Cas9 mediated genome editing. The mouse Atp6v1h gene (Transcript: ENSMUSG00000033793) is located on mouse chromosome 1 and has 14 exons with an ATG start codon in exon 2 and a TGA stop codon in exon 14. Atp6v1h knockout mice were generated using the CRISPR/Cas9 system. Exon 2 of the Atp6v1h gene was chosen to introduce frame shift changes and create an early stop codon. The sequence of gRNA (5'-GAT GGA CAT TCG AGG TGC TGT GG-3') was selected based on the ZiFiT tool. The gRNA expression vectors were constructed using a gRNA Cloning Vector from the Church Lab [39] by Gibson Assembly and a diagram of the resulting gRNA expression vector was produced. Atp6v1h-specific gRNA and Cas9 mRNA were generated according to mMESSAGE mMACHINE T7 Kit (Life Technologies), and then co-injected into fertilized eggs of C57BL/6 mice to create knockout mice (Cyagen, China). PCR genotyping and sequencing confirmation were performed to verify the frame shift changes for Atp6v1h knockout mouse generation. PCR primers for mouse tail genotyping were: gRNA1-Forward: 5'-AGTCAATAGAAAGCCT TTTATGGGG-3', gRNA1-Reverse: 5'-CTGAGGGCA CTCGCATTCAGGT-3'. The F0 founders had a 5bp (CGAGG) deletion and one base replacement (T>G) by the gRNA1/Cas9 procedure. All animals were treated in accordance with the ethical guidelines of the School of Stomatology, the Fourth Military Medical University (Xi'an, China).
Micro CT scanning and analysis. A micro CT scanner and Inveon Research Workplace (Siemens Inveon, Germany) were used for live animal scanning, tissue scanning and three-dimensional reconstruction. The femurs and vertebrae of mice were scanned with the parameters 80kV and 500μA, a 500 ms exposure time, and a scan angle of 360°. Three-dimensional images were reconstructed using standard methods. The bone density and other bone parameters were calculated and compared in the metaphyseal regions of distal femurs. The relative bone density, bone volume (BV/TV), bone surface (BS/TV), trabecular thickness (Tb.Th), trabecular number (Tb.N), trabecular separation (Tb.Sp), and cortical thickness (Ct. Th) of the femoral mid-diaphysis were calculated. The above trabecular morphometry measurements were also calculated in the third lumbar vertebrae of mice. Each group contained 8 mice.
Histology and histomorphometry. For dynamic histomorphometric analysis, 3-week-old littermate male wild-type and Atp6v1h +/mice were used. Calcein at 10 mg/kg body weight and xylenol orange at 90 mg/kg were injected on days 1 and 8, respectively. The mice were sacrificed 2 days after the second injection. The femurs were removed and fixed in 70% ethanol at 4°C, then subjected to standard tissue calcification procedures. Plastic sections were cut at a thickness of 50 μm with a Leica SP2600 (Leica Microsystems, Germany). Images were captured by FluoView™ FV1000 Confocal Microscope (Olympus Corporation, Japan), and the mineral apposition rate (MAR, mcm/d), the mineralizing surface (MS/BS, %) and bone formation rate (BFR/BS, mcm 3 /mcm 2 /d) were measured. Each group contained 6 to 8 mice.
Plastic sections were stained using the modified ponceau red method. Trabecular morphometry parameters, such as bone volume (BV/TV) and Tb.Th, were calculated and compared in the distal region of the femur. In each group, five 4-week-old littermate male wild-type and Atp6v1h +/mice were used. Each group contained 6 to 8 mice.
Femur samples from 3-months -old male wild-type and Atp6v1h +/mice were fixed with 4% paraformaldehyde and demineralized with 10% EDTA in PBS (pH 7.2~7.4) for 20~30 days at 37℃, followed by gradient dehydration with ethanol, embedded in paraffin and sectioned at a thickness of 5~7 μm. The slides went through with hematoxylin & eosin staining and modified Ponceau red staining and were stained for ALP using a commercial kit (Jiancheng Bioengineering Institute, Nanjing, China). The distribution and intensity of ALP staining in the distal region of the femur were measured and compared. Each group contained 8 mice. ALP+ osteoblasts adjacent to bone surface (N.Ob/mm 2 ) and the osteoblast surface (Ob.S/BS, %) were calculated following the rules recommended by the Nomenclature Committee of the American Society of Bone and Mineral Research [40].
Immunohistochemistry staining. The sABC immunohistochemistry staining kit (Boster, Wuhan, China) was used for immunohistochemical staining of ATP6V1H following the manufacturer's protocols. Rabbit polyclonal IgG anti-ATP6V1H (Origene, China was used in this experiment at a dilution of 1:400).
In vitro assay of osteoclasts. Bone marrow (BM) cells were separated from the femur and tibia of 6-weeks-old male wild-type and Atp6v1h +/mice and were quickly treated with Red Blood Cell Lysis Buffer (150 mM NH 4 Cl, 10mM KNCO 3 , 0.1 mM EDTA, pH 7.4) for 5min. The cells were cultured in the presence of M-CSF (50ng/ml) (R&D, USA) for 24 hours. After floating cells had been removed, the attached cells were used as osteoclast precursors and treated with M-CSF (50ng/ml) and RANKL (100ng/ml) (R&D, USA) for 7 days for the TRAP staining. The TRAP-positive multinucleated cells (≥5 nuclei per cells; >100µm in diameter) were counted. The cellular numbers (N.OC/well) and average areas (OC.S) were compared in the heterozygous and control groups.
Cells were fixed with 4% paraformaldehyde in PBS buffer and stained with rhodamine-phalloidin (Molecular Probes) for F-actin rings. In some experiments, TRAP staining was performed before F-actin ring staining.
BM cells were seeded on cattle bone slices (4mm×4mm×100µm) for 10-14 days with the RANKL and M-CSF treatment. The bone slices were fixed with 2.5% glutaraldehyde, stained with 0.1% toluidine blue and observed under a Leica DFC 550 Stereomicroscope (Leica Microsystems, Germany). The bone slices were also observed by scanning electron microscopy (Hitachi S-4800 high-resolution SEM, Japan). The number of resorptive lacunae of each bone slice was counted and the areas of the lacunae were measured using ImageJ software (NIH, USA) for at least 50 lacunae in each group.
The BM cells were also seeded on osteoassay surface polystyrene microplates (Corning, US) for 10 days with the RANKL and M-CSF treatment. The images of the wells were captured using a Leica DMI6000B microscope (Leica Microsystems, Germany). The number and size of osteoclasts areas, as well as the resorptive area (%) per well and osteoclast number normalized resorptive area (%) were measured using NIH Image J software (NIH, USA). siRNA transfection. Atp6v1h siRNA against the nucleotide 324 position of mouse Atp6v1h (NM_133826.4) was designed and synthesized (GenePharma, China). The sequence of Atp6v1h siRNA is: sense (5' to 3'): AAG CAA AUG CUU CAA ATT; antisense (5' to 3'): UUU GAA GCA UCU CUU GCU UTT. The negative control sequence is: sense (5' to 3'): UUC UUC GAA CGU GUC ACG UTT; antisense (5' to 3'): ACG UGA CAC GUU CGG AGA ATT. The BM cells from 6-week -old male wild-type and Atp6v1h +/mice were grown to 50% to 60% confluence and then transfected with siRNA using lipofectamine 2000 according to the manufacturer's protocols. The number of TRAP+ osteoclasts was determined and Q-PCR was performed.
Osteoblast functional assay. Calvarial osteoblasts (1 × 10 4 cells/well) were cultured in high glucose DMEM medium in the presence of 10% FBS (Gibco, US). For the osteogenic induction experiments, calvarial osteoblasts were grown in DMEM containing 10% FBS, 10mM β-glycerophosphate, 50 µg/ml ascorbic acid, and 10 -7 mM dexamethasone for 7 to 21 days. Osteogenic activity and mineralization nodule formation were assessed by staining with an ALP kit (Jiancheng Bioengineering Institute, Nanjing, China) and an alizarin red S kit (10 mg/ml, Sigma, US), respectively. Cell culture plates were photographed under a phase-contrast microscope using a Canon EF 100mm f/2.8 USM camera lens.
The third passage of wild-type osteoblasts as well as MC3T3-E1 cells were treated with Atp6v1h siRNA (same sequence as above), respectively. After 72 hours, the mRNA was collected and Q-PCR was performed as above. The mRNA of calvarial bone tissues from the newborn wild-type and heterozygous mice were also collected and Q-PCR was performed as above.
Co-culture of osteoblast and osteoclasts. Calvarial osteoblasts (1 × 10 3 cells/well) from wild-type and Atp6v1h +/mice were cultured in 96 well plate high glucose DMEM medium in the presence of 10% FBS (Gibco, US). The monocytes from wild-type and Atp6v1h +/mice were seeded on the osteoblasts with a ratio of 1:200 in the presence of 10 -8 M 1α,25(OH) 2 D 3 and 10 -6 M PGE2 (Sigma, US). After 9 days, TRAP staining was performed and TRAP+ osteoclasts were calculated.
Analysis of intracellular pH. BM cells from wild-type and Atp6v1h +/mice were treated with M-CSF and RANKL for 7 days, then stained with lysoSensor™ DND-160 (Invitrogen, US), and a pH calibration curve was performed at the end of each experiment to calibrate fluorescence intensity to actual pH values according to modified methods [43]. In each group, at least 50 cells were calculated and the experiment was repeated at least three times.
Biochemical and ELISA assay. Sera were collected from 3-month-old or 4-week-old male wild-type and Atp6v1h +/mice after overnight fasting.