Population density modulates insect progenitive plasticity through the regulation of dopamine biosynthesis

Insect fecundity is a quantitative phenotype strongly affected by genotypes and the environment. However, interactions between genotypes and environmental factors in modulating insect fecundity remain largely unknown. This study investigated the impact of population density on the fecundity of Nilaparvata lugens (brown planthopper; BPH) carrying homozygous high‐ (HFG) or low‐ (LFG) fecundity homozygous genotypes. Under low population densities, the fecundity and population growth rate of both genotypes showed similar increasing trends across generations, while the trends between HFG and LFG under high population densities were opposite. Through a combination of temporal analysis and weighted gene co‐expression network analyses on RNA‐seq data of HFG and LFG under low and high population densities in the 1st, 3rd, and 5th generations, we identified 2 gene modules that were associated with these density‐dependent progenitive phenotypes. Four pathways related to the neural system were simultaneously enriched by the 2 gene modules. Furthermore, Nlpale, which encodes a tyrosine hydroxylase, was identified as a key gene. The RNA interference of this gene and manipulation of its downstream product dopamine significantly affected the basic and density‐dependent progenitive phenotypes of BPH. These findings indicated that dopamine biosynthesis is the key regulatory factor that determines fecundity in response to density changes in different BPH genotypes. Thus, this study provides insights into the interaction of a typical environmental factor and insect genotype during the process of population regulation.


Introduction
Insects are the most diverse class of animals because of their excellent adaptability and high reproductive capacity (Grimaldi & Engel, 2005). The reproductive fitness of insects is a quantitative phenotype determined by both genotypes and environmental factors. Genotypic factors usually refer to genetic variants that change the function of fecundity-related proteins (Watt, 1992;Edward et al., 2014). Environmental factors associated with insect fecundity usually include predators (Rosenheim, 1998), temperature (Geister et al., 2008), xenobiotics (Ge et al., 2011), nutrition (Winkler et al., 2006), and insect population density (Stiling, 1988). However, how genotypes interact with the environment to modulate fecundity in a specific insect is largely unknown because of the lack of appropriate experimental materials.
Population density, which is a dynamic factor closely related to insect diet and nutrition availability, has significant impacts on the growth of insect populations (Hassell, 1975;Stiling, 1988). Population density is critical to the survival and development of larvae and consequently impacts the fecundity of emergent adults in various insect species. Aedes aegypti reared under higher larval density laid more eggs (Silva et al., 2020), whereas Chilo partellus moths at a lower pairing density had higher fecundity (Hari et al., 2008). An increased initial density results in prolonged nymphal development duration, lower emergence rate, and decreased longevity, and female fecundity in the rice planthopper (Mori & Nakasuji, 1991;Horgan et al., 2016). In addition, population density can influence many other processes, such as interactions between predators and preys, parasites and hosts, disease transmission, and competition (Arcese & Smith, 1988;Pacala & Hassell, 1991;Berg et al., 1992). Thus, density stress is ideal for studying interactions between insect genotypes and the environment during the regulation of population growth.
The population dynamics of insects can be influenced by individual variations (Hanski & Saccheri, 2006;Sun et al., 2015). However, whether the same population density has equal or diverse impacts on insect populations carrying different fitness-determined genotypes has remained unclear. To address this issue, we selected 2 genotypes of Nilaparvata lugens (brown planthopper; BPH), a serious rice insect pest. Based on our previous studies, we obtained 2 homozygous BPH populations (a highfecundity genotype [HFG] and a low-fecundity genotype [LFG] population) from the same initial population, with the ACE -862 and VgR -816 genotypes CCGG and AAAA, respectively (Sun et al., 2015;Liu et al., 2020). The ACE -862 (angiotensin-converting enzyme gene) and VgR -816 (vitellogenin receptor gene) genotypes were proved to be highly correlated with the fecundity of BPH (Sun et al., 2015), and the HFG and LFG populations homozygous in alleles ACE -862 and VgR -816 have steadily different fecundity phenotypes in normal environments . Here, we determined the population dynamics of the 2 BPH genotypes at different initial nymph densities through a density-dependent experiment. We then applied the RNA-seq and a series of bioinformatic analyses to unveil key factors involved in the density-dependent progenitive phenotypes in this insect. Overall, we aimed to uncover the potential mechanisms underlying the insect genotype and population density interactions during the modulation of insect fecundity.

