Transcriptome analysis of female western flower thrips, Frankliniella occidentalis, exhibiting neo-panoistic ovarian development

The western flower thrips, Frankliniella occidentalis, is one of the most devastating insect pests with explosive reproductive potential. However, its reproductive physiological processes are not well understood. This study reports the ovarian development and associated transcriptomes of F. occidentalis. Each ovary consisted of four ovarioles, each of which contained a maximum of nine follicles in the vitellarium. The germarium consisted of several dividing cells forming a germ cell cluster, presumably consisting of oocytes and nurse cells. The nurse cells were restricted to the germarium while the subsequent follicles did not possess nurse cells or a nutritive cord, supporting the neo-panoistic ovariole usually found in thysanopteran insects. Oocyte development was completed 72 h after adult emergence (AAE). Transcriptome analysis was performed at mid (36 h AAE) and late (60 h AAE) ovarian developmental stages using RNA sequencing (RNASeq) technology. More than 120 million reads per replication were matched to ≈ 15,000 F. occidentalis genes. Almost 500 genes were differentially expressed at each of the mid and late ovarian developmental stages. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that these differentially expressed genes (DEGs) were associated with metabolic pathways along with protein and nucleic acid biosynthesis. In both ovarian developmental stages, vitellogenin, mucin, and chorion genes were highly (> 8-fold) expressed. Endocrine signals associated with ovarian development were further investigated from the DEGs. Insulin and juvenile hormone signals were upregulated only at 36 h AAE, whereas the ecdysteroid signal was highly maintained at 60 h AAE. This study reports the transcriptome associated with the ovarian development of F. occidentalis, which possesses a neo-panoistic ovariole.


Introduction
The western flower thrip, Frankliniella occidentalis (Pergande) (Thysanoptera: Thripidae), is one of the most devastating insect pests to many horticultural crops, especially those in greenhouses [1]. Both the larval and adult stages cause damage to plants by directly feeding on leaves or flowers. Especially, adults transmit plant viruses including tomato spotted wilt virus (TSWV) [2]. TSWV infection becomes serious and causes massive economic loss in hot pepper production in Korea [3]. This pest, originally native to North America, has spread to more than 60 countries since the late 1970s, including Canada, Australia, the United Kingdom, and far East Asian countries [4]. Various techniques such as chemical insecticides, entomopathogens, and pheromone traps have been applied to control F. occidentalis without satisfactory efficacy due to the insect's specific hiding behavior and insecticide resistance [5]. The thrips exhibits arrhenotokous parthenogenesis, with females developing from fertilized eggs and males from unfertilized eggs [6]. A brief immature period less than 10 days along with this various reproductive modes allow the thrips to rapidly build up the field populations during crop cultivating periods and so frequently leads to outbreaks beyond economic injury level [6]. Along with high reproductive potential, this type of mating behavior contributes to a rapid population increase and the development of insecticide resistance [7]. However, the molecular processes underlying reproduction and its regulation in this species remain unclear.
To investigate the physiological processes of the ovarian development of F. occidentalis, transcriptome analysis is useful for understanding the expressed genes associated with reproduction. A draft genome (415.8 Mb) of F. occidentalis was sequenced and its 16,859 genes were annotated into different functional categories including chemosensory receptors, detoxification, salivary gland, immunity, and development [8]. This suggests that RNA sequencing (RNASeq) analysis would be highly validated by this genomic information.
To identify the genes associated with ovarian development, this study investigated the ovarian development of F. occidentalis after adult emergence. After determining the mid and late ovarian developmental stages, transcriptomes were assessed using the NovaSeq 6000 platform. Subsequent differentially expressed gene (DEG) analysis in different developmental stages of female adults predicted the genes associated with ovarian development.

Thrip rearing
F. occidentalis adults were obtained from Bio Utility, Inc. (Andong, Korea) and reared in a laboratory under conditions of 27 ± 1˚C constant temperature, a 16:8 h (light:dark) photoperiod, and relative humidity of 60 ± 5%. The insects were reared on sprouted bean seed kernels.

