2021 Volume 58 Issue 2 Pages 79-87
Skeletal muscle myoblasts are myogenic precursor cells that generate myofibers during muscle development and growth. We recently reported that broiler myoblasts, compared to layer myoblasts, proliferate and differentiate more actively and promptly into myocytes, which corresponds well with the muscle phenotype of broilers. Furthermore, RNA sequencing (RNA-seq) revealed that numerous genes are differentially expressed between layer and broiler myoblasts during myogenic differentiation. Based on the RNA-seq data, we herein report that chicken myoblasts transcribe endogenous retrovirus group K member (ERVK) genes. In total, 16 ERVKs were highly expressed in layer myoblasts and two (termed BrK1 and BrK2) were significantly induced in broiler myoblasts. These transcribed ERVKs had a total of 182 neighboring genes within ±100 kb on the chromosomes, of which 40% were concentrated within ±10 kb of the ERVKs. We further investigated whether the transcription of ERVKs affects the expression of their neighboring genes. BrK1 had two neighboring genes; LOC107052719 was overlapping with BrK1 and downregulated in the broiler myoblasts, and FAM19A2 was upregulated in the broiler myoblasts as well as BrK1. BrK2 had 14 neighboring genes, and only one gene, LOC772243, was differentially expressed between layer and broiler myoblasts. LOC772243 was overlapping with BrK2 and suppressed in the broiler myoblasts. These data indicate that the transcription of ERVKs may impact the expression of their neighboring genes in chicken myoblasts.
Broiler chickens are the breeds genetically selected for the rapid development of skeletal muscle tissue with a high feed efficiency in order to improve meat production. The skeletal muscle of the broilers is composed of a higher number of myofibers presenting larger diameters and faster growth compared to other breeds such as layers (Scheuermann et al., 2004). Myofibers are multinuclear contractile myocytes and have dozens of stem cells (named satellite cells) between the basal lamina and plasma membrane. In the process of muscle development and growth, satellite cells are transformed into the myogenic precursor cells termed myoblasts. Myoblasts proliferate through several rounds of cell divisions and then differentiate into mononuclear myocytes. Finally, myocytes fuse to form nascent myotubes, or fuse to preexisting myofibers (Dumont et al., 2015). Thus, the characteristics of myoblasts directly reflect the muscle phenotype of the animals.
We recently reported that myoblasts isolated from UK Chunky (UKC) broilers actively proliferate and promptly differentiate compared to those of White Leghorn (WL) layers. Comprehensive gene expression analyses by RNA sequencing (RNA-seq) revealed that numerous genes are differentially expressed between WL and UKC myoblasts during myogenic differentiation. UKC myoblasts highly express cell cycle genes in growing condition and muscular component genes in differentiation stages, which corresponds well with the cellular properties of UKC myoblasts and eventually with the muscle phenotype of broilers (Nihashi et al., 2019b). However, major myogenic regulatory factors (MRFs) such as MyoD and myogenin were not upregulated in UKC myoblasts as opposed to WL myoblasts. It is not fully understood why UKC myoblasts eventually achieve higher myogenic ability. Therefore, the precise understanding of the regulatory mechanisms of gene expression in myoblasts will contribute to improvements in chicken breeding and meat production.
Through the analyses of RNA-seq data, we found that chicken myoblasts transcribe endogenous retrovirus group K members (ERVKs). The RNA genome of the retrovirus, which commonly consists of a 5′ long terminal repeat (LTR), gag (group-specific antigen), pro (protease), pol (polymerase), env (envelope), and 3′ LTR, is reverse-transcribed and integrated into the genome of a host cell as a provirus. The provirus inserted into a germ cell is vertically inherited in that line. Such heritable proviruses are endogenous retroviruses (ERVs), which are able to be transcribed to make copies unless silenced (Gifford et al., 2018). Typical vertebrate genomes have thousands of ERVs. In chickens, eleven ERV families occupy > 21 Mb, which is ∼2% of the genome (Huda et al., 2008). Chicken ERVs mainly consist of ERVK and ERVL (Barr et al., 2005). Although ERVKs use a lysine (K) tRNA to bind to a viral primer binding site for reverse transcription, this does not imply that each ERVK member is phylogenetically related.
In humans, the mRNAs or the proteins derived from ERVKs play etiological roles in neurologic diseases (Manghera et al., 2014), cancer (Downey et al., 2015), and amyotrophic lateral sclerosis (Manghera et al., 2016). Accumulating evidence suggests that ERVK transcription potentially influences the physiological homeostasis of the host animal. As one of the mechanisms, ERVK is thought to regulate the expression of its neighboring genes. The LTR of ERVK serves as the promoter harboring multiple binding sites for transcription factors such as YY1 (Knössl et al., 1999), Sp1/3 (Fuchs et al., 2011), and nuclear factor-κB (NF-κB) (Manghera and Douville, 2013). The transcription factors recruited on the LTR of ERVK can modulate expression of not only ERVK itself but also its neighboring genes.
The human expressed sequence tag database revealed that skeletal muscles express ERVK but not ERVE, ERVH, or ERVW (Stauffer et al., 2004), which suggests that ERVK promoters are responsive in myogenic cells. In this study, we surveyed loci and neighboring genes of the ERVKs transcribed in chicken myoblasts. RNA-seq data demonstrated that the transcription of ERVKs may affect the expression of their neighboring genes in chicken myoblasts.
All experimental procedures were conducted in accordance with the Regulations for Animal Experimentation of Shinshu University. The animal experimentation protocol was approved by the Committee for Animal Experiments of Shinshu University.
Myoblasts from hindlimb muscle of WL or UKC chickens at embryonic day 10 were isolated, primary-cultured, and induced to differentiate as previously described (Takaya et al., 2017; Nihashi et al., 2019a, b). Myoblasts were cultured on dishes coated with collagen type I-C (Cellmatrix; Nitta Gelatin, Osaka, Japan) at 37°C with 5% CO2 throughout the experiments. Undifferentiated myoblasts were maintained in growth medium (GM) consisting of RPMI1640 (Nacalai, Osaka, Japan), 20% fetal bovine serum (FBS) (HyClone; GE Healthcare, UT, USA), 1% non-essential amino acids (Wako, Osaka, Japan), 1% chicken embryo extract (US Biological, MA, USA), 2 ng/mL basic fibroblast growth factor (Wako), and a mixed solution of 100 units/mL penicillin and 100 µg/mL streptomycin (PS) (Nacalai). To initiate myogenic differentiation, 2.0×105 myoblasts in GM were seeded on 60 mm dishes. Next day (defined as day 0), the myoblasts were induced to differentiate by replacing GM with differentiation medium consisting of DMEM (Nacalai), 2% FBS, and PS. Myoblasts from each breed were subjected to experiments three times independently. Total RNA was isolated from the myoblasts using NucleoSpin RNA Plus (Macherey-Nagel, Düren, Germany) on days 0, 1, and 2.
RNA-seq DataThe RNA-seq data used in this study have been reported (Nihashi et al., 2019b) and deposited into the DDBJ Sequence Read Archive (DRA; Research Organization of Information and Systems, National Institute of Genetics, Mishima, Japan) with the accession number DRA007964. Briefly, myoblasts from each breed were subjected to RNAseq in triplicate on differentiation days 0, 1, and 2. Expression levels of 26,640 transcripts were quantified as reads per kilobase per million reads (RPKM) or as transcripts per kilobase per million (TPM). The false discovery rate (FDR) was employed to correct the p values.
Survey for ERVK-neighboring GenesLoci of the ERVKs transcribed in chicken myoblasts were confirmed using Gene Tools (https://www.ncbi.nlm.nih.gov/gene/) (NCBI; National Center for Biotechnology Information, MD, USA). Gene annotations were performed using NCBI Gallus gallus Annotation Release 104 as a reference genome (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/Gallus_gallus/104/). The ERVK-neighboring gene was defined as the gene that has a start or end position located within ±100 kb of the ERVKs.
Statistical AnalysisThe data was presented as mean±standard error. Statistical differences between two groups were evaluated by unpaired two-tailed Student's t test. Statistical differences between more than three groups were evaluated by FDR p value. Statistical significance was set at p value <0.05.
RNA-seq detected the transcription of 44 ERVs in WL or UKC myoblasts on differentiation days 0, 1, or 2. One of them, LOC100857682, was an ERVFC1 encoding envelope polyprotein. The other 43 ERVs were ERVKs. We termed these ERVKs transcribed in chicken myoblasts “chMB-ERVKs”. As shown in Table 1, four chMB-ERVKs were Gag polyprotein genes and the remaining were Pol protein genes. Twenty-two chMB-ERVKs showed significant transcription levels (RPKM > 0.1). Sixteen of these were transcribed to a strikingly higher degree in layer WL myoblasts (named LyK1–LyK16) and two were significantly upregulated in broiler UKC myoblasts (named BrK1 and BrK2) on at least one of the differentiation days (FDR p value <0.05) (Fig. 1A). This clearly indicates that some chMB-ERVKs were differentially transcribed between WL and UKC myoblasts throughout myogenic differentiation.
Chr1 | Gene | Product | |
---|---|---|---|
1 | LOC107050425 | ERVK-25 Pol protein-like | |
1 | LOC107050867 | ERVK-18 Pol protein-like | |
1 | LOC107051640 | ERVK-113 Pol protein-like | |
1 | LOC107051645 | ERVK-18 Pol protein-like | |
1 | LOC107052718 | (BrK1) | ERVK-11 Pol protein |
1 | LOC107054241 | (LyK3) | ERVK-25 Pol protein-like |
1 | LOC107054837 | (LyK8) | ERVK-18 Pol protein-like |
1 | LOC107055731 | ERVK-25 Pol protein-like | |
2 | LOC770574 | (LyK5) | ERVK-25 Pol protein |
2 | LOC107051673 | (LyK15) | ERVK-18 Pol protein-like |
2 | LOC107052504 | (LyK10) | ERVK-18 Pol protein-like |
2 | LOC107052568 | (LyK2) | ERVK-113 Gag polyprotein-like |
2 | LOC107052661 | (LyK7) | ERVK-18 Pol protein-like |
2 | LOC107052685 | ERVK-11 Pol protein-like | |
2 | LOC107052829 | (LyK14) | ERVK-113 Pol protein-like |
2 | LOC107052830 | ERVK-113 Pol protein-like | |
3 | LOC107053122 | ERVK-18 Pol protein-like | |
4 | LOC107053216 | ERVK-113 Pol protein-like | |
4 | LOC107056268 | (LyK4) | ERVK-18 Pol protein-like |
5 | LOC107053591 | (LyK6) | ERVK-18 Pol protein-like |
5 | LOC107056419 | (LyK12) | ERVK-18 Pol protein-like |
8 | LOC107054036 | (BrK2) | ERVK-9 Pol protein-like |
18 | LOC107057172 | ERVK-18 Pol protein-like | |
18 | LOC107057173 | ERVK-8 Gag polyprotein-like | |
28 | LOC107055329 | (LyK11) | ERVK-18 Pol protein-like |
33 | LOC107049264 | (LyK9) | ERVK-8 Gag polyprotein-like |
W | LOC107049272 | ERVK-8 Gag polyprotein-like | |
W | LOC107049273 | ERVK-9 Pol protein-like | |
W | LOC107049274 | ERVK-18 Pol protein-like | |
W | LOC107049276 | ERVK-18 Pol protein-like | |
W | LOC107049277 | ERVK-9 Pol protein-like | |
W | LOC107055436 | (LyK16) | ERVK-8 Pol protein-like |
W | LOC107055437 | ERVK-18 Pol protein-like | |
W | LOC107055442 | ERVK-18 Pol protein-like | |
W | LOC107055443 | ERVK-25 Pol protein-like | |
W | LOC107055445 | ERVK-25 Pol protein-like | |
Z | LOC107049326 | ERVK-18 Pol protein-like | |
Z | LOC107049332 | (LyK13) | ERVK-113 Pol protein-like |
Z | LOC107049350 | ERVK-25 Pol protein-like | |
Z | LOC107051888 | ERVK-25 Pol protein-like | |
Z | LOC107051891 | ERVK-8 Pol protein-like | |
Z | LOC107052380 | (LyK1) | ERVK-25 Pol protein-like |
ND2 | LOC107049474 | ERVK-9 Pol protein-like |
Transcription of ERVKs and neighboring genes in chicken myoblasts. (A) Transcription levels (RPKM) of the 22 chMB-ERVKs with RPKM > 0.1 in WL and UKC myoblasts on differentiation days 0, 1, and 2. *FDR p<0.05 vs UKC. † FDR p<0.05 vs WL. (B) Total number of neighboring genes and distances from the 43 chMB-ERVKs (16 LyKs, two BrKs, and others) on the chicken chromosomes.
It has been reported that an intensive gene expression spreads into its physically neighboring genes to activate their expression. The radius of such an intergenic transcriptional ripple effect is estimated to be around 100 kb on the chromosome (Ebisuya et al., 2008). To assess the influence of chMB-ERVK on its neighbor transcription, the loci of the 43 chMB-ERVKs and their neighboring genes on the Gallus gallus genome were carefully surveyed. One chMB-ERVK, LOC107055731, had no neighboring genes. However, the remaining 42 chMB-ERVKs had a total of 182 genes within ±100 kb. Remarkably, 72 of the 182 genes (39.6%) were intensively located within ±10 kb of chMB-ERVKs (Fig. 1B). In other words, 37/43 (86.0%) of chMB-ERVKs had been integrated within ±10 kb of the genes. In humans, ERVKs are considered to be randomly and independently inserted into the genome (Kurdyukov et al., 2001). In chickens, ERVs are enriched in gene-sparse chromosomes and intergenic regions. Indeed, only 11.2% of chicken ERVKs are present in active transcription units (Barr et al., 2005). Contrary to such propensities of ERVs, most of the chMB-ERVKs had neighboring genes. It is thought that the ERVKs integrated near the genes have remained to be transcribed in chicken myoblasts.
A total of 182 chMB-ERVK-neighboring genes comprised 158 unique genes (Supplemental Table), of which 96 (60.8%) encoded proteins, and 44 (27.8%) were non-coding (Table 2). Intriguingly, the 96 coding genes included 15 ERVKs and four envelope genes, indicating that chMB-ERVKs form clusters with ERVs on the chromosomes. For instance, LyK2 and LyK10 were overlapping each other. Another major class of chMB-ERVK-neighboring genes was olfactory receptors but almost all of them (20/22) were LyK4-neighboring genes. The chicken genome has at least 554 olfactory receptor genes (Niimura and Nei, 2005), and they form genomic clusters on chromosomes as a result of duplication (International Chicken Genome Sequencing Consortium, 2004). LyK4 was perhaps accidentally integrated into the loci associated with an olfactory receptor gene cluster. This would explain why the olfactory receptor became one of the main components of chMB-ERVK-neighboring genes.
Class | No. of genes | (%) |
---|---|---|
Coding gene | ||
Characterized endogenous gene (except for olfactory receptor) | 28 | 17.7% |
Olfactory receptor | 22 | 13.9% |
ERVK | 15 | 9.5% |
Envelope glycoprotein gp95-like | 4 | 2.5% |
Predicted ORF1 | 27 | 17.1% |
Non-coding gene | ||
lincRNA2 | 38 | 24.1% |
miRNA3 | 5 | 3.2% |
snRNA4 | 1 | 0.6% |
Uncharacterized gene | 18 | 11.4% |
It was subsequently investigated whether chMB-ERVKs affect the transcription of their neighboring genes. RNA-seq data showed that the expression levels of all 79 neighboring genes of 16 LyKs were not significantly different between WL and UKC myoblasts except for a combination of LyK2 and LyK10 as the neighboring genes of each other. LyKs that highly expressed in WL myoblasts are considered not to influence their neighboring gene transcription.
Next, the expression levels of a total of 16 neighboring genes of BrK1 and BrK2 were dissected (Table 3). BrK1 had two neighboring genes, LOC107052719 and FAM19A2 (Fig. 2A). LOC107052719 expression was significantly suppressed in UKC myoblasts compared to those in WL myoblasts at early stages of myogenic differentiation both in RPKM and TPM values (Fig. 2C). It is of interest that BrK1 and LOC107052719 were overlapping on the identical strand and showed opposite transcription patterns. Conversely, FAM19A2 was strikingly upregulated in UKC myoblasts throughout differentiation (Fig. 2D), which is in accordance with the BrK1 transcription.
Gene | Product | Distance (bp) |
---|---|---|
BrK1: LOC107052718 | ||
LOC107052719 | Uncharacterized gene | Overlap |
FAM19A2 | Family with sequence similarity 19 member A2, C-C motif chemokine-like | 93,542 |
BrK2: LOC107054036 | ||
SLC1A7 | Solute carrier family 1 member 7 | −92,653 |
PODN | Podcan | −61,675 |
SCP2 | Sterolcarrier protein 2 | −49,303 |
LOC107054037 | Uncharacterized gene | −46,114 |
ECHDC2 | Enoyl-CoA hydratase domain containing 2 | −43,020 |
COA7 | Cytochrome c oxidase assembly factor 7 | −18,667 |
MIR6549 | miRNA | −1,396 |
LOC772243 | Envelope glycoprotein gp95-like | Overlap |
ZYG11B | Zyg-11 family member B, cell cycle regulator | Overlap |
FAM159A | Family with sequence similarity 159 member A | 2,405 |
GPX7 | Glutathione peroxidase 7 | 27,744 |
ZCCHC11 | Zinc finger CCHC-type containing 11 | 39,929 |
PRPF38A | Pre-mRNA processing factor 38A | 84,605 |
ORC1 | Origin recognition complex subunit 1 | 89,080 |
Loci and transcription levels of BrK1- and BrK2-neighboring genes. (A and B) Physical maps of BrK1, BrK2, and their neighboring genes. Arrows indicate sequence direction (5′–3′). (A) BrK1 and its neighboring genes on chromosome 1. (B) BrK2 and its neighboring genes on chromosome 8. (C–E) The transcription levels defined by RPKM and TPM of LOC107052719 (C), FAM19A2 (D), and LOC772243 (E) in WL and UKC myoblasts on differentiation days 0, 1, and 2. * p<0.05, ** p<0.01 vs WL (Student's t test at each day).
BrK2 had 14 neighboring genes, of which seven were located upstream, five were located downstream, and two were overlapping with BrK2 (Fig. 2B). RNA-seq data showed that expression levels of BrK2-neighboring genes except for LOC772243 were not different between WL and UKC myoblasts. Only the expression of LOC772243 was significantly downregulated in UKC myoblasts as opposed to BrK2 (Fig. 2E). LOC772243 encodes for a murine leukemia virus (MLV) -related proviral envelope polyprotein-like protein. MLV originally exists as an exogenous virus and serves as an ERV like ERVK (Stocking and Kozak, 2008). The BrK2 and LOC772243, which overlap on the same strand and negatively regulate the transcription of one another, might be an important instance of retrovirus-mediated gene expression in host cells.
This is the first report of the ERVKs that are transcribed in chicken skeletal muscle myoblasts. Generally in host cells, most ERVs are silenced by chromatin modifications such as hypermethylation (Schulz et al., 2006). It is predicted that the 43 chMB-ERVKs identified in this study are still epigenetically active in chicken myoblasts, which is supported by the fact that almost all chMB-ERVKs are locating near genes. Eighteen chMB-ERVKs were found to be differentially transcribed between WL and UKC myoblasts during differentiation. Previous studies have reported that the myoblast genome is epigenetically modified during myogenic differentiation (Blum, 2014). We recently showed that proliferating and differentiating abilities are different between WL and UKC myoblasts (Nihashi et al., 2019b), suggesting that the epigenome linking myoblast properties may also differ among the chickens. Genome-wide epigenetic comparison among the breeds will give an insight into the relationship between chromatin state and chMB-ERVK transcription.
chMB-ERVKs were physically associated with endogenous genes on chicken chromosomes despite the fact that ERVs tend to be integrated into intergenic regions (Barr et al., 2005). This suggests that the loci having chMB-ERVKs are active in myoblasts. However, there is no information about the roles of chMB-ERVK-neighboring genes in the MRF-driven myogenic program. At present, the transcription of chMB-ERVKs and their neighboring genes seems not to be directly involved in the principal myogenic process regulated by MRFs. The functions of chMB-ERVK-neighboring genes during myoblast differentiation should be clarified in the future.
RNA-seq data showed that the transcription levels of LyK-neighboring genes were not different between WL and UKC myoblast while both BrK1-neighboring genes (LOC107052719 and FAM19A2) and one BrK2-neighboring gene (LOC772243) were differentially expressed among the breeds. BrK1 and BrK2 displayed the first and second highest expression levels among chMB-ERVKs. Transcriptional ripple effect occurs at the locus intensively transcribed (Ebisuya et al., 2008). This could be one of the reasons why BrK1 and BrK2 seemed to affect the expression of their neighboring genes.
FAM19A2 was strongly upregulated in UKC myoblasts in parallel with BrK1. FAM19A2 (TAFA2 in human) encodes for the chemoattractant neurokine potently expressed in the brain (Tom Tang et al., 2004). A recent study reported that FAM19A2 protein enhances the migration, recruitment, and proliferation of human mesenchymal stem cells without affecting differentiation into osteocytes and adipocytes (Jafari et al., 2019). Mesenchymal stem cells can also differentiate into muscle cells. Therefore, it is possible that the myoblast-secreting FAM19A2 protein recruits mesenchymal stem cells to muscle tissue during muscle development and growth in broilers. RNA-seq data obtained in this study is still insufficient to conclude that BrK1 transcription actually enhances FAM19A2 expression. Although the transcriptional regulatory role of FAM19A2 has not been studied, FAM19A2 locus was reported to be clinically associated with various pathophysiological conditions such as intellectual disability (Borg et al., 2005), systemic lupus erythematosus (Lessard et al., 2012), vital capacity (Parker et al., 2014), and insulin sensitivity (Walford et al., 2016) in humans. These studies suggest that FAM19A2 expression can be altered by a variety of biological factors and processes. To clarify the relationship between BrK1 and FAM19A2 transcription, the role of FAM19A2 in myoblasts need to be further investigated.
The transcription levels of LOC107052719 and LOC772243 were significantly suppressed in UKC myoblasts in contrast to those of BrK1 and BrK2. LOC107052719 and LOC772243 were overlapping with BrK1 and BrK2 on the identical strands, respectively. Almost of such paired overlapping genes are expressed in parallel. However, a part of overlapping gene pairs exhibit the expression pattern that one is induced and another is suppressed in specific tissues (Chen et al., 2019). BrK1/LOC107052719 and BrK2/LOC772243 may also be the gene pairs that are differentially expressed in UKC myoblasts. The transcription of another BrK2-overlapping gene, ZYG11B, was not altered between WL and UKC myoblasts. ZYG11B was encoded on the antisense strand of BrK2, and its overlapping region occupied only a small portion of the ZYG11B gene structure. The transcription level of BrK2, which was similar to those of LyKs might be inadequate to modulate the expression of its neighboring genes including ZYG11B.
In any case, the causality of transcription levels of chMB-ERVKs and their neighboring genes so far is ambiguous. To address this issue, the transcription of chMB-ERVKs in myoblasts need to be carefully investigated in various situations because transcriptional regulation in myoblasts can be different depending on developmental origin or phase. For example in mice, regulatory systems of MRFs are distinct between facial and trunk muscles (Lu et al., 2002; Moncaut et al., 2012). It is possible that the transcription of ERVKs in myoblasts may also differ between breeds, growth stages, or muscle tissues. The myoblasts used in this study were initially isolated from hindlimb muscles of chicken embryos. It has not been strictly compared yet whether cellular properties of chicken myoblasts are distinctive for muscle tissues. The broiler breast muscle, which is mostly composed of fast-type myofibers frequently shows acute and chronic inflammatory and myopathic changes (Kuttappan et al., 2013; Hosotani et al., 2020). Intriguingly, ERVK transcription is known to be activated by pathogenic transcription factors such as NF-κB and interferon regulatory factors (Manghera and Douville, 2013). These findings suggest that ERVKs in myoblasts would be induced by extracellular factors or the environment. The transcription of chMB-ERVKs and their neighboring genes in vivo should be elucidated in further studies.
This study demonstrated that chicken myoblasts transcribe ERVKs, that transcription levels of these chMB-ERVKs are different among the breeds, and that chMB-ERVK transcription may be linked to the expression of their neighboring genes. These data suggest that chMB-ERVK-mediated gene expression may be involved in the molecular characteristics of myoblasts.
This study was supported in part by Grants-in-Aid from The Ito Foundation and from The Shinshu Foundation for Promotion of Agricultural and Forest Science to T. T., a Research Fellowship for Young Scientists from The Japan Society for the Promotion of Science (19J20888) and a Grant-in-Aid from The Fund of Nagano Prefecture to Promote Scientific Activity (H-29-3-12) to Y. N., and Grants-in-Aid from The Japan Society for the Promotion of Science to T. T. (19K05948), T. O. (16K07986), and H. K. (18K05939). RNA-seq data analysis was supported by the Cooperative Research Grant of the Genome Research for BioResource, NODAI Genome Research Center, Tokyo University of Agriculture.
The authors declare no conflict of interest.