\documentclass{article}
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
\usepackage{url}
\usepackage{hyperref}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}

\onecolumn
\firstpage{1}

\title[Non-cyanobacterial soil crust diazotrophs]{Non-cyanobacterial
diazotrophs dominate dinitrogen fixation in biological soil crusts during early
crust formation.} 
\author[\firstAuthorLast]{\Authors}
\address{}
\correspondance{}
\extraAuth{}% If there are more than 1 corresponding author, comment this line and uncomment the next one.
%\extraAuth{corresponding Author2 \\ Laboratory X2, Institute X2, Department X2, Organization X2, Street X2, City X2 , State XX2 (only USA, Canada and Australia), Zip Code2, X2 Country X2, email2@uni2.edu}
\topic{}% If your article is part of a Research Topic, please indicate here which.

\thispagestyle{empty}
\maketitle
\blfootnote{\scriptsize Running title: Non-cyanobacterial soil crust diazotrophs}
\blfootnote{\scriptsize Subject category: Microbial ecology and functional diversity of natural habitats}

%\bibliographystyle{plainnat}



\begin{document}

\title{Non-cyanobacterial diazotrophs mediate dinitrogen fixation in biological soil
crusts during early crust formation}


\author{Charles Pepe-Ranney}
\affil{Affiliation not available}

  
\author{Chantal Koechli}
\affil{Affiliation not available}
  
\author{Erin Eggleston}
\affil{Affiliation not available}
  
\author{Dan Buckley}
\affil{Affiliation not available}
  



\date{\today}

\maketitle




\section{Introduction} Biological soil crusts (BSC) are specialized microbial
communities that form at the soil surface in arid environments and they fill
a variety of important ecological functions. BSCs occupy plant interspaces and
cover a wide, global geographic range \citep{garcia2003estimates}\,. For
example, in some regions on the Colorado Plateau BSCs cover 80\% of the ground
\citep{karnieli2001}\,. The global biomass of BSC cyanobacteria alone is
estimated at 54 x 10$^{12}$ g C \citep{garcia2003estimates}\,. BSD nitrogen
fixation (N$_{2}$-fixation) is responsible for significant input of nitrogen
(N) to arid environments \citep{Evans_1999,belnap2003}\,. Interestingly, much
of this fixed N is exported from the crusts in dissolved form through
percolation or runoff and little is lost to volatilization
\citep{Johnson_2007}\,. The presence of BSC is positively correlated with
vascular plant survival due in part to N inputs from BSC (for review of
BSC-vascular plant interactions see \citealt{BelnapVascular}). These microbial
ecosystems are not immune to climate change and changes in precipitation and
temperature could alter BSC microbial community structure/membership and
possibly BSC diazotroph diversity and N$_{2}$-fixation
\citep{Garcia_Pichel_2013}\,.

BSC are highly susceptible to natural and anthropogenic disturbance
\citep{Garcia_Pichel_2013}. Succession in BSC communities is characterized by
transition from early successional "light" crusts to mature "dark" crusts
\citep{Yeager_2004,Belnap_2002}\,. Motile non-heterocystous
cyanobacteria(e.g.\textit{Microcoleus vaginatus} or \textit{M. steenstrupii}),
which cannot fix N$_{2}$ are pioneer colonizers of early successional crusts
and are abundant in all types of BSCs \citep{Garcia_Pichel_2013}\,.
Successional development of mature crust is accompanied by a change in color
produced by secondary colonization with non-motil N$_{2}$-fixing heterocystous
cyanobacteria which produce sunscreen compounds that reduce soil albedo
\Citep{Belnap_2002,Yeager_2004}\,. These heterocystous cyanobacteria (e.g.
\textit{Scytonema}, \textit{Spirirestis}, and \textit{Nostoc}) increase in
abundance during crust development and are more abundant in mature crusts
\citep{Yeager,Yeager_2012}. Heterocystous cyanobacteria are numerically
dominant in surveys of BSC \textit{nifH} gene diversity \citep{Yeager,
Yeager_2004, Yeager_2012}\,. For example, 89 percent of 693 \textit{nifH}
sequences derived from Colorado Plateau and New Mexico BSC were attributed to
heterocystous cyanobacteria \citep{Yeager}\,. Other BSC \textit{nifH} sequences
are attributed to \textit{Alpha-}, \textit{Beta-}, and
\textit{Gammaproteobacteria}, as well as a \textit{nifH} clade (\textit{nifH}
cluster III) that includes diverse anaerobes such as clostridia, sulfate
reducing bacteria, and anoxygenic phototrophs \citep{Steppe_1996, Yeager}\,.

Two lines of evidence suggest that nitrogen fixers other than phototrophs are
important in early-successional crusts.  First, the contributions of early
successional BSC to N$_{2}$-fixation in arid ecosystems may have been
systematically underestimated. The high abundance of heterocystous
cyanobacteria at the surface of mature crusts, where acetylene reduction assay
rates are often maximal, is generally taken as evidence that BSC
N$_{2}$-fixation occurs primarily in mature crusts and is dominated by
heterocystous cyanobacteria. However, rates of BSC N$_{2}$-fixation are
typically determined by areal measurements made at the crust surface with the
acetylene reduction assay and vary significantly across samples and studies
\citep{Evans_2001}\,. The reasons for inter-site and inter-study variability
are complex and likely include the spatial heterogeneity of BSC
\citep{Evans_2001}\,.  The acetylene reduction assay is also subject to
methodological artifacts that can complicate comparisons between samples that
differ in their physical and biological characteristics (see
\citealt{Belnap_2001} for review). In particular, N$_{2}$-fixation in early
successional BSC is maximal below the crust surface \citep{15643930} and hence
diffusional limitation (of both acetylene and ethylene) across the crust
surface can cause severe underestimates if they do not allow for sufficiently
long incubation times \citep{15643930}\,.  If BSC N$_{2}$-fixation is instead
estimated by integrating rates across a depth profile (which eliminates
constraints from diffusional limitation), then total rates of N$_{2}$-fixation
do not differ significantly between early successional and mature BSC
\citep{15643930}\,. This result suggests that diazotrophs other than
heterocystous cyanobacteria may be important contributors to N$_{2}$-fixation
in early successional BSC communities as early successional BSC possess few
heterocystous cyanobacteria and these are present near the crust surface.
Second, the bare soils that are colonized during the process of early crust
formation are unconsolidated and oligotrophic in many respects, with much lower
N content than adjacent crusts \citep{BERALDI_CAMPESI_2009}\,, and the
cyanobacteria that are typical colonization pioneers (\textit{Microcoleus}
spp., \citealt{Garcia_Pichel_2009})\,, are unable to fix nitrogen as they lack
that genetic capacity \citep{Starkenburg_2011, Rajeev_2013}\,. 