Dissection of ovaries and microscopic observation
Different ages of female western flower thrips were dissected in 1 × phosphate-buffered saline (PBS) under a stereomicroscope at 30× magnification. PBS was prepared with 100 mM phosphate buffer containing 0.7% NaCl (pH 7.4). The ovaries were pulled from the abdominal tip and fixed with 3.7% paraformaldehyde in a wet chamber under darkness at room temperature (RT) for 60 min. After washing three times with 1 × PBS, the cells in the ovarioles were permeabilized with 0.2% Triton X-100 in 1 × PBS at RT for 20 min. The cells were then washed three times and blocked with 5% skim milk (MB cell, Seoul, Korea) in 1 × PBS at RT for 60 min. After washing three times, the cells were incubated with DAPI (4 0 ,6-diamidino-2-phenylindole, 1 mg/mL) diluted 1,000 times in PBS at RT for 2 min for nuclear staining. After washing three times, the ovarian cells were observed under a fluorescence microscope (DM2500, Leica, Wetzlar, Germany) at 200× magnification.

RNA extraction and RNASeq analysis
Total RNAs were extracted from the whole bodies of female F. occidentalis at different ages (0, 36, and 60 h after adult emergence). Three independent samples were used for three replications at each age. Each sample consisted of 50 females. RNA extraction was performed using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. An RNA library was generated using TruSeq Stranded Total RNA with Ribo-Zero H/M/ R_Gold (Illumina, San Diego, CA, USA). RNA sequencing was performed on the NovaSeq 6000 platform (Illumina) from Macrogen (Seoul, Korea). The RNA sequence was trimmed using CLC Workbench (QIAgen, Hilden, Germany). To calculate relative transcript accumulation, reads per kilobase per million (RPKM) mappable reads of the F. occidentalis genome (GenBank accession number: GCF_000697945.2) were estimated using CLC Workbench based on a template of Focc_2.1 version with a trimmed sequence of more than 50 bp.

