Diatoms serve as the major link between the marine carbon (C) and silicon (Si) biogeochemical cycles through their contributions to primary productivity and requirement for Si during cell wall formation. Although several culture-based studies have investigated the molecular response of diatoms to Si and nitrogen (N) starvation and replenishment, diatom silicon metabolism has been understudied in natural populations. A series of deckboard Si-amendment incubations were conducted using surface water collected in the California Upwelling Zone near Monterey Bay. Steep concentration gradients in macronutrients in the surface ocean coupled with substantial N and Si utilization led to communities with distinctly different macronutrient states: replete (‘healthy’), low N (‘N-stressed’), and low N and Si (‘N- and Si-stressed’). Biogeochemical measurements of Si uptake combined with metatranscriptomic analysis of communities incubated with and without added Si were used to explore the underlying molecular response of diatom communities to different macronutrient availability. Metatranscriptomic analysis revealed that N-stressed communities exhibited dynamic shifts in N and C transcriptional patterns suggestive of compromised metabolism. Expression patterns in communities experiencing both N and Si stress imply that the presence of Si stress may partially ameliorate N stress and dampen the impact on organic matter metabolism. This response builds upon previous observations that the regulation of C and N metabolism is decoupled from Si limitation status, where Si stress allows the cell to optimize the metabolic machinery necessary to respond to episodic pulses of nutrients. Several well-characterized Si-metabolism associated genes were found to be poor molecular markers of Si physiological status; however, several uncharacterized Si-responsive genes were revealed to be potential indicators of Si stress or silica production. Methods Metatranscriptome analysis Seawater samples for metatranscriptome analysis were collected by filtration at <5 psi onto 47 mm, 1.2 µm pore-size polycarbonate membranes, flash frozen in liquid N2, and stored at -80°C. Total RNA was extracted from filters using the Trizol-RNeasy method with an additional bead beating step. In brief, filters were transferred to tubes containing RNase/DNase free 100 µm zirconia/silica beads with 1 mL of Trizol reagent, vortexed for 2 min, incubated at room temperature for 5 min, vortexed again for 2 min before following the standard Trizol RNA extraction protocol. The upper aqueous phase was transferred to a clean 1.5 mL tube with an equal volume of 70% alcohol, mixed, and added to QIAgen RNAeasy column. The standard RNeasy protocol with the optional on-column DNase treatment was followed. RNA was quantified using a Qubit® RNA HS Assay, and RNA integrity was determined using the Agilent Bioanalyzer RNA 6000 Pico kit eukaryotic assay. Total RNA (500 ng) was used for library prep with Illumina TruSeq RNA Library Prep kit v2 with mRNA purified using Oligo-dT bead capture polyA tails. Library quality was assessed on an Agilent Bioanalyzer and quantified via Qubit before being pooled in equal quantities. The pooled libraries were sequenced at UC-Davis Genome Center on HiSeq 4000 in lanes of paired-end-150bp flow cell. Primer and adapter sequences were removed from raw sequencing reads with trimmomatic v0.38 (Bolger et al., 2014). Paired-end reads were merged with interleave-reads.py from the khmer package v2.1.2 before removal of rRNA reads with sortmeRNA v2.0 utilizing the built in SILVA 16s, 18s, 23s, and 28s databases (Kopylova et al., 2012). Remaining reads were separated (using the deinterleave_fastq.sh https://gist.github.com/nathanhaigh/3521724) and assembled into contiguous sequences (contigs) with megahit v1.1.3 using ‘meta-large’ setting (Li et al., 2015). Open reading frames (ORFs) were predicted from assembled contigs using FragGeneScan v1.31 (Rho et al., 2010), and orfs shorter than 150 bp were removed. Reads counts for remaining orfs were estimated using salmon v0.6.0 `quant` quasi-mapping with seqBias and gcBias features. Orfs were annotated based on best homology (lowest E-value) using BLASTP v.2.5.0+ with an E-value threshold of <10-3. Taxonomic and functional annotations were assigned using the MarineRefII (roseobase.org/data), a custom databased maintained by the Moran Lab at the University of Georgia that includes sequences from PhyloDB v.1.076 which contains 24,509,327 peptides from 19,962 viral, 230 archaeal, 4,910 bacterial and 894 eukaryotic taxa, including peptides from KEGG, GenBank, JGI, ENSEMBL, CAMERA and 410 taxa from the Marine Microbial Eukaryotic Transcriptome Sequencing Project (MMETSP), as well as the associated KEGG functional and NCBI taxonomic annotations through a related SQL database. NCBI taxonomic annotations were further curated (Cohen et al., 2017); (https://github.com/marchettilab/metatranscriptomicsPipeline) to ensure the use of consistent taxonomic ranks. KEGG Ortholog (KO) classifications of diatom urea transporters, nitrite reductase, ammonium transporters were manually verified against known gene phylogenies and edited accordingly (Smith et al., 2019). Silaffin-like genes were identified in Pseudo-nitzschia multiseries (CLN-47 v1) genome, as described previously (Scheffel et al., 2011), using an initial signal peptide screen of translated nucleotide sequences and subsequence sliding window search for amino acid (AA) sequences containing a 100-2000 long segment with ≥10% lysine and ≥18% serine residues. Full length P. multiseries putative silaffin amino acid sequences were subsequently screened for the presence the characteristic silaffin-like motifs using the FIMO program within the Multiple EM for Motif Elicitation (MEME suite v.5.4.1) motif scanning software package (Grant et al., 2011). Additional functional annotation was performed using BLASTP with the NCBI nr database (E-value <10-10) and KEGG GhostKOALA (score <40). This process identified 157 putative unique P. multiseries silaffin-like genes (PNSLs) which were used with BLASTP+ v2.5.0 (E-value <10-5) to identify homologous sequences within our metatranscriptome assembly. To identify the distribution of these genes in diatoms and other microalgae, BLASTP was used against a database containing the 67 translated transcriptomes from the MMETSP (Keeling et al., 2014). PNSLswith no significant hits (E-value <10-10) to non-diatom genes were considered “diatom-specific”. For genes of interest that did not have an associated KO number, such as silicon transporters (SITs) and silicanin-1 (Sin-1), manual annotation was performed as previously described (Durkin et al., 2016; Brembu et al., 2017a; Kotzsch et al., 2017) using BLASTP+ v2.5.0 best hit with e-value cutoff of 10-5. Putative SIT sequences were further classified to clade designations, using a maximum-likelihood tree constructed from a reference alignment (Durkin et al. 2016) with RAxML version 8.2.12 – PROTGAMMAWAGF substitution model and 100 bootstrap replicates (Stamatakis, 2014). A HMM-profile was constructed from the reference alignment (using hmmbuild 3.2.1) and used with hmmalign to align amino acid sequences corresponding to SIT orfs. Alignment and tree files were packaged using taxtastic v.0.8.3, and SIT orfs were placed on the reference tree using pplacer v.1.1.alpha19 with posterior probability calculated (Matsen et al., 2010). The most closely related reference sequence was assigned to each SIT orf using guppy v.1.1.alpha19 (Matsen et al., 2010). Prior to differential abundance analysis, raw counts were aggregated within diatom genera by KO number. For genes lacking a KO number, e.g. SITs, ISIPs, Sin1, etc., aggregation was done based on gene assignment through KEGG annotation or BLASTX query of supplemental databases. Genus-specific aggregation of functionally annotated reads reduces redundancy and allows the use of tools originally designed for single organism RNAseq analysis (e.g. DESeq2, edgeR). This is necessary for microbial community transcriptomic analysis due to methodological and computational difficulties in resolving species-level, differential transcript abundance (Toseland et al., 2014; Alexander et al., 2015; Kopf et al., 2015; Cohen et al., 2017, 2021; Hu et al., 2018; Lampe et al., 2018). However, this approach does not resolve species-level contig abundance, nor does it a priori resolve clade-specific abundance patterns for genes that belong to multigene families. In the case of SITs, individual clades were identified through phylogenetic analysis as described above, and manually annotated to allow interrogation of SIT transcripts at the clade level. Additionally, annotations for urea transporters 1 and 2 were grouped together because of their close phylogenic relationship and shared abundance pattern in response to nitrogen supply (Smith et al., 2019). A similar approach for other multigene families, such as ammonium transporters and nitrate transporters, was not used because abundance of these genes in response to nitrogen supply does not appear to be clade-specific (Smith et al., 2019). Metabolic pathway groupings Genes indicative of nitrogen limitation were selected based on culture-based studies characterizing the response of the diatom Phaeodactylum tricornutum to depletion and replenishment of nitrate, nitrite, and ammonia (Smith et al., 2019). Using independent transcriptomic studies of Si-depleted and replenished Thalassiosira pseudonana cultures (Shrestha et al., 2012; Smith et al., 2016a; Brembu et al., 2017a), we identified a subset of genes conserved across those studies in either their response to silicon starvation or association with silica production and denote these here as SiLA (Si-limitation associated) or SiPA (silica production associated), respectively. SiLA genes were identified as those with a >2-fold increase in transcript

We consider a neighborhood random walk on a quadrant {(X1(t),X2(t),φ(t)):t≥0} with environment phase variable φ(t) modeled by a continuous-time Markov chain with φ(t)∈Snm when X1(t) = n, X2(t) = m. We describe this random walk using a two-dimensional level-dependent Quasi-Birth-and-Death process (2D-LD-QBD) with phase variable φ(t) and level variables X1(t),X2(t)∈{0,1,2,…} which change in a skip-free manner at the moments of jump in the process. We transform this random walk into a one-dimensional LD-QBD {(Z(t),χ(t)):t≥0} with level variable Z(t)∈{0,1,2,…} recording the maximum of the two level variables and phase variable χ(t)=(χ1(t),χ2(t),φ(t)) recording the remaining information about the random walk. Using this transformation, we perform transient and stationary analysis of the random walk, including first hitting times for various sample paths, using matrix-analytic methods. We also construct a sequence of neighborhood random walks, represented as two-dimensional QBDs ({(X1(k)(t),X2(k)(t),φ(t)):t≥0})k=1,2,…, converging in distribution to a two-dimensional stochastic fluid model (SFM) {(Y1(t),Y2(t),φ(t)):t≥0}, which describes a movement on a quadrant in which the position changes in a continuous manner according to rates dY1(t)/dt=c1,φ(t) and dY2(t)/dt=c2,φ(t) modulated by the underlying phase process {φ(t):t≥0}. Numerical examples are provided to illustrate the application of the methodology.

CI = Confidence Interval CXR = Chest radiograph.*Where the design effect was ≤1 CI’s were not adjusted for cluster design but binomial exact CI presented.†AUC = Area under the receiver operating characteristic curve.§Denominator for HIV-positive n = 52, HIV-negative n = 49, HIV-unknown n = 48.**3 cases did not have a CXR, 1HIV+, 1HIV-, 1HIVunknown, so denominators are 120, 51, 48 and 21 respectively. For specificity: 1347 missing records.