To determine the agency of nitrogen fixation in early developmental
crusts, we conducted $^{15}$N$_{2}$ DNA stable isotope probing (DNA-SIP)
experiments with early successional Colorado Plateau BSC conspicuously devoid
of significant surface populations of heterocystous cyanobacteria.  DNA-SIP
with $^{15}$N$_{2}$ has not been previously attempted with BSC. DNA-SIP
provides an accounting of \textit{active} diazotrophs on the basis of
$^{15}$N$_{2}$ assimilation into DNA whereas \textit{nifH} clone libraries
merely account for microbes with the genomic potential for N$_{2}$-fixation.
Further, we investigate the distribution of these active diazotrophs in surveys
of microbial diversity conducted on BSC over a range of spatial scales and soil
types \citep{Garcia_Pichel_2013, Steven_2013}\,. 


\section{Results} 
\subsection{DNA buoyant density changes in response to $^{15}$N$_{2}$} BSCs
were wetted and incubated for 4 days in transparent chambers with headspace
containing N$_{2}$ either from air or from 100 percent atom enriched
$^{15}$N$_{2}$. The chambers were illuminated with 16 h on / 8 h off cycles at
an intensity of 200 $\mu$mol photons m$^{-2}$ s$^{-1}$, which is the equivalent
of an overcast/rainy day. N$_{2}$-fixation as measured by acetylene reduction
increased from 4.8 $\mu$moles m$^{-2}$ d$^{-1}$ on day 2 to 14.5 $\mu$moles
m$^{-2}$ d$^{-1}$ on day 4. Amplicon sequences from $^{15}$N$_{2}$-labeled
samples and their corresponding unlabeled controls diverged specifically in
heavy gradient fractions (Figure~\ref*{fig:ord_heavy} and
Figure~\ref*{fig:ordination}) as assessed by Bray-Curtis dissimilarity
\citep{Bray_1957}\,, and this result was significant (Adonis test
\citep{Anderson_2001}\,; p-value: 0.001, r$^{2}$: 0.18). 

\subsection{OTUs responsive to $^{15}$N$_{2}$ are primarily
\textit{Proteobacteria} and \textit{Clostridiaceae}} 
OTUs that incorporated $^{15}$N into their DNA were detected by a differential
change in their abundance within heavy gradient fractions of
$^{15}$N$_{2}$-labeled samples relative to corresponding controls. A total
of 2,127 and 2,160 OTUs were detected in days 2 and 4, respectively, and
these OTUs were interrogated for evidence of $^{15}$N$_{2}$-labelling. Of
these OTUs, only 499 and 563, respectively, passed a sparsity threshold
applied to filter out OTUs with insufficient data for statistical analysis
(see \citet{Love_2014} for discussion of independent filtering). Of OTUs
passing the sparsity criterion, 34 were enriched significantly in heavy
fractions relative to control and this result is specifically expected for
OTUs that have $^{15}$N-labeled DNA (\textit{i.e.} $^{15}$N$_{2}$
``responders''). Of these, 19 are annotated as \textit{Firmicutes}, 12
as \textit{Proteobacteria}, 2 as \textit{Actinobacteria} and 1 as
\textit{Gemmatimonadetes} (Figure~\ref*{fig:l2fc},
Figure~\ref*{fig:trees}). If the responder OTUs are ranked by descending
enrichment in heavy gradient fractions versus control, 8 the top 10 responders
(\textit{i.e.} those most enriched in the heavy fractions of labeled gradients)
are either \textit{Firmicutes} (3 OTUs) or \textit{Proteobacteria} (5 OTUs)
(Figure~\ref*{fig:scatter_heavy}). Centroids (seed sequences) for strongly
responding \textit{Proteobacteria} OTUs all share high SSU rRNA gene sequence
identity ($>$98.48\%, Table~\ref*{tab:LTP_blast}) with isolates from genera
known to possess diazotrophs including \textit{Pseudomonas},
\textit{Klebsiella}, \textit{Shigella}, and \textit{Ideonella}. None of the
\textit{Firmicutes} OTU centroids in the top 10 responders share greater than
97\% SSU rRNA gene sequence identity with sequences in the Living Tree Project
(LTP) database of 16S rRNA gene sequences from type strains (release 115) (see
Table~\ref*{tab:LTP_blast}). OTUs that passed the sparsity threshold but were
not classified as $^{15}$N-responsive were subsequently tested with the null
hypothesis that the OTU fold enrichment in labeled gradient heavy fractions
versus control was above the selected threshold. Rejecting the second null
indicates that an OTU did \textit{not} incorporate $^{15}$N into biomass. There
were 86 and 89 ``non-responders'' at days 2 and 4, respectively. The $^{15}$N
labelling of OTUs that did not pass sparsity or could not be classified as
either a responder or non-responder cannot be determined conclusively.

Although we did not take replicate samples within a time point we can assess
the consistency of each OTUs response across the two time points. OTU fold
enrichment at day 2 and 4 was consistent (Figure~\ref*{fig:2v4}). There was
a significant correlation between OTU fold enrichment at day 2 versus day
4 (P-value 4.35e$^{-8}$). When the enrichment at day 2 is compared to day 4 via
an interaction term (day $*$ label/control, see methods), we found only two
OTUs had significantly different enrichment between time points ("OTU.227" and
"OTU.4037", Table~\ref*{tab:LTP_blast}). In addition, when day
2 and day 4 samples are treated as replicates (see methods) only five of the
OTUs we identified as responders OTUs were not significantly enriched in
labeled gradient heavy fractions versus control ("OTU.140", "OTU.4037",
"OTU.227", "OTU.137", and "OTU.263", Table~\ref*{tab:LTP_blast}). The labeling
of these OTUs should be interpreted with caution. None of the top 10 strongest
responding OTUs showed inconsistent enrichment across time points
based on the above analyses. Further, confidence in enrichment (i.e. lowest
enrichment P-values between Day 2 and Day 4) appears to be correlated with
consistency in response across both days (Figure~\ref*{fig:2v4}).

