Geographic distribution of nematodes in the Atacama is associated with elevation, climate gradients and parthenogenesis

Domain: biology
Published: 2026-01-09
Editorial outcome: minor_revision

Abstract
Soilecosystemsarecrucialforsupportinglife,yetlittleisknownabouttheir
biodiversityanditsdistributioninextremeenvironments.TheAtacama
Desert,thedriestnon-polardesertonEarth,hasscarcewater,highsalinity,and
metal-richwaterbodies,creatingchallengingconditionsformostorganisms.
Above-groundlifehasbeenpartiallydocumented,butitssoilsremainpoorly
studied.Hereweshowthatsoilnematodes,anabundantanddiversegroupof
invertebrates,displaydistinctbiodiversitypatternsacrosstheAtacamaat
multiplebiologicalscales,includinggenetic,taxonomic,community,andlife-
cyclelevels.Surveysacrossdunesystems,high-altitudemountains,saline
lakes,rivervalleys,andfogoasesrevealuniqueassemblagesineachhabitat.
Wefindthatasexualtaxaaremorecommonathigheraltitudes,consistentwith
patternsofgeographicalparthenogenesis.Genusrichnessfollowsalatitudinal
gradientandincreaseswithprecipitation.Theseresultsdemonstratethateven
inoneofthemostextremeterrestrialenvironments,stablesoilcommunities
canpersist.However,evidenceofsimplifiedsoilfoodwebssuggestsvulner-
abilitytofurtherenvironmentalchange.Ourfindingsprovidenewinsightsinto
themechanismsshapingbiodiversityinaridecosystemsandcaninformpre-
dictionsaboutsoilresilienceunderglobalclimate-drivenaridification.

Introduction
Understandingwhichfactorsshapethediversityanddistributionof
speciesisamajorfocusofecologicalresearch.Studiesexaminingthe
interactionoforganismswiththeirimmediateenvironmentfocuson
theinfluenceofecological,geological,andgeographicalparametersin
diverse systems such as aquatic and soil habitats1–5. From nutrient
cycling,mineralization,andpestcontrol6tocarbonsequestrationand
storage,soilsandtheirprocessesplayakeyroleinfoodsecurity,water
supplyandquality,andinshapingbiodiversity7.Globally,atleast58%ofspeciesinhabitsoils8,includingmicroorganismsanddifferentsize
classesoffauna9.Accordingly,soilsreceiveextensiveresearchatten-
tion, particularlyin arableareas.However, geographicalpatterns of
soil biota in relation to environmental parameters remain poorly
understood,largelyduetolimiteddataonthetruediversityanddis-
tributionofmanysoilorganisms8,10.
Evenlessinformationexistsforsoilsystemsinextremeenviron-
ments such as deserts. Abiotic factors such as strong temperature
fluctuations,minimalannualprecipitation,andhighevaporationrates
drivedesertaridityandposesignificantchallengesforlife11.Despite
theseconstraints,desertsandtheirsoilscontainspatiallydistributed
microhabitats and niches, including limestone areas, dune systems,
andsalinepatches,thatsupportdiverselifeforms11.Withtheglobal
trend of aridification12–14, studying biodiversity in desert soils is
essentialforpredictingorganismaldistributionsunderglobalchange.
TheAtacamaDesert,thedriestnon-polardesertonEarth,receives
lessthan200mmofrainfallannually.Despitethelong-termstabilityof
its geological system, regional environmental patterns are
apparent15,16. Long-term aridity shapes the distribution of sulfates,
chlorides,andnitratesinAtacamasoils,whichinturnaffectcommu-
nitycompositionacrosstaxa17–20.DiversemicrohabitatsintheAtacama
include alkaline soils, dune systems, fog oases, and soils with high
salinityandarseniccontent.Thesesupportspecializedtaxasuchas
microbiallifeforms,plants,andalimitednumberofvertebrates21,21–31.
InvertebratebiodiversityremainsunderstudiedintheAtacama,with
research focusing mainly on above-ground fauna32–34, leaving soil
invertebrates poorly sampled and their distribution patterns
unclear8,29,35–37. Nematodes, one of the most diverse and abundant
metazoan groups, contribute to soil nutrient turnover, carbon
sequestration,andregulationofbacterialpopulations,makingthem
important indicators of soil condition36,38–40. Their use as functional
indicatorsisbasedonlifecyclecharacteristicssuchasfeedingtype
and reproductive strategy38,41. Members of Nematoda adapt to
extremeconditions,includinglowoxygenindeep-seaenvironments42,
daily freezing in Antarctica43,44, and kilometer-deep mines45. Certain
families are prevalent and diverse in deserts, surviving low water
availability in systems suchas the Mojave46, Chihuahuan47, Monte48,
and Namib deserts49. Many nematode genera also evolve partheno-
geneticlineages,hypothesizedtohaveafounderadvantageinextreme
environments (geographical parthenogenesis)50,51. In the Atacama,
parthenogenesis may therefore be more frequent than sexual
reproduction.
Published records of nematodes from the Atacama are usually
based on very few or single individuals35. Consequently, it remains
unclear whether soil invertebrate biodiversity patterns align with
regionalandlocalecology,adaptivestrategiessuchasfeedingtypeor
reproductivemode,orbroaderglobaltrendssuchasenvironmental
gradients in species richness3. In this study, we analyze the factors
shapinginvertebratebiodiversityinAtacamasoilsusingnematodesas
amodelsystem.Wesamplesixlocalitiesacrossthedesert,including
sanddunes,salinelakeshores,andhigh-altitudemountainregions,to
characterizebiogeographicpatternsingeneticdiversity,lifehistory,
andcommunitycomposition.Despiteextremeandpersistentaridity,
wefindstronglocalpatterningacrossnichessuchassanddunes,saline
lakes, fog oases, and river valleys. Using genus richness and repro-
ductivemodemodelingapproaches,weidentifyannualprecipitation
and temperature heterogeneity as the main factors associated to
genusrichness,andfoundelevationtobethestrongestpredictorof
reproductivemode.Weshowthatevenunderextremeenvironmental
conditions, soil niches can sustain stable and diverse invertebrate
communities,althoughsomeareasshowsignsofsimplifiedfoodwebs.

Results
Inthisstudy,weaimedtoanalyzedifferentregionsoftheAtacamaand
thus established localities and sampling sites in major areas of the
desert: (1) the Altiplano area (Tarapacá region) is more humid and
knownforsupportingmorevegetationthanthehyperaridcoreofthe
Atacama52; (2) the area around the Quebrada de Aroma (Tarapacá
region)hasmoderatelyalkalinesoilswithincreasinggypsumasalti-
tudedecreasesandharborsseveralmicroorganismtaxasuchasPro-
teobacteriaandActinobacteria21,(3)theSanFranciscoarea-furtheron
referred to as Eagle Point (Antofagasta region)—is home to several
bacterialspecieswithpredominantFirmicutestaxa24,(4)theSalarde
Huasco(Tarapacáregion)experienceshighUVradiation,salinityand
arsenicconcentrationssupportingseveralmicrobialtaxa25,53aswellas
threeflamingospecies,anearlythreatenedpumaspeciesandseveral
othervertebrates23,(5)thePaposoarea(Antofagastaregion)receives
waterintheformoffog(calledcamanchaca),itharborsplantspecies-
richfogoases(knownaslomaformations)54,55andisanareaofhigh
endemism25,56,57andthe(6)Totoralarea(Atacamaregion)iscomposed
ofsandduneswhereendemicplantsoftheApocynaceaefamily22,as
wellasreptilesandbirdscanbefound.Inordertoanalyzepatternsof
biodiversityinthehyper-aridAtacamaDesertusingnematodesasa
model for soil organisms, we sampled, identified, sequenced, and
studieddiversityatdifferentlevels.Weanalyzeddiversityatthehap-
lotype, genus, family, and community level, including life-history
characteristics such as feeding types and reproductive modes. We
usedlinearmodelsaswellasarandomforestapproachtoassessthe
factorsthatcontributethemosttothegenusrichnessandreproduc-
tivemodepatternsacrossthedesert.

ArichdiversityofroundwormspersistsintheAtacama
Wecollected393nematodemorphotypesfrom112samplescollected
along sixlocalities throughout the Atacama Desert(Fig. 1B, supple-
mentaryFig.S1).Collectednematodeswerethensequenced(”Meth-
ods”)withapercentageofhigh qualitybasepairsaboveatleast20(HQ
20%),386sequencesfromdifferentmorphotypeswerekeptforfur-
theranalysis.Morespecifically,weanalyzedatotalof39nematodes
fromtheAltiplanolocality,37fromtheAromalocality,163fromthe
EaglePointlocality,44fromthePaposolocality,26fromtheSalars
locality(Salars),and77fromTotoralduneslocality.Oursamplingalso
included sites whereno nematodes couldbe found (Fig. 1B). These
sitesareusedasareferenceforenvironmentalconditionsunderwhich
nematodesdonotseemtopersist.Weclassifiedthenematodesinto21
familiesand36genera(listoflineagespersamplecanbefoundon
10.5281/zenodo.15517342—FamilyGeneraperSample.csv), illustrating a
rich diversity of roundworms flourishing in the Atacama. Seven
familieswerefoundintheAltiplanoaswellastheAromalocality,nine
familieswerefoundinEaglePoint,12intheSalars,eightinTotoral
Dunes,andteninPaposo(Fig.2B).Weobservedthatsomegeneralist
familieswereubiquitouslyfound,whereasothersweregeographically
restricted.Cephalobidaewerepresentinalllocalities,whereasrepre-
sentatives of families Panagrolaimidae and Qudsianematidae were
found in five of the localities (all except for the Salars). The family
AnguinidaewasuniquelyfoundintheEaglePointlocality.Thefamilies
TrichistomatidaeandMermithidaewereuniquetotheSalarslocality.
Tylencholaimidae was uniquely found in the Paposo locality, and
AcuariidaeintheAltiplano.