Insect populations and density-dependent bioassays
The BPH populations carrying HFG and LFG genotypes were obtained in our laboratory (Zhai et al., 2013;Sun et al., 2015). Briefly, a heterozygous initial BPH population collected in the field was used to select the high-and low-fecundity individuals in a greenhouse for more than 3 years. Then, the expressions of fecundityrelated genes were determined and the SNPs in corresponding genes were detected in insects with high and low fecundity. Finally, high-fecundity individuals with homozygous ACE -862 and VgR -816 genotypes CCGG were selected as the HFG population, and low-fecundity individuals with homozygous ACE -862 and VgR -816 genotypes AAAA were selected as the LFG population. All insects were reared on susceptible Fenghuazhan rice in a greenhouse at 26 ± 2°C, 70%-90% relative humidity, and a light/dark cycle of 16/8 h.
Groups of 10, 50, 90, 130, and 170 newly (<12 h) hatched nymphs from the 2 BPH populations were reared in 45-d-old rice plants that were transferred into pots covered with insect-proof nets, respectively. The density unit is nymphs/plant. Each group comprised at least 10 replicates. The rice plants were replaced frequently to ensure a constant supply of fresh food throughout the study. The development and survival of nymphs were recorded daily until they all reached adulthood. Two-d-old female adults from each group were weighted then immediately placed in liquid nitrogen and stored at −80°C for subsequent transcriptome sequencing and detection. All sample collections were conducted at 15:00-17:00 (light cycle) to exclude the potential influences of circadian rhythm. Adults remaining from each replicate were randomly mated, and female adults were counted. When nymphs emerged, we counted the offspring from each replicate until no nymphs hatched for 3 consecutive days, and the average number of offspring in corresponding replicate was determined as the total number of offspring divided by the number of female adults. Newly emerged offspring from each group were randomly selected, reared at the same density as the next generation, and processed as described above (Fig. 1). Density-dependent bioassays were conducted for six generations. The population growth rate (P) was calculated as follows: P = S × R × F × E (S: nymph survival rate, R: female ratio, F: fecundity for per female, E: egg hatchability) (Yu et al., 2013).

