Assessment of RACGAP1 as a Prognostic and Immunological Biomarker in Multiple Human Tumors: A Multiomics Analysis

Several recent studies have pointed out that arc GTPase activating protein 1 (RACGAP1) is a putative oncogene in many human tumors. However, to date, no pan-cancer analysis has been performed to study the different aspects of this gene expression and behavior in tumor tissues. Here, we applied several bioinformatics tools to perform a comprehensive analysis for RACGAP1. First, we assessed the expression of RACGAP1 in several types of human tumors and tried to correlate that with the stage of the tumors analyzed. We then performed a survival analysis to study the correlation between RACGAP1 upregulation in tumors and the clinical outcome. Additionally, we investigated the mutation forms, the correlation with several immune cell infiltration, the phosphorylation status of the interested protein in normal and tumor tissues, and the potential molecular mechanisms of RACGAP1 in cancerous tissue. The results demonstrated that RACGAP1, a highly expressed gene across several types of tumors, correlated with a poor prognosis in several types of human cancers. Moreover, it was found that RACGAP1 affects the tumor immune microenvironment by influencing the infiltration level of several immune cells. Collectively, the current study provides a comprehensive overview of the oncogenic roles of RACGAP1, where our results nominate it as a potential prognostic biomarker and a target for antitumor therapy development.

Introduction
The process of tumorigenesis has proven to be complex, involving a series of interactions between various genes that consequently transfer the cells from the normal state to the cancerous condition [1,2]. Hence, a deep study of different oncogenes is an essential process to explore the molecular mechanisms of several genes in cancer development [3]. A pan-cancer analysis that employs available databases such as The Cancer Genome Atlas (TCGA) gave us a great opportunity to analyze the expression and the behavior of a specific

Results
The abbreviations and the full names of analyzed tumors in the current study are shown in Supplementary Table S1.

RACGAP1 Increased Expression in Many Tumor Types
Initially, TIMER2 was used to examine the differential expression of our target gene between malignant and neighboring healthy tissues. RACGAP1 was observed to be significantly overexpressed in several cancers ( Figure 1A), including BLCA, BRCA, CHOL, COAD, ESCA, GBM, HNSC, KICH, KIRC, LIHC, LUAD, LUSC, PRAD, READ, STAD, as well as UCEC (p < 0.001), CESC, KIRP, and PCPG (p < 0.01). Due to the inaccessibility of normal tissue expressions to be used as a reference for some malignancies in TIMER2, we accessed the Interactive Gene Expression Profiling Analysis (GEPIA) database to acquire differential expression in normal and malignant tissues for these cancers. Most of the malignancies examined, including ACC, DLBC, LGG, O.V., SARC, SKCM, THYM, and UCS, exhibited a statistically significant upregulation of RACGAP1 in cancerous versus normal tissues ( Figure 1B). On the contrary, only two cancers, TGCT and LAML, did not exhibit significant differences and upregulation of RACGAP1 in healthy tissues compared to malignant tissue, respectively. the tumors that lack normal tissue for comparison in TIMER2.0 database and analyzed in the GEPIA database; (C) tumors experienced a positive correlation between RACGAP1 expression and tumor stage when analyzed with the TISDIB web server; (D) tumors experienced a consistent positive correlation between RACGAP1 expression and tissue type (normal-tumor-metastatic); * means statistically significant p > 0.05, ** means high statistically significant p > 0.001 and *** means p > 0.0001. Following examining the differential expression of RACGAP1 in normal and malignant tissues, we intended to investigate the association between RACGAP1 expression and cancer staging. Using the TISDIB web server, we determined that RACGAP1 expression was significantly associated with the disease stage of ACC, KIRC, KIRP, LIHC, LUAD (p < 0.001), BRCA, KICH, UCEC (p < 0.01), and LUSC (p < 0.05) ( Figure 1C). The "compare tumor, normal, and metastasis" module of the TNMplot online server was then utilized to correlate the RACGAP1 mRNA expression levels with the metastases. RACGAP1 expression was significantly elevated in the breast, kidney, and liver when normal and tumor tissues were compared, and this pattern was also maintained between metastatic and cancerous tissues ( Figure 1D).

Differential Protein Expression
Following investigating RACGAP1 at the transcriptional levels, we analyzed its protein level using the National Cancer Institute's CPTAC dataset of large-scale proteome data. Figure 2A-G demonstrate that the expression of the RACGAP1 protein was significantly elevated in the colon, HNSC, hepatocellular carcinoma, LUAD, PAAD, UCEC, and colon cancer compared to normal tissues (p < 0.001). Subsequently, we retrieved the IHC figures for the normal and malignant tissues to validate our earlier results. Normal tissues of the colon, nasopharynx, liver, lung, ovary, pancreas, and endometrium did not exhibit any staining, so the findings were consistent. Similarly, it was moderate to high in the comparable malignant samples (

Increased RACGAP1 Level Estimates Poor Clinical Outcomes
To investigate the connection between RACGAP1 expression patterns and survival rates, we utilized GEPIA as well as Kaplan-Meier (KM) plotter databases. From GEPIA, we determined that the overexpression of our target gene is associated with poor disease-free survival (DFS) prognosis for ACC (p < 0.001), LGG, PAAD, PRAD, UVM (p < 0.01), KIRP, LIHC, MESO, and SARC, as well as SKCM (p < 0.05) ( Figure 3). Furthermore, analysis of patients' overall survival (OS) revealed that in addition to ACC, KIRP, LGG, LIHC, MESO, PAAD, SARC, and SKCM patients with a poor prognosis (similar to the disease-free survival model), LUAD (p < 0.001) patients also had a poor survival with the overexpression of RACGAP1 ( Figure S1). Analysis results from the second server showed that the expression of RACGAP1 negatively influenced OS, distant metastasis-free survival (DMFS), and relapse-free survival (RFS) (p < 0.001) ( Figure 4A) in breast cancer patients. Regarding ovarian cancer, RACGAP1 was predicted to correlate with poor OS (p < 0.01) and progression-free survival (PFS) (p < 0.05), but not post-progression survival (PPS) ( Figure 4B). Moreover, our analyzed gene was related to poor OS, first progression (F.P.) (p < 0.001), and PPS (p < 0.01) in lung cancer ( Figure 4C). Gastric cancer was also affected by RACGAP1 expression in terms of patients' survival in relation to OS, PPS (p < 0.001), and F.P. (p < 0.01) ( Figure 4D). Finally, liver cancer demonstrated poor OS, disease-specific survival (DSS), PFS, and RFS (p < 0.001) ( Figure 4E) with RACGAP1 expression. The results of both databases indicate that an overall poor clinical outcome is expected with RACGAP1 expression in cancer patients.

RACGAP1 Mutation Analysis
In the present study, we investigated the copy number alteration (CNA) of RACGAP1 in clinical specimens. Based on the results of the cBioPortal web server, the endometrial cancer subtype had the highest frequency of RACGAP1 mutations, approximately 6%, where the mutation was the predominant form of genetic alterations. With the exception of glioma and leukemia, which exhibited deep deletion as the predominant form of genetic alterations, and pancreatic cancer, which revealed two different forms of alterations (mutation as well as deep deletion) with a comparable percentage of occurrence, it is notable that the majority of the investigated cancers exhibited either "amplification" or "mutation" as the predominant form of genetic alteration ( Figure 5A). Subsequently, we investigated the sites and types of RACGAP1 mutations ( Figure 5B), where we found 106 total mutations with missense mutations in the first place with 87 recorded samples. Regarding gene location, site R427 was reported to be the most altered site with three mutations (sample with breast invasive ductal carcinoma, second with COAD, and last with LUSC). Finally, we analyzed the correlation between RACGAP1 alteration and patients' survival, where genetic alteration of RACGAP1 showed no significant difference between altered and unaltered groups in terms of DFS, PFS, OS, and DSS ( Figure 5C).

