Transcriptional Analysis of Coccidioides immitis Mycelia and Spherules by RNA Sequencing

Coccidioides immitis and C. posadasii are dimorphic fungi that transform from mycelia with internal arthroconidia in the soil to a tissue form known as a spherule in mammals. This process can be recapitulated in vitro by increasing the temperature, CO2 and changing other culture conditions. In this study, we have analyzed changes in gene expression in mycelia and young and mature spherules. Genes that were highly upregulated in young spherules include a spherule surface protein and iron and copper membrane transporters. Genes that are unique to Coccidioides spp. are also overrepresented in this group, suggesting that they may be important for spherule differentiation. Enriched GO terms in young spherule upregulated genes include oxidation-reduction, response to stress and membrane proteins. Downregulated genes are enriched for transcription factors, especially helix–loop–helix and C2H2 type zinc finger domain-containing proteins, which is consistent with the dramatic change in transcriptional profile. Almost all genes that are upregulated in young spherules remain upregulated in mature spherules, but a small number of genes are differentially expressed in those two stages of spherule development. Mature spherules express more Hsp31 and amylase and less tyrosinase than young spherules. Some expression of transposons was detected and most of the differentially expressed transposons were upregulated in spherules.


Introduction
Coccidioides immitis and C. posadasii are primary pathogenic fungi that are endemic in the desert regions of the Western United States, Mexico, and Central and South America [1]. They cause pulmonary infections that range from asymptomatic infections to severe pneumonias and can disseminate beyond the lung [2,3]. The organisms grows as a mold in the soil and produce asexual spores, termed arthroconidia, within the mycelium. When the soil is disturbed, the mycelia can rupture and the arthroconidia are released. If inhaled by a susceptible host, the arthroconidium differentiates into a very different and unique form that is known as a spherule. In tissue, the spherule enlarges and can form many reproductive endospores. Mature spherules rupture and release endospores, which can then differentiate into the next generation of spherules. This form, the spherule, is the disease-associated form of the organism. If the disease is self-limited, the reproduction of spherules is limited, but if the disease is severe, then spherules continue to grow, rupture and give rise to new spherules, while eliciting inflammatory and immune responses.
This process goes on in many mammalian hosts, including desert rodents, which may be important in the ecology of the fungus [4].
The morphological transition between mycelia and spherule forms of C. immitis is dependent on sensing the host environment, and this transition can be recapitulated in the laboratory by changing the temperature and other growth conditions. C. immitis grow in saprobic form at 22-30 • C; culturing at 37-42 • C with 5-20% CO 2 in Converse media is used to convert arthroconidia to spherules [5,6]. Utilizing these conditions, whole genome-level transcriptional profiling studies of saprobic and parasitic forms have been performed [7]. That work reported that about 1300 genes are upregulated in mycelia, about 1900 genes are upregulated in spherules, and expression of known virulence genes is upregulated in the spherules of both C. immitis and C. posadasii, linking morphology to the virulence traits. For example, the spherule outer wall glycoprotein, which is the outermost layer on the spherule, is expressed only in spherules [8,9]. Another study compared mycelia to young spherules and mature spherules to mycelia and to each other using microarray technology [10]. In the current study, we analyze the RNA obtained from frozen samples of spherules and mycelia obtained from the previous microarray study by strand-specific RNA sequencing. Our results show that there are more than 1500 genes that are differentially regulated in spherules and mycelial phases of C. immitis. We have compared the up-and downregulated genes to previous studies and, where possible, analyzed the function of differentially expressed genes.

Culture Conditions
C. immitis RS strain was grown as mycelia or spherules as previously described [10]. To grow the mycelia, 2 × 10 6 arthroconidia/mL were incubated in 250 mL flat-bottom Erlenmeyer flasks (Corning, Corning, NY, USA) in 50 mL of GYE media. The flasks were cultured in a 30 • C incubator without shaking for 5 days. To grow the spherules, arthroconidia were washed 2 times in modified Converse media. The spores were inoculated at 4 × 10 6 arthroconidia/mL into a 250 mL baffled Erlenmeyer flask containing 50 mL of modified Converse media. The flasks were set up and grown on a shaker at 160 rpm, in 14% CO 2 at 42 • C. Four flasks were harvested 2 days after inoculation and the remaining four flasks after 8 days. Fresh Converse media was not added. The spherules did not rupture and release endospores within that time in this culture system.