Bioinformatics
DEGs were selected based on a fold change of � 2.0 and a P-value of < 0.05 with three biological replicates by comparing RPKM values at 36 h or 60 h after adult emergence (AAE) to those at 0 h AAE. KEGG pathway analysis was performed to test the statistical enrichment of DEGs using the KEGG mapper (https://www.kegg.jp/kegg/) by converting the National Center for Biotechnology Information (NCBI) Gene ID to a KEGG ID through the convert ID tool of KEGG mapper.

RT-qPCR
After RNA extraction, RNA was resuspended in nuclease-free water and quantified using a spectrophotometer (NanoDrop, Thermo Fisher Scientific, Wilmington, DE, USA). RNA (500 ng) was used for cDNA synthesis with RT PreMix (Intron Biotechnology, Seoul, Korea) containing oligo dT primers according to the manufacturer's instructions. All gene expression levels in this study were determined using a real-time polymerase chain reaction (PCR) machine (Step One Plus Real-Time PCR System, Applied Biosystem, Singapore) under the guidelines of [9]. Real-time PCR was conducted in a reaction volume of 20 μL containing 10 μL of Power SYBR Green PCR Master Mix (Thermo Scientific Korea), 3 μL of cDNA template (200 ng), and each 1 μL (10 pmol) of forward and reverse primers (S1 Table). After initial heat treatment at 95˚C for 2 min, qPCR was performed with 40 cycles of denaturation at 95˚C for 30 sec, annealing at 53~55˚C for 30 sec, and extension at 72˚C for 30 sec. The expression level of elongation factor-1 (EF-1, S1 Table) was used as a reference to normalize the target gene expression levels under different treatments. Quantitative analysis was performed using the comparative CT (2 -ΔΔCT ) method [10]. All experiments were independently replicated three times.

Statistical analysis
All the continuous variable data were subjected to a one-way analysis of variance using PROC GLM in the SAS program [11]. Means were compared with Duncan's multiple range test (DMRT) at type I error = 0.05.

Ovarian development in F. occidentalis
Dissection of 3-day-old female adults showed that a pair of ovaries contained eight ovarioles ( Fig 1A). The four ovarioles in each ovary were joined to a lateral oviduct and the two lateral oviducts were combined to the common oviduct. Each ovariole contained a string of follicles and was divided into a germarium and a vitellarium depending upon the presence of matured oocytes (Fig 1B). In the germarium, small cells were closely attached and formed a germ cell Ovariole structure divided into the germarium and the vitellarium. A total of 9 follicles are numbered from proximally to distally in the differential interference contrast ('DIC') picture. The white follicle represents chorionated oocytes. In the DAPI picture, each oocyte in the follicle is surrounded by follicular epithelium ('FE'). In the germarium, the germ cell cluster ('GCC') is located near to the oocyte ('OC'). (C) cluster. The vitellarium contained nine follicles, each of which was composed of oocytes and follicular epithelium without any nurse cells. The terminal follicles were usually chorionated.
Ovarian development occurred after adult emergence ( Fig 1C). The ovariole length began to extend just after adult emergence and reached a maximal size three days after emergence. Chorionated oocytes were visible after two days in some ovarioles. After three days, all the ovarioles had terminal chorionated oocytes, but some of the ovarioles lost them due to oviposition. This ovarian development pattern allowed us to determine three stages: early at 0 h, mid at 36 h, and late at 60 h after adult emergence.

Changes in the transcriptomes of female F. occidentalis adults during ovarian development
Total RNA was sequenced at three ovarian developmental stages of F. occidentalis. Each of the nine samples (= 3 stages × 3 replications) was sequenced from 10~13 Gb (Table 1). After trimming, 102~131 million reads in each sample were used to map to the F. occidentalis genome. With 71~82% mapping rates, the reads in each sample were matched to 14,042~14,521 genes among 16,859 predicted F. occidentalis genes. All the nine transcriptomes were deposited to GenBank with accession numbers of PRJNA833754.
When the transcriptomes of the three developmental stages were compared, they shared more than 95% (= 14,370/15,055) of the transcripts (Fig 2A). Forty-nine unique genes expressed at 36 h AAE were classified according to structure, gene regulation, and cell cycle along with several uncharacterized genes (S2 Table). However, their expression levels were extremely low at 0.0007~0.0468 RPKM. Forty-seven unique genes expressed at 60 h AAE were classified according to structure, protein-processing, and gene regulation along with several uncharacterized genes (S3 Table). Their expression levels were relatively high at 1.0016~1.6941 RPKM. However, this unique gene analysis did not identify genes apparently associated with oogenesis.
Ovarian development with female age. The ovariole length represents the germarium and vitellarium. Chorionated oocytes were counted per ovary. Each measurement used individual thrips and was replicated 5 times.
https://doi.org/10.1371/journal.pone.0272399.g001  4 Trimmed by CLC Workbench (QIAgen, Hilden, Germany). 5 Mapping to the F. occidentalis genome (GenBank accession number: GCF_000697945.2). 6 Total number of annotated genes was 16,859.  Table 1 To identify F. occidentalis genes associated with ovarian development, a differentially expressed gene (DEG) analysis was performed (Fig 2B). At 36 h AAE, 473 transcripts showed more than 2-fold increases in gene expression levels compared to those at 0 h AAE. At 60 h AAE, 443 transcripts showed more than 2-fold increases in gene expression levels compared to those at 0 h AAE. Only 19 transcripts showed changes between 36 h and 60 h AAE.
To characterize the DEGs at the mid (36 h AAE) and late (60 h AAE) ovarian developmental stages, their gene functions were predicted using KEGG analysis (Fig 3). Both DEGs were assigned to 52 KEGG functional categories but they did not exactly overlap. The DEGs of each developmental stage were assigned to 49 KEGG categories with three different missing categories. DEGs at 36 h AAE did not include the three categories of #21 (FoxO signaling pathway), #42 (protein export), and #50 (N-glycan biosynthesis), whereas the DEGs at 60 h AAE did not include the three categories of #17 (fatty acid biosynthesis), #30 (lysine degradation), and #39 (pentose and glucuronate interconversion). Among 46 common KEGG categories, most DEGs were classified into the metabolic pathway category (#33) in both developmental stages. The other major (> 10 DEGs) categories, biosynthesis of cofactors (#8) and lysosome (#31), were common in both developmental stages. However, the 36 h AAE samples had more DEGs in the inositol phosphate metabolism category (#28) than the 60 h AAE samples. In contrast, the glycan biosynthesis (#35) and purine metabolism (#44) categories had more DEGs at 60 h AAE than at 36 h AAE. The KEGG analysis suggested the upregulation of metabolic pathways associated with nucleic acid biosynthesis during ovarian development.

Expression profiles of egg proteins during ovarian development
To identify the specific genes associated with oogenesis in F. occidentalis, we selected genes that were highly expressed more than 8-fold at 36 h AAE or 60 h AAE compared to expression levels at 0 h AAE ( Table 2). The selected 99 genes were subdivided into structure, protein processing, lipid metabolism, gene regulation, and others. The structural protein category included typical egg proteins such as vitellogenin, chorion protein, mucin, and yellow melanization protein. In contrast, the highly suppressed genes at these stages included 37 genes (S4 Table). Especially, larval cuticular protein genes were included in the suppressed gene category.
The expression patterns of representative egg proteins during adult development were further analyzed (Fig 4). As expected, RNASeq analysis found that vitellogenin, chorion protein, mucin, and yellow genes were highly expressed in the mid and late ovarian developmental stages ( Fig 4A). However, there was little or no difference in the expression levels between the two developmental stages (36 h AAE and 60 h AAE). These transcript level profiles were further assessed by RT-qPCR with additional development stages (Fig 4B). The expression levels measured by RT-qPCR were mostly consistent with the expression profiles measured by RNA-Seq. However, RT-qPCR analysis indicated that the mucin and yellow genes were induced earlier than 36 h AAE. It also showed their expression patterns in late ovarian development after 60 h AAE, in which vitellogenin and mucin maintained the induced levels, whereas the expression levels of chorion protein and yellow significantly decreased.

Expression profiles of genes associated with endocrine signals during ovarian development
Juvenile hormone (JH), ecdysteroid, and insulin-like peptide (ILP) are well-known endocrine mediators during insect reproduction [12]. Genes associated with these endocrine signals were selected from the transcriptomes (Table 3). JH acid methyltransferase (JHAMT), JH esterase/ JH epoxide hydrolase, and Met involved in JH synthesis, JH degradation, and the JH receptor,   respectively, were identified. The RNASeq analysis showed that JHAMT was highly upregulated at 36 h AAE, whereas Met expression levels did not change during adult development ( Fig 5A). As ecdysteroid signaling genes, Shade and EcR were found in the adult transcriptomes. The expression levels of Shade were upregulated at mid and late ovarian development, during which EcR expression was slightly decreased. Insulin-like peptides and receptors were included in the adult transcriptomes. ILP was highly upregulated at 36 h AAE, whereas InR expression levels decreased during adult development. The transcriptome profiles related to endocrine signals were further supported by RT-qPCR analysis (Fig 5B). The RT-qPCR analysis assessed at several time points after adult emergence showed the upregulation of these endocrine signals during ovarian development. Especially, JH and insulin signals were more upregulated at 36 h AAE than at 60 h AAE, whereas the ecdysteroid signal was more upregulated at the late developmental stage.

Discussion
Various reproductive modes are observed in Thysanoptera. In F. occidentalis, a female produces progeny in bisexual or asexual arrhenotokous reproduction [13]. In arrhenotokous reproduction in the situation of few males (e.g., overwintering population), virgin females produce only male offspring [6]. When their sons are sexually mature, the females undergo bisexual reproduction with their sons and produce female-biased offspring. All thrips including F. occidentalis exhibit oviparous reproduction, in which females in either the asexual or bisexual mode undergo oogenesis in their ovaries. However, little is known about ovarian development in F. occidentalis. The current study analyzed the internal reproductive organ structure of F. occidentalis to assess the developing oocytes. Based on the temporal developmental pattern, RNASeq analysis was performed at early, mid, and late stages to determine the specific genes associated with ovarian development. The ovaries and associated internal reproductive organs of F. occidentalis were observed. Eight ovarioles from a pair of ovaries contained follicles, the end of which contained chorionated oocytes. Each follicle was composed of oocytes and follicular epithelium. The distal germarium contained a germ cell cluster, presumably consisting of interconnected nurse cells and oocytes. The ovarioles of insects are categorized into panoistic and meroistic types, in which the latter type is subdivided into polytrophic and telotrophic groups [14]. The panoistic type is considered to be the most ancestral because of its deficiency in transforming oogonia into nurse cells [15]. The polytrophic meroistic ovary has evolved from the panoistic ovary through the differentiation of nurse cells, and finally, the telotrophic meroistic type is derived from the polytrophic meroistic type by the restriction of nurse cells to the germarium. A deviation from the typical ovariole types is observed in Thysanoptera, in which the germ cell cluster is formed as seen in the ovary of the terebrantian thrip, Purthenothrips drucenae [16], suggesting that the panoistic follicles resulted from the secondary loss of nurse cells from the germ cell cluster. Stys et al. [14] called this type of ovary "neo-panoistic." Thus, the thysanopteran ovary provided new insight into the evolution of insect ovaries. Later, tubuliferan thrips also showed germ cell clusters, indicating that the neo-panoistic ovariole-type prevailed in Thysanoptera [17]. This was supported by our current study using F. occidentalis ovarioles, which did not have nurse cells in the follicles, while a germ cell cluster was found in the germarium. The ovaries of F. occidentalis grew just after adult emergence. During this period, the ovariole size increased along with oocyte development, and the final follicles in each ovariole had chorionated oocytes. Oogenesis is a sequential process consisting of previtellogenic development, vitellogenesis, and choriogenesis [18]. Previtellogenic development occurs in the germarium at the distal part of each ovariole and forms oocytes from the oogonial stem cells by mitosis and meiosis. Vitellogenesis is the process of accumulating vitellogenin (Vg) and other biomaterials into growing oocytes. After oocytes are fully grown, they are coated with chorion proteins secreted from the follicular epithelium to become eggs at the proximal part of the ovarioles. These eggs are then ovulated to the oviducts and fertilized just before oviposition. This ovarian developmental scenario suggests a sequence of oogenesis events in the neopanoistic ovariole of F. occidentalis. First, primary oocytes may be produced from germ cell clusters in the germarium. Second, the oocytes grow in size by accumulating Vg during vitellogenesis. Finally, the follicular epithelium forms the chorion of the fully grown oocytes during choriogenesis. Our transcriptome analysis supported the oogenesis processes by providing expression profiles of Vg, chorion proteins, mucin, and yellow genes during oogenesis.
RNASeq analysis used the NovaSeq platform, which sequenced more than 100 million reads per sample and resulted in more than 80% gene mapping rates. The first draft genome of F. occidentalis was reported and 16,859 genes were annotated [8]. Our RNASeq data from three ovarian developmental stages were mapped to 15,083 genes. Most of the mapped genes were shared among the three ovarian developmental stages. DEG analysis identified 473 and 443 DEGs in the mid and late stages, respectively. These DEGs were associated with metabolic pathways, especially related to nucleic acid biosynthesis and cofactor and amino acid biosynthetic pathways. The findings suggest that ovarian development in F. occidentalis requires a massive supply of raw materials such as nucleic acids and proteins.
Egg proteins were selected from the transcriptomes of the ovarian developmental stages. These genes represented highly expressed genes because they were increased more than 8-fold compared to the early ovarian developmental stage. The genes included mucin and yellow in addition to the well-known Vg and chorion protein genes. Mucins are high molecular and heavily glycosylated proteins with tandem repeats of identical or highly similar sequences rich in Ser, Thr, and Pro [19]. In insects, intestinal mucin is a major protein in the midgut peritrophic membrane, which facilitates the digestive process as well as protects the gut epithelium from microbial infections [20]. Salivary gland mucin might modulate the lubrication of insect mouthparts or defend plant attachment by inducing plant cell death through the formation of salivary sheaths [21,22]. A specific mucin protein is known to be associated with the formation of eggshells in Nilaparvapa lugens and Spodoptera exigua [23,24]. In our current study, the mucin gene highly expressed during F. occidentalis ovarian development suggests its function in chorion formation. Yellow and related major royal jelly protein (MRJP)-like proteins are widely found in insect genomes and these genes are classified into ten clades including Yellowb, -c, -d/e3, -e, -f, -g/g2, -h, -y, -x and MRJP-like protein [25]. In the Aedes albopictus mosquito, Yellow-g and Yellow-g2 are localized in the exochorion and outer endochorion, where they mediate darkening processes to physically strengthen the chorions [26]. Thus, the high expression of the yellow gene during F. occidentalis ovarian development suggests its function in chorion formation.
Transcriptome analysis also showed the expression profiles of genes associated with endocrine signals during F. occidentalis ovarian development. In insects, different endocrine signals are associated with ovarian development. JH is a sesquiterpenoid that mediates a status quo effect during the immature stage to prevent precocious metamorphosis. However, in adults, it stimulates ovarian development as a gonadotropin in various insects [27]. JH directly stimulates Vg biosynthesis in some insects and facilitates Vg uptake by growing oocytes by inducing follicular patency [28]. In mosquito females, 20-hydroxyecdysone acts as a gonadotropin [29]. ILPs are known to mediate ovarian development by stimulating oogonial proliferation to produce oocytes in the stem cell niche located in the germarium of the distal ovariole [30]. In F. occidentalis, JH and ecdysteroid play crucial roles in mediating metamorphosis. Krüppel homolog 1 (Kr-h1) and Broad (Br) are transcription factors leading to larval and pupal characteristics under JH and ecdysteroid hormones, respectively [31,32]. In F. occidentalis, Kr-h1 mRNA levels were high in the embryonic stage, remained at a moderate level in the larval and prepupal stages, and were low in the pupal stage. In contrast, Br mRNA levels were moderate in the embryonic stage and high at the larva-pupa transition stage. Except for Br expression in the embryonic stage, these two gene expression patterns followed the corresponding profiles of holometamorphic insects [33]. Furthermore, the adult specifier, E93, expression increased during immature development and its inhibition prevented adult metamorphosis [34]. However, little is known about JH and ecdysteroid mediation in oogenesis in F. occidentalis. Our current transcriptome analysis during ovarian development suggests that these endocrine signals play crucial roles in mediating F. occidentalis oogenesis based on their expression profiles. Increases in ILP and JHAMT expression at the mid ovarian developmental stage suggest their mediation of previtellogenesis by providing new oocytes from stem cells and vitellogenesis by stimulating Vg synthesis and uptake. Maintaining high levels of Shade expression suggest a high level of ecdysteroid during ovarian development, which may stimulate metabolic pathways, especially protein and nucleic acid biosynthesis, in addition to stimulating Vg synthesis with the cooperation of JH.
This study reports the comparative transcriptomes of F. occidentalis during different stages of ovarian development. Although the transcriptome analyses do not completely represent the protein expression profiles, they gave us valuable insights on the thrips reproduction. The results suggest an increase in metabolic pathways along with protein and nucleic acid biosynthesis. The high upregulation of egg proteins such as Vg, chorion protein, and sclerotizing agents during choriogenesis was also found. Finally, JH, ecdysteroid, and insulin signals may play crucial roles in mediating F. occidentalis oogenesis. A recent study showed that prostaglandin mediates oocyte devilment in early and late stages in addition to the endocrine signals [35]. This suggests the oogenesis of F. occidentalis would be a model system for an integrative analysis of endocrine signals mediating different reproductive processes of previtellogenesis, vitellogenesis, and choriogenesis.