RACGAP1 Phosphorylation Analysis
We analyzed the connection between RACGAP1 phosphorylation levels and tumorigenesis using phosphorylation analysis. Three phosphorylation sites, including S203, T342, and T580, were considerably elevated in breast cancer compared to normal tissues ( Figure 6A). Furthermore, the positions S257, T338, and T580 of RACGAP1 were significantly phosphorylated in tissues of hepatocellular carcinoma, LUAD, and ovarian cancer relative to healthy tissue ( Figure 6B-D). In contrast to normal tissues, our investigated locations of RACGAP1 exhibited higher phosphorylation in the investigated malignancies.

RACGAP1 Correlates with Immune Infiltration in Several Tumor Types
First, we tried to find the correlation between RACGAP1 expression and the infiltration of two types of immune cells with opposing roles against tumor growth. Analysis of MDSC infiltration, a cell with immunosuppressive functions in tumor [17], in the panel of TCGA tumors showed that more than 80% of the studied tumors experienced a positive correlation between RACGAP1 expression and MDSC infiltration. Notably, no cancer from the analyzed panel witnessed a negative correlation between RACGAP1 expression and MDSC infiltration ( Figure 7A). On the other hand, analysis of NKT cell infiltration, a cell that has a strong antitumor action and was selected as a target for cancer immunotherapy development [18], demonstrated that more than 70% of the analyzed tumors experienced a negative correlation between RACGAP1 expression and NKT cell infiltration ( Figure 7B). Again, no cancers in the panel analyzed witnessed a positive correlation between RACGAP1 expression and NKT infiltration. After filtering the results, we found 18 tumors, namely BLCA, BRCA, CESC, COAD, GBM, HNSC, KICH, KIRP, LGG, LUAD, LUSC, PAAD, PRAD, SARC, SKCM, STAD, THYM, and UVM, experiencing a positive correlation between the expression of RACGAP1 and MDSC in addition to a negative correlation between the same gene expression and NKT cell infiltration. The scatter plots that demonstrate the correlation between the expression of RACGAP1 and the infiltration level of MDSC in these 18 filtered tumors are shown in ( Figure 7C). Following that, the SangerBox online server was employed to find the correlation of RACGAP1 expression with the immune checkpoint, MSI, and TMB. Regarding the immune checkpoint, expressions of RACGAP1 in THCA, KIRC, LIHC, and PRAD were positively correlated with several immune checkpoint genes, while tumors UCS and KICH experienced no significant correlation between our target gene expression and most of the immune checkpoint genes ( Figure 8A). In GBM, UCEC, COAD, STAD, and SARC, the expression of RACGAP1 was discovered to be significantly positively associated with MSI. Furthermore, only one tumor, DLBC, revealed a significant negative association between RACGAP1 expression and MSI ( Figure 8B). Finally, our target gene expression analysis and TMB showed a significantly positive correlation in ACC, PCPG, GBM, LUAD, PRAD, UCEC, COAD, STAD, SKCM, KIRC, READ, KICH, and LGG ( Figure 8C).

Analysis of Interacting and Correlated Proteins to RACGAP1
According to the aforementioned findings, RACGAP1 has a significant correlation with the survival of patients with cancer and influences the immune cells in the tumor microenvironment. Therefore, it is essential to investigate the probable molecular pathways of this gene across several cancers. The top 50 experimentally confirmed RACGAP1-interacting proteins were extracted from the STRING database and displayed as a protein-protein interaction network ( Figure 9A). Furthermore, the GEPIA2 webserver was explored to obtain the 100 genes correlated to RACGAP1 in the panel of TCGA tumors. The "Correlation Analysis" module was employed to obtain the plots of the top five correlating genes. They were ordered as follows: KIF11 (R = 0.81), BUB1 (R = 0.81), CKAP2L (R = 0.80), KIF4A (R = 0.80), and NUSAP1 (R = 0.80) ( Figure 9D). Furthermore, a heatmap created by the "Gene Corr" module at TIMER validated the significant positive correlation between these five genes and RACGAP1 in the complete list of TCGA tumors (with the exception of TGCT, in which the connection with KIF4A was nonsignificant) ( Figure 9B). Seven genes, including KIF20B, SHCBP1, KIF14, PLK1, KIF20A, and PRC1, as well as KIF23, were discovered to be duplicated upon the comparison of the two lists created previously ( Figure 9C).
After the elimination of duplicates, a unique dataset was created by combining the two lists and uploaded to DAVID for Reactome and Gene Ontology (G.O.) enrichment analyses. Our gene list may be associated with DNA repair, chromatin remodeling, cellular response to DNA damage, and DNA replication, as indicated by the results of the biological process study. In terms of cellular components, the majority of genes were localized in the nucleus and nucleoplasm. Furthermore, the supplied gene list was enriched for protein, RNA, and DNA binding whenever its molecular function was evaluated ( Figure 9E). Ultimately, the Reactome pathway analysis revealed that RACGAP1 is strongly associated with the cell cycle and mitotic division ( Figure S2).