\subsection{$^{15}$N-responsive OTUs are found in low abundance in
    available environmental BSC SSU rRNA gene surveys} In total 13 of the
34 $^{15}$N-responsive OTUs have been observed previously in SSU rRNA gene
surveys of BSC communities (Figure~\ref*{fig:trees},
Figure~\ref*{fig:rspndr_dist}). Eleven of the 19 $^{15}$N-responsive
\textit{Firmicutes} OTUs are members of the \textit{Clostridiaceae}. Three
$^{15}$N-responsive \textit{Clostridiaceae} have been observed in previous
BSC SSU rRNA gene surveys.  Two $^{15}$N-responsive
\textit{Clostridiaceae} were found in "light" (\textit{i.e.} early
successional) crust during SSU rRNA gene sequence analysis of BSC
\citep{Garcia_Pichel_2013}\,, and one $^{15}$N-responsive
\textit{Clostridiaceae} OTU was found among the "below crust" BSC SSU rRNA
gene sequences described by \citet{Steven_2013} (Figure~\ref*{fig:trees}).
Five $^{15}$N-responsive proteobacterial OTUs (Table~\ref*{tab:LTP_blast})
were detected previously in BSC samples
\citep{Garcia_Pichel_2013,Steven_2013}.  The $^{15}$N-responsive
\textit{Gemmatimonadetes} OTU was observed in four
\citet{Steven_2013} samples and one $^{15}$N-responsive
\textit{Actinobacteria} OTU was found in three \citet{Steven_2013}
samples. \textit{Gemmatimonadetes} and \textit{Actinobacteria}
$^{15}$N-responsive OTUs were not observed in samples collected by
\citet{Garcia_Pichel_2013}.

\subsection{Comparison of SSU rRNA gene sequences from different BSC samples}
We compared the SSU rRNA gene sequences determined in this DNA-SIP experiment
with two previous surveys of SSU rRNA gene amplicons from BSC communities
\citep{Garcia_Pichel_2013, Steven_2013}\,. There were 3,079 OTUs
(209,354 total sequences after quality control) in the DNA-SIP data, 3,203 OTUs
(129,033 total sequences after quality control) in the
\citet{Garcia_Pichel_2013} study, and 2,481 OTUs (129,358 total sequences after
quality control) in the \citet{Steven_2013} study with a total of 4,340 OTUs
across all three datasets. Of the total 4,340 OTU centroids established
for this study, 445 have matches in the Living Tree Project (LTP) (a collection
of SSU rRNA gene sequences for all sequenced type strains \citep{Yarza_2008})
at or above a threshold of 97\% sequence identity (LTP version 115). That is,
445 of 4,340 OTUs are closely related to known isolates. The DNA-SIP data
shares 56\% OTUs with the \citet{Steven_2013} data and 46\% of OTUs with
the \citet{Garcia_Pichel_2013} data, while these latter two studies share
46\% of their OTUs. This result suggests that low frequency OTUs likely
remain undersampled in all datasets. 

Sequencing of DNA subjected to CsCl fractionation is expected to sample
a different subset of diversity than that sampled by sequencing of
unfractionated bulk DNA. For example, SIP enhances detection of OTUs that
incorporate $^{15}$N into their DNA, and these OTUs will be overrepresented in
the overall DNA-SIP sequence pool relative to their relative abundance in
unfractionated bulk community samples. In addition, the DNA-SIP sequencing
effort was directed at a relatively small number of "light" crust samples (n
= 4), while previous sequencing efforts \citep{Garcia_Pichel_2013, Steven_2013}
were spread across hundreds of samples from both "light" and "dark" crusts.
Hence, it is likely that the current study will be more likely to detect rare
OTUs present in early successional "light" crust communities, particularly
those that incorporate $^{15}$N into DNA. In all three BSC studies, most
sequences were annotated as either cyanobacteria or \textit{Proteobacteria},
though only in the DNA-SIP data did the sequences of \textit{Proteobacteria}
outnumber those of cyanobacteria. \textit{Proteobacteria} represented 29.8\% of
sequence annotations in DNA-SIP data as opposed to 17.8\% and 19.2\% for the
\citet{Garcia_Pichel_2013} and \citet{Steven_2013} data, respectively. In
addition, sequences annotated as \textit{Firmicutes} were more abundant in the
DNA-SIP data (19\%) than in the data from \citet{Steven_2013} and
\citet{Garcia_Pichel_2013} (0.21\% and 0.23\%, respectively)
(Figure~\ref*{fig:study_phy_dist}). Finally, and congruently with sampling
design sequences annotated to "Subsection IV" of cyanobacteria, which
encompasses the heterocystous cyanobacteria in the Silva taxonomic nomenclature
\citep{17947321}\,, comprised only 0.29\% of cyanobacteria sequences in the
DNA-SIP data while representing 15\%  and 23\% of cyanobacteria sequences from
the \citet{Steven_2013} and \citet{Garcia_Pichel_2013} data, respectively.


\section{Discussion}
% Fakesubsubsection:BSC N-fixation has long been attributed
BSC N-fixation has long been attributed to heterocystous
cyanobacteria and the preponderance of cyanobacterial \textit{nifH} genes
observed in molecular surveys of BSCs have generally supported this hypothesis
\citep{Yeager,Yeager_2004,Yeager_2012}\,. However, in this study
$^{15}$N$_{2}$-DNA-SIP reveals that non-cyanobacterial microorganisms fix
N$_{2}$ in early successional BSC samples. \textit{Proteobacteria} and
\textit{Clostridiaceae} were most abundant among $^{15}$N$_{2}$-responsive OTUs
as revealed by a robust statistical framework for quantifying and evaluating
differential OTU abundance in microbiome studies \citep{24699258,Love_2014}\,.
Many of these OTUs (about 40\%) have been observed previously in BSC
communities. Rarefaction curves of data from \citet{Steven_2013} and
\citet{Garcia_Pichel_2013} are still sharply increasing especially for
sub-crust samples (Figure~\ref*{fig:rarefaction}) suggesting the communities
remain undersampled. Parametric richness estimates of BSC diversity indicate
that the \citet{Steven_2013} and \citet{Garcia_Pichel_2013} sequencing efforts
recovered on average 40.5\% (s.d.  9.99\%) and 45.5\% (s.d. 11.6\%) of
predicted SSU rRNA gene OTUs from crust samples (inset
Figure~\ref*{fig:rarefaction}), respectively. Therefore, it would have been
surprising if all of the $^{15}$N-responsive OTUs had been observed in prior
environmental surveys of BSCs. Nitrogenase \textit{nifH} gene sequences related
to both \textit{Proteobacteria} and \textit{Clostridiaceae} have been
previously observed in BSC samples, though typically at relative abundance that
is much lower than \textit{nifH} gene sequences from heterocystous
cyanobacteria.