RNA extraction and sequencing
We collected samples of F1, F3, and F5 insect generations with HFG and LFG reared under low (10 nymphs per group) and high density (170 nymphs per group). Total RNA was extracted from each 5 pooled insect individuals using TRIzol Reagent according to the manufacturer, and used for tag library preparation using the Gene Expression Sample Preparation Kit (Illumina) as described by the manufacturer's instructions. Total RNA integrity and quality were assessed using the RNA 6000 Nano LabChip kit and Agilent 2100 Bioanalyzer (Agilent Technologies). The poly(A)-containing mRNA was purified using magnetic oligo(dT) beads and Oligotex mRNA kit (Qiagen). The transcriptome sequencing libraries were prepared according to a modified method (Su et al., 2012). The mRNA strands were cleaved into short fragments using fragmentation buffer and reverse transcribed for cDNA synthesis.
Two replicates of each treatment were used for RNA sequencing. The cDNA library (∼300-bp long) was sequenced using an Illumina HiSeq TM 4000 platform (Illumina). Adapter sequences, empty reads, and low-quality reads (i.e., reads with a ratio of unknown sequences of "N" > 5%) were removed to obtain clean reads that were mapped to the N. lugens genome (accession no. GCF_000757685.1) using Hisat2 software (Kim et al., 2015). Read counts for all transcripts were extracted using HTSeq-count (Anders et al., 2015), then transferred to Fragments Per Kilobase Million (FPKM) using GenomicFeatures in R (Bioconductor; https://git. bioconductor.org/packages/GenomicFeatures).

Temporal analysis
We generated an FPKM matrix of all transcripts in all samples. The samples were then divided into 4 groups: HFG under low density (HFG-LD), HFG under high density (HFG-HD), LFG under low density (LFG-LD), and LFG under high density (HFG-HD). The temporal expression profiles from F1 to F5 of all genes were calculated for each group with the corresponding FPKM matrix using Short Time-series Expression Miner (STEM) software (Ernst & Bar-Joseph, 2006). The cluster number was set to 12 for each analysis.

Identification and annotation of candidate co-expression gene modules
Genes with FPKM < 0.3 in a sample were considered to have "too low expression" Bai et al., 2015). We filtered the FPKM matrix by removing transcripts with low expression in >50% of the tested samples. Co-expression gene network modules were inferred from the filtered matrix using weighted gene coexpression network analysis (WGCNA) in R (Langfelder & Horvath, 2008). The parameters were power β = 6, minimum module size = 50, and ME miss thread = 0.25.
Candidate co-expression gene modules were selected based on the temporally analyzed gene list. The genes in the candidate modules were annotated to the KEGG pathway database, and pathway enrichment was analyzed using Fisher exact tests with Benjamini and Hochberg (BH) adjustment and the KEGG annotation of all genes in the BPH as background. Adjusted P < 0.05 was considered as significantly enriched.

Differential expression analysis and quantitative real-time PCR validation
Read counts generated by HTseq-count were normalized using the DESeq2 package in R (Love et al., 2014). The differentially expressed genes (DEGs) between HFG-HD and HFG-LD at F5 were estimated by DESeq2 according to the threshold of |log2 ratio| > 1 and adjusted P < 0.05.
We measured the expression of 20 randomly selected DEGs to support the analysis. Total RNA (1 μg) was reverse transcribed into first-strand cDNA using a PrimeScript TM RT reagent kit (Takara Bio). Quantitative real-time (qRT)-PCR was performed in 10-μL volumes containing 1 μL cDNA, 0.3 μL each of 10 μmol/L forward and reverse primers, and 5 μL SYBR ® FAST Universal qPCR mix (Kapa Biosystems) using the LightCy-cler 480 system (Roche Diagnostics). The cycling conditions were as follows: 5 min at 95°C, followed by 45 cycles at 95°C for 10 s, at 60°C for 20 s, and at 72°C for 20 s. The qPCR experiments were performed for each sample using 3 biological and 3 technical replicates. The expression of selected genes was normalized to that of N. lugens β-actin (Chen et al., 2013). Differential gene expression was calculated using the 2 − Ct method (Livak & Schmittgen, 2001).

RNA interference of Nlpale
The gene Nlpale that encodes tyrosine hydroxylase (TH), which is the rate-limiting enzyme in the dopamine production metabolic pathway, was selected as the key regulatory factor, and its function was validated by RNA interference (RNAi). We synthesized double-stranded RNA (dsRNA) by PCR amplification of a 466-bp fragment (Nlpale) using Nlpale cDNAs and the T7 Ribo-MAX Express RNAi System (Promega) as described by the manufacturer. All primers contained the T7 promoter sequence at each end (Table S1). The control was the dsRNA of Aequorea victoria green fluorescent protein (GFP, accession no. ACY56286) (Qiu et al., 2016;Pang et al., 2017). The dsRNA concentrations were quantified using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). The quality and size of the dsRNA were further verified by 1% agarose gel electrophoresis.
We injected the mesothorax of 1-d-old adult females from the HFG and LFG groups with approximately 250 ng dsRNA or control dsGFP as previously described (Chen et al., 2013). The RNAi efficiency was assessed 24 h later in 3 randomly collected females by qRT-PCR. The injected females were mated with male adults at a ratio of 1 : 1, and released to potted rice at the densities of 10 and 170 individuals, respectively. The methods of qRT-PCR and fecundity bioassay for the injected insects were as described above. The total number of offspring in each treatment was counted.