Discussion
RACGAP1 is an essential cellular protein that was found to be involved in many cellular processes, including cytokinesis, transformation, invasive migration, and metastasis [19]. Several studies have analyzed the role of RACGAP1 in different types of human cancers. For GIT-related cancers, there was a significant correlation of RACGAP1 protein expression at the invasive front of gastric cancer with older age, tumor size, lymph node metastasis, lymphatic invasion, vascular invasion, and advanced stage [20]. Moreover, higher RACGAP1 expression and Ki-67 index were correlated with unfavorable clinicopathological features in predicting poor outcomes of gastrointestinal stromal tumors [21]. Furthermore, RACGAP1 expression was found to increase malignant tumor potential and was used as a predictive biomarker for lymph node metastasis and poor prognosis in colorectal cancer [22]. Regarding its role in breast cancer progression, RACGAP1 was able to modulate ECT2-dependent mitochondrial quality control to drive breast cancer metastasis [23]. A study on 81 esophageal carcinomas (E.C.) patients showed that RACGAP1 could play a pivotal role in E.C. development, suggesting that it could be used as an indicator of prognosis in E.C. patients [24]. Moving to other types of human tumors, RACGAP1 was able to promote melanoma transendothelial migration through mediating adherens junction disassembly [25] and promoted proliferation and cell cycle progression by regulating CDC25C in cervical cancer cells [26]. RACGAP1 was also found to be a downstream effector of E2F7-dependent resistance to doxorubicin and a prognostic for OS in squamous cell carcinoma [27]. Considering hepatocellular carcinoma, RACGAP1 upregulation was significantly associated with the early recurrence of human hepatocellular carcinoma [28], as it activates the RACGAP1/Rho/ERK signaling axis as a competing endogenous RNA to promote early recurrence of hepatocellular carcinoma [29].
Although several studies have tried to analyze the oncogenic of RACGAP1 in several human cancers, there is a lack of a comprehensive study that can deal with the effect of RACGAP1 from many perspectives in a list of several human tumors. It is already established that the tumor microenvironment is complex, with several factors involved in tumor development, the immune response against this abnormal growth, patients' response to tumor therapy, and OS [30]. This complex status of the tumor confirms the requirement of a deep approach that can correlate a targeted gene with tumor progression through different analysis points. For this purpose, the current study applied a pan-cancer analysis for the oncogenic behavior of RACGAP1. An important characteristic of the oncogenic proteins is their upregulation in tumor tissue compared to the normal one. It is worth noting that the usage of bioinformatics tools and the exploration of tumor-related databases to study the roles and behavior of a specific gene in a list of human tumors or a targeted specific tumor have been applied in several previous studies. For example, IQGAP3 was investigated and expected to serve as an effective prognostic biomarker for pan-cancer immune-related therapy [31]. Additionally, multiomics studies have revealed CCNB1 and butyrophilins as potential prognostic biomarkers for ACC and breast cancer, respectively [32,33].
For this reason, our study started with the differential expression of RACGAP1 in a list of human tumors where it was found to be significantly Not only the tumor stage but also tumor metastasis showed a positive correlation with RACGAP1 expression in breast, kidney, and liver cancers. Our last differential comparison was based on the RACGAP1 protein level analysis in normal and tumor tissues. Again, the trend of elevated RACGAP1 expression in tumor tissues was observed in colon cancer, HNCS, LUAD, PAAD, UCEC, ovarian cancer, and hepatocellular carcinoma, where IHC staining, which was high for RACGAP1 in analyzed tumor tissues, confirmed our findings.
Survival analysis is a basic point of investigation for assessing disease progression and the patient's response to medical treatment [34]. Consequently, the current study aimed to find the correlation between RACGAP1 expression and patients' survival. The results from the GEPIA database demonstrated a positive correlation between RACGAP1 expression and the poor prognosis in ACC, KIRP, LGG, LIHC, PAAD, SARC, and SKCM in terms of DFS and OS. Furthermore, the output of the KM plot analysis confirmed this positive correlation in all models studied of breast, lung, gastric, and liver cancers, which recommends the use of RACGAP1 as a prognostic biomarker in the above-mentioned tumors. Mutations in several genes were found to be good prognostic markers for human cancer; examples include mutated KRAS that was correlated with a poor prognosis of pancreatic [35] and lung cancer [36], and mutated NRAS that was associated with a poor prognosis of metastatic melanoma [37]. Therefore, our next survival analysis step was to study whether the RACGAP1 genetic alteration could also affect patients' survival; we found that RACGAP1 genetic alteration did not significantly correlate with the patient's survival in four analyzed models of altered and unaltered groups.
The status of gene methylation has been extensively studied in several human cancers. Previous studies generally found that DNA hypermethylation was a major mechanism for silencing tumor suppressor genes [38]. On the other hand, the oncogenes experienced a hypomethylation state as a mechanism for their activation to induce tumor progression [39]; for example, the hypomethylation state for the oncogenes AQP1, LINE-1, and ELMO3 was reported in salivary gland adenoid cystic carcinoma [40], colorectal cancer [41], and lung cancer [42], respectively. From this point, we performed a methylation analysis for RACGAP1. As expected, several tumors, including UCEC, COAD, PRAD, BLCA, LIHC, HNSC, TGCT, BRCA, and THCA, showed hypomethylation in tumor samples versus the normal one. Additionally, CpG aggregated methylation data revealed that all of the significant results favored RACGAP1 hypomethylation in the tumor sample versus the normal one (except for CHOL).
Tumor immunotherapy, which witnessed a great evolution in the last few decades, became a well-established approach for fighting against cancer [43] where immune checkpoint inhibitors such as αPD-1 have been approved for the treatment of many types of human cancer, such as malignant melanoma, gastric carcinoma, and hepatocellular carcinoma [44]. For this purpose, it was important to study the correlation between elevated RACGAP1 expression in tumor tissue and the tumor infiltration of different types of immune cells. The first cell analyzed was MDSC, which was found to positively affect tumor cell survival and metastasis [45]. Additionally, it inhibits other cells with fighting ability against growing tumors (CD8 T cells and N.K. cells), supports tumor angiogenesis, and is involved in forming cancer stem cells [46]. Consequently, it was not surprising that the elevated level of MDSC infiltration was correlated with poor clinical outcomes for cancer patients [47]. The current study revealed that ACC, BLCA, BRCA, CESC, COAD, HNSC, HNSC-HPV-, HNSC-HPV+, KIRP, LIHC, LUAD, LUSC, O.V., READ, SARC, SKCM, STAD, TGCT, UCEC, and UCS experienced a positive correlation between RACGAP1 expression and MDSC infiltration. Notably, cytokines such as CCL2, CCL5, and CSF1 were found to be involved in the attraction of MDSCs to the tumor site [48]. However, our finding of a positive association between RACGAP1 expression and MDSC infiltration has not yet been fully investigated. The correlation between RACGAP1 upregulation and specific chemokine expression could be a possible mechanism that might explain this correlation. The second cell that was investigated for its correlation with RACGAP1 upregulation is the NKT cell. This kind of cell demonstrated important roles in fighting against early tumors, where it participates in cancer immune surveillance and secretes several effector molecules [49]. Due to its tumor suppressive roles, NKT cell abundance in the tumor tissue was a positive prognostic factor for patients' survival in several human cancers [50]. Our results revealed a significant negative correlation between RACGAP1 expression and NKT infiltration in BLCA, BRCA, CESC, COAD, HNSC, HNSC-HPV-, KIRC, KIRP, LUAD, MEOV, PRAD, SKCM, STAD, THYM, and UVM. Another interesting finding is that not even one tumor in our analyzed list experienced a positive between RACGAP1 expression and NKT infiltration. Putting the results of RACGAP1 expression and MDSC and NKT infiltration together, we can conclude that the upregulation of RACGAP1 expression could be used as a marker for a poor immune response against a growing tumor.
MSI and TMB are considered promising biomarkers for the patient's response to immunotherapy [51], where a robust antitumor effect of αPD1 treatment was observed with colorectal cancer patients with high microsatellite instability (MSI-H) [52]. Similarly, high TMB was positively correlated with a better clinical outcome across diverse tumors [53]. For this purpose, we tried to find a correlation between the upregulation of RACGAP1 in tumor tissue and those promising biomarkers. Our analysis revealed a positive correlation between MSI in READ, LUSC, UCEC, and BRCA and RACGAP1 expression. Additionally, ACC, PRAD, TGCT, LIHC, and READ experienced a positive correlation between RAC-GAP1 expression and TMB level. Collectively, our findings exposed a research question about the probability of relying on RACGAP1 expression in the above-mentioned tumor as a potential biomarker for patients' response to tumor immunotherapy.
As a final point of analysis, the current study aimed to investigate the molecular mechanism of RACGAP1 in tumor progression where the top 50 interacting proteins and top 100 correlated proteins to RACGAP1 in the tumor tissue were obtained from the STRING and GEPIA2 databases, respectively, and we interestingly found that PARP1 was a common protein in both of the generated datasets. This protein was a regulator for prostate cancer growth and progression through transcriptional regulatory functions [54]; moreover, it was highly expressed in SCLC, where its knockdown led to SCLC growth inhibition [55]. Moreover, PARP1 was a prognostic biomarker for poor clinical outcomes in breast cancer patients [56]; therefore, PARP1 inhibitors were extensively studied as a promising class of anticancer agents [57]. Analysis of RACGAP1 correlated proteins in the tumor tissue revealed that POLR3C, PRKAB2, SETDB1, GPATCH4, and MSTO1 were the top five. SETDB1 was implicated as an oncogene in several human tumors [58], where it was involved in tumor progression in HCC through the methylation of p53 [59]. Furthermore, SETDB1 has been involved in NSCLC progression through WNT-β-catenin pathway stimulation, and for these roles, it was nominated to be a therapeutic target to fight against numerous cancers [60]. It is noteworthy that the complex oncogenic functions of POLR3C, PRKAB2, SETDB1, and GPATCH4 have not been studied yet, and because of being from the top correlating proteins with RACGAP1 in tumor tissues, their potential oncogenic roles and interacting network should be further studied to present clues for novel tumor treatment strategies.