RNA Extraction, Purification and Sequencing
The mycelia and spherule samples were stored in QIAzol (Qiagen) at −70 • C for seven years. Samples were added to 2 mL ZR BashingBead lysis tubes with 0.5 mm beads (Zymo Research, Irvine, CA, USA) and the tubes were arranged in a pre-cooled Tissuelyzer II adapter (Qiagen, Germantown, MD, USA) and disrupted by shaking at 50 Hz for 25 min. Total RNA was purified from mycelia and spherule samples (2 replicates/condition) using chloroform extraction and isopropanol precipitation and quantified using a Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, CA, USA). For RNA-Seq, strand-specific, paired-end libraries were prepared from total RNA by ribosomal depletion using the Yeast Ribo-Zero rRNA Removal Kit (Illumina, San Diego, CA, USA) and then using the TruSeq Stranded total RNA-Seq kit (Illumina) according to manufacturer's instructions. Then, 100 bases were sequenced from both ends using an Novaseq 6000 instrument according to the manufacturer's instructions (Illumina, San Diego, CA, USA). A total of 10 million reads per sample were acquired.

Analysis of RNA-Seq Data
Duplicate sets of paired strand-specific reads were mapped to predicted mRNA sequences and quantified using Kallisto [11]. Read count tables from the Kallisto output were analyzed by DESeq2 in R. The results were filtered for a Benjamini-Hochberg adjusted p value < 0.05; upregulation and downregulation for the gene tables are defined as more than 2 log2 or less than −2 log2 spherules/mycelia fold change (FC). This cut off is chosen to minimize the false positive identification of differentially regulated genes. All FC in Tables  are expressed as log 2 values; FC changes in the text are expressed as arithmetic values.