Geneticdiversityandhaplotypeuniquenessdifferbetween
nematodetaxa
Given the diversity of retrieved nematode genera and the disparity
betweenubiquitousandregionallyrestrictedtaxa,wewantedtoana-
lyzepatternsofbiogeographyonthegeneticlevel.Weassessedhap-
lotype clustering and genetic diversity (π, θ) for seven taxonomic
groupsusing18SrRNA(Table1).Westudiedgeneticdiversityintwo
ways: per taxonomic group, including all habitats, and per habitat,
analyzingonetaxonomicgroupatatime.
Therewerebetween17to166segregatingsitesforthedifferent
soil nematodes in our networks (supplementary Table S8), and we
identified differences in haplotype structure and clustering at the
generaandfamilylevels(supplementaryFigs.S3–S7).Forthegenera
Acrobeloides,Acrobeles,andPanagrolaimuswefound24,22,and20
differenthaplotypes,respectively(Table1).Forthefamiliesincludedin
theorderDorylamidaandthesuperfamilyAphelenchoidea,16haplo-
typeswereidentified,respectively.ThefamilyPlectidaewascomposed
often different haplotypes. For Alaimidae, only one haplotype was
found(twoincludingthereferencesequence)andissharedinthree
habitats(Table1).Overall,wewerenotabletofindanydistinctgeo-
graphicalclusteringofhaplotypesbasedonthenetworks,i.e.,indivi-
dual haplotypes were present in geographically distant localities (a
complete list of all identified haplotypes can be found on 10.5281/
zenodo.15517342).Wefoundgeneticdiversityatthegenusleveltobe
thehighestforPanagrolaimus(π=0.05788,θ=0.04221,Fig.3A)and
thelowestforAcrobeloides(π=0.01278,θ=0.02192).Atthefamily
level,thehighestgeneticdiversitywasfoundforAphelenchoidea(π=
0.11199,θ=0.08095,Fig.3B),andthelowestwasfoundforPlectidae(π
=0.00633,θ=0.00993).
Tobetterunderstandthepossibleeffectofgeographyongenetic
distance,wecalculatedthegeneticdistanceasEuclideandistanceand
comparedittothegeographicdistanceofthedifferentlocalitiesusing
MantelstatisticsbasedonSpearman’srankcorrelation(ρ).ForPana-
grolaimus,asignificant(p<0.05)isolation-by-distancepatterncould
be found (r= 0.6,p = 0.041667, 119 permutations); no pattern was
found for other taxa analysed. With this approach, it appears that
geographically distanthaplotypesofPanagrolaimusare alsogeneti-
callymoredistantintheir18SrRNAregion.Specificallyforthisgenus,
weseethehighestdiversityintheEaglePoint(π=0.06548)andthe
lowest in Totoral Dunes, where only one haplotype was found for
severalindividuals(supplementaryTableS9).Theseresultsshowcase
thehighvariabilityofgeneticdiversitywithinsoilnematodetaxainthe
AtacamaDesertandrelationtotheirhabitat.

Thenematodecommunitycompositionvariesbetweenhabitats
Weassessedsoilnematodecommunitycompositionvariationinthe
sixlocalitiesemployingtheJaccardsimilarityindexatthefamilylevel.
WefoundtheAromaandEaglePointlocalitiestosharethehighest
proportion of families (J=0.78), followed by Altiplano and Aroma
(J=0.75).Conversely,weobservedtheleastamountofsharedfamilies
between Eagle Point and the Totoral Dunes (J=0.21) and the Salars
locality when compared to Totoral Dunes (J=0.25)(Fig. 2A). As
expectedundertheisolation-by-distancemodel,communitycompo-
sitionappearstobecomemoredissimilarwithincreasedgeographical
distance between the different localities (Mantel statistic based on
Spearman’srankcorrelationρ,r0.5153,p=0.048611,719permuta-
tions). We were interested to know if a co-occurrence analysis of
familieswouldidentifyaspecificgeographicalpattern,butusing82
familypairs(108pairswereremovedduetouniquepresenceinspe-
cificlocalities)couldnotretrieveasignal.

SoilsintheAtacamaaredominatedbydifferentcategoriesof
nematodes
We wanted to use nematodes as indicators of soil condition
throughout different habitats in the desert by assessing their
differentfunctionalgroups.Wecategorizednematodesintodiffer-
entfeedingtypes,aswellasintothecolonizer-persisterseriesbased
on common life-history characteristics reported in literature and
assigned using the NINJA application41. Nematodes that are
r-strategists (c-p1 colonizer) are considered indicators of resource
availability and unstable environments, whereas those that are
k-strategists(c-p5persister)areconsideredindicatorsoffoodweb
complexityandmorestablehabitats58.Wefoundthatthedifferent
habitats show a different structure regarding nematode functional
groups. Despite the presence of nematodes belonging to all
colonizer-persister categories (c-p) across the desert, r-strategists
(categoriesc-p1andc-p2)arepredominantintheAltiplanolocality
(57.2%offamiliesfoundinthelocality)aswellasintheAromaand
TotoralDuneslocalities.Whereask-strategists(categoriesc-p4and
c-p5)werethepredominanttypeoffamiliesinthePaposoandSalars
localities (supplementary Table S11). By assessing the different
feeding types of the nematodes analysed, we found that most
families in the different habitats are bacterivores and omnivores
(41.2% and 20.1% respectively). The remaining nematode families
were categorized as predators, fungivores, and herbivores (18.7%,
14.4%, and 5.6% respectively)(Fig 2B). Families with unicellular
eukaryoticfeedingtype(UEF)werenotfoundinthisstudy(supple-
mentary Table S12). A structuring of the different soil habitats
regardingfeedingtypesisapparent,astheSalarsonlyharborbac-
terivore, predatory and omnivore families whereas other areas as
Totoral Dunes have a more diverse feeding type assemblage with
herbivores,fungivores,bacterivores,predatorsandomnivorespre-
sent(Fig.2B).

ClimateandlatitudedeterminesoilbiodiversityintheAtacama
To understand the driving forces and constraints of nematode
taxonomicbiodiversityacrosstheAtacamaDesertmoreglobally,we
used a modeling approach incorporating environmental data. We
first aimed to predict genus richness and found that the best pre-
dictive model using a probabilistic approach incorporated soil
thicknessandrangeoftemperature.Themodelshowedthehighest
Akaike weight (0.4709556) and a statistically significant slope for
both environmental variables (supplementary Table S16), which
indicatesthatnematodegenerarichnessintheAtacamaDesertcould
berelatedtoreducedsoilthicknessandmoreclimateheterogeneity,
hererepresented asthe rangeoftemperature.Toaccountforvar-
iationacrosslocalitiesinourpredictions,wefurtherassessedgenus
richness using linear mixed-effect models specifying the sampling
locality as a random effect. Our analysis confirmed that sampling
localitieshada significantinfluenceonthemodelingprocess,sug-
gesting that spatial variation plays a crucial role in shaping the
observed patterns, potentially due to underlying environmental
heterogeneitybetweenthedifferentlocations.Incorporatinglocality
effectsintothemodelswasthereforeessentialtocapturetheeco-
logical drivers of genus richness given our dataset. Moreover, we
notedthatprecipitationpatternsarenotcorrelatedtothedistance
betweenlocalities(Supplementarynote1).
Withthisapproach,wethenfoundthatthebestpredictivemodel
incorporated mean annual precipitation and range of temperature,
representingtheseasonalthermalamplitudeanorganismexperiences
withinagivenarea(temperatureheterogeneity).Thepreferredmodel
showedthehighestAkaikeweight(0.5263228)andastatisticallysig-
nificantslopefortheinterceptandbothclimaticvariables(equation1)
(Table2).
Tofurther validatethe results ofthetested linear models, we
implemented a random forest approach. With this approach, we
found that the most important variables for genera richness pre-
diction in our analysis were elevation, mean annual precipitation,
latitude,andrangeoftemperature(supplementaryTableS17).This
corroboratesthatmeanannualprecipitationandrangeoftempera-
ture,incongruencewiththelinearmodels,haveastrong,consistent
influenceongenusrichnessprediction.Withthethreeapproaches,
wefoundthatfewergeneraarepredictedtooccurinthehyper-arid
coreofthedesert(e.gzerofortheregionofPampadelTamarugal),
andmoregeneraarepredictedtooccur,forinstance,intheAltiplano
region (Fig. 4A). Analyzing the raw data relationship between
environmentalvariables and genus richness,wesee thatthelinear
associationsweregenerallyweakornon-significant(supplementary
Fig.S11),therandomforestmodelcapturesmorecomplex,poten-
tially nonlinear effects, as shown in the partial dependence plots
(supplementary Fig. S12). This highlights that genus richness
respondstoenvironmentalgradientsinawaynotfullycapturedby
onlylinearmodels.Duetothedifferingassumptionsunderlyingeach
modelingapproach,asingledefinitivedistributionofgenusrichness
cannot be determined. However, temperature heterogeneity
consistentlyemergedasastrongpredictoracrossallmodels.Addi-
tionally, mean annual precipitation was identified as a significant
factor in both linear mixed-effects and random forest models,
underscoring the important role of climate in shaping genus rich-
ness.