Gene Expression Analysis
Initially, data collected from the Tumor Immune Estimation Resource, version 2 (TIMER2.0) were used to visualize the level of RACGAP1 expression in various cancers compared to normal tissue [56], where the list of tumors of the Cancer Genome Atlas (TCGA) (Supplementary Table S1) was employed for this comparison. Following TIMER2.0 analysis, it was discovered that several cancers lacked normal tissue for comparisons. Subsequently, RACGAP1 expression in healthy and malignant tissue was compared in several human cancers using the Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn/index.html (accessed on 20 September 2022)) [57]. We used the TISBID website to examine the connection between the expression of the RACGAP1 gene and the grading of the tumor [58]. Last but not least, the expression of RACGAP1 in the tumor, normal, and metastatic tissues was evaluated using the TNMplot (differential gene expression analysis in Tumor, Normal, and Metastatic tissues) online server and its Kruskal-Wallis test for significance testing [59].

Protein Expression and Immunohistochemistry (IHC) Staining
To examine the expression patterns of RACGAP1 protein between normal and malignant specimens, the UALCAN program, which conducts protein expression analysis using data collected from the Clinical Proteomic Tumor Analysis Consortium (CPTAC), was utilized [60]. In addition, we utilized HPA to track IHC images of RACGAP1 expression in normal and malignant tissues for tumors, in which UALCAN analysis revealed significant differences in order to validate our findings.

Survival Prognosis Analysis
The GEPIA2.0 webserver was initially utilized to evaluate survival outcomes. We tested RACGAP1 in the "survival analysis" section, chose the whole tumor list from the TCGA cohort, and acquired the heatmap for the two possible approaches from the server (OS and disease-free survival). Then, we analyzed the relationship between RACGAP1 expression as well as survival rates in five different cancers using the KM plotter [61]. (Breast, ovarian, lung, gastric, and liver cancer).

Gene Alteration Analysis
We utilized the cBioPortal web server to conduct an extensive and comprehensive investigation of RACGAP1 polymorphisms [62]. In order to perform this study, we chose "TCGA PanCancer Atlas Studies" as our data source and followed three key guidelines. The "cancer types summary" tab was used to retrieve the mutation frequency as well as the mutation type results. Furthermore, the "Mutations" tab was utilized to view RACGAP1 alteration sites. We visited the "Comparison/Survival" page to examine the relationship between RACGAP1 genetic alterations and survival rates.

RACGAP1 Phosphorylation Analysis
Numerous enzymes and receptors are activated/deactivated by phosphorylation and dephosphorylation processes, which makes protein phosphorylation a significant physiological regulation mechanism [63]. Using data retrieved from the Clinical Proteomic Tumor Study Consortium, we performed a protein phosphorylation investigation of RACGAP1 in malignant tissues versus normal tissues using the UALCAN software (CPTAC).

Immune Reactivity Assessment
Initially, we utilized the TIMER2 web server [64] to determine the association between RACGAP1 expression patterns and immune cells that could also play positive or negative roles in the development of TCGA cancers. We inserted the name of our target gene into the "gene" module within the "immune" partition. We chose two invading immune cell types, including myeloid-derived suppressor cells (MDSCs) as well as natural killer T cells (NKT), having opposing functions in tumor formation. Subsequently, we generated heatmaps and scatter plots to illustrate our investigated correlations. Then, we accessed the SangerBox web server to investigate the relationship between our target gene and three variables: immunological checkpoints, microsatellite instability (MSI), and tumor mutational burden (TMB).