Manipulation of dopamine in the BPH
The dopamine content in BPH was determined using Insect Dopamine ELASA kit (Boshen Biotechnology) following the manufacturer's instruction, which states that the kit does not cross-react with other soluble structural analogs (http://www.boswio.com/cp_view.asp? id=30921). Briefly, 5 1-d-old female adults comprised 1 repeat, and each treatment had 10 repeats. Homogenate, standard dopamine, and horse radish peroxidase (HRP)labeled dopamine antibody were incubated at 37°C with micropores that were precoated with the dopamine antibody, washed, and then stained with 3ʹ,3ʹ,5ʹ,5ʹ-tetramethyl benzidine (TMB) substrate for color rendering. The optical density (OD; absorbance) was measured at a wavelength of 450 nm to determine the dopamine content in the samples.
Dopamine hydrochloride (100 nL of 80 μg/μL; Sigma-Aldrich) and chlorpromazine (CPZ; 100 nL of 20 μg/μL), a dopamine receptor antagonist, were dissolved in Ringer's solution, respectively. Afterward, they were injected into the mesothorax of the BPH as previously described (Ma et al., 2011;Chen et al., 2013). The insects injected with Ringer's solution were used as control. The qRT-PCR and fecundity bioassay after injection was conducted as described above.

Measurements of 20-hydroxyecdysone content
The biosynthesis of 20-hydroxyecdysone (20E) in insects can be regulated by dopamine (Rauschenbach et al., 2007); therefore the content of 20E in BPH was determined using Insect Ecdysone ELISA Assay kit (Boshen Biotechnology) according to the manufacturer's instruction. Briefly, 5 1-d-old female adults comprised 1 repeat, and each treatment comprised 10 repeats. The homogenate, standard 20E, and HRP-labeled ecdysone antibody were sequentially added to the micropores that were precoated with the ecdysone antibody and incubated at 37°C. The incubations were washed, then 3ʹ,3ʹ,5ʹ,5ʹ-tetramethyl benzidine (TMB) substrate was added for color rendering. The absorbance (OD) was measured at a wavelength of 450 nm to determine the ecdysone content in the samples.

Statistical analysis
The impact of density, generation, genotype, and their combinations on insect fecundity and population growth rate was analyzed by multi-factor (multi-way) ANOVA and Tukey's test using SPSS 18.0 (SPSS). Differences in pairwise comparisons of gene expression and BPH offspring numbers were analyzed by Student's t-test using SPSS 18.0. Statistical differences were considered significant at P < 0.05. Differences among multiple comparisons were assessed using one-way ANOVA, followed by Duncan's multiple range test (P < 0.05). Differences in density-dependent responses between insects injected with dsGFP and dsNlpale were assessed using two-way ANOVA, followed by Tukey's multiple comparison test.

Population density modulates BPH fecundity
Population growth trends and dynamics between HFG and LFG at densities of 10, 50, 90, 130, and 170 nymphs were measured for 6 continuous generations. The density, generation, genotype, and combinations of these factors significantly impacted the number of offspring per female (Table S2). Both genotypes had lower fecundity at higher individual densities ( Fig. 2A), which might result from phenotypic plasticity in response to density changes. Additionally, some transgenerational effect was evident in most treatments. Under densities of 90, 130, and 170 nymphs, the number of HFG offspring tended to decrease across generations, but displayed the opposite at low densities (10 and 50 nymphs). The LFG showed higher fecundity as the generations increased at all tested densities except the highest one (170 nymphs). Most adjacent generations in the LFG did not significantly differ at density of 170 nymphs.
Differential transgenerational plasticity in response to the population density between these 2 genotypes of BPH was also evident from the results of the population growth rates (Table S2). The population growth rates of the HFG group significantly increased and decreased under low (10 and 50 nymphs) and high (90, 130, and 170 nymphs) population densities, respectively (Fig. 2B). The population growth rate of the LFG group similarly changed under low densities, but such changes were not obvious at high densities (Fig. 2B).