% Fakesubsubsection:Notably, the DNA-SIP analysis did not
We propose three mechanisms that could bias
\textit{nifH} clone libraries against heterotrophic diazotrophs. First, extreme
polyploidy in cyanobacteria (up to 58x ploidy in stationary
phase,\citep{Griese_2011}) can be expected to inflate the representation of
cyanobacteria \textit{nifH} gene sequences in community DNA relative to the
frequency of N$_{2}$-fixing heterocysts. Although, cyanobacteria often
have relatively large cells so ploidy per cell is probably greater than ploidy
per unit volume. Second, heterocysts make up a small fraction of total cells
along a trichome, though all cells in the trichome possess the \textit{nifH}
gene. As a result of polyploidy and heterocyst frequency in a cyanobacterial
filament, the ratio of cyanobacterial \textit{nifH} gene copies to
heterotrophic \textit{nifH} gene copies may be inflated as much as 10$^{3}$
times relative to the corresponding ratio of N$_{2}$-fixing cells (i.e.
the ratio of heterocyst number to the cell number of heterotrophic
diazotrophs). Third, \textit{nifH} PCR primers, which are highly degenerate,
could be biased against heterotrophic diazotrophs. For example, the
\textit{nifH} PCR  primers used in the second round of a widely used nested PCR
protocol \citep{Yeager,Yeager_2004,Yeager_2012} have fairly low coverage for
\textit{Proteobacteria} and \textit{Clostridiales} \citep{Gaby_2012}\,. Primer
``nifH11'' is biased against ``Cluster III'' \textit{nifH} gene sequences which
includes those of the \textit{Clostridiales} (50\% \textit{in silico} coverage
of reference \textit{nifH} sequences). In addition, primer ``nifH22'' has low
coverage of reference sequences from \textit{Proteobacteria}, cyanobacteria and
``Cluster III'' \textit{nifH} gene sequences (16\%, 23\% and 21\% \textit{in
silico} coverage, respectively) \citep{Gaby_2012}\,. Hence, it is reasonable to
assume that heterotrophic diazotrophs may have been underestimated in previous
analyses of early successional BSC communities. Our DNA-SIP results, which do
not require PCR of functional genes, suggest that BSC N-fixation in early
successional BSC may include a large non-cyanobacterial component. This is
consistent with small-scale, spatially resolved functional measurements of
nitrogen fixation in BSCs \citep{15643930} that show a subsurface maximum that
does not coincide spatially with maxima in chlorophyll \textit{a}  (a proxy for
phototrophic biomass) in early-successional crusts, and a surface maximum of
N$_{2}$-fixation in mature crust that coincides with the maximum in chlorophyll
\textit{a}.  

% Fakesubsubsection:We did not observe incorporation of
We did not observe incorporation of $^{15}$N$_{2}$ into the DNA of
heterocystous cyanobacteria in the early successional BSC samples used in this
study. It is possible that $^{15}$N$_{2}$-fixation by heterocystous
cyanobacteria could go undetected in DNA-SIP. One possible explanation for this
result is that the early successional BSC samples used in this study possessed
too few heterocystous cyanobacteria to statistically evaluate their
$^{15}$N-incorporation. Indeed, heterocystous cyanobacteria represented only
0.29\% of sequences from the DNA-SIP data (see results) as opposed to 15\% and
23\% of total sequences in the \citet{Steven_2013} and
\citet{Garcia_Pichel_2013} data, respectively. OTUs that correspond to
heterocystous cyanobacteria (as defined by \citet{Yeager}), all fall below the
sparsity threshold used in our analysis (see methods). Given the sparsity of
heterocystous cyanobacteria sequences in the light crust DNA-SIP data, it is
not possible to conclusively determine whether heterocystous cyanobacteria
incorporated $^{15}$N during the incubation. Our results show that
heterotrophic diazotrophs can contribute to N$_{2}$-fixation in early
successional BSC but they do not exclude the potential for fixation by
heterocystous cyanobacteria. Indeed, heterocystous cyanobacteria if present,
active, and limited for nitrogen would be expected to form heterocysts and fix
N$_{2}$. It is likely that scarcity limits their contribution to
N$_{2}$-fixation in early successional crusts. Heterocystous
cyanobacteria form sessile colonies and they require stabilization of the crust
environment before they can successfully colonize soil; and this stabilization
is performed by other pioneering members of the crust community
\citep{Castenholz_2002}\,.  N$_{2}$-DNA-SIP would also fail to identify
N$_{2}$-fixing bacteria if N$_{2}$-fixation were uncoupled from
DNA replication over the time frame of the experiment (i.e. 4 days), that is
$^{15}$N$_{2}$-DNA-SIP will not detect bacteria that fix N$_{2}$ but do
not incorporate the $^{15}$N-label into DNA. Therefore, the contribution of
heterocystous cyanobacteria (or any other microbe) to N$_{2}$ would be
underestimated if their cell division is uncoupled from N$_{2}$-fixation
at time frames of up to 4 days. We should also note that $^{15}$N can be
incorporated into biomass from trophic interactions although in this case the
$^{15}$N labeling would likely be weaker than that for a N$_{2}$-fixer as
a result of label dilution.
 
