{'fraction_non_alphanumeric': 0.033588985400110884, 'fraction_numerical': 0.01894289410460174, 'mean_word_length': 5.0125, 'pattern_counts': {'":': 1, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Loss of the HIF pathway in a widely distributed intertidal crustacean, the copepod Tigriopus californicus

Significance Oxygen availability is essential for development, growth, and viability of aerobic organisms. The genes in the hypoxia-inducible factor (HIF) pathway are considered master regulators of oxygen sensitivity and distribution inside cells, and they are hence highly conserved across animal groups. These genes are frequent targets of natural selection in organisms living in low-oxygen environments, such as high-altitude humans and birds. Here, we show that the abundant tidepool copepod Tigriopus californicus can withstand prolonged exposure to extreme oxygen deprivation, despite having secondarily lost key HIF-pathway members. Our results suggest the existence of alternative mechanisms of response to hypoxic stress in animals, and we show that genes involved in cuticle reorganization and ion transport may play a major role.


Copepod collection and maintenance
Copepods were collected from high intertidal rocky pools in Ocean Beach, San Diego, California (SD: 32.7333°N, 117.2500°W). Cultures were maintained at common garden conditions in incubators (20°C, 12hr:12hr light:dark cycle) across multiple 400-mL beakers with artificial sea water (ASW; Instant Ocean) at 35 parts-per-thousand salinity. Copepods were fed a mixture of the live microalgae Isochrysis galbana and Tetraselmis suecica and finely ground dry fish flakes (Tetramin). Once per month, water and food were changed. Populations were maintained at approximately 200-500 individuals per beaker, and beakers from each population were periodically mixed to maintain outbreeding. In order to reduce effects from recent environmental experiences and to quantify only heritable variation, copepods were maintained in laboratory cultures for a minimum of three generations before performing tolerance or RNA-seq experiments.

Anoxia tolerance experiments
We examined the effect of anoxia on developmental rate, viability, and survival in juvenile and adult copepods in multiple experiments. The effect of exposure to anoxia on developmental rate by picking well-developed egg sacs from females, placing each sac on a separate well of a 6-well plate, allowing them to hatch one day prior to treatment, and splitting each brood in two sets destined for normoxia and anoxia treatments (control normoxia: ≥ 6.5 mg O 2 /L, "anoxia": < 0.50 mg O 2 /L), allowing for a paired designed. On the day of the anoxia stress treatment, experimental conditions were achieved by bubbling nitrogen through filtered seawater in a 2-liter beaker until specific dissolved oxygen (DO) levels were attained and then filling several 20-mL glass scintillation vials with the water by submerging them while nitrogen was still being added. Oxygen levels were measured with a FiveGo Portable DO meter (Mettler-Toledo). Nauplii were then transferred quickly into the vials with minimal liquid transfer, and the vials were closed with cone urea caps, with no air gaps, and were para-filmed around the edge as well as over the top. Experimental vials were kept in incubator conditions as above for 24 hrs.
After 24 hrs, nauplii from each replicate were transferred to a well of 6-well culture plate containing normal ASW and counted under a stereo microscope. Cultures were fed live algae ad libitum and monitored daily for various stage milestones: date of first metamorphosis (date at which the first nauplius metamorphosed into copepodid juvenile), date of full metamorphosis (date at which all nauplii in the brood had metamorphosed into copepodis), date of maturity (date at which the first adult male was seen in each brood), and survival (fraction alive at day 20).
We also quantified survival in the adults aged 20 days by subjecting them to 24-hr treatments in the same way. However, since no mortality was observed (survivorship was 100% 3 days post treatment, Dataset S1A), a new experiment was conducted in which a time series of longer exposures was used to determine the time point at which mortality becomes significant. A new batch of 30 well-developed egg sacs were collected from females from stock cultures and hatched in a 400-mL beaker, and were cultured as usual for 20 days, which controls for age. On the day before treatment, 40 sets of 10 individuals, each containing 5 males and 5 females, were rinsed with fresh ASW in a petri dish and then transferred to 6-well dishes with fresh ASW and dry food. No live algae were used in order to minimize their carryover into experimental vials, which could affect oxygen content. The 40 replicates were then randomly assigned to either anoxic or control (normoxic) treatments and the experiment performed as before. At each of four time points (24,48,72, 96 hrs), copepods were removed from 5 anoxia and 5 control vials, and placed into 6-well dishes with fresh (normoxic) ASW and food, and their survival was monitored for 3 days. The fraction of individuals alive or dead in each replicate were regressed on exposure time using a probit function to estimate the median lethal time (LT 50 ).
The time series experiment above was also performed with 1-day old nauplii. Broods were picked from well-developed sacs and placed into a 24-well plate with normal ASW plus dry-food, and allowed to hatch overnight. Treatment was performed exactly as above, with 10-39 nauplii per replicate.

Search for HIF machinery
We queried the protein sequences from the T. californicus reference genome (1), as well as translated sequences from five independently assembled de novo transcriptomes from multiple populations (see Additional transcriptome assemblies below). We reasoned that, while our search via the genome sequence are still reliant on a possibly incomplete assembly, by querying multiple transcriptome assemblies that were sequenced and generated independently and from different initial tissue, we would increase the chances of detecting these two genes if they are expressed in T. californicus. Our analyses involved two approaches. In order to be classified as a "true" HIF-a, HIF-b/ARNT and EGLN, the corresponding proteins need to [1] have the known canonical domain architecture elements, and [2] group phylogenetically with previous identified gene family members. The method followed in this study is similar to what was used in Graham and Presnell (2) in regards to HIF transcription factor members.
Pfam Hidden Markov Models (HMMs) were used for identification rather than BLAST searches because HMMs are more flexible and provide full probabilistic models for detection of pattern similarities by utilizing multiple sequence alignments that can accommodate variable lengths with a focus on domain architecture (3). We used the hmmsearch command from the HMMER 3.0 program (3), along with perl scripts to identify proteins that contained a PAS domain (PF00989.24, present in HIF-a or HIF-b) or a P4HC (PF13640.5, present in EGLN).
Ultimately, we included hits that had an e-value < 0.01, which is considered a "conservative" cut- Similarly, multiple invertebrate genomes were included in the P4HC-containg protein set, but only those from humans were used for delineating gene family identities across this tree. A multiple sequence alignment for the full protein sequence for each gene set was built with MUSCLE (4) using the default parameters. Maximum Likelihood analyses were performed using RAxML (5) with bootstraps (100 replicates) using the model determined to be best fit for the alignment via ModelFinder in IQ-Tree (6, 7). The models used were VT+F+R7 for the PAS tree and PMB+F+R7 for P4HC tree. Trees were then visualized and formatted with iTOL v3 (8).
Ultimately, there were two phylogenetic trees -one with the PAS-containing sequences for identifying HIF homologs, and one with the P4HC domain-containing sequences for identifying EGLN homologs. Sequences from T. californicus that grouped with the human sequences or that formed unique clades were further examined for domain structure using InterproScan 5 (9).

Additional transcriptome assemblies
RNA-seq reads (100 bp paired-end) from three divergent populations were obtained from previous studies: Abalone Cove (AB), Santa Cruz (SCN), and San Diego (SD), all from California (10) (GenBank accessions: SRX1460514, SRX1460515, SRX1460520, SRX1460521, SRX1460522, and SRX1460523). Transcriptome assemblies for these populations were generated using CLC Genomics Workbench with parameters mismatch and insertion/deletion costs = 2, length fraction = 0.5, similarity=0.90, and conflicts resolved by majority vote. We also reads on half of a lane in a HiSeq3000 instrument. After adaptor trimming using cutadapt (11), a total of ~89 million paired reads for YA and ~84 million for SH were used for transcriptome assembly of each population using Trinity (12). Finally, we obtained predicted coding sequence and their peptide translations for all five transcriptome assemblies using TransDecoder (12), while incorporating homology BLASTP results against the UniProt/Swiss-prot database.

Transcriptional response to hypoxia
Gravid adult females (n = 24) from the SD population were picked at random from stock cultures and their egg sacs were removed and transferred to individual wells of a 6-well plate with ASW, where they hatched within 6 hrs. Each brood was raised separately and used as a biological replicate. The larvae were fed ad libitum and raised for 20 days, resulting in adult copepods of the same age across all replicates. One day prior to the anoxia treatment, the copepods were transferred to fresh ASW and fed on dry food only, rather than live algae, to minimize carryover of live algae to the experimental treatments. On the day of the exposure, anoxic and normoxic water and vials were prepared, and all individuals from each brood  adults) were then transferred quickly into the vials, as above. A total of 12 normoxia and 12 anoxia replicates were set up and placed in an incubator with typical culture conditions. Four replicates from each condition were terminated at three time points: 3 hours (3hr), 24 hours (24hr), and 24 hours after the 24-hr treatment (recovery). The anoxia time-points (3 hr, 24 hr) were chosen to reflect "acute" and "chronic" exposures that are relevant for comparisons to assays in other systems (13)(14)(15)(16), and our longer time-course experiment above revealed that 24 hrs of anoxia is still well within survival limits. Immediately after each treatment time point, all individuals from each replicate were transferred into 1.5 mL tubes and immersed in 350 µL of Trizol. For the recovery replicates, individuals were taken out of the vials after 24 hrs and transferred to normal culture plates with fresh (normoxic) ASW and dry food for another 24 hours, at which point they were sacrificed in Trizol. Reads were trimmed of adapter sequences, then filtered by length and quality using cutadapt (11), and only reads with length ≥ 50 nt and with an average Phred score ≥ 30 were retained for further analysis. Each library was mapped using Salmon (17) against the reference transcriptome generated from the SD genome (1). Differential expression was determined using edgeR (18), where libraries were normalized with trimmed mean of M-values method (TMM; (19) . Only transcripts with at least one mapped read per million in at least four libraries were included in the analysis. To identify genes differentially expressed in each time point, we compared anoxia samples for each stage to their respective normoxia control group. Transcripts with an FDR-adjusted P-value ≤ 0.10 were considered to be differentially expressed. Specific commands and scripts utilized in these analyses are available at https://github.com/amgraham07. Gene Ontology (GO) annotations were retrieved using topGO, and enrichment of GO terms within significantly upregulated genes and significantly downregulated genes, relative to the rest of the genes, was assessed with a threshold at FDR of 0.05. Summary of sequencing effort and mapping results can be found on Dataset S2A.

Quantitative PCR validation
We used quantitative reverse transcription PCR (RT-qPCR) to validate differential expression observed in our RNA-seq results. We selected 16 genes that encompass the main patterns detected in our results (Fig. 3) and designed primers based on the genomic transcripts that were used as mapping reference in the RNA-seq analyses (Table S4).
For RT-qPCR validation, the anoxia experiment was repeated in its entirety with new individuals, which were raised together to age 20 days, and then subjected to anoxia or normoxia treatments in the same way as above (see Transcriptional response to hypoxia). We sampled copepods at the 24-hr anoxia/normoxia and recovery time points, but not at the 3-hr time point, since the vast majority of differential expression in our RNA-seq results occurred between the two latter time points. Each treatment had five replicates each containing 15-18 adult copepods, for a total of 20 replicates. RNA isolation, DNase treatment, and RNA purification were performed as above. First-strand cDNA was synthesized from each replicate with 300 ng of total RNA in 20-µL reactions using the High-Capacity RNA-to-cDNA kit (Applied Biosystems).
Before quantifying gene expression, we first determined the optimal annealing temperature for each primer pair by performing gradient PCR reactions with a single sample, and then estimated primer efficiencies by performing qPCR reactions of a dilution series of cDNA, in duplicate.
Efficiencies were determined to be within the range recommended for quantitative PCR (90-110%, Table S4; (20,21). In addition to the 16 target genes, we optimized PCR conditions of two stable reference genes (RpS21: ribosomal protein S21, and ETIF3: eukaryotic translation initiation factor 3) to be used for normalization. For each sample, technical duplicates for each target gene were averaged and the geometric mean of RpS21 and ETIF3 was used for normalization. The normalized expression for each target gene was calculated as 2 -∆Ct , where ∆Ct = Ct target gene -Ct reference genes . Treatment groups were then compared via t-tests, and fold-change estimated relative to the mean of normalized expression of the respective control group (Fig. S2).

Identification of putative hypoxia response elements
Promoter regions (1 kb upstream) of every gene in the T. californicus genome were extracted and searched for the presence of an HRE, based on its motif: (22). This search was performed using a combination of Bioconductor R packages (BSgenome, Biostrings, GenomicFeatures) and bedtools (23). Scripts used for this analysis are available on GitHub (https://github.com/amgraham07). Among all genes annotated with GO terms (n = 10,727), we tested whether GO terms associated with oxygen homeostasis ("heme-", "hemato-", "angio-", "oxygen" and "hypoxia) were significantly overrepresented in genes containing a putative HRE. In addition, potential enrichment of HREs was assessed in the promoters of the set of genes differentially expressed at the 24-hr time point relative to the frequency of HREs among genes that did not change significantly during that exposure time.

TFBS motif discovery and enrichment
In lieu of the presence and activity of the HIF pathway transcription factors in T. californicus, the possible involvement of other transcription factors was examined by assessing the presence of potential of transcription factor binding site (TFBS) motifs, using programs designed for motif discovery and enrichment. Using the MEME suite (24), we searched for potential TFBS motifs in the promoter regions of genes that were significantly up-or downregulated in the 24-hr treatment. First, the presence of known motifs from the JASPAR database (25) (core_insect_nonredundant) were queried against the promoter regions using AME (Analysis of Motif Enrichment). Motifs that were considered significantly enriched were ones that passed an optimal enrichment P-value of the motif according to a Fisher's exact test, adjusted for multiple tests using a Bonferroni correction. Second, the presence of "novel" motifs were assayed using MEME, with the potential identity of the motifs evaluated using Tomtom (26).

Assessment of co-expression network analysis
Analysis of gene co-expression networks is a powerful approach for examining high-level patterns of putative gene interactions, and it has become useful in revealing patterns of coexpression across physiological treatments using RNA-seq data (27)(28)(29). Estimating coexpression networks during hypoxia in T. californicus is likely a promising approach to gain initial insights into what 'hub' genes may be involved in response, given the absence of the master regulatory machinery in this species.
For our current data set, however, this type of analysis is not appropriate because of the non-ideal experimental size and small sample size. The primary goal of a co-expression analysis is to detect strong correlations in expression among genes, without requiring that these genes are differentially expressed. However, the robustness of these correlations can be sensitive to samples size. In other words, the composition of each of gene network module can change drastically (i.e. be less robust) when sample size is low. For example, Iancu, et al. (30) showed empirically that a minimum of 35 replicates are necessary to maintain strong correlation among genes; decreasing sample size much lower than that resulted in large reorganization of the gene associations, which indicate spurious correlations.
We assessed the use of the Weighted Gene Coexpression Network Analysis (WGCNA) model (31) to explore our current dataset. As input, we used the variance-stabilized gene expression levels produced with edgeR during our main analysis and we followed the WGCNA package instructions closely. Briefly, we found that ~10,000 genes in our data set were distributed among 13 network modules (Fig. S6). When tested for association between module expression and our hypoxia treatments, we found no significant correlation between any module and the change in hypoxia treatment (Fig. S7). We also briefly assessed the robustness of these module groupings by downsampling our data set from 24 replicates to 20 replicates, by randomly selecting 4 replicates to be excluded. We created two different sets of downsampled data, and repeated the co-expression analysis with each set. We then compared module membership between the full data set above and each of the two downsampled data sets. We found that, even though the number of modules detected was similar among all data sets (13 or 14 modules), the co-occurrence (or co-assignment) of genes across data sets changed drastically, often by more than 50-70%, meaning that the groups of genes that were detected to be co-expressed in one data set were not detected to be co-expressed in each subsequent downsampled data set. Instead, the co-expression modules were "shuffled" with regards to their gene content. This suggests that our full data set of 24 replicates is still too low to produce robust meaningful patterns.
The authors of this model recommend, in general, a minimum of 20 samples, and that this minimum should be increased in data that are from highly heterogeneous samples, such as our outbred, mixed tissue, RNA samples. Most published studies that show strong statistical power in these analyses follow these recommendations and use N >30-50 replicates, and many with N >100 replicates (32,33), especially studies of samples from natural populations (Rose et al. 2015). With regards to hypoxia response, a useful recent reference that shows robust patterns of co-expression across hypoxia timepoints is Rand, Mossman, Zhu, Biancani and Ge (27). In that study, a total of 35 RNA-seq replicates were used, with 24 total hypoxia replicates (12 per hypoxia time point).
We intend to follow statistical recommendations for increasing our sample size and experimental design in future experiments aimed specifically at quantifying co-expression networks associated with hypoxia response in T. californicus.     Figure S2. RT-qPCR validation of RNA-seq results. Shown are relative expression levels of 16 genes after 24 hrs of anoxia (A) and after 24 hrs of recovery post anoxia (B). This validation experiment used an entirely new set of individuals and was repeated identically to the original RNA-seq experiment. Fold change was calculated as the ratio of each anoxia replicate (n = 5) to the mean of the respective normoxia control group (n = 5). Expression levels for RNA-seq replicates were estimated from the number of Counts Per Million mapped reads (CPM), while those for qPCR were based on threshold Ct values normalized to two internal stable reference genes (RpS21 and ETIF3). The list of gene names and abbreviations can be found on Table S10. *Statistical tests for change in expression between anoxia and control samples were quasilikelihood F-tests for RNA-seq and t-tests for qPCR.  Figure S3. Motif discovery (MEME) results for upstream regulatory regions of significantly differentially expressed genes at 24-hrs. The motif associated with hb (top) represents the best match as determined by Tomtom against a motif found to be enriched from the MEME analysis (bottom).   Figure S7. Tests of correlation between hypoxia treatment levels (oxygen level and time point) and each co-expression module. In each cell, upper value shows the correlation coefficient r and in parenthesis is the P value for the correlation test.