Temporal expression trends revealed by RNA-seq
We explored key factors responsible for the different responses to density changes between HFP and LFP. Insects reared at densities of 10 (low density, LD) and 170 (high density, HD) nymphs were selected, and samples at F1, F3, and F5 were assessed using RNA-seq, which produced an average of 33 320 239 clean paired-end reads per sample with a Q30 value >92% (Table S3). The average mapping rate of all samples to the reference genome of BPH was 54.43% (Table S3), and an FPKM matrix of 24 306 transcripts in all samples was generated based on mapping results (Table S4). After filtering out all transcripts that had very low expression levels, 15 147 transcripts were used for subsequent analysis.
We independently applied temporal analysis with the filtered FPKM matrix to the HFG-LD, HFG-HD, LFG-LD, and LFG-HD groups. According to the fecundity phenotype, the key factors in both HFG and LFG should show a linear correlation with fecundity under low density. Therefore, transcripts with continuous increasing or decreasing trends alongside generations in the HFG-LD and LFG-LD groups were selected (Fig. 3A). The key factors should have differential expression profiles between high and low densities in both genotypes, and the expression patterns of key factors between HFG and LFG at high density should also be inconsistent. By applying these filter parameters, we identified 31 candidate genes (Fig. 3B). The expression of these genes demonstrated a significantly high correlation between HFG-LD and LFG-LD, but not between HFG-LD and HFG-HD, LFG-LD and LFG-HD, and HFG-HD and LFG-HD (Fig. 3C).

Identification of the key regulatory factor for density-dependent progenitive phenotypes
We conducted WGCNA of all filtered transcripts to determine the key factors for density-dependent progenitive phenotypes. The analysis resulted in 16 distinct gene modules (Fig. 4A, Table S5), and 31 candidate genes were distributed into 8 modules. Turquoise and dark red  modules accounted for most candidate genes (21/31, Table S6). The KEGG pathway enrichment of the genes in these 2 modules was analyzed. The genes in the turquoise module were significantly enriched in 69 pathways, while the number of enriched pathways for the genes in the dark red module was 25 (Table S7). Four pathways were enriched in both modules, which included endocrine and other factor-regulated calcium reabsorption, axon regeneration, dopaminergic synapse, and regulation of actin cytoskeleton (Fig. 4B). These findings suggest that these pathways might play important roles in regulating fecundity between different population densities.
According to the phenotypic differences, the key factors should be significantly differentially expressed between HFG-HD and HFG-LD at F5. Differential expression analysis of these datasets identified 464 DEGs in HFG-HD compared with HFG-LD, including 311 upregulated and 153 downregulated genes (Fig. 5A, Table S8). qRT-PCR analysis of the randomly selected DEGs validated the reliability of the RNA-seq results (Fig. S1). Among the DEGs, only 6 were involved in the 4 pathways mentioned above. Importantly, only 1 of the 6 DEGs was among the candidate transcripts obtained from the temporal analysis. This candidate gene, LOC111049383 (Nlpale), encodes tyrosine hydroxylase, which is the key rate-limiting enzyme for dopamine synthesis in organisms (Morgan et al., 1996;Ma et al., 2011). The decreased expression of Nlpale at a low population density was irrelevant to the fecundity genotypes of BPH (Fig. 5B). However, the expression of Nlpale was continuously elevated in HFG under high population density with increasing generations, but a similar trend was not observed in LFG under the same population density. Overall, the changes in Nlpale expression were negatively correlated with the phenotypic changes in fecundity in both BPH genotypes.