% Fakesubsubsection:The OTUs with significant evidence of
The OTUs with significant evidence of $^{15}$N-incorporation during the
incubation were predominantly \textit{Proteobacteria} and \textit{Firmicutes}.
The \textit{Proteobacteria} OTUs with the strongest signal of
$^{15}$N-incorporation all shared high sequence identity ($>$98.5\%) with SSU
rRNA gene sequences from genera known to contain diazotrophs
(Table~\ref*{tab:LTP_blast}). In contrast the \textit{Firmicutes} that
displayed signal for $^{15}$N-incorporation (predominantly
\textit{Clostridiaceae}) were not closely related to any known cultivars
(Table~\ref*{tab:LTP_blast}). Hence, we have little knowledge of the ecology of
these organisms. Assessing the physiological characteristics of these
diazotrophic \textit{Clostridiaceae} may be useful for predicting how
environmental change will affect the development and stability of BSC. Prior
intense cultivation efforts from these crusts in separate studies did not yield
any members of the \textit{Clostridiaceae} \citep{Gundlapally_2006}\,. Although
under sampled in environmental data sets, $^{15}$N-responsive OTUs were indeed
more abundant in sub-crust or in early successional BSC samples relative to
crust surface or mature crust samples (Figure~\ref*{fig:trees} and
Figure~\ref*{fig:rspndr_dist}). While members of \textit{Clostridiaceae} have
been found in low abundance in molecular surveys of BSC, most surveys are
carried out on desicated crust samples, where thick-walled spores would
predominate relative to vegetative cells, thus increasing the likelihood for
their underrepresentation in DNA surveys. It should also be noted that crusts
were incubated in an atmosphere of He and N$_{2}$ rather than O$_{2}$ and
N$_{2}$. While cyanobacteria in the presence of light rapidly produce oxygen
super saturation in BSC relative to air \citep{JPY:JPY774}\,, and whereas
heterotrophic N$_{2}$-fixation by many microorganisms is inhibited in the
presence of atmospheric levels of O$_{2}$, it remains possible that the
conditions present in microcosm are not representative of field conditions and
may have favored N$_{2}$-fixation by crust organisms that are less active in
situ. Further experiments will need to be performed to verify that these
heterotrophic diazotrophs are contributing to the N budgets of early
successional crusts in the field.  

% Fakesubsubsection:Our results generate more refined hypotheses
Our results generate more refined hypotheses pertaining to the contribution of
diazotrophs during the development of BSC communities.  Specifically,
N$_{2}$-fixation in BSC may not be tied solely to the climax of
heterocystous cyanobacteria in mature crusts. Rather, N$_{2}$-fixation
may occur throughout crust development with the transition between early
successional and mature crusts marked by a transition between heterotrophic and
phototrophic N$_{2}$-fixation in the crust community. Therefore,
sub-biocrust soil may contribute significantly to the arid ecosystem N budget
and may be of considerable importance in the early phases of BSC establishment.
We propose that interactions between fast-growing heterotrophic diazotrophs
such as members of the \textit{Clostridiaceae} and filamentous
(non-heterocystous) cyanobacteria are important in the early establishment of
BSC communities.  During progressive desication, cyanobacteria, such as
\textit{M. vaginatus}, accumulate compatible solutes such as trehalose and
sucrose \citep{Rajeev_2013}\,. Upon wetting, microorganisms rapidly excrete
compatible solutes to prevent cell lysis due to osmotic shock
\citep{Poolman_1998}\,. Among them are dihexoses (such as sucrose and
trehalose), which are observed in natural crusts upon wetting and then are
rapidly depleted in the soil solution \citep{baran}\,.  Many
\textit{Clostridiaceae} have a saccharolytic metabolism with the potential for
rapid growth rates on substrates such as trehalose and/or sucrose
\citep{Wiegel_2006}\,. Wetting of crust may allow for rapid germination and
growth of these organisms as the time required for germination of clostridial
spores can be less than 30 minutes \citep{Stringer_2005}\,. Indeed, intense
blooms of clostridia have been detected in crusts within tens of hours of
wetting \citep{karaoz}\,. N$_{2}$-fixing clostridia are common in soils
\citep{Wiegel_2006} and it is notable that \textit{C. pasteurianum}, isolated
from soil, was the first N$_{2}$ fixing bacterium ever described
\citep{wino_1985}\,. \textit{C. pasteurianum}, though an anaerobe, grows
readily in the presence of oxygen when co-cultured with aerobic organisms that
reduce oxygen tension \citep{chester}\,. We propose that during a typical
precipitation event, water saturation and heterotrophic activity rapidly render
the interior of the crusts anoxic \citep{JPY:JPY774}, presenting optimal
conditions for growth of anaerobic, dihexose-fermenting, N$_{2}$ fixing
clostridia. Clostridial organic nitrogen would then become available to other
members of the community, including the primary producers, when carbon
limitation induces sporulation and mother cell lysis. Mother cell lysis, the
last step in sporulation, releases rich sources of P and N into the environment
in the form of nucleotides and peptides \citep{Hoch_2002}\,. 

\subsection{Conclusion} The abundance of $^{15}$N-responsive OTUs from
\textit{Clostrideaceae} and \textit{Proteobacteria} found in this study, the
\textit{nifH} gene sequences of \textit{Clostrideaceae} and
\textit{Proteobacteria} observed previously in BSC \citep{Steppe_1996}\,, and
the evidence for subsurface N$_{2}$-fixation in early successional BSC
\citep{15643930}\,, taken together, suggest that heterotrophic diazotrophs may
be important contributors to N$_{2}$-fixation in the subsurface of early
successional BSC. Heterocystous cyanobacteria are also key contributors to the
BSC N-budget, however and it is clear that heterocystous cyanobacteria increase
in abundance with BSC age \citep{Yeager_2004}\,. It is less clear if the
transition to mature crust is marked mainly by a change in the abundance and
activity of heterocystous cyanobacteria, or rather represents a succession
within the diazotroph community from early crusts where N$_{2}$-fixation
is dominated by \textit{Clostridiaceae} and \textit{Proteobacteria} to mature
crusts where it is dominated by heterocystous cyanobacteria. Predicting the
ecological response of BSC to climate change, altered precipitation regimes,
and physical disturbance requires an understanding of crust establishment,
stability, and succession. Diazotrophs are critical contributors to all of
these phenomena and their activities make critical contributions to the
N-budget of arid ecosystems worldwide.