Functional Analysis of Genes
All differentially expressed genes were evaluated using data from FungiDB version 51 (https://fungidb.org/fungidb/app, accessed on 4 January 2021) [12]. Orthology searches were done using the orthology tool at FungiDB. The FungiDB tool compares two sets of proteins by reciprocal BLASTP and computes the percent match length with a threshold for orthology, where a blast match e value is < 10 −5 with a percent match length of ≥ 50%. Finally, the FungiDB tool obtains paralogs, orthologs and ortholog groups using OrthoMCL Pairs release 6.4 (https://orthomcl.org, accessed on 4 January 2021) [13]. C. immitis genes without a BLASTP match with an e value < 10 −8 to all fungal species other than Coccidioides spp. and less than ten orthologs were considered genus-specific. Multiple alignment of transcription factors was done with Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/, accessed on 4 January 2021) [14]. GO enrichment analysis was done using FungiDB (https://fungidb.org/fungidb/app, accessed on 4 January 2021) and FungiFun (https://elbe.hki-jena.de/fungifun/, accessed on 4 January 2021) [15]. GO term enrichment using the FungiFun tool was deemed significant if the Benjamini-Hochberg adjusted p value was < 0.05). GO enrichment analysis using the FungiDB tool was considered significant if the p value was < 0.05. Enrichment was carried determined using a Fisher's Exact test with the background defined as all genes from the organism being queried. p-values corrected for multiple testing are provided using both the Benjamini -Hochberg false discovery rate method and the Bonferroni method. GO term enrichment are presented in Supplemental Table S2.
The FungiDB pathway tool was used to generate data in Supplemental Table S3. The pathway enrichment analysis of differentially expressed genes was conducted using the tool at FungiDB, which searches for annotation at MetaCyC (https://metacyc.org/, accessed on 4 January 2021) [16] and KEGG (https://www.genome.jp/kegg/, accessed on 4 January 2021) [17].

Results and Discussion
The RNA-Seq gene expression from the three conditions tested were compared to each other and the results are shown in Figure 1 and Supplemental Table S1. Comparing young and mature spherules to mycelia, at least 8% of genes were upregulated and 16-18% were downregulated (FC > 2 log2 or < −2 log2 , adjusted p value less than 0.05). The number of differentially expressed genes comparing young spherules to mature spherules was much smaller (1.4-2%).

Differential Gene Expression in Young Spherules vs. Mycelia
The 20 most differentially expressed genes in young spherules are shown in Table 1. The most highly upregulated gene is the spherule outer wall glycoprotein (SOWgp). This protein is expressed in very large amounts on the external surface of spherules but not mycelia and is involved in pathogenicity [8,9]. Two other spherule antigenic proteins, expression-library immunization protein-1 and parasitic-phase-specific protein PSP-1, [18,19], are also very highly upregulated, as are three transporters, including a copper transporter. There are only three genes that are predicted to be copper transporters in C. immitis, so this is a substantial increase in transporter expression, which suggests that spherules may have a much higher need for copper than mycelia do. Recent experiments in Paracoccidioides brasiliensis showed that copper deprivation has a major effect on fungal metabolism, but similar experiments have not been done in Coccidioides spp. [20]. In addition, a comparison of gene expression in C. posadasii spherules and an engineered chitinase deletion that does not endosporulate found that a major difference was the upregulation of iron and copper uptake genes, which is further evidence for the importance of iron and copper uptake for spherule and endospore formation [21]. Three of the twenty very highly expressed genes are genus specific. Although this group of genes was highly upregulated, this does not imply that they are functionally related.

Differential Gene Expression in Young Spherules vs. Mycelia
The 20 most differentially expressed genes in young spherules are shown in Table 1.   Examining all of the 408 upregulated genes in young spherules, there is no correlation between gene length and FC. Those that are genus specific are more common in differentially expressed genes than in the whole genome. A total of 24% of all C. immitis genes are genus specific compared to 36% of differentially expressed genes. This suggests that some genus-specific genes may be important for spherule differentiation ( Table 2).  Out of all differentially expressed genus-specific genes, 11% have signal peptides (compared to 5% of all C. immitis genes), 7% have a transmembrane domain (compared to 16% of all genes) and 1.6% have a predicted Pfam domain (compared to 40% of all genes). Differentially expressed genus-specific genes were much shorter than conserved genes ( Figure 2). copper uptake for spherule and endospore formation [21]. Three of the twenty very highly expressed genes are genus specific. Although this group of genes was highly upregulated, this does not imply that they are functionally related. Examining all of the 408 upregulated genes in young spherules, there is no correlation between gene length and FC. Those that are genus specific are more common in differentially expressed genes than in the whole genome. A total of 24% of all C. immitis genes are genus specific compared to 36% of differentially expressed genes. This suggests that some genus-specific genes may be important for spherule differentiation ( Table 2).  Table 2: Comparison of the number of differentially expressed (young spherule/mycelia) genus-specific genes to all genus-specific genes. a) p value < 0.05 (Chi-square test) comparing the fraction of genus-specific genes in up-or downregulated genes to the fraction in all genes.
Out of all differentially expressed genus-specific genes, 11% have signal peptides (compared to 5% of all C. immitis genes), 7% have a transmembrane domain (compared to 16% of all genes) and 1.6% have a predicted Pfam domain (compared to 40% of all genes). Differentially expressed genus-specific genes were much shorter than conserved genes (Figure 2). . Figure 2. Length of genus-specific and conserved genes. Legend Figure 2: Lengths (Y-axis) of conserved and differentially expressed genus-specific proteins.
Enriched GO terms in the upregulated genes (young spherule/mycelia) included oxidation-reduction processes, integral membrane components, transmembrane transport and response to stress (Supplemental Table S2). One of the oxidation/reduction genes is superoxide dismutase, which has been shown to be important for pathogenicity in Histoplasma capsulatum [22]. Changes in expression of oxidation/reduction genes are to be expected since the atmosphere for mycelial culture is air with 0.4% CO2 compared to 14% CO2 in spherule culture conditions. Several metal transporters, including copper and iron transporters are also upregulated. A cluster of genes, including an iron siderophore, are Enriched GO terms in the upregulated genes (young spherule/mycelia) included oxidation-reduction processes, integral membrane components, transmembrane transport and response to stress (Supplemental Table S2). One of the oxidation/reduction genes is superoxide dismutase, which has been shown to be important for pathogenicity in Histoplasma capsulatum [22]. Changes in expression of oxidation/reduction genes are to be expected since the atmosphere for mycelial culture is air with 0.4% CO 2 compared to 14% CO 2 in spherule culture conditions. Several metal transporters, including copper and iron transporters are also upregulated. A cluster of genes, including an iron siderophore, are induced by iron-deprivation and are required for pathogenicity in H. capsulatum [23]. Five of the six iron-related genes (including the siderophore biosynthesis cluster) have homologs in C. immitis that are upregulated in young spherules. These genes are tightly clustered (<25 kb) on contig 1 and have the upstream regulatory sites for the GATA transcription factor Sre1 that is identified in identified in H. capsulatum. The Blastomyces dermatitidis homolog of Sre1 has been knocked out and the resulting mutant is unable to differentiate from mold to yeast [24]. However, the C. immitis homolog of Sre1, SreP, is downregulated in young spherules. Taken together, these results show that the organization of iron-related genes in C. immitis and H. capsulatum is very similar and that most of the C. immitis iron acquisition genes are upregulated in both of these dimorphic fungi. Although this may imply the upregulation of spherules in vivo, in this experiment, mycelia and spherules were grown in different media, which may have played a role in the upregulation of iron-regulated genes in this experiment.
Cellular component GO terms associated with cell membranes are dramatically overrepresented in upregulated genes (young spherule/mycelia) (p < 10 −4 ). In addition, 22% of the up-regulated genes have at least one predicted transmembrane domain, compared to 17% of total genes (p < 0.05). However, there is no difference in the proportion of up-regulated genes compared to all genes with a predicted signal peptide.
A pathway enrichment analysis was also conducted (Supplemental Table S3). An analysis of upregulated genes (young spherule/mycelia) found that 24 pathways were enriched with a p value of less 0.01 but none had a Benjamini-Hochberg adjusted p value of less than 0.05. The most highly enriched pathways in upregulated genes were pentose and glucuronate interconversion (see Supplemental Table S3). Other enriched pathways included the degradation of proponate, furfural and benzoate. Two pathways for polyketide synthase were also enriched.
Many genes that have previously been identified to be upregulated in spherules were also found to be upregulated in this study. In addition to SOWgp [8,9] and the parasitic phase specific gene [19], expression of the urease and the ureidoglycolate hydrolase genes was also upregulated [25,26]. One gene that was found to be upregulated in the yeast (or spherule) phase of all dimorphic fungi is 4-hydroxyphenylpyruvate dioxygenase (4-HPPD or HpdA) [27]. This gene is involved in tyrosine catabolism, which plays a role in the synthesis of melanin [28]. Chemical inhibition of 4-HPPD blocks the formation of yeast in P. brasiliensis and deletion of the gene blocks' differentiation to yeast [29]. There are two genes coding for 4-HPPD in C. immitis-one is upregulated and the other is downregulated. Boyce described a cluster of genes involved in tyrosine catabolism in Pennicillium marneffei and other pathogenic dimorphic fungi. The expression of this cluster of genes is upregulated when tyrosine is the only nitrogen source [28]. However, C. immitis spherules are grown in media containing ammonium salts as the primary nitrogen source.
There are twice as many downregulated genes in young spherules than are upregulated, so many genes that are expressed well in mycelia are expressed poorly in spherules. Downregulated genes are enriched for the cytochrome p450 superfamily, transcription factors and genes with oxidoreductase GO terms (Table 3). There are two homologs of the cytochrome p450 ERG11 gene in C. immitis: CIMG_00573 and CIMG_07469. ERG11 is involved in azole anti-fungal activity in other fungi, and both of these are downregulated in young spherules.
The Stu1 transcription factor was downregulated 5-fold in spherules. This transcription factor was found to be required for optimal hyphal growth in H. capsulatum [30], so the difference between expression of this gene in the mycelial and parasitic form is shared by these two pathogenic, dimorphic fungi. However, direct experiments about the role of Stu1 have not been done in Coccidioides spp. Eleven out of 29 C2H2 type zinc finger domain-encoding proteins were also downregulated more than 2-fold and only four were upregulated. The downregulated genes are clustered together within a phylogenetic tree of the C. immitis C2H2 transcription factors suggesting that they may share functions (Figure 3).  The expression of HLH transcription factors was also downregulated. Both C2H2 and HLH transcription factors influence growth rate and differentiation in Neurospora crassa [31]. Most transcription factors containing the fungal Zn(2)-Cys(6) binuclear cluster domain, the most common class of zinc finger protein in the C. immitis., are not differentially expressed. In H. capsulatum, four transcription factors (Ryp1, Ryp2, Ryp3 and Ryp4) are required for differentiation into yeast [32]. These transcription factors are needed for changes in the transcriptional responses to an increase in temperature that triggers yeast formation [28]. The C. immitis homologs of these proteins, Ryp2 and Ryp4 (also known as FacB), are upregulated 2.7-and 4.72-fold in spherules, but the other transcription factor is not. These factors are critical regulators of yeast morphology and fungal virulence in H. capsulatum [33]. The upregulation of Ryp2 and Ryp4 in C. immitis suggests that this regulator may play a similar role in the Coccidioides spp. as it does in H. capsulatum.
Thirty-four pathways were enriched in downregulated genes (p < 0.01); five had a Benjamini-Hochberg p adjusted value < 0.05. Several pathways involved in detoxification of compounds by glutathione are highly enriched. Other enriched pathways included metabolism of xenobiotics by cytochrome p450, metabolism of chloroalkane and methane. This suggests that the mycelial form may encounter more toxic compounds than spherules. Another interesting observation is that the ergosterol pathway is enriched in downregulated genes, which is consistent with the GO term enrichment results and suggests that mycelia and spherules may have somewhat different susceptibilities to azole antifungal drugs.

Comparison to of Differential Gene Expression in Young Spherules to Previous Studies in Coccidioides Immitis
The RNA used in this study was derived from mycelia and spherules obtained from frozen samples of spherules and mycelia obtained from the previous microarray analysis [10]. The microarray study found that 2.5% of the total genes were upregulated more than 4-fold in young spherules; the number of differentially expressed genes in the current RNA-Seq analysis was significantly larger (8%). However, 69% of the differentially expressed genes identified in the microarray study were also differentially expressed in this RNA-Seq study. These results suggest that the RNA-Seq is a more sensitive technique for determining differential expression, but differences observed in microarray tend to be found in RNA-Seq too.
GO enrichment analysis of the up-and downregulated genes were similar in the two studies. Upregulated genes in the microarray study were enriched for oxidationreduction processes and terms for sulfate and sulfite biosynthesis. Enrichment of metal transporter and homeostasis genes was not recognized in the microarray study. However, GO analysis of the down-regulated genes in the microarray study revealed enrichment of transcription factors. One finding of the microarray study was that 25 protein kinase genes were downregulated in spherules. In the current RNA-Seq study, 16 of these genes were also found to be downregulated.
Whiston has previously published a study comparing C. immitis spherules (Day 4 maturity) to mycelia [7]. Reanalyzing Whiston's data using our Kallisto/DESeq2 pipeline, we found that 902 genes were differentially expressed by at least 2 spherule/mycelia FC in both studies. The FC values for differentially expressed genes in the Whiston study and young spherules/mycelia in this study are compared in Figure 4. There is a positive correlation between the two studies (73% of the genes are in the same quadrants) but there are also obvious exceptions. The difference in spherule maturity may account for some of these disparities.

Differential Gene Expression in Mature Spherules vs. Mycelia
Gene expression in mature spherules was also compared to the gene expression in mycelia (Supplemental Table S1). There were many more genes upregulated in mature spherules (960 genes) than in young ones (408 genes), but the great majority of the genes upregulated in young spherules were also upregulated in mature organisms. In contrast, there were many genes that were downregulated in young spherules but not mature spherules and vice versa ( Figure 5). Genes that were differentially expressed in both young and mature spherules had very similar FC values (Supplemental Figure S1). This suggests that upregulation of some genes may be needed for all maturities of the spherule phase.  The enrichment of GO terms in the upregulated genes were similar to the results in young spherules; oxidation-reduction, transmembrane transport and integral membrane component terms were highly enriched (Supplemental Table S2). In contrast, one of the The enrichment of GO terms in the upregulated genes were similar to the results in young spherules; oxidation-reduction, transmembrane transport and integral membrane component terms were highly enriched (Supplemental Table S2). In contrast, one of the most significantly enriched GO terms in the genes that were downregulated in mature spherules were associated with microtube activity and kinesin. The functional consequences of downregulation of microtubule genes are very difficult to predict, since they play so many different roles in cell biology. The cellular component GO term for the cell wall was also highly enriched in the downregulated genes.

Mature Spherules Compared to Young Spherules
Gene expression in mature spherules was also compared to that in young spherules. A relatively small number of genes were differentially expressed. The most differentially expressed genes are in Table 4.  The most highly upregulated gene is a homolog of Hsp31, which is a methylglyoxalase, a chaperone stress-response gene in yeast [34]. The expression of this gene increases in response to DNA replication stress [34]. A ferritin-like protein, involved in iron regulation and oxidation reactions is also upregulated, as are genes involved in mRNA splicing, transcription, and meiotic recombination. Expression of alpha-amylase-1 gene was also up-regulated; this gene is required for pathogenicity in H. capsulatum, but its role in the pathogenicity of C. immitis is unknown [35]. The most dramatically downregulated gene in mature spherules was tyrosinase, a gene that plays a central role in the synthesis of melanin. The beta-glucan synthesis-associated protein Kre6 and beta-glucosidase, two genes that influence beta-glucan metabolism were also downregulated, as were two other genes influencing DNA replication (cell division cycle protein Cdc20 and DNA replication helicase Dna2). The downregulation of these genes suggests that remodeling of the cell wall and DNA synthesis decrease in mature spherules.

Expression of Transposable Elements in Young Spherules Compared to Mycelia
We have previously identified 1309 Gypsy, Copia and TcMar transposons in the C. immitis RS genome and found that proximity to a transposable element (TE) was associated with a lower level of protein-encoding gene expression in C. immitis mycelia [31]. In this study, we compared expression of TE mRNA in mycelia and young spherules. Only 350 TE genes had an adjusted p < 0.05, and the median baseMean of these genes was 20.4, indicating most TEs were poorly expressed. A total of 230 were upregulated in young spherules and only 15 were downregulated. Gypsy TE were the most common upregulated transposons. Analysis of upregulated predicted protein-encoding genes in young spherules showed that a few of them are transposon proteins. A previous proteomic study found one transposon protein in C. posadasii spherules [32]. There was a total of 77 genes that were within 1 kb up-or downstream of the upregulated TE. The median FC of these genes was 3.90, so genes near the up-regulated TE were somewhat overexpressed in spherules.

Summary
The transition from the environmental arthroconidia/mycelial form to the spherule is required for the pathogenesis of Coccidioides spp. The morphologic changes are dramaticthe organism changes from a mycelium with internal spores to a round structure that enlarges circumferentially and divides internally to form a large number of endospores [36]; these are released and can differentiate into mature spherules. The importance of spherule maturation and spherule release is demonstrated by the observation that an engineered chitinase deletion mutant that does not form spherules is avirulent [37]. For these reasons, understanding the transcription program of spherules as they mature is clearly important for comprehending the biology of this dimorphic fungus.
Our findings in this study confirm previous experiments in showing that there are substantial differences in gene expression in spherules and mycelia [7,10]. One of the most interesting findings of this study is that differentially expressed genes are more likely to be genus-specific than those that are not. Although the function of these genes in unknown, and they likely have evolved fairly recently, the results in our study suggest that they play an important role in mycelial/spherule differentiation.
The expression of copper and iron transporters is upregulated in spherules, which is consistent with previous studies in C. posadasii [21] and H. capsulatum [23]. Two of the four Ryp transcription factor genes that have been found to be important for yeast formation and virulence in H. capsulatum are also upregulated in C. immitis, suggesting that these two organisms may use similar transcription programs for phase differentiation.
Another interesting group of genes that is downregulated in young spherules are the transcription factors. The APSES family transcription factor Stu1, several HLH and many C2H2 transcription factors are downregulated. Identification of the genetic programs driven by these transcription factors in C. immitis should be useful for understanding the biological significance of these changes.
Analysis of mature spherules shows that almost all the upregulated genes in young spherules remain upregulated. This suggests that this group of genes may be required for life in the spherule phase; examining the transcriptome of spherules of intermediate maturities would address this question. Furthermore, there were a very small number of differentially expressed genes between the two spherule maturities, which also suggests that young and mature spherules have very similar transcription programs. This study confirms that many genes are differentially expressed in mycelia and spherules and suggests that further investigation of the expression profiles of these morphologic forms may be useful for understanding the biology of this organism.

Institutional Review Board Statement:
This study did not involve experimental animals.

Informed Consent Statement:
This study did not involve human subjects.

Data Availability Statement:
The data for this manuscript are available at NCBI GEO, accession number GSE171286.

Conflicts of Interest:
The authors declare no conflict of interest.