Protein-Protein Interaction and Enrichment Analysis
A total of two web servers were assessed in order to investigate the proteins that may have significant interactions or correlations with RACGAP1. To generate the network of RACGAP1-interacting proteins, we subsequently consulted the STRING database [65]. Setting "Experiments" as the active interaction source and "low confidence" as the interaction score yielded the top 50 interacting proteins. The "Correlation Analysis" module on the same database and the "Gene Corr" module in TIMER were employed to generate the correlation curves and the heatmap for the top five associated genes. Then, we accessed the online server (http://bioinformatics.psb.ugent.be/webtools/Venn/ (accessed on 20 September 2022)) in order to identify the proteins shared by both the "RACGAP1 interacting" and "RACGAP1 correlating" lists. After integrating both lists and deleting duplicate records, the resulting set of data was uploaded to the Database of Annotation, Visualization, and Integrated Discovery (DAVID) [66] for functional enrichment analysis, with the findings shown utilizing "ggplot2" package of R. (4.2.0).

Conclusions
RACGAP1 is an oncogene that was shown to be overexpressed as mRNA and protein in a significant number of human malignancies, and its overexpression was discovered to be correlated with poor health outcomes. RACGAP1 expression was positively correlated with immunosuppressive cell (MDSC) infiltration. On the contrary, the gene was adversely associated with the invasion of tumor-fighting cells (NKT cells). Furthermore, TMB and MSI were connected with RACGAP1 expression in various human malignancies; hence, these results could designate RACGAP1 as a reliable prognostic biomarker, a marker for patients' responsiveness to immunotherapy, and a promising target for cancer therapeutics. Moreover, the current study explored the interaction network of RACGAP1 with other proteins in the tumor microenvironment, where future deep investigation of our primary data can explain the mechanism of RACGAP1 in the induction of tumor progression and open the avenue for the development of novel antitumor treatment.

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


extfieldsofstudy: Medicine Biology

s2fieldsofstudy:
Biology
Environmental Science

text:
Caenorhabditis nematodes colonize ephemeral resource patches in neotropical forests

Abstract Factors shaping the distribution and abundance of species include life‐history traits, population structure, and stochastic colonization–extinction dynamics. Field studies of model species groups help reveal the roles of these factors. Species of Caenorhabditis nematodes are highly divergent at the sequence level but exhibit highly conserved morphology, and many of these species live in sympatry on microbe‐rich patches of rotten material. Here, we use field experiments and large‐scale opportunistic collections to investigate species composition, abundance, and colonization efficiency of Caenorhabditis species in two of the world's best‐studied lowland tropical field sites: Barro Colorado Island in Panamá and La Selva in Sarapiquí, Costa Rica. We observed seven species of Caenorhabditis, four of them known only from these collections. We formally describe two species and place them within the Caenorhabditis phylogeny. While these localities contain species from many parts of the phylogeny, both localities were dominated by globally distributed androdiecious species. We found that Caenorhabditis individuals were able to colonize baits accessible only through phoresy and preferentially colonized baits that were in direct contact with the ground. We estimate the number of colonization events per patch to be low.

among species raises questions about their long-term phenotypic stasis, species coexistence, and the niches they occupy. Previous studies of wild populations of Caenorhabditis find that they live on microbe-rich patches of decaying fruit and vegetable matter (Crombie et al., 2019;Félix et al., 2013;Félix & Duveau, 2012;Ferrari et al., 2017;Frézal & Félix, 2015;Schulenburg & Félix, 2017), the stages on which niche partitioning and interspecific competition play out. Stochastic colonization and extinction rates on these ephemeral resources are key parameters in understanding the local coexistence of species (Dubart et al., 2019).
Some Caenorhabditis species occupy a substrate niche as specialists (Dayi et al., 2021;Kanzaki et al., 2018;Li et al., 2014). The majority, however, have no obvious substrate preference. Studies of Caenorhabditis microbiomes in both the laboratory (Berg et al., 2016) and the field (Dirksen et al., 2016;Zhang et al., 2017) suggest that animals regulate the composition of their gut flora on substrates with differing microbial composition. From these data, one could hypothesize that species with overlapping ranges specialize by occupying niches defined by what they eat. However, it is unclear which microbes are the primary food source of worms in the wild (Schulenburg & Félix, 2017). Beyond food, other factors including predators and pathogens along with nonbiological sources of variation like humidity and temperature may play a role in determining where Caenorhabditis species colonize and proliferate (Crombie et al., 2019;Félix & Duveau, 2012). Still missing is substantial evidence that these Caenorhabditis species preferentially colonize substrates like fruits versus flowers. However, one field study found the degree to which a patch is rotting may influence the incidence of species found on those patches, suggesting that priority effects and ecological succession may also be involved in species coexistence (Ferrari et al., 2017).
Equally critical to understanding Caenorhabditis species adaptations to ephemeral resource patches are determining modes of dispersal. Two models described by Slatkin (1977) represent the extremes of a theoretical spectrum. In the propagule pool model, all colonists are derived from a single patch, whereas in the migrant pool model, colonizers come from every patch in the metapopulation. In Caenorhabditis, these models encompass potential modes of dispersal, either by a phoretic host (Kiontke, 1997;Sudhaus et al., 2011;Woodruff & Phillips, 2018;Yoshiga et al., 2013) or by a semimobile "seed bank" of dauer larvae crawling towards or waiting for a fresh patch (Cutter, 2015). These contrasting modes of dispersal may have profound effects on the level of inbreeding and genetic diversity (Li et al., 2014). In addition, propagule size may contribute to the evolution of a female-biased sex ratio and the evolution of self-fertile hermaphroditism as a means of generating a population growth advantage and reproductive assurance, respectively (Cutter et al., 2019;Hamilton, 1967;Lo et al., 2021;Theologidis et al., 2014).
To better understand Caenorhabditis diversity and the factors that influence it, we performed field surveys and experiments in two of the most intensively studied lowland tropical forests on Earth: Barro Colorado Island (BCI), Panamá, and La Selva, Costa Rica. Barro Colorado Island lies in the center of the Panama Canal in the man-made Lake Gatún. Shortly after its formation, the island was designated a protected nature reserve and has been hosting field research for the last 100 years (Leigh, 1999). Likewise, La Selva Biological Field Station in Sarapiquí, Costa Rica, has been a protected research forest for nearly 70 years (McDade et al., 1994).
We focused our collection efforts on these two localities as they are relatively undisturbed by human activity, and their histories of intensive research provide a rich source of information about the local ecology. One nematode metagenetic study previously found Caenorhabditis DNA in a soil and leaf litter sample at BCI, but the species were not identified (Porazinka et al., 2010). In contrast to the majority of previous work on Caenorhabditis in the tropics, which involved transporting substrates out of the country and isolating animals from nematode growth medium plates days or weeks later, we isolated and cultured all animals immediately in the field, potentially reducing sampling biases that favor species that survive transport and grow well on nematode growth medium. One other study used a combination of these approaches (Félix et al., 2013).
In total, we collected seven species of Caenorhabditis, four of them known only from these collections. We formally describe two new species, C. krikudae sp. nov. and C. agridulce sp. nov., and we place them within the Caenorhabditis phylogeny. Each locality was dominated by globally distributed self-fertile species. We assayed several ecological features related to patch accessibility, patch specificity, and co-occurrence of species. Using baits that vary in their accessibility, we demonstrate that Caenorhabditis are able to colonize baits that are only accessible by phoresy. Further, the colonization rate varied significantly with accessibility where baits making direct contact with the ground were preferentially colonized. We found that individual species tended to occur in habitat patches close to other patches of conspecifics, and we use the frequency of uncolonized patches to estimate the number of colonization events per patch.
Taken together, our data support a model that many Caenorhabditis species are habitat generalists whose population biology is strongly influenced by metapopulation dynamics.

| Collections
We collected nematodes on BCI in May 2012 (wet season), March 2015 (dry season), and August 2018 (wet season). Schemes for sampling varied within and among sampling sessions and included opportunistic sampling and the use of baits as described in the results.
In all cases, worms were isolated from substrates and transferred to Nematode Growth Medium (NGM) plates at the BCI field station and identified as Caenorhabditis by morphology under a stereomicroscope. In 2012 and 2015, material from the forest (e.g., rotting fruits and flowers) was placed directly onto NGM plates, and Caenorhabditis worms were picked to new plates to establish cultures (Barrière & Félix, 2005). These individual patches of organic material are defined as samples in our dataset and were evaluated for the presence of nematodes. For the majority of samples collected in 2018, worms were isolated by the Baermann funnel technique (Baermann, 1917;Tintori et al., 2022) and subsequently cultured on NGM plates. These cultures were transported to New York for species determination as described below.
We collected nematodes at La Selva, Costa Rica, in July 2019, by Baermann funnel. We used two methods to identify Caenorhabditis to species. Individual Caenorhabditis worms were chopped with razor blades, transferred to Whatman paper as described by Marek et al., 2014, and transported to New York. There, the stored nematode DNA was used to identify the specimen to species by ITS2 sequencing. Separately, we established isofemale cultures on NGM plates. These plates were stored at La Selva for six months prior to their transport to New York, where surviving cultures were revived and species identified as described below. Complete collection data are reported in Supporting Information File 1.

| Species identification
Species were identified by sequencing a fragment of rDNA to derive a prediction. Subsequently, experimental crosses were performed with isolates of known species identity to establish a biological species assignment (Félix et al., 2014;Ferrari et al., 2017;Kiontke et al., 2011;Stevens et al., 2019).
Mating tests were performed with worms of known species identity. Cross plates were monitored for the presence of viable progeny. For isolates of androdioecious species, hermaphrodites were crossed to males from strains of C. briggsae and C. tropicalis. In some cases, we used males of wild-type strains AF16 and JU1373, respectively, and monitored plates for male progeny. In other cases, we used males of strains QG2801, an AF16 derivative carrying GFPexpressing transgene syIs803 (Inoue et al., 2007), and QG3501, a derivative of C. tropicalis NIC58 carrying mCherry-expressing transgene qgIs5 (Noble et al., 2021), and we monitored for wild-type green-or red-fluorescent offspring.

| Sequencing and assembling the transcriptome of Caenorhabditis krikudae n. sp.
We generated the C. krikudae n. sp. inbred line QG3077 by 28 generations of full-sibling mating from isofemale line QG3050. We generated RNA-seq mRNA transcriptome data using a pool of five mixed-stage populations of QG3077, with each population being subjected to a different condition. All worms were grown at 25°C on 10 cm NGMA plates (for 1 L: 3 g NaCl, 5 g bacto-peptone, 10 g agar, 7 g agarose, 1 ml cholesterol 5 mg/ml in ethanol, 1 ml CaCl 2 1 M, 1 ml MgSO 4 1 M, 25 ml KPO 4 1 M). One population was fed with CemBio strains (Dirksen et al., 2020), and the other four were fed with E. coli OP50. The conditions for OP50 populations consisted of (1) mixed-stage, (2) starved, (3) heat-stressed, and (4) cold-stressed.
Temperature stress consisted of exposing the worms to either 35°C or 4°C for 2 h followed by a 2-hour recovery prior to RNA extraction.
Total RNA was isolated using TriZol following the protocol described in Green and Sambrook (2020). The mRNA library was constructed using the Illumina Stranded mRNA Prep Ligation protocol. The library was sequenced using a NextSeq 500 MidOutput 2X150 for 300 cycles. Paired-end sequences were trimmed with Trim Galore (https://github.com/Felix Krueg er/TrimG alore). Trimmed sequences were assembled into a transcriptome using Trinity  also running default parameters for paired-end reads.
We then generated the longest predicted ORFs using TransDecoder (https://github.com/Trans Decod er/Trans Decoder) for use in phylogenetic analyses.

| Sequencing and assembling the genome and transcriptome of Caenorhabditis agridulce n. sp.
Isofemale strain QG555 was grown on 9 cm NGMA plates. We harvested nematodes just after starvation and washed using M9 several times to remove E. coli. For genomic DNA extraction, the nematode pellets were suspended in 600 μl of Cell Lysis Solution (Qiagen) with 5 μl of proteinase K (20 μg/μl) and incubated overnight at 56°C with shaking. The following day, the lysate was incubated for one hour at 37°C with 10 μl of RNAse A (20 μg/μl) and the proteins were precipitated with 200 μl of protein precipitation solution (Qiagen). After centrifugation, we collected the supernatant in a clean tube and precipitated the genomic DNA using 600 μl of isopropanol. The DNA pellets were washed in 70% ethanol and dried for one hour before being resuspended in 50 μl of DNAse-free water. For RNA extraction, we resuspended 100 μl of nematode pellet in 500 μl of Trizol (5 volumes of Trizol per volume of pelleted nematodes). The Trizol suspension was frozen in liquid nitrogen and then transferred to a 37°C water bath to be thawed completely. This freezing/thawing process was repeated four to five times and the suspension was vortexed for 30 s and let rest for 30 s (five cycles). A total of 100 μl chloroform was added and the tubes were shaken vigorously by hand for 15 sec and incubated for 2-3 min at room temperature. After centrifugation (15 min at 13,000 rpm and 4°C), the aqueous (upper) phase containing the RNA was transferred to a new tube and precipitated with 250 μl of isopropanol. The pellets were washed in 70% ethanol and dried for 15-20 min before being resuspended with 50-100 μl of RNAse-free water. An aliquot of each DNA and RNA preparation was run on agarose gel to check their quality and quantitated with Qubit (Thermo Scientific). Two short-insert (insert sizes of 300 and 600 bp, respectively) genomic libraries and a single short-insert We performed quality control of our genomic and transcriptomic read sets using FastQC (v0.11.9; Andrews, 2010) and used fastp (0.20.1; Chen et al., 2018;--length_required 50) to remove low-quality bases and Illumina adapter sequence. We generated a preliminary genome assembly using SPAdes (v3.14.1; Bankevich et al., 2012;--only-assembler --isolate -k 21,33,55,77) and iden- We also mapped the genomic reads to the genome assembly using bwa mem (0.7.17-r1188; Li, 2013). We provided the assembly, the BAM file, and the BLAST and Diamond files to blobtools (1.1.1; Laetsch & Blaxter, 2017) to generate taxon-annotated, GC-coverage plots, which we used to identify contaminant contigs. Any read pair that mapped to the contaminant contigs was discarded. Using this filtered read set, we generated a final assembly using SPAdes (--isolate -k 21,33,55,77,99). We also generated a transcriptome assembly using Trinity (Trinity-v2.8.5; Haas et al., 2013), which we then used to scaffold the genome assembly using SCUBAT2 (available at https://github.com/GDKO/SCUBAT2). We used numerical metrics and BUSCO (v4.1.4;Simão et al., 2015;−l nematoda_odb10 -m genome) to assess assembly quality and biological completeness, respectively. Prior to gene prediction, we generated a species-specific repeat library using RepeatModeler (2.0.1; Smit & Hubley, 2010;-engine ncbi), and combined this library with known Rhabditid repeats from RepBase (Jurka et al., 2005). This repeat library was then used to soft-mask the genome using RepeatMasker (open-4.0.9; Smit et al., 1996;−xsmall). We predicted genes in the genome by aligning trimmed transcriptomic data to the genome using STAR (2.7.3a; Dobin et al., 2013;-twopassMode Basic) and providing the resulting BAM file to BRAKER2 for gene prediction (2.1.5; Brůna et al., 2021;--softmasking). We used BUSCO (−l nem-atoda_odb10 -m proteins) to assess gene set completeness.

| Phylogenetic analysis
We identified a set of orthologous proteins by running BUSCO (Seppey et al., 2019) using the nematode_odb10 dataset on each of the 36 nematode genomes found in Table S1. Multisequence fasta files for each ortholog were extracted using busco2fasta (https:// github.com/lstev ens17/ busco 2fasta) with the setting -p 0.8, meaning each ortholog was required to be in 80% or 28 of the 36 species. Orthologous sequences were then aligned with MAFFT (Katoh & Standley, 2013) and ML gene trees estimated using IQ-TREE (Nguyen et al., 2015), both on default settings. Newick trees were concatenated into a single file and a species tree was estimated using ASTRAL-III (C. Zhang et al., 2018), which uses a coalescent framework. We also generated a species tree using a supermatrix of all concatenated orthologs. To generate the supermatrix we used TrimAl (Capella-Gutiérrez et al., 2009) to remove poorly aligned regions using the settings -gt 0.8 -st 0.001 -resoverlap 0.75 -seqoverlap 80. Sequences were subsequently concatenated using catfasta2phyml (https://github.com/nylan der/catfa sta2p hyml). A tree was then inferred with IQ-TREE using the LG substitution model (Le & Gascuel, 2008), modeling the rate variation among sites using a Discrete Gamma model (Yang, 1994) with 4 categories. Support was estimated using 1000 ultrafast bootstrap replicates (Hoang et al., 2018). We then estimated ASTRAL-III tree branch lengths in units of replacements per site rather than coalescent units using IQ-TREE with the same parameters as the supermatrix analysis while fixing the tree by the output of the ASTRAL-III analysis using the -te setting. All newick trees were visualized using the ITOL web browser (Letunic & Bork, 2019).

| The Caenorhabditis faunas of BCI and La Selva
We recovered Caenorhabditis nematodes from 225 samples col- To estimate the frequency of Caenorhabditis across samples while minimizing variation due to differences in sampling technique, we consider a dataset of 177 samples collected and processed by a single investigator in August 2018. These samples included a range of rotten fruits, flowers, stems, fungi, and leaf litter.
To assess the completeness of our survey, we used rarefaction of the chao2 incidence-based estimator (Chao et al., 2014;Hsieh et al., 2016), which generated an estimated species richness of 6 ± 0.34 (95% CI) ( Figure 1). These data suggest that we have recovered the maximum number of species at BCI, conditional on our sampling strategy. The two most abundant species, C. briggsae and C. tropicalis, are androdioecious (males and self-fertile hermaphrodites), and their geographic distributions are cosmopolitan and pantropical, respectively. The other species are gonochoristic (males and females). One of these species, C. agridulce n. sp., which we formally describe in the Appendix to this paper, has also been found in French Guiana (Ferrari et al., 2017), Mexico, and Southern California (Appendix 1).
We successfully recovered Caenorhabditis nematodes from 77 samples at La Selva, Costa Rica ( Figure 1; Supporting Information File 1). These derive from an opportunistic sampling of rotten fruits, flowers, mushrooms, and litter in 2019. These samples yielded only 3 different species, one of which is known only from our collections at La Selva (C. sp. 60). La Selva differed from BCI in that C. tropicalis was most prevalent (present in 55 samples), followed by C. briggsae (32 samples). Gonochoristic C. sp. 60 was isolated from a single substrate, which contained thousands of individuals. The rarefaction of the chao2 incidence-based estimator generates a species richness of 3 ± 0.48 (95% CI) ( Figure 1). This suggests that the lower number of observed species at La Selva is not due to inadequate sampling given our sampling strategy. We measured substrate temperature for 22 samples that contained F I G U R E 1 Collection sites for Caenorhabditis species used in this study. Caenorhabditis were collected at two localities: Barro Colorado Island, Panamá, and La Selva, Sarapiquí, Costa Rica. (a-c) Distribution of species collected from opportunistic sampling from each locality by year. Each marker represents a patch positive for that species. Patches may be plotted multiple times if species co-occurred on the same patch. Patches are jittered to prevent overpotting. (d) A field of rotting Spondias mombin substrates (e,f). Rarefaction curve of the chao2 incidence-based estimator for both localities. The solid line represents the predicted species richness the dotted line represents an extrapolation of species richness. The gray area is the 95% confidence interval. Caenorhabditis; these ranged from 24.1 to 28.4°C. There was no difference in temperature between patches containing different species, with a mean temperature of 26°C for each species (Supporting Information File 1).
To understand the phylogenetic positions of the undescribed species, we sequenced and assembled transcriptome for C. krikudae n. sp. and a genome for C. agridulce n. sp. Using these assemblies and the assemblies of 34 additional Caenorhabditis species, we identified 1931 single-copy orthologs that were represented in at least 28 of the 36 species. We reconstructed the Caenorhabditis phylogeny using two approaches. First, we used a coalescentbased approach with individual gene trees as input. Second, we used a maximum likelihood approach using a concatenated alignment of all orthologues as input. The resulting phylogenies We name this clade the Auriculariae group, defined as species more closely related to C. auriculariae than to C. elegans. Based on data from ITS2 sequence only, C. sp. 60 is sister to C. macrosperma within the Japonica group (NCBI Accession: OL960095).
Overall, the species found at BCI and La Selva span the

| Substrate specificities
To test whether the common Caenorhabditis species show substrate specificity, we analyzed the dataset of 177 samples processed by a single investigator in 2018 (Table 1, Supporting Information File 1).
Each of the four most common species was collected from multiple types of fruit and flower. Classifications of the substrates, at high levels (fruit vs. other) or lower levels (fruit type), revealed no significant association between Caenorhabditis generally or any of the common species specifically and any substrate (logistic regression, p > .05; see Supporting Information File 1 for specific values).
Acknowledging the very limited statistical power of these tests, we interpret this as evidence that the common species are substrate generalists, colonizing and proliferating in any available habitat patch.

| The spatial patterning of patch occupancy
To understand the spatial patterning of Caenorhabditis among habitat patches, we performed a hierarchical spatial sampling of a single substrate type, rotten flowers of Gustavia superba, in May 2012. We selected four G. superba trees spread across the island. At each tree, we established three well-separated 1 m 2 quadrats. Within each quadrat, we sampled four rotten flowers, each at least 10 cm apart.
From each flower that yielded Caenorhabditis, we established isofemale or isohermaphrodite lines from four or more randomly selected worms from each flower. At one tree, only two quadrats were sampled. In total this sampling scheme involved 44 samples of G. superba flowers.
Thirty-six of 44 G. superba flowers (82%) contained There is no evidence that the presence of one species affects the probability of observing a second species within a sample. For example, C. briggsae and C. tropicalis are present in 66% and 9% of the 44 samples; the expected co-occurrence under independence is 2.6/44 and we observe co-occurrence of 3/44 samples.

| Colonization patterns among classes of bait
To test how substrate type and accessibility affect rates of colonization by Caenorhabditis, we set up arrays consisting of several bait types. At each of the seven sites on BCI, we set up a 7-by-7-meter field site with five arrays of baits (four in the corners, one in the center). Each bait array consisted of six agar baits, each bait of a different type (Figure 4), arranged 3 × 2 with 30 cm spacing between the 6-cm diameter baits. Our experiment as a whole therefore included 210 baits in total.
Two of the bait types consisted of 1% agar supplemented with 5% peptone and 1% glucose. In one such bait type, the agar was

Frequency
Colonization Events briggsae specifically (p < 10 −4 ). The worms preferentially colonized baits that made direct contact with the ground over baits that were isolated from the ground by plastic. In both of those classes of bait, the worms preferentially colonized those with peptone enrichments over those with heat-defaunated Gustavia superba flower slurry. And among Gustavia plugs, they preferentially colonized those that were not supplemented with raw Gustavia slurry. Caenorhabditis showed a bait-type distribution that does not differ significantly from the distribution of baits colonized only by non-Caenorhabditis nematodes (Fisher's exact test, p = .92), though the power of this test is limited by the small size of the dataset. Another way to state this is that the probability of Caenorhabditis colonizing a bait type is correlated with the probability of only non-Caenorhabditis worms colonizing that same bait type (r 2 = .98, p < .001). This means that Caenorhabditis and non-Caenorhabditis nematodes preferred the same baits and colonized each bait type with similar proportions.

| Test of colonization by phoresy
We used size-selective exclosures to determine whether coloniza- The lids had a circular opening with 3.1-cm diameter, which was either totally open or covered with a nylon mesh to restrict access by animals larger than the mesh size. The mesh openings restricted passage to animals smaller than 4, 1, 0.064, or 0.01 mm.
After 5 days in the field, we collected the slurry samples and transferred a small volume (approximately 1 cm 3 ) to NGM plates. If worms emerged, we attempted to establish cultures. Surviving cultures were cryopreserved in New York, and species were identified by sequencing and mating tests.
One bait was lost, and of the 143 baits that we recovered, we found nematodes in 30, including three species of Caenorhabditis and at least ten additional species ( Figure 5; Supporting Information File 1). Because some baits were colonized by multiple species, we count 34 species observations overall.
C. briggsae and C. tropicalis both colonized baits inside plastic cups, demonstrating that these animals can colonize new substrates by phoresy on other animals. Conversely, Oscheius tipulae, which colonized seven baits, only colonized baits that were accessible directly from the soil or leaf litter. Among the animals found in the plastic cups with Caenorhabditis were mites, dipterans, hemipterans, coleopterans, and hymenopterans; fly larvae and pupae were common.
We observed substantial heterogeneity among the plots ( Figure 5).
Bait accessibility significantly affected colonization rates by nematodes generally (p = 5.4 × 10 −8 ; analysis of deviance from logistic regression) and by Caenorhabditis specifically (p = .007). Caenorhabditis colonized only the three most accessible classes of bait, suggesting that their phoretic hosts did not pass through mesh with pores of a millimeter or smaller.

| DISCUSS ION
Over the past twenty years, a community effort to study Caenorhabditis elegans and its relatives in their natural context has been fruitful. The catalogue of Caenorhabditis species and wild isolates has increased dramatically and along with it the ability to apply population, quantitative, and comparative genomic methods (Andersen & Rockman, 2022;Cook et al., 2017;Stevens et al., 2019). Despite these advances, a well-supported model of any Caenorhabditis species' population biology is still missing. Here, we present a deep sampling of Caenorhabditis natural diversity in two of the most extensively studied neotropical field sites, along with a collection of experiments aimed at understanding Caenorhabditis species ecology and metapopulation structure. In total we collected seven species, four of which were only found in these collections (BCI: C. becei, C. panamensis, and C. krikudae n. sp.; La Selva: C. sp. 60). We estimate that we recovered the total number of species at both field sites accessible to our sampling scheme, which was limited by various factors like time of year, selection of visibly rotting material, nematode isolation method, and proximity of sampling localities to trails.
Species from four major clades within Caenorhabditis were found in these forests, including representatives of the Elegans, Japonica, Angaria, and Auriculariae groups. Our findings comport with biogeographic hypotheses about the history of Caenorhabditis diversity (Cutter, 2015). In particular, we find three species (C. becei, C. panamensis, and C. sp. 60) that are part of a neotropical-endemic clade within the Japonica group. Species in this group can be locally abundant in neotropical forests, but their geographic ranges appear to be quite narrow. Each species is known only from a single region, with no overlap among the species in this group found in La Selva, BCI, French Guiana, or Dominica (Marie-Anne Félix, personal communication; Stevens et al., 2019). Most parts of the neotropics have not yet been surveyed for Caenorhabditis, and we infer that many Japonica group species remain to be discovered there. Conversely, Elegans group species are represented exclusively by two widely F I G U R E 3 Species are patchily distributed among rotting Gustavia superba flowers. (a) 10 × 10 meter plots were systematically sampled at each of four focal trees. At each plot, four flowers were collected from two or three 1-meter quadrats. Each box represents a flower; each color represents the species present on that flower (b). The distribution of C. briggsae colonization events per flower under a simple Poisson model (mean = 1.08).   substrates. We found that patterns of substrate colonization were highly correlated between Caenorhabditis and non-Caenorhabditis nematodes, suggesting that these Caenorhabditis communities are substrate generalists. This conclusion is consistent with our opportunistic sampling data, which found no associations between any substrate type and incidence of Caenorhabditis.

F I G U R E 4 Colonization rates vary in
Our data suggest that Caenorhabditis species on BCI disperse by phoresy. In our exclosure experiment, Caenorhabditis colonized baits that were directly accessible from the ground, isolated from the ground in a cup, and isolated in a cup and further blocked by mesh with openings of 4 mm or greater. Baits isolated by mesh with openings of 1 mm or smaller were not colonized. These data reinforce the idea of a propagule pool model of dispersal in which individuals migrate from a single patch. However, the design of this experiment provides no quantitative measure of the proportion of phoretic versus soil colonizations. In contrast to Caenorhabditis species, Oscheius tipulae only colonized baits making direct contact with the ground suggesting that they were not colonizing baits using phoresy. While O. tipulae is commonly found in soil (Félix, 2006), it was originally isolated from the cadavers of larvae of Tipula paludosa, a marsh cranefly (Lam & Webster, 1971 season in 2015, collections at that same tree yielded only C. tropicalis. One model is that habitat patches are colonized randomly from the local species pool, as suggested by the patchy species distribution of G. superba flower occupancy. An alternative is that species differences among years illustrate ecological succession at larger scales than the level of an individual substrate and its lifespan. Félix and Duveau (2012) more systematically describe seasonal shifts in the abundance of C. briggsae and C. elegans in a French orchard, paralleling their finding that C. briggsae outcompetes C. elegans at higher temperatures in the lab. In the neotropics, changes between wet and dry seasons impact the availability of fruit and flower patches, and the availability of phoretic vectors (Leck, 1972).
Species in our spatial sampling dataset appeared to differ in their distributions across sampling sites and quadrats. C. briggsae was present in every quadrat at every focal tree sampling site, while other species had a patchier distribution over a scale of meters and at scales between focal tree sampling sites. These patterns could indicate differences in colonization efficiency and differences in the scale of dispersal between species, which might be picked up by a larger dataset. Under the assumption that animals colonize patches independently and randomly, we estimated that about 44% of patches occupied by C. briggsae had multiple colonizations. Richaud et al. (2018) modeled C. elegans founder number using a Poisson distribution given the proportion of genotypes they observed at a given distance between two patches. They varied how they modeled local haplotype frequencies to account for the unknown proportions of said haplotypes in the source population and came to a mean number of 3-10 founders. Our estimate adds growing support to the hypothesis that colonization event numbers are low for many species across Caenorhabditis and that their population biology is affected by living in an ephemeral metapopulation structure. The estimates in our study and in Richaud et al. (2018) are based on androdiecious species. For gonochoristic species at least one individual of each sex must reliably colonize a patch to found a new subpopulation, assuming dispersal is achieved by prereproductive dauer individuals. Founder numbers may be higher for these species while phoresy may ensure that multiple individuals colonize a patch simultaneously. Anecdotally, however, we have on several occasions isolated unmated adult female Caenorhabditis from samples that contain no Caenorhabditis males, and males from samples that contain no females. Analogously, it has been suggested that the colonization of Réunion island by exclusively hermaphroditic Pristionchus species is a likely product of reproductive assurance (Herrmann et al., 2010).
Our data join with comparable field studies in tropical lowland sites in French Guiana and Hawaii to suggest that androdioecious species not only have larger global ranges than dioecious relatives but are also locally dominant (Table 2). Our collection efforts identified C. briggsae as the predominant species at BCI followed by C. tropicalis, as in lowland Hawaii (Crombie et al., 2019). At La Selva, C. tropicalis was the most abundant and the sole dioecious isolate was C. sp. 60. This contrasts with the findings at Nouragues, French Guiana, where C. tropicalis predominates among the androdioecious species, but the gonochoristic C. nouraguensis is the most abundant overall (Ferrari et al., 2017). Taken together, this suggests that the hypothesized benefits of self-fertile hermaphroditism, including reproductive assurance, population growth advantages, and resistance to Medea elements (Cutter et al., 2019;Noble et al., 2021), are adaptive at multiple spatial scales.

| OUTLO O K
These experiments help inform projects which could more systematically build a model of Caenorhabditis species ecology and metapopulation dynamics which includes species co-occurrence and competition, dispersal dynamics, founding numbers, and the ef-

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information can be found online in the Supporting Information section at the end of this article.  Etymology: This name is derived from kri-kudé, the word for "tree branch" in the native Panamanian language Guaymí (Alphonse, 1956).
We here describe C. krikudae n. sp. based on its distinct DNA barcode sequences, reproductive incompatibility with closely related Instead, it shows a ventrally attached tip in the light microscope.
In strain QG3050, GP1 was frequently thin or missing on one side.
In one male, GP1 was missing and GP3 was thin. Phasmids are in- side is weakly cuticularized ("velum"). The proximal third of the blade displays a conspicuous "hole." The spicule head is small and slightly offset ( Figure A2). Comparison with related species: Based on the molecular phylogeny (Figure 2), C. krikudae n. sp. is the sister species of a clade consisting of C. monodelphis and C. auriculariae. In addition, phylogenetic analyses of rRNA gene sequences place C. sonorae (Kiontke, 1997) as the sister species of C. monodelphis (Dayi et al., 2021;own unpublished data). A GP6 with a tip to the outside is also seen in the fan-less C. parvicauda Stevens and Félix (2019). However, light and SEM observations F I G U R E A 2 Comparison of stoma, male tail and spicules in species related to C. krikudae n. sp. Stoma and male tail morphology differ profoundly in these species. The phylogenetic relationships shown on the bottom are based on this study and Dayi et al. (2021). The position of C. sonorae is based on SSU rRNA sequence data only.
C. sonorae C. auriculariae C. krikudae C. monodelphis (Kiontke, 1997) showed that the tip of this GP is embedded in the bursal velum in C. sonorae, and this derived trait characterizes all other Caenorhabditis species, which form the sister group to C. parvicauda. Note that the ribosomal DNA sequences may vary within the species.

Caenorhabditis agridulce
From the available sequence data, the named species most closely related to C. agridulce n. sp. are C. quiockensis Stevens and Félix (2019) and C. dolens (Figure 2). In reciprocal laboratory crosses be-

F I G U R E A 4
Spicules (drawings) and stomata of species in the Angaria group in comparison, and phylogenetic relationships based on this study and an analysis of partial sequences of 15 protein-coding genes and rRNA genes (Karin Kiontke, unpublished analysis). The spicule tip is more similar in the respective sister species: it is slightly enlarged in C. angaria and C. castelli, narrow in C. dolens and C. quiockensis, and paw-shaped in C. agridulce n. sp. and C. sp. 8. The stoma is relatively short in all species, but only in C. angaria and C. castelli is the stegostom unusually short. In all other species, the stegostom and gymnostom contribute at least equally to the stomatal tube; the gymnostom is shorter than the stegostom in C. sp. 8. All species display flaps at the lips. The metastegostom carries a projection with a bifid tip (compare Baldwin et al., 1997) in all species. In the light microscope, this structure appears shorter in C. angaria and C. castelli than in the other species.
Gubernaculum narrow, in lateral view with a very slight bulge in the distal third; distal end is curved ventrally, distal and proximal end rounded. Median sperm area is 48 μm 2 (measured as in Vielle et al. (2016)). Their overall shape is similar in all species, but the shape of the spicule tip varies ( Figure A2): It is rounded in C. angaria and C. castelli, narrower in C. quiockensis and C. dolens but widened and pawshaped in C. agridulce n. sp. and C. sp. 8. This character distribution is consistent with the hypothesized relationships. In addition, the lips carry flaps that project into the center of the mouth. The stoma is shorter than in most other Caenorhabditis species. However, only in C. angaria and C. castelli is the mesostegostom exceptionally short.

Comparison with previously described species in the
In all other species, the stegostom and gymnostom are of about the same length, or the gymnostom is even shorter than the stegostom as in C. agridulce n. sp. All species have a bifid tooth on each side of the metastegostom. This structure is small in C. angaria and C. castelli and slightly larger in all other species.