\section{Materials and Methods}
\subsection{BSC sampling and incubation conditions}
BSC samples were taken from the Green Butte site near Moab, Utah as previously
described (site "CP3"; latitude N 38$^{\circ}$44'55.1", longitude
W 109$^{\circ}$44'37.1"; \citealt{BERALDI_CAMPESI_2009}). All samples were from
early successional 'light' crusts as described by \citep{15643930}\,. Early
successional BSC samples (37.5 cm$^{2}$, average mass 35 g) were incubated in
sealed chambers under controlled atmosphere and in 16 h light / 8 h dark for
4 days. Crusts were sampled and transported while dry and wetted at initiation
of the experiment. Water was added to each sample to fully saturate the soil,
but avoid visible ponding. The samples were then placed in air-tight sealed
incubation containers for the rest of the experiment, so that soil and
atmosphere remained saturated through the incubation period. The water was
amended with calcium bicarbonate to yield a final concentration of 3 mM, so
that autotrophy could proceed unimpeded. The control treatment received
a headspace of air and the experimental treatment received a headspace
containing $^{15}$N$_{2}$ ($>$98\% atom $^{15}$N$_{2}$). $^{15}$N$_{2}$
(100\%) gas was purchased from Sigma-Aldrich (St. Louis, MO). We used
a composition of 75\% $^{15}$N$_{2}$ in helium for the initial incubation
headspace. Four crust samples were treated and incubated (two control and two
experimental). One control/experimental crust pair was collected at day 2 and
the other at day 4. Acetylene reduction rates were measured daily. Acetylene
reduction rates increased over the course of the experiment (0.8, 4.8, 8.8,
and 14.5 $\mu$moles m$^{-2}$ hr$^{-1}$ ethylene for days 1 through 4,
respectively). 

\subsection{DNA extraction}
DNA was extracted for DNA-SIP at 2 and 4 days. DNA was extracted from 1 g of
BSC. DNA from each sample was extracted using a MoBio (Carlsbad, CA) UltraClean
Mega Soil DNA Isolation Kit (following manufacturer's protocol, but lysis was
done as previously described \citep{Strauss2011}), and then gel purified to
select high molecular weight DNA ($>$4 kb) using a 1\% low melt agarose gel and
$\beta$-agarase I for digestion (manufacturer's protocol, New England Biolabs,
M0392S). Extracts were quantified using PicoGreen nucleic acid quantification
dyes (Molecular Probes).

\subsection{Formation of CsCl equilibrium density gradients} 
CsCl gradient fractionation was used to separate the DNA into 36 gradient
fractions on the basis of buoyant density. CsCl density gradients were formed
in 4.7 mL polyallomer centrifuge tubes filled with gradient buffer (15 mM
Tris-HCl, pH 8; 15 mM EDTA; 15 mM KCl) which contained 1.725 g mL$^{-1}$
CsCl. CsCl density was checked with a digital refractometer as described below.
A total of 2.5-5.0 $\mu$g of DNA was added to each tube, and the tubes mixed,
prior to centrifugation. Centrifugation was performed in a TLA-110 fixed angle
rotor (Beckman Coulter) at 20$^{\circ}$C for 67 hours at 55,000~rpm.
\citep{17369332}\,. Centrifuged gradients were fractionated from bottom to top
in 36 equal fractions of 100 $\mu$L, using a syringe pump as described
previously \citep{17369332}\,. The density of each fraction was determined
using an AR200 refractometer modified to accommodate 5 $\mu$L samples as
described previously \citep{17369332}\,. DNA in each fraction was desalted on
a filter plate (PALL, AcroPrep Advance 96 Filter Plate, Product Number 8035),
using four washes with 300 $\mu$L TE per fraction. After each wash, the filter
plate was centrifuged at 500~x~g for 10 minutes, with a final spin of 20
minutes. Purified DNA from each fraction was resuspended in 50 $\mu$L of TE
buffer.

\subsection{PCR, library normalization and DNA sequencing} 
To characterize the distribution of SSU rRNA genes across density gradients,
SSU rRNA gene amplicons were generated from 20 gradient fractions per gradient
for both unlabeled controls and $^{15}$N$_{2}$ labeled samples. The 20
fractions analyzed are those expected to contain DNA (both labeled and
unlabeled) having buoyant density in the range of 1.66 g mL$^{-1}$ to 1.77
g mL$^{-1}$. Barcoded PCR of bacterial and archaeal SSU rRNA genes was carried
out using primer set 515F/806R \citep{21349862} (primers purchased from
Integrated DNA Technologies).  The primer 806R contained an 8 bp barcode
sequence, a "TC" linker, and a Roche 454 B sequencing adapter, while the primer
515F contained the Roche 454 A sequencing adapter. Each 25 $\mu$L reaction
contained 1x PCR Gold Buffer (Roche), 2.5 mM MgCl$_{2}$, 200 $\mu$M of each of
the four dNTPs (Promega), 0.5 mg mL$^{-1}$ BSA (New England Biolabs), 0.3
$\mu$M of each primers,
1.25~U of Amplitaq Gold (Roche), and 8 $\mu$L of template. Each sample was
amplified in triplicate. Thermal cycling occurred with an initial denaturation
step of 5 minutes at 95$^{\circ}$C, followed by 40 cycles of amplification (20 s at 95$^{\circ}$,
20 s at 53$^{\circ}$, 30 s at 72$^{\circ}$), and a final extension step of 5 min at 72\selectlanguage{ngerman}°C.
Triplicate amplicons were pooled and purified using Agencourt AMPure PCR
purification beads, following manufacturer's protocol. Once purified,
amplicons were quantified using PicoGreen nucleic acid quantification dyes
(Molecular Probes) and pooled together in equimolar amounts. Samples were sent
to the Environmental Genomics Core Facility at the University of South
Carolina (now Selah Genomics) where they were run on a Roche FLX 454
pyrosequencing machine (FLX-Titanium platform). 

\subsection{Data analysis} 
\subsubsection{Sequence quality control}
Sequences were initially screened by maximum expected errors at a specific read
length threshold \citep{23955772} and this has been shown to be as
effective as denoising with respect to removing pyrosequencing errors.
Specifically, reads were first truncated to 230 nucleotides (nt) (all reads
shorter than 230 nt were discarded) and any read that exceeded a maximum
expected error threshold of 1.0 was removed. After truncation and max expected
error trimming, 91\% of original reads remained. Forward primer and barcode
were then removed from the high quality, truncated reads.  Remaining reads were
taxonomically annotated using the "UClust" taxonomic annotation framework in
the QIIME software package \citep{20383131, 20709691} with cluster seeds from
Silva SSU rRNA database \citep{17947321} 97\% sequence identity OTUs as
reference (release SSU Ref
111). Reads annotated as "Chloroplast", "Eukaryota", "Archaea", "Unassigned" or
"mitochondria" were removed from the dataset. Finally, reads were aligned
to the Silva reference alignment provided by the Mothur software package
\citep{19801464} using the Mothur NAST aligner \citep{16845035}\,. All
reads that did not align to the expected amplicon region of the SSU rRNA
gene were discarded. Quality control parameters removed 34,716 of 258,763 raw
reads. Raw sequences have been uploaded to MG-RAST (MG-RAST ID 4603397.3).

\subsubsection{Sequence clustering}
Sequences were distributed into OTUs using the UPARSE methodology
\citep{23955772}\,. Specifically, OTU centroids (i.e. seeds) were identified
using USEARCH on non-redundant reads sorted by count. The sequence
identity threshold for establishing a new OTU centroid was 97\%. After initial
OTU centroid selection, select SSU rRNA gene sequences from \citet{Yeager}
were added to the centroid collection. Specifically, \citet{Yeager} Colorado
Plateau or Moab, Utah sequences were added which included the SSU rRNA gene
sequences for \textit{Calothrix} MCC-3A (accession DQ531700.1), \textit{Nostoc
commune} MCT-1 (accession DQ531903), \textit{Nostoc commune} MFG-1 (accession
DQ531699.1), \textit{Scytonema hyalinum} DC-A (accession DQ531701.1),
\textit{Scytonema hyalinum} FGP-7A (accession DQ531697.1), \textit{Spirirestis
rafaelensis} LQ-10 (accession DQ531696.1). Original centroid sequences that
matched selected \citet{Yeager} (above) sequences with greater than to 97\%
sequence identity were subsequently removed from the centroid collection. With
USEARCH/UPARSE, potential chimeras are identified during OTU centroid selection
and are not allowed to become cluster centroids effectively removing chimeras
from the read pool. All quality controlled reads were then mapped to cluster
centroids at an identity threshold of 97\% again using USEARCH. A total of
95.6\% of quality controlled reads could be mapped to centroids. Unmapped reads
do not count towards sample counts and were removed from downstream analyses.
The USEARCH software version for cluster generation was 7.0.1090.
\citet{Garcia_Pichel_2013} and \citet{Steven_2013} sequences were quality
screened by alignment coordinates (described above) and included as input to
USEARCH for OTU centroid selection and subsequent mapping to OTU centroids. 

\subsubsection{Phylogenetic analysis}
Alignment of SSU rRNA genes was done with SSU-Align which is based on Infernal
\citep{Nawrocki_2013, 19307242}\,. Columns in the alignment that were not included in
the SSU-Align covariance models or were aligned with poor confidence (less than
95\% of characters in a position had posterior probability alignment scores of
at least 95\%) were masked for phylogenetic reconstruction. Additionally, the
alignment was trimmed to coordinates such that all sequences in the alignment
began and ended at the same positions. FastTree \citep{20224823} was used to
build the tree.

\subsubsection{Identifying OTUs that incorporated $^{15}$N into their DNA}
DNA-SIP is a culture-independent approach towards defining identity-function
connections in microbial communities \citep{Buckley_2011, 17446886,
Radajewski_2001}\,. Microbes are identified on the basis of
isotope assimilation into DNA. As the buoyant density of a macromolecule is
dependent on many factors in addition to stable isotope incorporation (e.g.
G+C-content in nucleic acids \citep{25139123})\,, labeled nucleic acids from
one microbial population may have the same buoyant density as unlabeled nucleic
acids from another. Therefore, it is imperative to compare results of isotopic
labelling to results obtained with unlabeled controls where everything mimics
the experimental conditions except that unlabeled substrates are used. By
contrasting heavy gradient fractions from isotopically labeled samples relative
to corresponding fractions from controls, the identities of microbes with
labeled nucleic acids can be
determined 

We used an RNA-Seq differential expression statistical framework
\citep{Love_2014} to find OTUs enriched in heavy fractions of labeled
gradients relative to corresponding density fractions in control gradients
(for review of RNA-Seq differential expression statistics applied to
microbiome OTU count data see \citealt{24699258}). We use the term
\selectlanguage{english}"differential abundance\selectlanguage{english}" (coined by \citealt{24699258}) to denote OTUs that
have different proportion means across sample classes (in this case the only
sample class is labeled:control).  CsCl gradient fractions were categorized
as "heavy" or "light". The heavy category denotes fractions with density
values above 1.725 g mL$^{-1}$. Since we are only interested in enriched OTUs
(labeled versus control), we used a one-sided Wald-test to test the statistical
significance of regression coefficients (the null hypothesis is that the
labeled:control fold enrichment for an OTU is less than a selected threshold).
We independently filtered out sparse OTUs prior to P-value correction for
multiple comparisons. The sparsity threshold was set to the value which
maximized the number of p-values under a false discovery rate (FDR) the
specific sparsity threshold was 0.3 meaning that an OTU not found in at least
30\% of heavy fractions (control and labeled gradients) in a given day were not
considered further and not included in P-value adjustment for multiple
comparisons. P-values were corrected with the Benjamini-Hochberg method
\citep{citeulike:1042553}\, and a FDR of 0.10 was applied (this rate is the
typical FDR threshold adopted during RNASeq analysis). We selected a log$_{2}$
fold change null threshold of 0.25 (or a labeled:control fold enrichment of
1.19). DESeq2 was used to calculate the moderated log$_{2}$ fold change of
labeled:control proportion means and corresponding standard errors for the
Wald-test (above). Fold change moderation allows for reliable ranking such that
high variance and likely statistically insignificant fold changes are
appropriately shrunk and subsequently ranked lower than they would be as
unmoderated values. Those OTUs that exhibit a statistically significant
increase in proportion in heavy fractions from $^{15}$N$_{2}$-labeled samples
relative to corresponding controls have increased significantly in buoyant
density in response to $^{15}$N$_{2}$ treatment; a response that is expected
for N$_{2}$-fixing organisms.

We also assessed the consistency of enrichment between time points by including
the interaction of day and label:control in a DESeq2 generalized linear model.
The interpretation of the interaction coefficient is the change in OTU
enrichment per unit time. P-values for the interaction coefficient were
adjusted for all OTUs that passed the sparsity threshold in the label versus
control comparison (above) and we used the default null model that the
coefficient equaled zero. Additionally, we assessed fold change between labeled
and control gradient heavy fractions after pooling day 2 and day 4 data when
treating the different time points as replicates. The same null model as the
label versus control comparison (above) was used in this replicate analysis
(log$_{2}$ fold change in abundance between label and control is less than or
equal to
0.25). We included all OTUs that passed sparsity based independent filtering at
either day (above) for p-value adjustment in the replicate analysis.

\subsubsection{Community and Sequence Analysis}
BLAST searches were done with the "blastn" program from BLAST+ toolkit
\citep{20003500} version 2.2.29+. Default parameters were always employed and
the BioPython \citep{19304878} BLAST+ wrapper was used to invoke the blastn
program. Pandas \citep{citeulike:11241428} and dplyr \citep{dplyr} were used to
parse and manipulate BLAST output tables.  

Principal coordinate ordinations depict the relationship between
samples at each time point (day 2 and 4). Bray-Curtis distances were
used as the sample distance metric for ordination. The Phyloseq
\citep{24699258} wrapper for Vegan \citep{vegan} (both R packages) was
used to compute sample values along principal coordinate axes. GGplot2
\citep{ggplot2} was used to display sample points along the first and
second principal axes. Adonis tests \citep{Anderson_2001} were done with default
number of permutations (1000).

Rarefaction curves were created using
bioinformatics modules in the PyCogent Python package \citep{Knight_2007}\,.
Parametric richness estimates were made with CatchAll using only the best model
for total OTU estimates \citep{BUNGE_2010}\,.

All code to take raw sequencing data through the presented figures (including
download and processing of literature environmental datasets) can be
found at:

\url{https://github.com/chuckpr/NSIP\_data\_analysis}



\section{Acknowledgements} We would like to thank T Whitman, CHD Williamson, AN
Campbell and EK Hall for helpful comments in the preparation of this
manuscript. This material is based upon work supported by the Department of
Energy Office of Science, Office of Biological \& Environmental Research Genomic
Science Program under Award Numbers DE-SC0004486 and DE-SC0010558. This project
was also supported by Agriculture and Food Research Initiative Competitive
Grant no. 2007-35107-18299 from the USDA National Institute of Food and
Agriculture. This report was prepared as an account of work sponsored by an
agency of the United States Government. Neither the United States Government
nor any agency thereof, nor any of their employees, makes any warranty, express
or implied, or assumes any legal liability or responsibility for the accuracy,
completeness, or usefulness of any information, apparatus, product, or process
disclosed, or represents that its use would not infringe privately owned
rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise does not
necessarily constitute or imply its endorsement, recommendation, or favoring by
the United States Government or any agency thereof. The views and opinions of
authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.91\columnwidth]{figures/ordination-heavy-day-facet/ordination-heavy-day-facet}
\caption{{Ordination of heavy gradient fractions by Bray-Curtis distances on the basis ofOTU content. Each point represents a gradient fraction OTU profile. Pointscloser together have more similar OTU content than those further apart.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.42\columnwidth]{figures/study-phylum-dist/study-phylum-dist}
\caption{{Distribution of sequences into top 9 phyla (phyla ranked by sum of all sequence annotations).%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/ordination-all-day-facet/ordination-all-day-facet}
\caption{{Ordination of Bray-Curtis sample pairwise distances for eachincubation time. Point area is proportional to the density of the CsClgradient fraction for each sequence library, and color/shape reflects control(red triangles) or labeled (blue circles) treatment. Each point represents theOTU profile for a single gradient fraction. Points closer together are moresimilar in OTU content than those further apart.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.56\columnwidth]{figures/l2fc1/l2fc}
\caption{{Moderated log$_{2}$ fold change of OTUs proportions for labeled versus
control gradients (heavy fractions only, densities $>$1.725 g/mL). All OTUs passing the 
sparsity treshold (see methods) at a specific incubation day are
shown. Red color denotes a proportion fold change that has a corresponding
adjusted p-value below a false discovery rate of 10\% (ratio is significantly
greater than 0.25, black line.)%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/scatter-heavy-topN2/scatter-heavy-topN}
\caption{{Relative abundance values in heavy fractions (density greater orequal to 1.725 g/mL) for the top 10 $^{15}$N "responders" (putativediazotrophs, see results for selection criteria of top 10) at eachincubation day. Each point is a relative abundance value for the indicatedOTU in a CsCl gradient fraction SSU rRNA gene collection. SeeTable~\ref*{tab:LTP_blast} for BLAST results against the LTP database (release115). Point area is proportional to CsCl gradient fraction density, and colorsignifies control (red) or labeled (blue) treatment.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.84\columnwidth]{figures/rarefacation-curves1/raref-and-boxplot}
\caption{{Rarefaction curves for all samples presented by \protect\citet{Garcia_Pichel_2013} and \protect\citet{Steven_2013}. Inset is boxplot of estimated sampling effort for all samples in \protect\citet{Garcia_Pichel_2013} and \protect\citet{Steven_2013} (number of observed OTUs divided by number of CatchAll \protect\cite{BUNGE_2010} estimated total OTUs).%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/rspndr-dist/rspndr-dist}
\caption{{Counts of "responder" OTU occurrences in samples from \protect\citet{Steven_2013} and
 \protect\citet{Garcia_Pichel_2013}. \protect\citet{Steven_2013} collected BSC samples (25
 samples total) and samples from soil beneath BSC (17 samples total, "below"
 column in figure). \protect\citet{Garcia_Pichel_2013} collected samples from "dark" (9
 samples total) and "light" (12 samples total) crusts in addition to "Lichen"
 (2 samples total) dominated crusts.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.56\columnwidth]{figures/trees/trees}
\caption{{Phylogenetic trees of OTUs passing sparsity threshold for
\textit{Proteobacteria} \textbf{A} and 
\textit{Firmicutes} \textbf{B}. $^{15}$N-responders are identified by dots present in
column \textbf{i}. Log$_{2}$ of OTU proportion fold change (labeled:control
samples) for each OTU are presented as a heatmap in column \textbf{ii} with
results from days 2 and 4 on the left and right sides of the column
respectively. High fold change values indicate $^{15}$N incorporation.
Presence/absence of OTUs (black indicates presence) in lichen, light, or dark
environmental samples \protect\citep{Garcia_Pichel_2013} is shown in column
\textbf{iii}.  Presence/absence of OTUs (black indicates presence) in crust and
below crust samples \protect\citep{Steven_2013} is shown in column \textbf{iv}.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/2v4l2fc/2v4l2fc}
\caption{{Scatter plot of fold enrichment values (label versus control heavy fractions) for OTUs passing sparsity criteria in day 2 and 4. P-value is the minimum value from day 2 or 4. Blue line has slope of one.%
}}
\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\bibliographystyle{plain}
\bibliography{bibliography/converted_to_latex.bib%
}

\end{CJK}\end{document}