generarichness = −3.723925+0.015028*precipitation (1)
+0.688582*temperaturerange

Reproductivemodeassociationtoenvironmentalvariables
The theory of geographical parthenogenesis predicts an excess of
asexualtaxainmarginalhabitats.Givenourfindingsoftaxonomically
structurednematodecommunitiesindifferentregionsoftheAtacama
Desert,wewantedtotestforthespecificlife-historytraitofparthe-
nogenesis(asexualreproduction)usingourmodelingapproach.We
established single-worm cultures for a representative subset of our
samples. Ten out of the eighteen cultures showed nematodes with
parthenogenetic reproduction, i.e., no presence of males and suc-
cessfulpropagationfromasingleunfertilizedfemale.Theremaining
seven cultures showed sexual reproduction, i.e. presence of males
(supplementary Fig. S2G) and no propagation from single females.
Sexualstrainsshowedapropagationrateof0%fromsingleunferti-
lized female cultures, whereas the asexual strains’ propagation rate
was83.3% (± 10.48%).We initiallytested variables ofelevation and
latitudeaspredictorsofreproductivemodealongwithmeanannual
precipitation and range of temperature, the former two to test for
commonvariablesassociatedwithgeographicalparthenogenesis,the
lattertwogiventheireffectongenerarichnesspredictionfromour
linearmixed-effectsandrandomforestmodels.Wetestedgeneralized
linearmodelsaswellasgeneralizedmixed-effectsmodelstoaccount
fortheeffectofsamplinglocality(supplementaryTablesS20andS21).
Inaddition,wetestedamodelincludingvegetation(supplementary
Fig. S14, supplementary Table S20). The best predictive model for
reproductivemodeincorporatedonlyelevation,showingthehighest
Akaike weight and a statistically significant slope (supplementary
Table S16), while, for example, plant association did not show any
effect on the prediction of reproductive mode. We also tested the
effectofsamplinglocalityonthepredictionofreproductivemodeand
found, in congruence with the previous approach, that the model
incorporatingonlyelevationwaspreferred(Table3).Ourvalidationof
these linear models with a random forest approach supports that
increasingelevationisthemaindriverforasexualreproductivemode
(representedas1)insoilnematodesintheAtacama(supplementary
Table S18). Our models show that the likelihood of nematode taxa
beingasexualisgreaterathigherelevationsthaninlowerones(e.g.
coastalarea)(Fig.4).Toaccountforvariationacrosslocalitiesinour
predictions, we further assessed genus richness using generalized
linearmixedeffectmodelsspecifyingthesamplinglocalityasaran-
domeffect.Inthiscase,wefoundthatthesamplinglocalitieshaveno
significanteffectonthelikelihoodofareproductivemode.

Methods
Sampling,extractionandmorphologicalidentification
All research conducted in this study complies with relevant ethical
regulations. Fieldwork and sampling were carried out under the
requiredpermitsand permissions,including authorization fromthe
CorporaciónNacionalForestal(CONAF;Permit09/2023)forworkin
Salarde Huasco. The study wasconducted infullcooperation with
localcommunities,includingtheComunidadLickanAntaydeToconao
andtheComunidadIndígenaAymaraLagunadelHuasco,whogranted
accesstotheirterritories.Localresearcherswereactivelyinvolvedin
theresearchprocess.AfieldschoolwasconductedinCollaboration
withChileanresearchers,universityandschoolstudentswheretrain-
ingonsamplingmethods,molecularworkanddataanalysiswascon-
ducted.Additionally,reportsontheresultsobtainedwereprovidedto
theCorporaciónNacionalForestalinSpanishassoonastheresults
were obtained. The research was designed to be locally relevant,
addressingenvironmentalandconservationquestionsofimportance
to the communities. Local and regional literature and prior studies
conductedintheareawerecarefullyconsideredandcitedtoensure
contextualrelevance.
Samplesweretaken in sixdifferentlocalities acrossthedesert.
Herewedefinealocalityasanenvironmentallydistinctarea,identified
based on visible ecological and landscape features. We determined
thembasedondirectfieldobservationandpreviouslydefinedpriority
sites in ref. 117, considering factors such as dominant vegetation,
presence of water bodies and landforms while considering habitat
suitabilityfornematodes.WenamedthelocalitiesasAltiplano(ALT),
Aroma(ARO),EaglePoint(EPT),SalardeHuascoandLagunaGrande
(Salars),Paposo(PAP)andTotoralDunes(TDT).Samplinginthepro-
tectedSalardeHuascowasconductedundertheauthorization(09/
2023) of the Corporación Nacional Forestal (CONAF). In this study,
7 samples from Altiplano, 11 samples from Salars, 15 samples from
Paposo,20samplesfromTotoralDunes,23samplesfromAromaand
21 samples from the Eagle point were analyzed (supplementary
TablesS1–S6).Otherareassampledwherenonematodeswerefound
andarenotpartofthedefinedlocalitiesarelateronincludedinthe
modelingprocess(supplementaryTableS7).
Soilsamplesofapproximately500gweretakeninthedifferent
locationsusingashovel,fromtheupper0–10cmandupto30cmof
eachsiteandstoredinzip-lockbags.Samplestakenneartoplants
corresponded to sediment ~10cm away from the plant stem, or
fromplantlitterinthecaseofPaposo.Theywereinitiallychecked
for the presence/absence of nematodes on-site by extracting a
portion from the soil samples using seeding trays with one layer
Kimberly-Clark Kimtech science precision wipes 7551. Trays were
floodedfromthebottomandsampleswereallowedtore-hydrate
overnight. The water was then filtered through sieves of 80 and
120μm. Nematodes were collected on petri dishes and examined
usingaZeissStemi2000microscope;thisprocedurewasalsofol-
lowed in a laboratory setting for further morphological and
molecularwork.
Adult nematodes weremorphologically identified using a Zeiss
Axioplan2 lightmicroscope and further examined at 10X,20X and
40Xfollowingmultiplediagnosticfeaturesasshowninsupplementary
Fig. S2. Examined individuals were categorized as morphotypes to
accountformorphologicalvariationwithinnematodefamilies.Mor-
photyperepresentativesforeachfamilyfoundonthesampleswere
storedin10μlofnuclease-freewaterandstoredat-20degreesCenti-
gradeforfurtherprocessing.

DNAextraction,barcodingandsequenceprocessing
DNAextractionwasperformedonwholesingleindividualnematodes
usingamodifiedversionoftheHotSHOTDNAextractionmethod118.
PCR was performed using the Cytiva PuReTaq Ready-To-Go™ PCR
BeadsorREDTaq®ReadyMix™PCR-mix.Thegeneticmarker18SrRNA
was amplified using the forward (CGCGAATRGCTCATTACAACAGC)
andreverse(GGGCGGTATCTGATCGCC)primersdescribedinref.119,
targeting the V2, V3 and V4 regions. The thermocycling conditions
usedwere:denaturationat94°Cfor5min,35cyclesofamplification
(94°C for 30s; 54°C for 30 s; 72 °C for 70s), followed by a final
extensionat72°Cfor5min.Gelelectrophoresiswasperformed(1.5%
gel)toassesstheamplicons.PCRproductswerepurifiedusingtheExo-
CIP™RapidPCRCleanupKit(NEB-E1050L).DNAsamplesweresub-
mittedforsequencingtoEurofinsGenomicsandGenewiz.
Sequences were analyzed using Geneious Prime® (v. 2024.0.5).
Onlysequenceswithatleast20%high-qualitybases(HQ%)werekept
forfurtherprocessing.Sequencesweretrimmedattheends(regions
withmorethan1%chanceofanerrorperbase)andsequencesimilarity
wasanalyzedusingBLAST+(v.2.12.0)implementingthecurated18S-
NemaBase120 using the blastn function incorporated in Geneious
Prime.BesthitswereassignedaccordingtothelowestE-valueanda
highpercentageofidentity(listwithaccessionnumbersandsample
informationcanbefoundon10.5281/zenodo.15517342).