We investigated the basic expression level of gene
Nlpale and the dopamine content of the BPH genotypes with HFG and LFG of BPH under low population fecundity (10 nymphs). The qRT-PCR results showed that Nlpale mRNA expression was more abundant in LFG than in HFG (Fig. 6A). Accordingly, the LFG had higher dopamine content than the HFG (Fig. 6B). These results corresponded to the transcriptomic analysis that identified a negative correlation between Nlpale gene expression and BPH fecundity. This relationship was further confirmed by RNAi of Nlpale. The expression of Nlpale in both genotypes significantly decreased from 24 to 72 h after dsRNA injection, compared to the GFP control (Fig. 6C). The expression of NlVg, an important molecular marker of insect fecundity, was significantly increased in both genotypes at 48 and 72 h after RNAi of Nlpale (Fig. 6D). Correspondingly, the number of HFG and LFG offspring was significantly elevated after RNAi of Nlpale (Fig. 6E). In addition, when female adults were supplemented with extra dopamine hydrochloride by microinjection, the number of their offspring significantly decreased. The injection of CPZ, an inhibitor of dopamine receptors, resulted in more offspring in the corresponding adult females (Fig. 6F). The negative correlation between dopamine content and fecundity did not differ between insects of HFG and LFG, indicating the intrinsic regulatory role of dopamine biosynthesis in BPH fecundity.
We further investigated whether Nlpale and dopamine were key factors in determining the difference in response to population density changes between HFG and LFG insects. A significant negative correlation between the dopamine content and fecundity was observed in insects with both genotypes under different population densities across 6 generations (R = −0.88 and P = 1.52 × 10 −6 for HFG, and R = −0.90 and P = 3.45 × 10 −7 for LFG; Fig. S1). We then compared the fecundity of insects of both genotypes that had RNAi of Nlpale with that of the GFP control. In terms of control traits, the average fecundity of both genotypes significantly decreased as the population density increased (Fig. 7A). When Nlpale was knocked down, the density-dependent change in LFG fecundity was significantly weakened (P < 0.001, Fig. 7B, Table S9). However, the density effect was absent only in the HFG under high population densities (130 and 170 nymphs). We subsequently modulated the dopamine content in the HFG and LFG insects. The injected dopamine hydrochloride and CPZ had opposing effects on the overall fecundity of BPH, but the density effect on fecundity was weakened in both groups ( Fig. 7C and 7D). These results suggested that disrupting of dopamine biosynthesis impacted on the density-dependent progenitive plasticity of BPH significantly.
Dopamine can regulate the biosynthesis of various insect hormones (Gruntenko et al., 2005a(Gruntenko et al., , 2005bRauschenbach et al., 2007). Thus, we assessed the relationship between 20E and dopamine at different densities. The 20E content was negatively correlated with dopamine content but was positively correlated with BPH fecundity under different densities across different generations (Fig. S2). Additionally, the RNAi of Nlpale led to an increase in 20E in both genotypes (Fig. S3), indicating an intrinsic regulatory relationship between dopamine and 20E in this insect species.   . 6 Role of Nlpale mRNA level and dopamine content in BPH fecundity. The Nlpale mRNA levels (A) and dopamine contents (B) in 2 BPH genotypes are shown. The mRNA levels of Nlpale (C) and NlVg (D) in the 2 BPH genotypes were measured at 24, 48, and 72 h after RNA interference of Nlpale. Insects injected with dsGFP were used as control. (E) Fecundity of HFG and LFG after injection with Nlpale dsRNA. (F) The fecundity of HFG and LFG after injection with dopamine receptor antagonist chlorpromazine, dopamine hydrochloride, or water control. All data are shown as means ± SE (n = 3 for panel A, C, and D; n = 10 for panel B, E, and F). *P < 0.05; ** P < 0.01 (Student's t-test). Values sharing different letters are significantly different at P < 0.05 (one-way ANOVA and Duncan's multiple range tests).

Discussion
In the present study, we validated the progenitive plasticity in BPH populations with two different fecundityrelated genotypes. Under continuous low-density conditions, the BPH insect population expanded with elevated individual fecundity, independent of the genotype. The individual fecundity population growth rate of the HFG insects was continuously reduced under continuous high density, but this decreasing trend was not found in LFG (Fig. 2).
Our findings can be largely explained by the traditional theory of density-dependent population regulation. Nicholson (1933) proposed that a population could adjust density according to its own natural and environmental conditions. All species would balance their population density with the resources necessary to live. Accordingly, a population density that exceeds the balance can be restricted by a series of morphological, physiological, behavioral, and genetic factors to prevent infinite growth in the population (Chitty, 1960). In turn, the factors suppressing population growth will be weakened if the population density is lower than the balance (Murdoch, 1994;Turchin, 1999). The changes in the population growth rate of the HFG-HD, LFG-LD, and HFG-LD groups were basically consistent with this hypothesis. The phenotypic plasticity of these insects might have caused these changes in response to crowding or dispersing, and the transgenerational effect was probably present in these samples since trends across subsequent generations were significant (Fig. 2). The fecundity of the LFG group increased under a high density, but the overall population growth rate was stable across different generations. Therefore, it is possible that thresholds for population balance under different densities might exist in this species. At high density, the threshold was advantageous for the LFG group. Thus, the LFG under higher density rapidly reached equilibrium, as indicated by the stabilized population growth rates.
RNA-seq and a series of bioinformatic analyses were utilized to reveal the determinant for the phenotypic differences between HFG and LFG under high population density. The analyzed results centralized at 4 pathways ( Fig. 4B), all of which were closely associated with the neural system. The regulation of actin cytoskeleton is essential during axonogenesis, because neuroblast growth must undergo engorgement, which involves actin depolymerization (Bradke & Dotti, 1999) followed by the penetration of the microtubules into the central and peripheral domains (Kunda et al., 2001). Endocrine and other factorregulated calcium reabsorption is a pathway associated with calcium regulation in the organism, while calcium ion is involved in synaptic signaling (vesicle docking to the presynaptic membrane via synaptogamin) (Li et al., 1995;Tucker et al., 2004). Importantly, the pathways of axon regeneration and dopaminergic synapse are components of the neural system (Neve et al., 2004;Hammarlund & Jin, 2014).
Our subsequent investigation revealed that dopamine was the key factor in modulating BPH fecundity (Figs. 5 and 6). Dopamine is a catecholamine neurotransmitter that is synthesized in neuronal and nonneuronal tissues. In insects, tyrosine was transferred to L-DOPA (L-3,4dihydroxyphenylalanine) by tyrosine hydroxylase (TH), and L-DOPA was enzymatically catalyzed into dopamine by DOPA decarboxylase (DDC) (Noh et al., 2016). Then, the dopamine was transported from cytoplasm to synaptic vesicles via the assistance of vesicular monoamine transporter (VMAT) (Yaffe et al., 2018). After pumping out from presynaptic terminal to synaptic cleft through exocytosis, the dopamine may specifically bind to the dopamine receptors (DARs) located in postsynaptic cells to active the signal transmission, or be reuptaken into presynaptic cells via the dopamine transporter (DAT) to prevent overreaction . Here, the expression of genes encoding TH, DDC, VMAT, and DAT were determined in BPH (Fig. 5B), and the TH coding gene (Nlpale) showed negative correlation with the fecundity changes in both LFG and HFG insects. Since TH is the rate-limiting enzyme for dopamine synthesis, the dopamine content in BPH was also observed to be negatively correlated with the fecundity phenotypes (Fig. S2).
Notably, in many insect species, dopamine can promote ovarian development as a gonadotropin and thus have positive effects on fecundity (Neckameyer, 1998;Boulay et al., 2001), which were opposite to our findings for BPH. Actually, some studies have showed the gene expression of biogenic amine receptors in the ovaries of insects (Vergoz et al., 2012;Okada et al., 2015), which revealed that biogenic amines are likely to directly influence insect fecundity. However, different biogenic amines may pose different impacts on insect fecundity. For example, the reproduction of the honey bee could be promoted by dopamine but inhibited by octopamine (Sasaki & Nagao, 2001;Salomon et al., 2012). The tyramine can motivate the ovarian development of honey bees, while octopamine has no such effect (Salomon et al., 2012). Further, the same biogenic amine may also have different effects on different insect species. For instance, serotonin plays a key role in the mating of Drosophila, but has no such function in bruchid (Yamane & Miyatake, 2010;Becnel et al., 2011). Therefore, the negative effect on fecundity might be a specific function of dopamine in BPH. A hypothesis for this observation might fall in the regulation of dopamine on insect hormones. It has been reported that changes in dopamine content can cause changes in other hormones such as 20E and juvenile hormone, indicating its potential role in regulating these hormones (Gruntenko et al., 2005a(Gruntenko et al., , 2005bRauschenbach et al., 2007). Our data also suggested that dopamine negatively regulates NlVg expression and 20E content (Fig. 6D, Fig. S2). In BPH, the 20E content in adult female BPH is significantly associated with their progenitive capacity (Yu et al., 2014;Zhou et al., 2020). Thus, the way that dopamine biosynthesis regulates BPH fecundity might proceed via the 20E signaling pathway, and then act on the oogenesis process. Besides, dopamine and 20E might antagonistically bind to the dopamine/ecdysteroid Receptor (DopEcR) to activate the phosphoinositol 3 kinase pathway or the MAP kinase pathway (Srivastava et al., 2005;Evans et al., 2014). The involvement of DopEcR in the reproductive behavior has been reported (Abrieux et al., 2013;Crickmore & Vosshall, 2013;Abrieux et al., 2014), thus in our case, the dopamine and 20E might also oppositely regulate the reproduction of BPH via activation of different pathways in the downstream of DopEcR. Further research is needed to investigate the regulatory approach of dopamine and its relevance with 20E in the fecundity of BPH.
Furthermore, the present study showed that dopamine contributes to the density-dependent progenitive plasticity of BPH. The RNAi of the TH coding gene (Nlpale) and dopamine content modulation in insects significantly disturbed the density effect on BPH fecundity (Fig. 7). In Locusta migratoria, the transitions between solitary and gregarious phases in response to population density changes are closely related to dopamine biosynthesis (Ma et al., 2011;Zhang et al., 2020), and gregarious locusts tended to contain higher content of dopamine (Rogers et al., 2004). A similar finding was also observed in crowded crickets, which showed significantly higher amounts of dopamine than the isolated ones (Iba et al., 1995). In the present study, both LFG and HFG insects under higher population density showed relatively higher dopamine content than those insects under lower population density (Fig. S2A), which was to some extent in accordance with the previous reports. We assumed that higher population densities would release stronger signals to the brains of insects, which may upregulate the dopamine synthesis pathway and cause corresponding changes in female fecundity. The transgenerational effect was determined by the population genotype, which is associated with proximity to population balance under a specific density.
In summary, we showed that dopamine biosynthesis is a key regulatory factor for density-dependent progenitive plasticity in BPH. Our findings provide insights into the interactions between genotypes and environmental factors regulating insect population dynamics, and facilitates the understanding of the mechanisms underlying population regulation.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.  Table S1 The primers used in this study.

Table S2
Analysis of ecology fitness of high-and lowfecundity genotypes of the BPH under different population densities across 6 generations.

Table S3
The statistics of all RNA-seq datasets. Table S4 The FPKM matrix of all transcripts in the BPH. Table S5 List of the gene modules obtained by WGCNA. Table S6 The module information of candidate genes.

Table S7
Results of the KEGG pathway enrichment analysis for the turquoise and dark red gene modules. Table S8 The information of differentially expressed genes in HFG-HD compared to HFG-LD at F5. Table S9 ANOVA of the effect of dopamine biosynthesis in the density-dependent effects on the fecundity of BPH.