Geneticdiversityandhaplotypeuniqueness
Nematodefamiliesrepresentedbyfiveormore18srRNAsequences
from different samples and present in more than one locality were
usedtogeneratehaplotypenetworks.Sequencesweregroupedbased
ontheirsequencesimilaritytoreferencesequencesofthecurated18S-
NemaBasedatabase120,identifiedtofamilyandifpossiblegenuslevel.
Duetolackofresolutioninthesequences,somehaplotypenetworks
werecreatedbasedonfamilylevelinsteadofgenuslevel.Thegrouped
sequenceswerealigned(multiplealignments)usingClustalOmega(v.
1.2.3). Short sequences (<300 bp) or with lack of overlap were
removed. Alignments were manually trimmed in Geneious prime
accordingtosequenceoverlap.Trimmedandcuratedalignmentswere
usedasinputforDnaSP(v.6.12.03)121tocalculatenucleotidediversity
(π)andtheWatterson’sestimator(θ).Haplotypedatafileswerecre-
ated in DnaSP and according to this trait matrices defining the
occurrencesofthedifferenthaplotypesweregeneratedandthenused
asinputforPopART(v.1.7)122.Haplotypenetworkswereconstructed
using the TCS inference method123. Networks were finalized using
Inkscape(v.1.3.2).SequencesetsweredefinedinDnaSPtoestimateπ
foreachlocationwhereenoughsequenceswerepresent(alignments,
haplotypematrixandfinaloutputcanbefoundonhttps://doi.org/10.
5281/zenodo.15517342).

CommunitycompositionthroughouttheAtacama
We estimated the Jaccard disimilarity index (β diversity) using pre-
sence/absencedataoffamiliesfoundin6localities(Altiplano,Aroma,
Eagle Point, Salars, Paposo and Totoral Dunes), here we analyzed
aggregatesofcommunitiesperlocalityanddonotdiscriminatebysoil
sample(bag).Weusedtheveganpackage(v2.6-4)124andvisualizedthe
resultsinverted(Jaccardsimilarity)asaheatmapusingpheatmap(v.
1.0.12)125inR.Wefurtheronanalyzedfamilyco-occurrenceusingthe
cooccur function of the cooccur package (v. 1.3)126. To test for
isolation-by-distance patterns, we estimated distance between all
combinationsofsamplinglocalitiesoncentroidsusingthesf(v.1.0-
16)127,128andgeodist(v.0.1.0)129packages.Following,amanteltestwas
performedcomparing Jaccarddissimilaritytogeographicaldistance
usingtheMantelfunctionoftheveganpackagewiththeSpearman’s
rankcorrelationρ.Thenucleotidediversityforeachlocationofthe
different groups was also tested for isolation-by-distance patterns
using the same approach. Nucleotide diversity was transformed to
geneticdistancebyestimatingtheEuclideandistanceusingthebaseR
packagestats130.

Nematodesassoilindicators
Presence-absenceinformationattheFamilylevelforeachlocalitywas
used as input for the application NINJA: Nematode INdicator Joint
Analysis41.Weassessedthecoloniser-persisterstructureassemblageof
soil nematode families found in the Atacama Desert as well as the
feedingtypecompositionineachofthesamplinglocalities.

Modelingecologicalconstraintsofnematodebiodiversity
We evaluated climatic (environmental, e.g. precipitation) and geo-
morphological(edaphicandgeographic,e,g.topographiccomplexity
andlatituderespectively)variablesaspredictorsforgenerarichnessin
the Atacama Desert. We tested 57 probabilistic models using gen-
eralized least squares linear regression from the package nlme (v.
3.1)131,132onRand57linearmixedeffectmodelsusinglmefromthe
packagenlme,usingthesamplinglocalitiesasarandomeffect.Dataon
genera presence in the Atacama Desert was obtained as described
through the morphological and genetic identification of the nema-
todes. The geomorphological variables evaluated were elevation,
topographiccomplexity,curvature,soilthickness,rocktypeandlati-
tude(supplementaryTableS13).Weobtainedelevationdatafroma
digital elevation model at 15 arc sec resolution133 and calculated a)
topographic complexity as the standard deviation of elevation in a
coarser0.01degreeresolutiongrid134andb)curvatureasthesecond
derivativeofelevationwithrespecttotheaspectina3×3grid135as
implemented in the RichDEM package (v. 0.0.3)136 in Python. Soil
thickness and rock-type data were obtained from publicly available
datasets compiled respectively by Pelletier et al.,137 and Hartmann
etal.,138.Theevaluatedclimaticvariablesweremeanannualtempera-
ture,meanannualprecipitation,temperaturerange,andprecipitation
range. Here, ’range’ is defined as the difference between seasonal
(ratherthandaily)minimumandmaximumvalues,servingasaproxy
forclimaticheterogeneity.Temperatureandprecipitationdatawere
obtained from the openly accessible Chelsea version 2 climate
dataset139,140.Acorrelationmatrixplotwasobtainedusingthecorand
corrplot functions from the corrplot package (v. 0.92)141 to assess
variablecorrelations(>∣0.65∣)(SupplementaryFig.S9),onlyarepre-
sentative variable from highly correlated pairs was included in the
analysistoavoidredundancy.Wefurhterconductedaprinciplecom-
ponent analysis using the prcomp function of the stats package (v.
4.3.2)130, the PCA was visualized using the factoextra package (v.
1.0.7)142 in R. We log-transformed soil thickness (supplementary
Fig.S8)tofulfillmodelassumptions.Bivariaterelationshipsbetween
genus richness and each numeric environmental predictor were
explored with faceted scatter-plots to obtain the coefficient of
determination (R2) using broom (v. 1.0.7)143. Bioinformatic pipelines
wereoptimizedusingChatGPT.
Weevaluatedatotalof114models,includingnullmodels(Sup-
plementaryTablesS14andS15),whichincorporatedvariouscombi-
nationsofselectedpredictivevariables.Halfofthetestedmodelswere
least squares models (GLS), while the other half accounted for the
samevariablesbutincludedlocalitiesasarandomeffecttocontrolfor
influence of the localities (linear mixed effect models–LME). The
variableswereselectedbasedontwocriteria:1)variablesshouldnot
displaystrong correlation with each other,and 2) the variable con-
tributiontoexplainthevariabilityalongthefirstandsecondprincipal
componentsmustbehighandamongthetopvariables(supplemen-
taryFig.S10 and S13).Furthermore, the residuals of each model fit
werevisuallyevaluatedandavariancestructureadjustmentwasper-
formedwherenecessaryfollowing144.Thepreferredpredictivemodels
perapproach(GLSandLME)werethenselectedbasedontheBayesian
informationcriterion(BIC)andfurtherestimationofAkaikeweights
(wm)145,asthelatterprovidestheprobabilitythatamodelisthebest
model according to the set of models tested and the experimental
data146.Thepreferredpredictivemodelswerethenusedtoobtaina
maprepresentingthepossiblenematodediversityintermsofgenera
richness.Wefurtherontestedvariableimportanceofthepreviously
usedvariablesetusingthemachinelearningalgorithmRandomForest
usingthepackagerandomForest(v.4.7-1.1)147tovalidatetheresults
obtainedinthelinearmodels.Trainingandtestingsetwerecreated
usingthersample(v.1.2.1).Themodelwasthenevaluatedusingcaret
(v.6.0-94).Partialdependenceplots(PDPs)werecreatedtovisualize
the modeled effects of key environmental variables on predicted
genus richnessusingthe fitted randomforestmodelusing thepdp
package(v.0.8.2)148.

Reproductivemodedeterminationandassociationtoenviron-
mentalvariables
Laboratorycultureswereestablishedfromnematodesextractedfrom
the soil samples. Worms were propagated on nematode growth
medium(NGM)platescontainingE.colistrainOP50asafoodsource.
Single-strain established cultures were analyzed morphologically to
identifymalesandfemalesofastrain.Cultureswithonlyfemaleswere
categorizedaspossibleparthenogens,andcultureswithbothmales
and females present were categorized as possible gonochoristic. In
addition,singlejuvenilestagenematodecultureswerepropagated(to
ensurenomating)todeterminethereproductivemodeofthestrains
asdescribedinref.149.Between5to30singlejuvenilecultureswere
setupforreproductivemodecorroborationperstrain.Weconsider
asexualcultureswherelarvaegrewintoadultsandlaideggswithout
fertilizationfromamale.
Models testing elevation, precipitation and mean annual pre-
cipitation as predictors for reproductive mode (sexual or asexual).
Environmental data was collected as in the previous section. Two
modeltypeswereemployed:(a)alogisticregressionanalysisusingthe
glmfunction fromthestatspackage in R(v. 4.3.2),with abinomial
family distribution (family = binomial), and (b) a generalized linear
mixed-effectsmodel,implementedwithglmerfromthelme4package
(v.1.1-36),incorporatinglocalitiesasarandomeffecttoaccountfor
spatial variabilityinthe samples.Following,additionalmodelswere
tested including vegetation coverage and percentage of phaner-
ophytesestimatedforthesampleswherereproductivemodecouldbe
assessed.Plantidentificationwasdonethroughmorphologicaliden-
tification.Wefurthercheckedwhichvariablesweremostimportantby
usingaRandomForestmodelwiththerandomForestpackage(v.4.7-
1.1)147,toconfirmthefindingsfromthelogisticregressionandlinear
mixed-effectmodel.Thedatawassplitintotrainingandtestingsets
withrsample(v.1.2.1),andthemodel’sperformancewasthentested
usingcaret(v.6.0-94).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors conducted a survey of nematode communities in 231 soils from six different habitats from the Atacama Desert in Chile. They employed several techniques to gain ecological insights from this survey. First, they characterized the diversity (morphospecies/family/genus) of the communities through morphological analyses, augmented by partial 18S barcoding (haplotypes). Second, they modelled the structure of these communities incorporating several variables, including latitude, elevation, precipitation and temperature. Finally, they studied the reproductive habits of nematode that they were able to culture and correlated that to environmental characteristics.\nKey results:\n1. Characterization of the nematode taxa present in Atacama Desert soil – this is new information as no work has been done on these nematodes that I know of for over 50 years.\n2. Assessment of B-diversity (Jaccard Index) through incorporation of different sites with different characteristics (latitude, elevation, MAP)\n3. Finding a relationship between geographic and genetic distance for the genus Panagrolaimus through evaluation of biogeography and 18S genetic diversity of select taxa.\n4. Identification of positive effect of a latitudinal diversity gradient and increasing precipitation.\n5. A preliminary finding that asexual reproduction is positively correlated to elevation.\nValidity:\nIn so far as key results 1 & 2, the expertise of the team and the combination of morphological and barcoding approaches to identify the nematodes means the data is valid.\nI think that use of the Bongers Index to characterize “poor soil food webs” that are “sensitive to disturbance” is a big stretch. The Index was developed mainly for agricultural systems and to investigate the impact of disturbance in those systems. I think, in the context of this study, the Index has little utility. It’s trying to impose what is known about the ecology of nematodes from eutrophic systems to the oligotrophic desert, and I don’t believe it is valid or necessary. Are any of these sites being actively disturbed in some manner? If not, saying the nematode communities are “sensitive to disturbance” isn’t really valid. While lower cp nematodes might be over-represented in a disturbed agricultural soil, their presence in a desert soils is due to different phenomena.\nRegarding key result 2, I think the modelling effort is interesting, but I found that some of the environmental data used was strange (see below). Discussion of how the soil habitat changes with the variables that were important was absent (eg, how do MAP, latitude, and elevation actually influence soils as a habitat for nematodes).\nRegarding 3: The sequencing effort seems valid and sufficient to be able to examine the haplotype diversity.\nRegarding 4: I think this is very, very interesting data, and I appreciate that the authors were candid in the paper regarding the limited sample size. I think saying in the abstract that the distribution of “asexual taxa appears to be in the marginal ranges of high altitudes” is not valid. As far as I can tell, you found a positive correlation with elevation for the small number of cultures established. The data does not show that asexual taxa are not found at lower elevations.\nSignificance:\nVery few studies have attempted to examine nematode genetic diversity within taxa across sites, as the authors have done here. As well, very little work has been done to connect nematode reproductive strategy to the environment. So, these are good questions, and this is novel work of significance to the soil nematode community.\nData and methodology:\nI think everything seems clear and appropriate.\nAnalytical approach:\nRegarding the modelling effort to investigate environmental variables and nematode diversity: I think some edaphic and vegetation-related variables would have been useful in this effort, and I wondered if some of this data was collected. For example, soil properties (organic C is very important to nematodes) or NDVI. “Plant association” was mentioned at Line 228 as a variable, vegetation coverage and other variables are in Table S11. For me, the inclusion of Rock type and Soil thickness were not justified beyond being publicly available data (and indeed, these variables did not prove important in any of the models – I would have not put them in the paper at all – maybe leave them out of the figure, as done with topographic complexity and curvature, Line 416). To me, I think the significance of latitude, MAP and elevation on nematode communities is likely due to impacts on vegetation and thus soil organic matter (overall soil habitability), and I think these mechanisms should be at least discussed in the paper. I appreciate not wanting to include codependent variables in the modeling effort, however (but aren’t latitude and precipitation codependent? Surprising not, per S9).\nSuggested improvements:\nWhat is a “transect” in this study? It seems that samples were collected haphazardly within sites, but “transect” implies some scheme that was not described in the paper. I would perhaps not describe the sampling as along a transect.\nAssuming the number was different for each site, can you provide a range for the number of samples taken? (Line 345) Also, can you clarify the degree of distance from plants that the samples were collected? Were you under the canopy? Was their litter?\nI think that addition of a table in the paper or supplement that is a full list of the nematode family/genera and number of haplotypes found at each site is essential and would increase the utility of this work.\nI would leave off the Bongers analysis – I don’t think it is informative.\nClarity and context:\nI have made some comments above as to where I think more context would be helpful (interpretation of modelling results).\nIn Figure A: What is “temperature range”? Is that yearly? This is just not a variable I am familiar with.\nLines 153-154: I have read this many times and I cannot make sense of this text, but I feel like maybe it’s important.\nFigure 3: Remove UEF from legend – presence is confusing and unnecessary.\nReferences\nThe paper cites the literature appropriately.\nMy expertise: I feel comfortable reviewing all aspects of the work.", "reviewer_2": "This manuscript presents a comprehensive investigation of nematode diversity in various subregions of the Atacama. Diversity (genetic, taxonomic, reproductive strategies, functional groups) are investigated and statistical approaches used to determine best predictors (elevation, MAP) of those patterns. This is a unique environment, with limited investigations of the soil biodiversity, so this work presents an important contribution to our survey of life in this region. The authors use existing ecological frameworks (e.g. isolation by distance, geographical parthenogenesis) to test specific hypotheses about the drivers of biodiversity of nematodes. The use of modelling approaches could be useful to help identify areas of this region that would be good to focus on for future soil biodiversity surveys.\nThe title is fairly broad, I would suggest including nematodes in the title to better convey the focus of this study.\nA stated goal in the abstract was to study resilience in a unique ecosystem (desert), but I didn’t see how “resiliency” was measured or assessed in this study.\nIs the relationship between diversity and precipitation driven more by the fact that you had clusters of samples from 6 different sites (which also had different soil/mineral types, vegetation etc.)? How much of these other predictor variables were accounted for in the model? Is it possible to conclude that MAP is the driver, or is it simply due to spatial distances? Some focused within-transect modeling/testing may help strengthen and support the conclusion being drawn.\nSimilarly, with the relationship to latitude: I’m not an expert on this region, but the distribution of sample sites shows that the northern sites have more inland representation while the southern sites are closer to the coast – could that be affecting patterns by latitude?\nHow temporally stable are these patterns of diversity? This datasets represents a single time point; how much variation would you expect if you returned to those sites to sample?\nHow many samples were taken from each site? How much local spatial variability is expected at each location?\nIn general, the description of how the data was handled with respect to within- and between-transect analyses needs some more detail to clarify what was done (see specific comments for the methods section)\nSpecific suggestions by line number:\n30 (abstract) - “desert” - lowercase; “indicators” instead of “indicatives”?\n107-108 – capitalize taxa names\n159 – not entirely clear what is meant by “we calculated the latter”.\n159-167 – it wasn’t clear to me if the genetic distance results were within a geographical region or between regions. Unclear what “randomly” means here.\n168 – subtitle of this section indicates random (no driving pattern), but then the paragraph goes on to suggest there is a relationship with geographic distance.\n178 – no comma after “interested”\n324 – “taxa are at an…”\n349 – How many soil samples were taken at each site?\n393 – Need more detail here on how the data was analyzed. How were the 6 transects treated in this analysis? Were the data from each site along the transect combined to get an aggregate community for the transect? Or means of each taxa calculated?\n406 – Same comment as previous - what does “each transect” mean – a combination of the data from each site along a single transect? Mean? Other?\n411 – which edaphic properties were included?\nFig S3 legend seems to be missing a reference", "reviewer_3": "The authors present findings on ecological extremes and the impact on biodiversity, using nematodes as bioindicators.\nThe authors must be commended as this is a nicely written paper on a difficult but important theme and, from an extremely difficult terrain, logistically.\nThe results are original and solid for publication. There are limits to extrapolation from a relatively small sample size but the authors have not overdone any conclusions.\nThe methodology is sound however the M&M section would benefit with added detail on i) L351 how were the initially assessments (for nematodes on-site) done; ii) how many samples were taken at each of the six locations; iii) whether any systematic process (or not) was used to select individual sampling locations (bias?); and iv) it would be great to have a table with all 194 sampling GPS coordinates (if available).\nOne bit of advice is to re-read and be sure that it is nice and clear that this is a paper on the relationships between habitats and biodiversity, not a nematode paper. At times it strays a little too much about the nematode. But this is not a blocker.\nThere are a few typos so a couple of re-reads are required (e.g. L186, 236 and 248)"}, "author_response": "Reply: We thank the reviewer for raising an important point in regard to “sensitive areas” in the Atacama Desert, as these are not imminently obvious to the reader. However, there are a multitude of threats to the ecosystem(s) there, such as for example the activity of mining or farming. For sites in our study, these threats were present in particular where more “sensitive” taxa were found. Consequently, we now mention mining, plant extraction, and uncontrolled tourism in the text for a better explanation (lines 311-316).\nFurthermore, we are aware that the Bongers index was initially developed in an agricultural context. Our main goal in using the index is to provide a base line for further studies as this can impact conservation efforts in the area. We believe the index is indeed meaningful in this, as several publications (e.g. https://doi.org/10.3897/mbmg.8.111307, https://doi.org/10.3389/fevo.2016.00084, https://doi.org/10.1002/ece3.8061) have shown its utility in natural systems. The advantage of using an established index is that the interpretation of the value ranges is known, making it possible to understand the findings in the context of better-studied systems. Nevertheless, we agree that the oligotrophic conditions of the desert environment can have an impact on the index and for this we refer to a rather relative interpretation of our findings between sites within our study (lines 309-317).\nRegarding key result 2, I think the modelling effort is interesting, but I found that some of the environmental data used was strange (see below). Discussion of how the soil habitat changes with the variables that were important was absent (eg, how do MAP, latitude, and elevation actually influence soils as a habitat for nematodes).\nReply: We thank the reviewer for pointing this out. We have now included more information about the effect of the main drivers on soil communities to complement the previous text. In the case of genera richness (taxonomic diversity per site), we highlight that precipitation is used as a proxy for water availability as done by other researchers previously, and how higher precipitation improves nutrient availability in soils, promotes higher bacterial densities (food source for nematodes) and increases soil moisture, which is necessary for nematodes to move as they live in the water film between soil particles (lines 335-338, 342-347). We also mention the effect of of temperature variation (lines 348-354) and explain the theory of latitudinal diversity gradient, where underlying mechanisms for the pattern remain unclear (lines 355-369). Particularly for elevation,\nRegarding key result 3: The sequencing effort seems valid and sufficient to be able to examine the haplotype diversity.\nReply: We want to mention here that some of the sequences in the analysis have now been removed due to insufficient quality and resolution for uploading them to Genbank. This reduced the number from 393 to 386. Hence, we repeated the calculations, but the conclusions remain the same.\nRegarding 4: I think this is very, very interesting data, and I appreciate that the authors were candid in the paper regarding the limited sample size. I think saying in the abstract that the distribution of “asexual taxa appears to be in the marginal ranges of high altitudes” is not valid. As far as I can tell, you found a positive correlation with elevation for the small number of cultures established. The data does not show that asexual taxa are not found at lower elevations.\nReply: We agree that the statement was not a true representation of the results and have corrected it in the text. This section now states, “We also find that distribution of asexual taxa is more likely to occur at higher altitudes” (line 30 -abstract), and we highlight in the discussion that there is a higher probability of asexual taxa to occur at higher elevations, this however does not imply that asexual taxa can’t be found at lower elevations. This can also be corroborated on 4 (D-E) where the probability of being asexual is closer to 1 (light orange to beige) with increasing elevation.\nSignificance:\nVery few studies have attempted to examine nematode genetic diversity within taxa across sites, as the authors have done here. As well, very little work has been done to connect nematode reproductive strategy to the environment. So, these are good questions, and this is novel work of significance to the soil nematode community.\nData and methodology:\nI think everything seems clear and appropriate.\nAnalytical approach:\nRegarding the modelling effort to investigate environmental variables and nematode diversity: I think some edaphic and vegetation-related variables would have been useful in this effort, and I wondered if some of this data was collected. For example, soil properties (organic C is very important to nematodes) or NDVI. “Plant association” was mentioned at Line 228 as a variable, vegetation coverage and other variables are in Table S11. For me, the inclusion of Rock type and Soil thickness were not justified beyond being publicly available data (and indeed, these variables did not prove important in any of the models – I would have not put them in the paper at all – maybe leave them out of the figure, as done with topographic complexity and curvature, Line 416). To me, I think the significance of latitude, MAP and elevation on nematode communities is likely due to impacts on vegetation and thus soil organic matter (overall soil habitability), and I think these mechanisms should be at least discussed in the paper. I appreciate not wanting to include codependent variables in the modeling effort, however (but aren’t latitude and precipitation codependent? Surprising not, per S9).\nReply: We agree that plants can play an important role in the determination of nematode’s communities, and this is why we tried to include this information at a finer scale in the reproductive mode models from on-site measurements. For these models, we recorded the major plant groups present in the sampling localities, percentage of vegetation cover and type of coverage (shrubs, cacti and herbs). However, for the genera richness determination, we couldn’t define such fine scales (overall for the environmental variables tested) given that the reanalysis data implemented doesn’t not provide resolution at the sampling spot but rather on a grid. We discussed implementing NDVI information, but unfortunately, this satellite-based measures seem to have a large discrepancy between the real vegetation and the recorded vegetation, for instance, shrubs which are a major component of the areas we analysed are not recorded with NDVI. We thus rely on our more precise records, based on repeated extensive sampling efforts.\nWith regard to the choice of environmental variables, we have included more information in the discussion, why we selected rock type and soil thickness as variables that could influence nematode diversity (lines 355-359). We have also included a statement acknowledging that plant presence is an important factor for nematode diversity but clarify that we haven’t looked into this for the genera richness aspect (lines 359-363). In short, we do consider plants to be an important factor in determining community diversity for nematodes, but have not included this in the modeling approach for the above mentioned reasons.\nWhile latitude and precipitation show a week positive correlation (0.43), we didn’t consider the exclusion of this variables since the correlation was bellow |0.65| For instance, mean annual temperature and elevation had a negative correlation of 0.79, hence only elevation was considerd. This is now clarified in the methods (lines 517-519)\nSuggested improvements:\nWhat is a “transect” in this study? It seems that samples were collected haphazardly within sites, but “transect” implies some scheme that was not described in the paper. I would perhaps not describe the sampling as along a transect.\nReply: We thank the reviewer for pointing this out. It is true that we selected environmentally distinct areas, identified based on ecological and landscape features (e.g. dominant vegetation, presence of water bodies, different and landforms), rather than following strict transects (apart from Eagle Point, which is a transect on a mountain slope). We have been using the term transect as our sampling localities align with defined transects that have been established for the collaborative research cluster our project is part of (doi.org/10.1016/j.gloplacha.2020.103275). We have now changed the text to use the term “localities” and we specify criteria for selecting these as mentioned above.\nAssuming the number was different for each site, can you provide a range for the number of samples taken? (Line 345) Also, can you clarify the degree of distance from plants that the samples were collected? Were you under the canopy? Was their litter?\nReply: We have included more details about sample collection in the Methods section (lines 412-425), i.e. samples taken near to plants corresponded to sediment approximately 10 cm away from the plant stem, or from plant litter in the case of Paposo. Plant litter are drying and dried leaves, twigs, potentially dried fruit found under a plant (explained in lines 434-435 in the text now). We state the number of samples per locality: 7 samples from Altiplano, 11 samples from Salars, 15 samples from Paposo, 20 samples from Totoral Dunes, 23 samples from Aroma and 21 samples from the Eagle point (lines 428-429) and provide sampling coordinates of all sampling spots (Supplementary tables S1 - S7)\nI think that addition of a table in the paper or supplement that is a full list of the nematode family/genera and number of haplotypes found at each site is essential and would increase the utility of this work.\nReply: This information is now included in the supplementary files as table S10. The table, sorted by locality, contains the number haplotypes per taxonomic group. One sequence (accession number PQ587770.1) was initially identified as Acrobeloides but after curation corrected to Acrobeles Hap_12. Due to uncertainties in the differentiation between Aporcelaimidae, Qudsianematidae and Pararhyssocolpidae based on the current database only the order level was assigned for most identifications.\nI would leave off the Bongers analysis – I don’t think it is informative.\nReply: See our comments above for a justification of the use of this index. Nevertheless, we have reduced the emphasis on these results.\nClarity and context:\nI have made some comments above as to where I think more context would be helpful (interpretation of modelling results).\nIn Figure A: What is “temperature range”? Is that yearly? This is just not a variable I am familiar with.\nReply: We have included an explanation of how the temperature range was defined: as the difference between seasonal (rather than daily) minimum and maximum values, serving as a proxy for climatic heterogeneity.\nLines 153-154: I have read this many times and I cannot make sense of this text, but I feel like maybe it’s important.\nReply: This statement has been changed to make it clearer “Overall, we were not able to find any distinct geographical clustering of haplotypes based on the networks, i.e. individual haplotypes were present in geographically distant localities” (lines 154-155).\nFigure 3: Remove UEF from legend – presence is confusing and unnecessary.\nReply: UEF has been removed from the legend.\n2. Reviewer #2 (Remarks to the Author):\nThis manuscript presents a comprehensive investigation of nematode diversity in various subregions of the Atacama. Diversity (genetic, taxonomic, reproductive strategies, functional groups) are investigated and statistical approaches used to determine best predictors (elevation, MAP) of those patterns. This is a unique environment, with limited investigations of the soil biodiversity, so this work presents an important contribution to our survey of life in this region. The authors use existing ecological frameworks (e.g. isolation by distance, geographical parthenogenesis) to test specific hypotheses about the drivers of biodiversity of nematodes. The use of modelling approaches could be useful to help identify areas of this region that would be good to focus on for future soil biodiversity surveys. The title is fairly broad, I would suggest including nematodes in the title to better convey the focus of this study.\nReply: We have included the term “Nematode” in the title for clarity. However, we would like to point to the comments by reviewer 3 here, who asked for less focus on nematodes in the text in general. We have now aimed to balance this better, acknowledging that our study is based on nematodes, but inferences can potentially be extended to soil invertebrates.\nA stated goal in the abstract was to study resilience in a unique ecosystem (desert), but I didn’t see how “resiliency” was measured or assessed in this study. Is the relationship between diversity and precipitation driven more by the fact that you had clusters of samples from 6 different sites (which also had different soil/mineral types, vegetation etc.)? How much of these other predictor variables were accounted for in the model? Is it possible to conclude that MAP is the driver, or is it simply due to spatial distances? Some focused within-transect modeling/testing may help strengthen and support the conclusion being drawn.\nReply: We previously hadn’t consider within transect testing, and we thank the reviewers for pointing it out as an important concern, as there could be clustering from a sampling bias. In order to address this, we went back to our data and re-grouped the samples per locality (formerly named transects) resulting in a total of 112 samples. Specifically, we combined some subsamples of the Eagle Point Transect as these were taken in very close proximity with no distinctive ecological features (e.g. no proximity to different plant species, rock or different soil type). This ensures better comparability of this sampling point with all other samples in the analysis.\nWe then tested our dataset, as previously done in the first version of the manuscript, using a test of normality, PCAs, and variable correlation, before proceeding test the GLS and Random Forest models again. This test confirmed our previous conclusion that according to the random forest model the most important variables in determining genera richness are elevation, latitude and mean annual precipitation (with the addition of temperature range). We can accordingly exclude a clustering bias in our data.\nAs the GLS model seems to overestimate genera richness, not accounting for sampling bias and locality selection, we tested 57 more models using a linear mixed-effects model approach. Overall, these confirmed that climate is a predictor for genera richness in the desert given. See lines 524-527 (Methods), lines 219-229 (Results), and lines 335-354 (Discussion).\nRegarding spatial distances and their potential role as the primary driver of diversity instead of mean annual precipitation, we found a moderate isolation-by-distance pattern when analyzing community dissimilarity (Mantel test: ρ=0.5153ρ=0.5153, p=0.048611, 719 permutations) (lines 176–178). However, this pattern depends on the specific taxa inhabiting each area. Since our model focuses on the number of different taxa rather than specific genera or families, we believed this didn’t directly impact the results. We have now performed a correlation test specifically between geographic distance of the localities and the average precipitation per locality, the result showed that there is no significant correlation between these two characteristics (rho= 0.1321, p= 0.24861, 719 permutations), discarding that geographical distance is correlated to precipitation, geographical distance does not reflect the precipitation patterns. Moreover, thanks to reviewer comments with included a modeling approach accounting for sampling locality as depicted above.\nWe also tested the effect of the locality in the reproductive mode models, but this had no significant effect.\nSimilarly, with the relationship to latitude: I’m not an expert on this region, but the distribution of sample sites shows that the northern sites have more inland representation while the southern sites are closer to the coast – could that be affecting patterns by latitude? How temporally stable are these patterns of diversity? This datasets represents a single time point; how much variation would you expect if you returned to those sites to sample?How many samples were taken from each site? How much local spatial variability is expected at each location? In general, the description of how the data was handled with respect to within- and between-transect analyses needs some more detail to clarify what was done (see specific comments for the methods section)\nREPLY: We have performed the analysis including locality as a random effect as detailed above. Furthermore, it is important to consider that in Chile the distribution from the coast to the more in-land locations are well-described by the altitudinal variation. The coastal cordillera is marked by a large altitudinal increase closer to the coast and the altiplano then marks the plateau area at higher elevation before the uprise of the Andes. That way, the longitudinal distribution of localities is reflected by variation in altitude.\nWe have clarified how stable we expect this patterns to be: on shorter time-scales climatic variation is present (seasonal patterns), our results reveal nematode diversity across different levels but do not provide insights into abundance, which can be expected to be highly variable across seasons (e.g., during the Altiplano rainfall season). Nevertheless, we expect that taxonomic diversity remains relatively stable throughout the year for the localities and sampling points analyzed, given that several taxa are capable of undergoing anhydrobiosis (lines 400-405).\nFor better clarity, the number of samples per locality is now specified in the materials and methods (lines 428-431).\nSpecific suggestions by line number:\n30 (abstract) - “desert” - lowercase\nReply: Has been changed\n“indicators” instead of “indicatives”?\nReply: Has been changed\n107-108 – capitalize taxa names\nReply: Has been changed: Proteobacteria and Actinobacteria.\n159 – not entirely clear what is meant by “we calculated the latter”.\nReply: We changed it to Eucledian distance for clarity (lines 160-161).\n159-167 – it wasn’t clear to me if the genetic distance results were within a geographical region or between regions. Unclear what “randomly” means here.\n168 – subtitle of this section indicates random (no driving pattern), but then the paragraph goes on to suggest there is a relationship with geographic distance.\nReply: We thank the reviewer for pointing this out, we have rephrased the title for a better representation of moderate effect of geographical distance.\n178 – no comma after “interested”324 – “taxa are at an…”\nReply: Has been corrected in the text.\n349 – How many soil samples were taken at each site?\nReply: We have included the number of samples taken a teach spot in the materials and methods (lines 428-429).\n393 – Need more detail here on how the data was analyzed. How were the 6 transects treated in this analysis? Were the data from each site along the transect combined to get an aggregate community for the transect? Or means of each taxa calculated?\nReply: We acknowledge this wasn’t clear in the methods, we have included a statement clarifying that we analyzed the communities in the localities (see above for a description on why we use the term “locality” now) as aggregates since we only relied on presence absence data but not abundances. This information can be found on lines 480-482.\n406 – Same comment as previous - what does “each transect” mean – a combination of the data from each site along a single transect? Mean? Other?\nReply: We have now specified how we selected and defined the different sampling localities. The rationale is explained in a reply to ref 1 above and in the text in lines 421-425.\n411 – which edaphic properties were included?Fig S3 legend seems to be missing a reference\nReply: We analysed the soil thickness (depth to mother rock) as edaphic property and the type of mother rock (parental material). The reference in the supplementary figure has been included.\nReviewer #2 (Remarks on code availability):\nREADME file is blank\nReply: We apologize for this faulty file. The readme file has been updated a now includes all information where each of the analysis can be found.\n3. Reviewer #3 (Remarks to the Author):\nThe authors present findings on ecological extremes and the impact on biodiversity, using nematodes as bioindicators.\nThe authors must be commended as this is a nicely written paper on a difficult but important theme and, from an extremely difficult terrain, logistically.\nThe results are original and solid for publication. There are limits to extrapolation from a relatively small sample size but the authors have not overdone any conclusions.\nThe methodology is sound however the M&M section would benefit with added detail on\ni) L351 how were the initially assessments (for nematodes on-site) done\nReply The procedure of on-site identification of presence-absence of nematodes was the same as for the laboratory checking for further morphological and molecular work (we brought the Zeiss Stemi 2000 to the field). This is now also clarified in the text (line 435).\nii) ii) how many samples were taken at each of the six locations\nReply: We have included the number of samples taken a teach spot in the materials and methods (lines 428-429)\niii) whether any systematic process (or not) was used to select individual sampling locations (bias?)\nReply: We have now specified how we (systematically) selected and defined the different sampling localities. The rationale is explained in a reply to reviewer 1 above and in the text in lines 421-428. In general sampling was conducted in a systematic way by sampling in straight lines with evenly spaced out distances in at each of the localities.\niv) it would be great to have a table with all 194 sampling GPS coordinates (if available).\nReply: We have included 6 supplementary tables displaying the exact coordinates of each sample analysed from the different localities as well as one table with coordinates from samples outside of those defined areas (tables S1-S6 of supplementary material).\nReply: As we decided to group samples in localities and test for sampling bias using extra 57 models with a linear mixed-effects approach that account for the effect of the sampling locations, we noticed that subsets of samples from EPT should be grouped together: e.g. EPTA01, EPTA02, EPTA03 have been group as sample EPTA to account for 500g of sediment in specific coordinate set that shares all identical environmental measures. This samples were analyzed in subsets, hence we initially had them labelled differently, but noticed upon revision that this wouldn’t guarantee fully comparability with all other samples. In the analysis they are now grouped as explained above (grouped by close proximity where no distinctive features, e.g. plant nearby, were present).\nTherefore, the number of total samples is 112 as can be seen in the tables with coordinates. All analysis has been repeated where this could have an impact, and it has been specified in lines 524-526 (Methods), lines 219-229 (Results), and lines 342-354 (Discussion).\nOne bit of advice is to re-read and be sure that it is nice and clear that this is a paper on the relationships between habitats and biodiversity, not a nematode paper. At times it strays a little too much about the nematode. But this is not a blocker.\nReply: Thank you for the important suggestion, we appreciate it, as it is far too easy to drift off into describing these fascinating animals. We have gone through the text to shift the focus even more towards habitats and biodiversity. However, to make clear that we used nematodes as our study system, and based on the suggestion of reviewer 2, we added the term nematode to the title.\nThere are a few typos so a couple of re-reads are required (e.g. L186, 236 and 248)\nReply: We combed through the text to fix typos. Thank you for pointing these out."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I will not repeat my remarks from an earlier review. The authors have mainly addressed any suggestions or questions that I had.\nI did request that the authors provide a full list of nematode taxa found. In the response, they say this information is in STable10, but it is not. That table only lists taxa that were subjected to 18S analysis. The paper reports that 21 families (36 genera) were found. This is the list I was expecting to see.\nThe paper still needs to be carefully read for comma placement and punctuation. Things that should be capitalized are not (eg, Firmicutes). Capitalization in the references needs attention.", "reviewer_2": "The authors have done a thorough job of addressing reviewer comments in the revised manuscript. In particular, the additional modelling efforts to address within site variability and the attention given to other co-variates (e.g. vegetation, soils) in the discussion has strengthened the manuscript.", "reviewer_4": "This study investigated soil nematode diversity across multiple levels - genetic, taxonomic, community, and life-cycle characteristics - in various Atacama Desert habitats including dune systems, high-altitude mountains, saline lakes, river valleys, and fog oases. The results that genus-level richness of soil nematodes exhibited a distinct latitudinal diversity gradient, showing positive correlation with increasing precipitation levels. The annual precipitation and temperature heterogeneity are the major drivers of genera richness. These findings are not surprising. The most interesting results is that asexual taxa is more likely to occur at higher altitudes, suggesting that elevation serves as the strongest predictor of reproductive strategy.\nMy major concerns:\n1. There are too much supplementary tables and figures. It seemed that there were lots of data. However, the present tabular data and graphical representations prove challenging to decipher. The abundances of total nematodes, different trophic groups, family/genera/species, and the biodiversity of each sample are not enumerated, and the relationships between nematode variables and environmental variables are also unclear. Figure 4 and figure 1 are ambiguous. Figures of regression analysis between nematode and environmental variables are better than such figures.\n2. How the nematode species was determined as asexual taxa or parthenogenesis?\nDetailed comments:\nL94-103: There is no need to describe the finding here.\nFigure4: I am unfamiliar with this analysis. Which model makes the most accurate predictions?\nL415-416: Two paragraphs?\nL435: Was the top 0-10cm soil removed? Only the 10-30cm soil was collected?\nFig. S2: Panel D is missing."}, "author_response": "Point-by-point response to reviewers\nReviewer #1 (Remarks to the Author):\nI will not repeat my remarks from an earlier review. The authors have mainly addressed any\nsuggestions or questions that I had.\nI did request that the authors provide a full list of nematode taxa found. In the response, they\nsay this information is in STable10, but it is not. That table only lists taxa that were subjected\nto 18S analysis. The paper reports that 21 families (36 genera) were found. This is the list I was\nexpecting to see.\nThe paper still needs to be carefully read for comma placement and punctuation. Things that\nshould be capitalized are not (eg, Firmicutes). Capitalization in the references needs attention.\nWe thank the reviewer for their thoughtful comments and for highlighting the lack of clear access\nto the requested dataset. While the data were previously included as part of the model input in\nCode Ocean (10.24433/CO.6395549.v2), we acknowledge that this was not an intuitive or easily\naccessible location. To improve transparency and accessibility, we have now deposited the dataset\nin a Zenodo repository (https://doi.org/10.5281/zenodo.14930774), under the file\nname \"FamilyGeneraPerSample\" . This table contains sample coordinates, sample IDs, and the\nidentified families and genera. We have also revised the main text to clearly reference this dataset\n(line135).\nReviewer #2 (Remarks to the Author):\nThe authors have done a thorough job of addressing reviewer comments in the revised manuscript.\nIn particular, the additional modelling eIorts to address within site variability and the attention given\nto other co-variates (e.g. vegetation, soils) in the discussion has strengthened the manuscript.\nWe would like to thank again the reviewer #2, we also feel implementing his previous comments\nimproved the analysis.\nReviewer #4 (Remarks to the Author):\nThis study investigated soil nematode diversity across multiple levels - genetic, taxonomic,\ncommunity, and life-cycle characteristics - in various Atacama Desert habitats including\ndune systems, high-altitude mountains, saline lakes, river valleys, and fog oases. The results that genus-\nlevel richness of soil nematodes exhibited a distinct latitudinal diversity gradient, showing positive\ncorrelation with increasing precipitation levels. The annual precipitation and temperature\nheterogeneity are the major drivers of genera richness. These findings are not surprising. The most\ninteresting results is that asexual taxa is more likely to occur at higher altitudes, suggesting that\nelevation serves as the strongest predictor of reproductive strategy.\nMy major concerns:\n1. There are too much supplementary tables and figures. It seemed that there were lots\nof data. However, the present tabular data and graphical representations prove\nchallenging to decipher. The abundances of total nematodes, diIerent trophic groups,\nfamily/genera/species, and the biodiversity of each sample are not enumerated, and\nthe relationships between nematode variables and environmental variables are also\nunclear. Figure 4 and figure 1 are ambiguous. Figures of regression analysis between\nnematode and environmental variables are better than such figures.\nWe thank the reviewer for highlighting areas where the presentation of results could be\nimproved for clarity. As noted in our response to Reviewer #1, a table listing all identified\nfamilies and genera for each sampling location is now explicitly referenced in the main text\n(line 135) and is available via Zenodo (https://doi.org/10.5281/zenodo.14930774).\nRegarding nematode abundance, our study does not provide quantitative insights due to\nthe temporal variability in sampling (i.e., different seasons), which could introduce\ninconsistencies in abundance data. Instead, we focused on identifying the adult taxa\npresent in each sample. As discussed in the manuscript (lines 409–414), while we expect\ntaxonomic composition to remain relatively stable throughout the year, nematode\nabundance is likely to fluctuate substantially with seasonal dynamics, particularly during\nevents such as the Altiplano rainfall season.\nTo improve the clarity of Figures 1 and 4, we have revised their legends to better articulate\ntheir purposes: Figure 1 contextualizes the climatic and geological (geomorphological)\ncharacteristics of the Atacama Desert, while Figure 4 illustrates how predictions vary\nspatially depending on the modeling approach and variables used. Additionally, to address\nthe reviewer’s concern about data presentation, we have added the figures of the\nregression analysis between nematode diversity and environmental variables in the\nsupplementary material (Supplementary Figure S11). These plots highlight the generally\nweak linear associations when considering individual predictors, which is the main reason\nwe used modelling approaches less sensitive to these non-linear effects, such as random\nforest. To complement this, we include partial dependence plots (Supplementary Figure\nS12), which illustrate the non-linear effects of key predictors within the random forest\nmodel framework.\n2. How the nematode species was determined as asexual taxa or parthenogenesis?\nNematodes reproductive mode was assessed following the protocol from Lewis 209 and\nthis is specified in the text in lines 555-563.\nDetailed comments:\n3. L94-103: There is no need to describe the finding here.\nWe appreciate the reviewer’s suggestion. However, we would prefer to retain this section in the\nintroduction, as we believe it helps set the stage for the reader by providing context and a clear\nsense of the study's direction and relevance.\n4. Figure4: I am unfamiliar with this analysis. Which model makes the most accurate\npredictions?\nWe have included a statement acknowledging that, due to the di`ering assumptions underlying\neach modeling approach and the limited set of observations used, it is not possible to select a\nsingle “best” model at the moment. Therefore, we focus on identifying consistent trends across\nmodelling approaches. Notably, genus richness shows responses to environmental gradients\nthat are not fully captured by linear models alone. Thus, at this stage, random forest models\ncould be considered as the most accurate, due to the nonlinear e`ects observed in the\nobservations. However, temperature heterogeneity consistently emerged as a strong predictor\nacross all modeling approaches. Additionally, mean annual precipitation was a significant\nfactor in both linear mixed-e`ects and random forest models, highlighting the key role of\nclimate in shaping genus richness (lines 233–238). These preferred models are contingent on\nour set of observations, as new information about nematode diversity emerges from the region,\nthese models should be revised with the up-to-date observations.\nL415-416: Two paragraphs?\n5. L435: Was the top 0-10cm soil removed? Only the 10-30cm soil was collected?\nWe now clarify in the text (line 443) that samples were collected from depths ranging\nbetween the upper 10 cm (0–10 cm) and up to 30 cm, depending on soil characteristics.\nThe final depth was determined by factors such as soil coarseness and the point at which\nrock was encountered, limiting further sampling.\n6. Fig. S2: Panel D is missing.\nThe figure has been fixed and the letter labelling is correct now."}]

Decision letter