UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO POSGRADO EN CIENCIAS BIOLÓGICAS FACULTAD DE CIENCIAS Diferenciación de nicho a lo largo de un gradiente de estrés hídrico en un pastizal semiárido TESIS QUE PARA OPTAR POR EL GRADO DE: DOCTORA EN CIENCIAS PRESENTA: Alejandra Martínez Blancas TUTOR PRINCIPAL DE TESIS: Dr. Carlos Martorell Delgado Facultad de Ciencias, UNAM COMITÉ TUTOR: Dr. Jorge Arturo Meave del Castillo Facultad de Ciencias, UNAM Dra. Karina Boege Paré Instituto de Ecología, UNAM Ciudad Universitaria, CD. MX. 2023 UNAM – Dirección General de Bibliotecas Tesis Digitales Restricciones de uso DERECHOS RESERVADOS © PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL Todo el material contenido en esta tesis esta protegido por la Ley Federal del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). El uso de imágenes, fragmentos de videos, y demás material que sea objeto de protección de los derechos de autor, será exclusivamente para fines educativos e informativos y deberá citar la fuente donde la obtuvo mencionando el autor o autores. Cualquier uso distinto como el lucro, reproducción, edición o modificación, será perseguido y sancionado por el respectivo titular de los Derechos de Autor. UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO POSGRADO EN CIENCIAS BIOLÓGICAS FACULTAD DE CIENCIAS Diferenciación de nicho a lo largo de un gradiente de estrés hídrico en un pastizal semiárido TESIS QUE PARA OPTAR POR EL GRADO DE: DOCTORA EN CIENCIAS PRESENTA: Alejandra Martínez Blancas TUTOR PRINCIPAL DE TESIS: Dr. Carlos Martorell Delgado Facultad de Ciencias, UNAM COMITÉ TUTOR: Dr. Jorge Arturo Meave del Castillo Facultad de Ciencias, UNAM Dra. Karina Boege Paré Instituto de Ecología, UNAM Ciudad Universitaria, CD. MX. 2023 COORDINACIÓN DEL POSGRADO EN CIENCIAS BIOLÓGICAS FACULTAD DE CIENCIAS DIVISIÓN ACADÉMICA DE INVESTIGACIÓN Y POSGRADO OFICIO FCIE/DAIP/645/2022 ASUNTO: Oficio de Jurado M. en C. Ivonne Ramírez Wence Directora General de Administración Escolar, UNAM P r e s e n t e Me permito informar a usted que en la reunión ordinaria del Comité de Posgrado en Ciencias Biológicas, celebrada el día 15 de agosto de 2022 se aprobó el siguiente jurado para el examen de grado de DOCTORA EN CIENCIAS de la estudiante MARTINEZ BLANCAS ALEJANDRA con número de cuenta 304586916 con la tesis titulada: “Diferenciación de nicho a lo largo de un gradiente de estrés hídrico en un pastizal semiárido”, realizada bajo la dirección del DR. CARLOS MARTORELL DELGADO: Presidente: DRA. MARIA DEL CARMEN MANDUJANO SANCHEZ Vocal: DR. JUAN ENRIQUE FORNONI AGNELLI Vocal: DRA. JOSÉ ALBERTO BURQUEZ MONTIJO Vocal: DR. ANDRÉS LIRA NORIEGA Secretario: DR. KARINA BOEGE PARÉ Sin otro particular, me es grato enviarle un cordial saludo. A T E N T A M E N T E “POR MI RAZA HABLARÁ EL ESPÍRITU” Ciudad Universitaria, Cd. Mx., a 06 de diciembre de 2022 COORDINADOR DEL PROGRAMA Agradecimientos Al Posgrado en Ciencias Biológicas, UNAM. Al Consejo Nacional de Ciencia y Tecnología (CONACyT) por la beca otorgada durante los estudios de doctorado. Al programa de Apoyo de Investigación e Innovación Tecnológica (PAPIIT/DGAPA UNAM) por el financiamiento para la realización de este proyecto (PAPIIT IN-212618 a cargo del Dr. Carlos Martorell Delgado) Al Dr. Carlos Martorell por la dirección de la tesis. A los miembros el comité tutor, la Dra. Karina Boege Paré y al Dr. Jorge Meave del Castillo por la revisión y los comentarios y las observaciones hechas al trabajo. Agradecimientos a título personal Al jurado y a mi comité tutoral, Karina Boege y Jorge Meave, por leer esta tesis y por sus comentarios siempre acertados. A la comunidad de Concepción Buenavista por su apoyo y su amistad para desarrollar este proyecto de tesis. A Marco Antonio Romero por su apoyo en el laboratorio. A mis compañeros de laboratorio que pronto se convirtieron en buenos amigos. Gracias por los diálogos, las discusiones, el trabajo y la diversión en campo, la comprensión, las fiestas y las risas. Gracias por mil comidas sentados todos juntos en la pequeña mesa del laboratorio. Gracias en particular a Verónica Zepeda, Fernanda Herce, Marcela Martínez, Luis Boullosa, Fernando Pedraza, Carmen Vázquez, Ian Belaustegui, Lizeth Luna, Gonzalo Martínez, Emiliano Reyes, Sandra Piña, Cristina Alonso, Edgar González, Rodrigo Muñoz y Gerardo Cervantes. A mis amigos por acompañarme no solo a lo largo del doctorado si no toda la vida. Gracias por infinitas risas, diversión, complicidad y hasta por las lágrimas que compartimos. Agradezco a Alejandra Dávila, Bárbara Casillas, Daniel Méndez, Andrea Ibarra, Tatiana Brofft, Jessica Betancourt, Dulce Kadise, Sofía Solar, Andrea Arévalo, Jonathan Solórzano, Javier Noyola, Carlos Terrón e Isabel Rojas. A mi abuela Gloria por todo su amor, por ser mi mayor cómplice de niña y uno de mis más grandes apoyos de grande. A Diego García por ser mi principal compañero durante el doctorado, por su amor, comprensión y eterno apoyo durante las partes más complicadas por compartir sudor y lágrimas en el campo y durante la modelación y por siempre estar ahí cuando me caía y me tenía que levantar. A Carlos Martorell por ser un tutor excepcional, por siempre impulsarme a cumplir mis proyectos sin importar lo ambiciosos que puedan ser. Le agradezco infinitamente su paciencia, apoyo y su amistad. A mis papás, Rosario y Roberto, por su amor incondicional y por todo el apoyo que me han brindado no solo a lo largo del doctorado si no en cada momento de mi vida. Reconozco que sin su esfuerzo no estaría cumpliendo esta meta. Me proporcionaron un número enorme de herramientas para alcanzar mis sueños y me es imposible poner en palabras el agradecimiento que les tengo. Índice Resumen……………………………………………………… …………………………1 Abstract……………………………………………………………………………………3 Capítulo I. Introducción general………………...…………………. ……………………5 Teoría moderna de la coexistencia………………………………. ………………6 Implicaciones de las condiciones ambientales …………………. ……………….7 Agrupamiento de especies en las comunidades …………………………………10 Objetivos y preguntas de investigación …………………………………………………13 Estructura de la tesis …………………………………………………….. ……………16 Referencias bibliográficas ………………………………………………………… ……20 Capítulo II. Changes in niche differentiation and environmental filtering over a hydric stress gradient………………………………………………………………. …… ……..26 Capítulo III. Species alliances and hidden niche dimensions drive species clustering along a hydric gradient in a semiarid grassland ………………………..…………. …………49 Capítulo IV. Discusión general ……………………………………………….….……139 Conclusiones ………………………………………………………………….…..……147 Referencias bibliográficas ………………………………………………… …..………149 1 Resumen Entender cómo coexiste la diversidad de especies que observamos en la naturaleza ha sido uno de los principales retos de la ecología. La teoría moderna de la coexistencia establece que las especies pueden coexistir mediante dos mecanismos antagónicos, los mecanismos estabilizadores y los mecanismos igualadores. Cuando operan los mecanismos estabilizadores, las especies pueden coexistir si difieren en su nicho, de manera que no interfieren en la adquisición de recursos de otras. Cuando esto ocurre, la competencia intraespecífica será más fuerte que la competencia interespecífica. Los mecanismos igualadores hacen referencia a las capacidades competitivas de las especies. Se postula que la coexistencia solo es posible cuando las diferencias de nicho son suficientemente grandes como para superar las diferencias en las capacidades competitivas. Las condiciones ambientales estresantes pueden modificar a las interacciones intra- e interespecíficas. Sin embargo, no se ha estudiado cómo estas condiciones ambientales afectan a la coexistencia de especies y en particular a la diferenciación de nicho. La acción simultánea de los mecanismos estabilizadores e igualadores puede generar un patrón denominado neutralidad emergente, el cual consiste en la formación de grupos de especies similares en las comunidades. En este escenario, las especies podrían coexistir al interior de los grupos gracias a las similitudes en sus habilidades competitivas y en sus nichos, mientras que la coexistencia de los diferentes grupos sería posible por la diferenciación de nicho entre ellos. No obstante, se ha sugerido que existen nichos adicionales, no contemplados en los estudios de neutralidad emergente, que en realidad estabilizan la coexistencia al interior de los grupos. También se ha encontrado que se 2 pueden formar grupos como resultado de alianzas entre especies que se facilitan. Sin embargo, no se sabe si alguno de estos escenarios ocurre en las comunidades vegetales. En el Capítulo II investigamos cómo cambian la diferenciación de nicho y los filtros ambientales a lo largo de un gradiente hídrico. Utilizamos un experimento en campo donde se sembraron seis especies de hierbas a lo largo de un gradiente hídrico en un pastizal semiárido bajo tres tratamientos: sin interacciones, con interacciones intraespecíficas y con interacciones interespecíficas. Analizamos cómo cambió la proporción de competencia inter- e intraespecífica y, por lo tanto, la diferenciación de nicho a lo largo del gradiente. Todas las especies, excepto una, mostraron diferenciación de nicho en algún punto del gradiente y la diferenciación de nicho fue ligeramente mayor con estrés. En el Capítulo III investigamos si las especies de una comunidad vegetal forman grupos a lo largo del gradiente hídrico como resultado de la neutralidad emergente, la presencia de nichos ocultos o alianzas de facilitación. Utilizamos datos observacionales de 35 especies de hierbas para parametrizar modelos poblacionales, que usamos para simular a la comunidad en diferentes escenarios de interacciones. Encontramos que la competencia es más débil al interior de los grupos, evidencia de que otros ejes del nicho podrían contribuir a estabilizar la coexistencia al interior de los grupos. Asimismo, la facilitación indirecta fomentó la formación de grupos pues los miembros de los grupos incrementan la abundancia de otros. En el Capítulo IV se discuten las aportaciones de esta tesis al entendimiento de la coexistencia en comunidades vegetales. Nuestro trabajo se suma a estudios previos que han encontrado una gran diferenciación de nicho promoviendo la coexistencia. Queda pendiente definir cuáles son los ejes del nicho que permiten la coexistencia entre plantas. 3 Abstract Understanding how the diversity of species we observe in nature coexists has been one of the main challenges in ecology. Modern coexistence theory states that species may coexist stably in a community through a balance between stabilizing mechanisms and equalizing mechanisms. When stabilizing mechanisms come into play, species may coexist if they differ in their niches, and they do not interfere with other species’ resource acquisition. If this occurs, intraspecific competition is stronger than interspecific competition, and each species will regulate its population more strongly than it regulates others. Equalizing mechanisms refer to each species’ fitness. Thus, if there are fitness differences among species, they will not be able to coexist unless niche differences are large enough to overcome such fitness differences. Environmental conditions, such as stress, may modify intra- and interspecific interactions. However, we have limited information on how these environmental conditions may affect coexistence and particularly stabilizing mechanisms. It has been hypothesized that the simultaneous presence of equalizing and stabilizing mechanisms may generate a pattern known as emergent neutrality in communities, where clusters of similar species emerge. When this is the case, species within clusters would be able to coexist because of similarities in both their niches and their fitness. Species in different groups would be able to coexist through niche differences. Alternatively, it has been suggested that additional niche dimensions not contemplated in studies of emergent neutrality (referred to as hidden niches) could contribute to stabilizing coexistence within clusters. Furthermore, clusters have also been observed when facilitation alliances occur between species in plant communities. 4 In Chapter II we analyzed how niche differences and environmental filtering change along a hydric stress gradient. To this end, we sowed six species of annual herbs along a hydric gradient in a semiarid grassland. Species were sowed according to three treatments: without interactions, with intraspecific interactions, and with interspecific interactions. We then analyzed how the ratio of intra- and interspecific interactions and, in consequence niche differences, changed along the hydric stress gradient. All but one species had niche differentiation in at least some portion of the gradient and we found slightly greater niche differentiation with greater stress. In Chapter III we tested if the species of a plant community clustered along a hydric stress gradient because of species interactions. We then assessed if these clusters resulted from emergent neutrality, hidden niches, or facilitation alliances. We used data from 35 herb species to parameterize population models. We then used these models to simulate the community over the hydric gradient in different interaction scenarios. We found that species cluster along the gradient and that competition is weaker within clusters, evidence that additional hidden niches stabilize coexistence within clusters. Furthermore, indirect facilitation may be important for clusters formation, where members increase each other’s abundance when they occur together. In Chapter IV we discuss our findings and contributions to the understanding of species coexistence in plant communities. This dissertation adds to previous studies that have also found great niche differentiation between plant species that promote coexistence. Future directions include studying what niche axes are allowing for stabilizing niche differences. 5 Capítulo I. Introducción General En la naturaleza frecuentemente observamos un gran número de especies coexistiendo en un determinado sitio, en un cierto hábitat o una comunidad. Una pregunta importante en la ecología ha sido el porqué de esta gran riqueza de especies que viven juntas, ya que por lo general los individuos requieren de un número limitado de recursos. Por ejemplo, las plantas compiten por luz, agua y algunos cuantos nutrientes como fósforo, magnesio y nitrógeno (Barot y Gignoux 2004; Silvertown 2004). Por lo tanto, dada la gran cantidad de especies que frecuentemente encontramos viviendo en una misma comunidad, es difícil imaginarse cómo podrían coexistir los organismos si solo tomamos en cuenta a los recursos que utilizan: unas cuantas especies podrían acaparar la totalidad de los recursos y excluyendo competitivamente al resto. Una posibilidad es que la heterogeneidad de condiciones ambientales permita la diferenciación en la adquisición de recursos. Sin embargo, la importancia del ambiente en la coexistencia de especies a la luz de la teoría moderna de la coexistencia ha sido poco explorada hasta el momento. Incorporar el efecto del ambiente al estudio de la coexistencia de especies nos podría ayudar a entender cómo es que se mantiene la gran riqueza de especies que observamos en la naturaleza. Teoría moderna de la coexistencia En décadas previas, la ecología se debatió en torno a dos teorías para tratar de explicar la coexistencia de especies en las comunidades (Chesson 1991). Por un lado, el principio de exclusión competitiva establecieron que dos especies deberían tener nichos diferentes para coexistir (den Boer 1986) mientras que la teoría neutral sugería que dos especies idénticas 6 podrían coexistir por largos periodos de tiempo porque ninguna excluye competitivamente a la otra (Hubbell 2001). La teoría moderna de la coexistencia reconcilia estas dos teorías al reconocer que ambos componentes pueden operar al mismo tiempo (Chesson 2000; Adler, HilleRisLambers y Levine 2007). Según la teoría moderna de la coexistencia, la coexistencia estable de las especies dentro de las comunidades está determinada por los mecanismos estabilizadores y los mecanismos igualadores (Chesson 2000; Adler, HilleRisLambers y Levine 2007). Los mecanismos estabilizadores actúan cuando la competencia intraespecífica es más fuerte que la competencia interespecífica, lo que implica que cada especie se regula a sí misma más fuertemente de lo que regula a las demás (Broekman et al. 2019). Esto sucede cuando las especies difieren en su nicho de manera que no interfieren en la adquisición de recursos de otras especies (Chesson 2000; Broekman et al. 2019). Cuando hay diferenciación de nicho y las especies se vuelven raras en una comunidad, sus tamaños poblacionales podrán incrementar ya que sus individuos conespecíficos serán menos abundantes, aliviando la competencia. Es decir, la diferenciación de nicho es un mecanismo denso dependiente. En el contexto de la teoría moderna de la coexistencia el nicho se define como el efecto que tienen una especie sobre los recursos, por lo que repercute en las interacciones intra e interespecíficas. Por lo tanto, podríamos decir que se asemeja al nicho Eltoniano (Leibold 1995; Chesson 2000, Soberón 2007). La diferenciación de nicho entre las especies surge cuando hay diferencias en los recursos que más limitan su crecimiento, cuando hay diferencias en sus interacciones con depredadores especialistas o cuando hay diferencias espaciotemporales en los ambientes que más favorecen su desempeño (Turcotte y Levine 2016). 7 Sin embargo, la coexistencia no depende únicamente de los mecanismos estabilizadores y la diferenciación de nicho. Depende también de los mecanismos igualadores. Los mecanismos igualadores se refieren a las asimetrías en las habilidades competitivas de las especies pero que son independientes de la densidad (Chesson 2000). En la teoría moderna de la coexistencia, si hay asimetrías entre las habilidades competitivas de las especies, la coexistencia no es posible a menos de que las diferencias de nicho sean suficientemente grandes como para superar las diferencias en las capacidades competitivas (Chesson 2000). Algunos ejemplos de las diferencias en las capacidades competitivas de las especies son las diferencias en las habilidades para agotar un recurso compartido, las diferencias en la capacidad reproductivas cuando los recursos no están limitados o las diferencias en la tolerancia a un consumidor generalista (Turcotte y Levine 2016). La exclusión competitiva y la teoría neutral representan casos especiales extremos dentro de la teoría moderna de la coexistencia. En la exclusión competitiva se considera que sólo están operando los mecanismos estabilizadores y las especies podrían coexistir establemente siempre y cuando la competencia intraespecífica sea mayor que la competencia interespecífica, es decir, cuando ocupen nichos distintos (den Boer 1986). En la teoría neutral, las especies podrían coexistir en una comunidad cuando no haya diferencias en sus habilidades competitivas ni en sus nichos (Hubbell 2001). Implicaciones de las condiciones ambientales La importancia de los mecanismos estabilizadores y los mecanismos igualadores, así como de la competencia y la facilitación en las comunidades, dependen de las condiciones ambientales (Brooker et al. 2008). En particular, las diferencias en la humedad del suelo 8 pueden contribuir a la diferenciación del nicho de las plantas y pueden ser un filtro ambiental importante especialmente en condiciones de aridez (Silvertown, Araya y Gowing 2015). La intensidad y el signo de las interacciones dependen de las condiciones ambientales, especialmente del estrés (Brooker et al. 2008; Callaway 2007). Por ende, es razonable suponer que las condiciones ambientales como la disponibilidad de agua, tengan un efecto sobre la importancia relativa de los mecanismos de coexistencia y los filtros ambientales. El estrés se define como los factores externos a las plantas, como escasez o exceso de recurso o también condiciones abióticas extremas, que limitan su capacidad de producción (Grime 1977). Esta disminución en la capacidad de producción puede ocasionar que la competencia disminuya en condiciones estresantes (Grime 1977). Esto concuerda con la hipótesis del gradiente de estrés, la cual establece que en los ambientes benignos predomina la competencia y en los estresantes la facilitación (Bertness y Callaway, 1994; Brooker et al. 2008; Callaway 2007). Sin embargo, cuando el estrés es causado por la escasez de recursos tales como el agua, se ha observado el patrón contrario: la competencia por el recurso escaso se vuelve intensa en los sitios más estresantes donde escasea dicho recurso (Maestre et al. 2009; Tielbörger y Kadmon 2000; Maestre y Cortina 2004). Dada esta complejidad, es posible suponer que las condiciones estresantes deben tener un efecto sobre la importancia relativa de los mecanismos de coexistencia. En particular, el estrés podría determinar la intensidad de los mecanismos estabilizadores ya que éstos dependen de la intensidad de la competencia intraespecífica respecto a la interespecífica. 9 El ambiente también puede determinar la presencia de las especies en las comunidades. Es decir, el ambiente puede funcionar como un filtro para el establecimiento de las especies. Un filtro ambiental ocurre cuando las condiciones ambientales son demasiado severas (o los recursos muy escasos) y la mayoría de las especies no se pueden establecer en la comunidad (Krafft et al. 2015). En particular, la humedad del suelo puede ser importante no solo para las interacciones, sino que también puede ser un importante filtro ambiental (Silvertown et al. 2015). En el contexto de la teoría moderna de la coexistencia podríamos esperar que las diferencias de nicho no sean muy grandes en las comunidades de ambientes benignos. Un ambiente benigno es aquel en donde la mayoría de las especies pueden tolerar las condiciones ambientales (Chase 2007); por lo tanto, puede haber una gran riqueza de especies en estos ambientes. Para las comunidades muy diversas se ha hipotetizado que un aumento en el número de especies a lo largo de un gradiente ambiental finito podría implicar un mayor sobrelapamiento de nichos, dando como resultado especies ecológicamente equivalentes (Gravel et al. 2006; Chisholm y Pacala 2011). Además, podríamos esperar que en los vecindarios de individuos en las comunidades muy diversas las combinaciones de especies que se presentan sean una fracción muy pequeña de todas las posibles combinaciones de vecindarios de individuos que podrían existir (Lieberman y Lieberman 2007; Hubbell y Foster 1986; Wang et al. 2016). Por lo tanto, la frecuencia con la que se presenta una determinada combinación de especies alrededor de un individuos tomado al azar en una comunidad muy diversa es muy baja y en consecuencia no evolucionan estrategias competitivas específicas para lidiar con cada una de las especies (Hubbell y Foster 1986). Por lo tanto, en comunidades muy ricas en especies no se espera 10 un resultado determinístico de la competencia ni una gran diferenciación de nicho, por lo que las especies tenderían a la equivalencia ecológica (Hubbell y Foster 1986). Sin embargo, los trabajos que han analizado hasta el momento, el efecto del estrés sobre los mecanismos de coexistencia y los filtros ambientales no han encontrado un patrón consistente (Wainwright et al. 2018; Bimler et al. 2018; Germain, Mayfield y Gilbert 2018). En contraste, en ambientes estresantes como el desierto (Cody 1991) y la alta montaña (Halloy 1990), la diversidad de formas de crecimiento (Rowe y Speck, 2005) es de las más altas del planeta. Esto sugiere que existe una gran diferenciación de nicho entre las especies de los ambientes adversos. Sin embargo, si nos movemos al extremo del gradiente de ambientes sumamente estresantes, es probable que dominen los filtros ambientales y encontraríamos que sólo coexisten las especies capaces de establecerse en dichos ambientes (Chase 2007). Además, en el extremo del gradiente de estrés la competencia puede perder importancia (Grime 1977), y por lo tanto podrían también perder importancia los mecanismos estabilizadores que se rigen por el balance entre las competencias intra- e interespecífica (Broekman et al. 2019). Agrupamiento de especies en las comunidades La presencia simultánea de los mecanismos estabilizadores e igualadores en una comunidad puede generar un patrón denominado neutralidad emergente. Este patrón se ha observado en comunidades ricas en especies donde se generan grupos discretos de especies semejantes (Scheffer y van Nes 2006). La neutralidad emergente se observó por primera vez como resultado de simulaciones de comunidades usando el modelo de Lotka-Volterra 11 (Scheffer y Nes 2006). Al inicio de la simulación, las especies se colocaron aleatoriamente sobre el gradiente de un eje de nicho hipotético resultando en una competencia interespecífica proporcional a la cercanía entre ellas en el gradiente. En este escenario, se esperaba que las especies que estuvieran demasiado cerca entre sí se extinguieran y que al final de la simulación permanecieran solamente especies que fueran suficientemente distintas en sus nichos como para coexistir. Para sorpresa de los autores, se observó que surgieron grupos de especies similares a lo largo del gradiente del eje del nicho, lo que los llevó a concluir que las especies podían coexistir al interior de sus grupos por similitudes en sus habilidades competitivas y similitudes en sus nichos sensu Chesson. Es decir, los grupos se compartan como una comunidad neutral (Hubbel, 2001) donde no hay exclusión competitiva gracias a que las especies son equivalentes. Las especies que forman parte de distintos grupos pueden coexistir por diferenciación de nicho. El agrupamiento de especies se ha observado frecuentemente como resultado de la competencia en simulaciones computacionales (Doebeli y Dieckmann 2003; Bonsall 2004; Scheffer y van Nes 2006; Roelke y Eldridge 2008; Scranton y Vasseur 2016). En la naturaleza, este fenómeno se ha estudiado en menor grado pero se ha observado en comunidades de plancton (Vergnon, Dulvy, y Freckleton 2009; Segura et al. 2011), en comunidades de escarabajos (Scheffer et al. 2015) y en aves (Thibault et al. 2011). Sin embargo, las interacciones de las especies dentro y entre los grupos han sido poco estudiades (pero véase Segura et al. 2011). Para plantas solamente hay un trabajo donde se estudia el agrupamiento de las especies a lo largo de diversos gradientes que podrían ser ejes importantes del nicho (D’Andrea et al. 2020). El estudio se realizó en una selva tropical húmeda, donde se encontró que las especies se agrupan de acuerdo con los caracteres que 12 están asociados a la competencia por luz. Sin embargo, los autores señalan que las especies al interior de los grupos no necesariamente son similares en sus requerimientos lumínicos. Por lo tanto, los grupos podrían no ser resultado de la neutralidad emergente. Sin embargo, no se estudiaron las interacciones de las especies de esta comunidad ni los mecanismos de coexistencia. El agrupamiento de especies no sólo se puede explicar por la neutralidad emergente, ya que hay otros procesos que pueden generar agrupamiento en las comunidades. De hecho, existe una controversia respecto a la neutralidad de los grupos que se observan en las comunidades (Barabás et al. 2013; Vergnon, van Nes, y Scheffer 2013). Barabás et al. (2013) discuten que, en realidad, para que los grupos puedan permanecer en las comunidades se requiere de diferenciación de nicho dentro de los grupos y no solo entre los grupos. Estos autores argumentan también que los estudios que han detectado la presencia de neutralidad emergente pueden no estar considerando otros ejes del nicho, que pueden ser los responsables de estabilizar la coexistencia entre las especies dentro de un mismo grupo. A estos nichos no estudiados, que pueden estabilizar la coexistencia dentro de los grupos, se les llamó nichos ocultos. Otro proceso que puede generar agrupaciones de especies dentro de las comunidades puede ser la formación de alianzas entre especies a través de la facilitación indirecta. La facilitación indirecta puede ocurrir cuando una especie suprime a un competidor fuerte, aliviando así la presión competitiva sobre una tercera especie (Levine 1999; Davidson 1985). Es más probable que este escenario se presente cuando la primera y la tercera especie requieran de recursos diferentes. La facilitación indirecta entre plantas se ha observado en comunidades herbáceas (Levine 1999; Callaway y Pennings 2000). 13 También se ha observado que la facilitación indirecta sobre más de una especie puede generar agrupamientos en las comunidades (Wilson 1983). Este es el caso de las alianzas entre artrópodos carroñeros, en las que diferentes especies de escarabajos facilitan indirectamente a los ácaros carroñeros al alimentarse de las larvas de las moscas que compiten con los ácaros (Wilson 1983; Wilson y Knollenberg 1987). Un fenómeno similar se observó en una comunidad de erizos acuáticos, en la que una especie facilita indirectamente a otras tres al disuadir a los depredadores de erizos (Duggins 1981). Objetivos y preguntas de investigación El ambiente abiótico tiene un efecto sobre las interacciones entre las plantas. En particular, el estrés puede cambiar la intensidad y el signo de las interacciones (Brooker et al. 2008; Callaway 2007). Sin embargo, no conocemos las consecuencias que pueden tener estos cambios sobre la coexistencia de especies y en específico sobre las diferencias de nicho que pueden estabilizar la coexistencia en comunidades vegetales. Además, el ambiente abiótico también puede funcionar como un filtro, impidiendo el establecimiento de ciertas especies que no pueden tolerar condiciones ambientales extremas. Por lo tanto, entender el efecto del ambiente sobre la diferenciación de nichos y el establecimiento de las especies en las comunidades puede ayudarnos a entender mejor el mantenimiento de la diversidad de las comunidades. Por otro lado, a pesar de que los agrupamientos de especies a lo largo de los gradientes ambientales parecen ser un fenómeno relativamente común en la naturaleza, existe una controversia en cuanto a los mecanismos que los producen. De hecho, la gran mayoría de los estudios que encuentran grupos utilizan simulaciones 14 matemáticas puramente teóricas, por lo que, para entender mejor los efectos del ambiente sobre la coexistencia de especies y la estructura de las comunidades hacen falta estudios empíricos que midan las interacciones entre las especies que coexisten y el efecto del ambiente sobre las mismas. El objetivo general de esta tesis fue determinar el efecto de un gradiente hídrico sobre la diferenciación de nicho entre las plantas de en un pastizal semiárido. Para alcanzar este objetivo, utilizamos datos experimentales para evaluar el efecto del gradiente hídrico sobre las interacciones intra- e interespecíficas, así como las consecuencias de los cambios en estas sobre la diferenciación de nichos entre las especies. También utilizamos estos datos experimentales para determinar si el estrés hídrico representaba un filtro ambiental en la comunidad. Por otro lado, utilizamos datos observaciones para examinar si las especies se organizaban en grupos a lo largo del gradiente de estrés hídrico como resultado de las interacciones. Con estos datos parametrizamos modelos poblacionales que utilizamos para simular a la comunidad con y sin interacciones, a fin de establecer los procesos que generan estos grupos: neutralidad emergente, nichos ocultos o alianzas de facilitación. Debido a la importancia de la disponibilidad de agua para la estructuración de las comunidades vegetales, en esta tesis se hipotetizó que la profundidad de suelo, relacionada con la disponibilidad de agua, permite la diferenciación de nichos entre las plantas de la comunidad. 15 Objetivo general: Analizar el efecto de un gradiente de estrés hídrico sobre la coexistencia de especies vegetales, evaluando en particular la diferenciación de nicho entre las plantas de un pastizal semiárido. Objetivos particulares: 1) Evaluar experimentalmente la importancia de los mecanismos estabilizadores y de los filtros ambientales a lo largo de un gradiente de estrés hídrico. 2) Analizar si existen grupos de especies a lo largo de un gradiente hídrico como resultado de una de las tres hipótesis que predicen agrupamiento de especies en las comunidades: neutralidad emergente, nichos ocultos o alianzas entre especies. Sitio de estudio El sitio de estudio es el pastizal semiárido de Concepción Buenavista, Oaxaca, México. En este sistema la diversidad vegetal es muy alta. Hay más de 600 especies de plantas vasculares en este pastizal (Martorell et al. 2022) y de hecho comparte el récord mundial con las praderas de Estonia que tienen hasta 25 especies en 0.1 ´ 0.1 m (Martorell et al. 2017). Esta característica lo hace un sistema ideal para abordar el problema de la coexistencia de especies. La profundidad del suelo varía de 1 a 20 cm y la disponibilidad de agua aumenta con la profundidad (Villarreal-Barajas y Martorell 2009; Martorell y Martínez-López 2014). Estudios previos han determinado que a lo largo de este gradiente hay un recambio notable de especies (Martorell et al. 2015). 16 Estructura de la tesis Esta tesis está conformada por cuatro capítulos. El Capítulo I y el Capítulo IV corresponden, respectivamente, a la Introducción General y a la Discusión General, y constituyen al marco conceptual del trabajo donde se sintetiza el alcance de la investigación. A continuación, se detallan los objetivos de investigación abordados en los Capítulo II y III. Capítulo II - Changes in niche differentiation and environmental filtering over a hydric stress gradient. A. Martínez-Blancas y C. Martorell. 2020. Journal of Plant Ecology. 13: 185–194. DOI: 10.1093/jpe/rtz061 En este Capítulo se planteó la hipótesis particular de que las diferencias de nicho y los filtros ambientales se fortalecen conforme el estrés ambiental se agudiza (Figura 1.). Para analizar esto, sembramos individuos de seis especies de hierbas anuales en el campo a lo largo de un gradiente de estrés hídrico en tres tratamientos distintos: sin interacciones (sin otras plantas), con interacciones intraespecífica (en presencia de individuos conespecíficos) y con interacciones interespecíficas (en presencia de individuos heteroespecíficos). Al final del experimento evaluamos la biomasa final y registramos el número de inflorescencias producidas para evaluar la intensidad de las interacciones. Para interpretar los resultados, planteamos que si la competencia intraespecífica era más fuerte 17 que la interespecífica, esto significaría que hay evidencia de la diferenciación de los nichos de las especies. Esperábamos que este escenario se presentara predominantemente en el extremo estresante del gradiente. Por el contrario, si el estrés ambiental limitaba el establecimiento de las plantas, esto significaría que el estrés hídrico está funcionando como un filtro ambiental (Figura 1. panel c) Figura 1. Hipótesis puestas a prueba en el Capítulo II. En el panel a) se muestra cómo esperaríamos que cambiara la competencia intraespecífica e interespecífica con el estrés hídrico. La línea azul representa los tratamientos donde las especies se sembraron con individuos conespecíficos mientras que la línea verde representa los tratamientos donde las especies se sembraron con conespecíficos. Conforme aumenta el estrés hídrico, la competencia intraespecífica incrementa con el estrés más de lo que aumenta la competencia interespecífica. Esto resulta en lo que se observa en el panel b) donde la diferenciación de nicho es mayor conforme aumenta el estrés. El panel c) muestra un Estrés hídrico C o m p e te n c ia Competencia interespecífica D if e re n c ia c ió n d e n ic h o Estrés hídrico F il tr o a m b ie n ta l Competencia intraespecífica a) b) c) 18 incremento en el filtro ambiental conforme aumenta el estrés hídrico medido en ausencia de interacciones (tratamientos de un solo individuo). Capítulo III - Species alliances and hidden niche dimensions drive species clustering along a hydric gradient in a semiarid grassland. A. Martínez-Blancas, I. X. Belaustegui y C. Martorell. 2022. Ecology Letters. 25:2651- 2662. DOI: 10.1111/ele.14122 La hipótesis particular que se puso a prueba en este capítulo establece que las especies de la comunidad se agrupan en ciertas porciones del gradiente hídrico como resultado de las interacciones y como consecuencia de neutralidad emergente, nichos ocultos o alianzas de facilitación. Para analizar el efecto de las interacciones sobre la estructura de la comunidad, parametrizamos modelos de crecimiento poblacional para 35 especies con datos de campo para estimar la competencia y la facilitación entre pares de especies. Después usamos estos modelos para simular las dinámicas de la comunidad en tres escenarios distintos de interacciones y de miembros de la comunidad. Para determinar si los grupos surgen a partir de las interacciones, comparamos la formación de grupos en los escenarios con y sin interacciones. Para establecer si los grupos eran neutrales o si involucraban nichos ocultos, comparamos la intensidad de la competencia entre los grupos con la competencia al interior de los grupos. Esperábamos que la competencia entre grupos fuera más fuerte en un escenario de neutralidad emergente (Figura 2. panel b) pero lo opuesto en un escenario de nichos ocultos (Figura 2. panel c). 19 Para determinar si los grupos surgen por alianza de facilitación, analizamos si los grupos desaparecían al eliminar la facilitación del modelo y si la facilitación era más fuerte al interior de los grupos (Figura 2. panel d). Para probar si surgen alianzas de facilitación como resultado de interacciones indirectas, analizamos si las especies se veían favorecidas por la presencia en conjunto de las demás especies en su grupo. Figura 2. Los patrones esperados bajo los diferentes mecanismos que pueden propiciar la formación de grupos en las comunidades. Cada línea representa una especie. Si no hay grupos de especies en las comunidades se esperaríamos encontrar lo que se observa en el panel a). Bajo neutralidad emergente (panel b) hay un traslape completo de nichos y por lo tanto la competencia es fuerte al interior de los grupos y débil entre grupos. La coexistencia es posible debido a que las habilidades competitivas son iguales. En un escenario de nichos ocultos (panel c) la competencia es débil dentro de los grupos porque hay diferenciación de nicho en ejes del nicho no observados en el estudio. Las especies que tienen diferenciación de nicho en estos están representadas con colores distintos. Competencia débil Competencia débil Competencia fuerte Competencia fuerte Competencia fuerte Competencia fuerte Competencia fuerte Competencia débil Competencia débil Competencia débil Facilitación directa Facilitación directa Facilitación directta No hay grupos Neutralidad emergente Nichos ocultos Alianza de especies Eliminación de interacciones Eliminación de interacciones Elim inación de facilitación Elim inación de interacciones A b u n d a n ci a Disponibilidad de agua (profundidad de suelo) A b u n d a n ci a Disponibilidad de agua (profundidad de suelo) a) b) c) d) Competencia fuerte Competencia fuerte Competencia débil Competencia débil Competencia débil 20 Bajo la hipótesis de las alianzas de especies (panel d) las especies se agrupan porque se facilitan unas a otras directamente (morado) o indirectamente al suprimir a los competidores fuertes que tienen en común (negro). En todos los casos, la eliminación de las interacciones debilita la formación de grupos de especies (panel a). Si hay alianzas de especies a través de facilitación directa, la eliminación de la facilitación sería suficiente para ocasionar dicho debilitamiento. Por otro lado, si eliminamos a todas menos una de las especies de cada grupo, observaríamos que la abundancia de la especie incrementaría fuertemente en un escenario de neutralidad emergente, débilmente en un escenario de nichos ocultos y decrecería si hay alianzas de especies. Referencias Adler PB, HilleRislambers J, Levine JM (2007) A niche for neutrality. Ecol Lett 10:95– 104. Anderson HM, Hutson V, Law R (1992) On the conditions for permanence of species in ecological communities. Am Nat 139:663–668. Barabás G, D’Andrea R, Rael R, Meszéna G, Ostling A (2013) Emergent neutrality or hidden niches? Oikos 122:1565–1572. Barot S, Gignoux J (2004) Mechanisms promoting plant coexistence: Can all the proposed processes be reconciled? Oikos 106:185–192. Bertness MD, Callaway R (1994) Positive interactions in communities. Trends Ecol Evol 9:191–193. Bimler MD, Stouffer DB, Lai HR, Mayfield MM (2018) Accurate predictions of coexistence in natural systems require the inclusion of facilitative interactions and 21 environmental dependency. J Ecol 106:1839–1852. Bonsall MB (2004) Life History Trade-Offs Assemble Ecological Guilds. Science (80- ) 306:111–114. Broekman MJE, Muller-Landau HC, Visser MD, Jongejans E, Wright SJ, Kroon H de (2019) Signs of stabilisation and stable coexistence. Ecol Lett 22:1957–1975. Brooker RW, Maestre FT, Callaway RM, et al. (2008) Facilitation in plant communities: the past, the present, and the future. J Ecol 96:18–34. Callaway RM (2007) Positive Interactions and Interdependence in Plant Communities. 1st ed. Dordrecht, Netherlands: Springer. Callaway RM, Pennings SC (2000) Facilitation may buffer competitive effects: Indirect and diffuse interactions among salt marsh plants. Am Nat 156:416–424. Chase JM (2007) Drought mediates the importance of stochastic community assembly. Proc Natl Acad Sci U S A 104:17430–17434. Chesson P (1991) A need for niches?. Trends Ecol Evol 1:26–8. Chesson P (2000) Mechanisms of maintenance of species diversity. Annu Rev Ecol Syst 31:343–366. Chisholm RA, Pacala SW (2011) Theory predicts a rapid transition from niche-structured to neutral biodiversity patterns across a speciation-rate gradient. Theor Ecol 4:195– 200. Cody ML (1991) Niche theory and plant growth form. Vegetatio 97:39–55. D’Andrea R, Guittar J, O’Dwyer JP, et al. (2020) Counting niches: Abundance-by-trait patterns reveal niche partitioning in a Neotropical forest. Ecology 101:e03019. Davidson DW (1985) An Experimental Study of Diffuse Competition in Harvester Ants. 22 Am Nat 125:500–506. den Boer PJ (1986) The present status of the competitive exclusion principle. Trends Ecol Evol 1:25–8. Doebeli M, Dieckmann U (2003) Speciation along environmental gradients. Nature 421:259–264. Duggins DO (1981) Interspecific facilitation in a guild of benthic marine herbivores. Oecologia 48:157–163. Germain RM, Mayfield MM, Gilbert B (2018) The ‘filtering’ metaphor revisited: competition and environment jointly structure invasibility and coexistence. Biol Lett 14:20180460. Gravel D, Canham CD, Beaudet M, Messier C (2006) Reconciling niche and neutrality: the continuum hypothesis. Ecol Lett 9:399–409. Grime JP (1977) Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory. Am Nat 111:1169–1194. Halloy S (1990) A morphological classification of plants, with special reference to the New Zealand alpine flora. J Veg Sci 1:291–304. Hubbell SP (2001) The Unified Neutral Theory of Biodiversity and Biogeography. Princeton, New Jersey: Princeton University Press. Hubbell SP, Foster RB (1986) Biology, chance, and history and the structrue of tropical rain forest tree communities.In Diamond J, Case TJ (eds). Community Ecology. New York: Harper & Row, 314–330. Leibold MA (1995) The niche concept revisited: mechanistic models and community context. Ecology 76:1371–82. 23 Levine JM (1999) Indirect facilitation: Evidence and predictions from a riparian community. Ecology 80:1762–1769. Lieberman M, Lieberman D (2007) Nearest-neighbor tree species combinations in tropical forest: the role of chance, and some consequences of high diversity. Oikos 116:377–386. Maestre FT, Callaway RM, Valladares F, Lortie CJ (2009) Refining the stress-gradient hypothesis for competition and facilitation in plant communities. J Ecol 97:199–205. Maestre FT, Cortina J (2004) Do positive interactions increase with abiotic stress? A test from a semi-arid steppe. Proc R Soc B 271:S331–S333. Martorell C, Almanza-Celis CAI, Pérez-García EA, Sánchez-Ken JG, Perez-Garcia EA, Sanchez-Ken JG (2015) Co-existence in a species-rich grassland: Competition, facilitation and niche structure over a soil depth gradient. J Veg Sci 26:674–685. Martorell, C., García-Meza, D., Martínez-Blancas, A., Zepeda, V. & Vazquez-Ribera, C. (2022) Pastizales de la región Chocholteca: un despreciado récord mundial de diversidad vegetal. En: Cruz Angón, A., Nájero Cordero, K.C., Cruz Medina, J. & Melgarejo, E.D. (Eds.) La biodiversidad en Oaxaca estudio de estado. Ciudad de México: Comisión Nacional para el Conocimiento y Uso de la Biodiversidad, pp. 290–301. Martorell C, Martínez-López M (2014) Informed dispersal in plants: Heterosperma pinnatum (Asteraceae) adjusts its dispersal mode to escape from competition and water stress. Oikos 123:225–231. Martorell C, Zepeda V, Martínez-Blancas A, García-Meza D, Pedraza F (2017) A diversity world record in a grassland at Oaxaca, México. Bot Sci 95:1–7. 24 Roelke DL, Eldridge PM (2008) Appendix Mixing of Supersaturated Assemblages and the Precipitous Loss of Species. Am Nat 171:162–175. Rowe N, Speck T (2005) Plant growth forms: an ecological and evolutionary perspective. New Phytol 166:61–72. Scheffer M, Nes EH Van (2006) Species of groups of similar species. Proc Natl Acad Sci 103:6230–6235. Scheffer M, Nes EH van (2006) Self-organized similarity, the evolutionary emergence of groups of similar species. Proc Natl Acad Sci 103:6230–6235. Scheffer M, Vergnon R, Nes EH Van, et al. (2015) The Evolution of Functionally Redundant Species; Evidence from Beetles. PLoS One 10:e0137974. Scranton K, Vasseur DA (2016) Coexistence and emergent neutrality generate synchrony among competitors in fluctuating environments. Theor Ecol 9:353–363 10.1007/s12080-016-0294-z. Segura AM, Calliari D, Kruk C, Conde D, Bonilla S, Fort H (2011) Emergent neutrality drives phytoplankton species coexistence. Proc R Soc B Biol Sci 278:2355–2361 10.1098/rspb.2010.2464. Silvertown J (2004) Plant coexistence and the niche. Trends Ecol Evol 19:605–611. Silvertown J, Araya Y, Gowing D (2015) Hydrological niches in terrestrial plant communities: a review. J Ecol 103:93–108. Soberón J (2007) Grinnellian and Eltonian niches and geographic distributions of species. Ecol Lett 10:1115–1123. Thibault KM, White EP, Hurlbert AH, Ernest SKM (2011) Multimodality in the individual size distributions of bird communities. Glob Ecol Biogeogr 20:145–153. 25 Tielbörger K, Kadmon R (2000) Temporal environmental variation tips the balance between facilitation and interference in desert plants. Ecology 81:1544–1553. Turcotte MM, Levine JM (2016) Phenotypic plasticity and species coexistence. Trends Ecol Evol 31:803–813. Vergnon R, Dulvy NK, Freckleton RP (2009) Niches versus neutrality: uncovering the drivers of diversity in a species-rich community. Ecol Lett 12:1079–1090. Vergnon R, Nes EH van, Scheffer M (2013) Interpretation and predictions of the Emergent neutrality model: a reply to Barabás et al. Oikos 122:1573–1575. Villarreal-Barajas T, Martorell C (2009) Species-specific disturbance tolerance, competition and positive interactions along an anthropogenic disturbance gradient. J Veg Sci 20:1027–1040. Wainwright CE, HilleRisLambers J, Lai HR, Loy X, Mayfield MM (2018) Distinct responses of niche and fitness differences to water availability underlie variable coexistence outcomes in semi-arid annual plant communities. J Ecol 107:293–306. Wang X, Wiegand T, Kraft NJB, et al. (2016) Stochastic dilution effects weaken deterministic effects of niche-based processes in species rich forests. Ecology 97:347–360. Wilson DS (1983) The Effect of Population Structure on the Evolution of Mutualism: A Field Test Involving Burying Beetles and Their Phoretic Mites. Am Nat 121:851– 870. Wilson DS, Knollenberg WG (1987) Adaptive indirect effects: the fitness of burying beetles with and without their phoretic mites. Evol Ecol 1:139–159. 26 Capítulo II Changes in niche differentiation and environmental filtering over a hydric stress gradient A. Martínez-Blancas & C. Martorell. 2020. Journal of Plant Ecology 13: 185 – 194. Journal of Plant Ecology JOURNAL OF PLANT ECOLOGY | doi:10.1093/jpe/rtz061 185 © The Author(s) 2020. Published by Oxford University Press on behalf of the Institute of Botany, Chinese Academy of Sciences and the Botanical Society of China. All rights reserved. For permissions, please email: journals.permissions@oup.com Changes in niche differentiation and environmental filtering over a hydric stress gradient Alejandra Martínez-Blancas1,2,†, and Carlos Martorell1,*,†, 1Departamento de Ecología y Recursos Naturales, Facultad de Ciencias, Universidad Nacional Autónoma de México, Coyoacán, 04510 Ciudad de México, Mexico, 2Posgrado en Ciencias Biológicas, Universidad Nacional Autónoma de México, Av. Ciudad Universitaria 3000, C.P. 04510, Coyoacán, Ciudad de México, Mexico *Corresponding author. E-mail: martorell@ciencias.unam.mx †These authors contributed equally to this work. Received: 9 May 2019, Revised: 13 December 2019, Accepted: 27 December 2019, Advanced Access publication: 10 January 2020 Handling Editor: Shaopeng Wang Citation: Martínez-Blancas A, Martorell C (2020) Changes in niche differentiation and environmental filtering over a hydric stress gradient. J Plant Ecol 13:185–194. https://doi.org/10.1093/jpe/rtz061 Abstract Aims Diversity in communities is determined by species’ ability to coexist with each other and to overcome environmental stress that may act as an environmental filter. Niche differentiation (ND) results in stronger intra- than interspecific competition and promotes coexistence. Because stress affects interactions, the strength of ND may change along stress gradients. A greater diversity of plant growth forms has been observed in stressful habitats, such as deserts and alpine regions, suggesting greater ND when stress is strong. We tested the hypothesis that niche differences and environmental filters become stronger with stress. Methods In a semiarid grassland in southern Mexico, we sowed six annual species in the field along a hydric stress gradient. Plants were grown alone (without interactions), with conspecific neighbors (intraspecific interactions) or with heterospecific neighbors (interspecific interactions). We analyzed how the ratio of intra- to interspecific competition changed along the gradient to assess how water availability determines the strength of ND. We also determined if hydric stress represented an environmental filter. Important Findings We observed stronger intra- than interspecific competition, especially where hydric stress was greater. Thus, we found ND in at least some portion of the gradient for all but one species. Some species were hindered by stress, but others were favored by it perhaps because it eliminates soil pathogens. Although strong ND was slightly more frequent with stress, our species sample was small and there were exceptions to the general pattern, so further research is needed to establish if this is a widespread phenomenon in nature. Keywords: hydrological niche, species interactions, environmental constraints, species coexistence, stabilizing mechanisms, equalizing mechanisms 摘要：群落的多样性取决于物种相互共存的能力,以及克服环境胁迫的能力,从而可能扮演一个环境过滤器的角色。生态位分化导致了比种间 竞争更强的种内竞争,促进了物种共存。由于胁迫影响相互作用,生态位分化的强度可能沿胁迫作用的梯度发生变化。在有胁迫的栖息地,如沙 漠和高山地区,植物的生长形式更加多样化,这表明在环境胁迫下植物具有更强的生态位分化。我们验证了生态位差异和环境过滤作用随着环 境胁迫的增加而增强的假设。在墨西哥南部的半干旱草原上,我们沿着水分胁迫梯度在田间种植了6种一年生植物。植物进行三种种植处理： 单独种植(没有相互作用),与同种植物为邻(种内相互作用)或与异种植物为邻(种间相互作用)。我们分析了种内与种间竞争的比值是如何沿水 分胁迫梯度变化的,以评估水分如何决定生态位分化的强度。我们还测定了水分胁迫是否代表了一种环境过滤作用。我们发现了种内竞争比 种间竞争更强,尤其是在水分胁迫更大的情况下。因此,我们的结果表明,除了一个物种外,所有物种至少在部分水分胁迫梯度下存在生态位分 化。一些物种受到胁迫的阻碍,而另一些物种可能因为胁迫消除了土壤病原体而受到它的青睐。虽然高强度的生态位分化在环境胁迫下出现的 频率略高,但由于我们的物种样本很小,而且这种通用的模式也存在例外,因此需要进一步的研究来确定这是否是自然界普遍存在的现象。 关键词：水文生态位、物种相互作用、环境制约、物种共存、稳定机制、均衡机制 D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGY RESEARCH ARTICLE 186 JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 INTRODUCTION Species can only co-occur in a community for indefinitely long periods of time if there is a process that stabilizes their coexistence (Chesson 2000). Stabilizing mechanisms promote population growth when a species is rare, actively precluding the extinction that would follow if the population size kept diminishing (Chesson 2000). Niche differentiation (ND) can act as a stabilizing mechanism because, if species use different resources, competition is mostly intraspecific and thus diminishes with the species’ abundance. This allows the population to recover when rare, avoiding extinction (Anderson et al. 1992; Chesson 2000). Stronger intraspecific competition relative to interspecific competition is the signature of ND, and is a sufficient condition for stabilization to occur (Broekman et al. 2019; Levine and HilleRisLambers 2009). Other than interactions, the environment may also determine coexistence in a community. When environmental conditions are too severe, many species cannot get established (Kraft et al. 2015). Under this scenario, the environment acts as a filter. Spatial variations in the environment may allow the coexistence of species that would otherwise exclude each other competitively. This can happen, e.g. if one of the species is filtered-out from some areas, allowing its competitor to grow (Adler et al. 2013). The environment also affects biotic interactions and ND. Competition may decrease when stress limits plant’s ability to use and acquire resources (Grime 1977). Thus, competition is thought to lose importance with increasing stress (Grime 1977), even giving way to positive interactions, as stated in the stress gradient hypothesis (SGH; Brooker et al. 2008; Callaway 2007). However, it is still unclear how these changes in competition over stress gradients affect the relationship between intra- and interspecific competition that arises from ND (Hart and Marshall 2013). ND may not drive coexistence in benign (i.e. not stressful) environments. A benign environment has been defined as one where the majority of the species in the regional pool tolerate the physical conditions (Chase 2007), and thus would have high richness. It has been hypothesized that a large number of species combined with a finite variation in resources and conditions does not allow for large ND (Chisholm and Pacala 2011; Gravel et  al. 2006). In species-rich communities, the identity of neighbor plants is unpredictable (Hubbell and Foster 1986; Lieberman and Lieberman 2007; Wang et al. 2016). Thus, strategies that confer competitive superiority over particular species are unlikely to evolve, so species are expected to have equal competitive abilities (Hubbell and Foster 1986). In contrast, coexistence in stressful environments may be dominated by ND. For example, deserts (Cody 1991) and alpine regions (Halloy 1990) have a great diversity of plant growth forms, suggesting large ND. Diversity in water acquisition traits increases with hydric stress, also suggesting greater ND in stressful habitats (Bernard-Verdier et al. 2012). However, ND was found to diminish during an extreme drought event, suggesting that stabilization may be undermined under stressful conditions (Matías et al. 2018). Two other studies have found no trend in ND along stress gradients (Bimler et al. 2018; Wainwright et al. 2018), and Zambrano et al. (2017) report different relationships between ND (measured as trait differentiation) and stress in temperate and tropical forests. In this study, we tested the hypothesis that niche differences and environmental filters become stronger with stress. We conducted a field experiment using six herb species from a semiarid grassland grown over a stress gradient to address the following questions: (i) Does intraspecific competition become stronger than interspecific competition (the signature of ND) with stress? (ii) Does the performance (growth and fecundity) of species worsen with stress, perhaps even precluding them from occurring where stress is stronger? We used a soil depth gradient related to hydric stress because of the importance of hydrological niches in structuring plant communities (Silvertown et al. 2015), as it has been observed at the study site (Martorell et al. 2015). MATERIALS AND METHODS Study site Our study site is a semiarid grassland in Concepción Buenavista, Oaxaca, southern Mexico. It has an elevation of 2275 m a.s.l, an annual precipitation of 530.3 mm and mean temperature is 16.3°C. Most plants in this community are small (<5 cm tall) (Martorell et al. 2017). It shares the world record for the highest plant species richness with up to 25 species per dm2 (Martorell et al. 2017), making it ideal for studies on species coexistence. Soil depth changes rapidly over short distances (usually a few centimeters) and it is rarely deeper than 30 cm throughout the grassland (Martorell et al. 2015; Villarreal-Barajas and Martorell 2009). Water potential increases with soil depth: soil with a 2 cm depth has a potential three times more negative than 30 cm soil (Martorell and Martínez-López 2014; Villarreal-Barajas and Martorell 2009). Thus, shallow areas are dry and more stressful than deeper soils, and species performance responds to the soil depth gradient (Villarreal- Barajas and Martorell 2009). Species composition also changes with soil depth (Martorell et al. 2015). Deeper soils generally have more organic matter, but other abiotic soil variables are not correlated with depth (Martorell et al. 2015). Experimental setup Six annual species were selected for the experiment (Table 1). Because the experiment was conducted over a hydric stress gradient, we chose species with different drought stress tolerance (as measured by Martorell et al. 2015). We chose annuals in order to follow them throughout their entire life cycle. In our study site, annuals germinate at the beginning of the rainy season in July and set seed after the rainy season in November. Seeds for the six species were collected from >20 different mother plants in November 2016. In a 0.5 ha site protected from cattle, we selected 35 blocks differing in soil depth trying to cover the complete soil depth gradient of the study site (4–30 cm) (Fig. 1). Each block contained 13 0.15 m × 0.15 m plots. Soil depth was measured in each plot by inserting a metal rod into the soil until it reached parental material. Plots had a mean soil depth of 13.44 cm with a standard deviation of 7.66 cm. The deepest plot’s soil depth was 36 cm while the shallowest was 2 cm. All plants growing in the plots were manually removed. The plots were periodically hand-weeded to ensure that no other plants grew in them. Treatments were randomly assigned to the plots: (a) one individual plant (six treatments, one per species), (b) six individuals of the same species (again six treatments, one per species) or (c) six individuals, one of each of the selected species for the experiment (Fig. 1). In treatments (b) and (c) the number of individuals was the same and were sown close together (in a 3 cm diameter circle) to assure that plant would homogenously interact with each other. In treatment (c) a mixture of seeds from each species was sowed so that the Table 1: Species in the experiment, the family they belong to and they are optimum soil depth (Martorell et al. 2015) Species Family Niche soil depth position (cm) Aristida adscensionis L. Poaceae 29.61±2.64 Tagetes micrantha Cav. Asteraceae 24.24 ± 1.26 Plantago nivea Kunth Plantaginaceae 16.76 ± 0.54 Crusea diversifolia (Kunth) W.R. Anderson Rubiaceae 14.74 ± 0.58 Porophyllum linaria (Cav.) DC. Asteraceae 13.86 ± 1.09 Conyza filaginoides (DC.) Hieron. Asteraceae 9.8 D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGYRESEARCH ARTICLE 187JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 position of each species was different in each plot and no systematic bias would appear due to the distribution of plants in the plot. If conspecific and heterospecific individuals have the same effects on their neighbors, then individuals of a given species would have the same performance in treatments (b) and (c). Such result would correspond to a scenario with no ND. If there is ND we expected that performance would be better in the treatment with heterospecifics (Levine and HilleRisLambers 2009). Plants in plots with only one individual did not experience inter- or intraspecific competition, and thus informed us of the strength of environmental constraints at different soil depths. At the beginning of the rainy season in July, within a span of 3  days, the seeds were sown. We sowed many more seeds than strictly needed based on preliminary tests of seed germination in the field. In the case of Aristida adscensionis, Conyza filaginoides, Plantago nivea, Porophyllum linaria and Tagetes micrantha we sowed 3 seeds in the plots where 1 individual was needed, and 30 in the plots with 6 conspecifics. For Crusea divaricata, 4 and 40 seeds were sown, respectively. The extra individuals that germinated were removed in the first few days after germination to maintain the numbers required in each treatment. Individual plant survival was recorded monthly until the end of the plants’ life cycle, from July to November. A  total of 165 individuals of Aristida, 85 individuals of Conyza, 121 individuals of Crusea, 142 individuals of Plantago, 61 individuals of Porophyllum and 171 individuals of Tagetes survived until the end of the experiment. Individual plants that survived until the end of their life cycle were harvested, dried at air temperature, and weighed to obtain final aboveground biomass. We also recorded the number of inflorescences produced as a measure of fecundity (number of inflorescences correlates with number of seed produced, see online Appendix S1). No individual of Conyza reproduced in our experiments. Data analysis To analyze species’ response to environmental conditions and competition, we used linear mixed-effects models in the case of biomass and generalized mixed-effects models in the case of fecundity. Each individual plant in our treatments was modeled as a focal individual with biomass and fecundity as response variables to number of neighbors and soil depth. In such models, environmental constraints should be reflected in the effects of soil depth on biomass in the absence of neighbors, whereas competition can be appreciated as the effect of the number of neighbors. Despite the fact that at the beginning of the experiment the number of neighbors was the same in all experimental units, many of them died before the end of the experiment. Therefore, as a measure of the potential strength of Figure 1: Experimental setup within our 0.5 ha site and an example of our experimental blocks. Each block contained 13 plots with a randomly assigned treatment (a) one individual plant, (b) six individuals of the same species, or (c) six individuals one of each species. There were six variants of treatments (a) and (b), one for each of the six species in each block. D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGY RESEARCH ARTICLE 188 JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 competition in each plot, we used N =  sum of neighbors that were counted each month (for instance, a plant with five neighbors in September, four in October and four in November had N = 13. Note that we are assuming for simplicity that the effects of neighbors do not depend on the development of the focal individuals. In August all treatments with competition had six individuals, so the figure is not informative). Therefore, the models contained soil depth, N, and their interactions as explanatory fixed variables. We determined that squared soil depth and squared N (but not cubic terms, see online Appendix S2) should also be included in the model to allow for non- linear responses. Block was introduced as a random variable. Models were fit using lme4 (Bates et al. 2015) and glmmTMB (Brooks et al. 2017) for R (R Development Core Team 2018) As a measure of plant performance we used size and fecundity. Thus, the response variables were log biomass at the end of the experiment and the number of inflorescences produced. The error was normal for size, and Poisson or negative binomial (depending on the species; the selection of the type of error was done by comparing the Akaike Information Criterion (AIC) of models with each distribution) for reproduction. Survival could not be analyzed as a metric of performance. The number of competitors, N, was necessary for estimating competition, and thus crucial to test our hypotheses. However, the number of neighbors throughout the experiment was not fixed but instead depended on survival. Thus, quadrats with large N necessarily had high survival, leading to an artificial positive correlation between both variables. Two issues are of great importance for our study. First, the soil depth gradient is a potential axis over which niches differ (Martorell et al. 2015), so different species occur at different depths. Thus, coexistence may be possible at a spatial scale such that there is enough variation in soil depths. Thus we tested whether all species respond similarly to environmental constraints, i.e. soil depth. Second, niche differences depend on the relationship between intra and interspecific interactions, so we examined whether there is a different effect of conspecifics and heterospecifics on individual performance. Therefore, in a first phase of the analysis (Fig. 2) we constructed two different scenarios based on the general model described in the previous paragraph. In some scenarios we included focal species identity (the focal species is defined as the species on which the effect of the interaction are measured; Goldberg and Scheiner 2001) and its interactions with N and soil depth, so the responses to competition and depth differ between species (species-specific response models). Statistical interactions with focal species were omitted in other scenarios (same-response models; see online Appendix S3 for details on each one). Species-specific models were expected to perform better (have lower AICc) than same- response models if there was ND over the soil depth gradient. There are different ways in which the same scenario can be modeled, but still represent the same biological case. For instance, a model with the interactions focal species identity × soil depth + focal species identity × soil depth2, and a model only with focal species identity × soil depth, would both represent a case in which species niches are diferentiated along the stress gradient, and thus both are species-specific response models. We had no reason to choose a priori between these variants. Thus, we first fitted the most complex alternative for each of the two scenarios, and then used the dredge function of package MuMIn (Bartoń 2018) for R to obtain all the possible models nested in the complex ones (nested models have only a subset of the parameters in a more complex model; see online Appendix S3). Using AICc we compared the best model from each of the two scenarios. The model with the lowest AICc came from the species-specific response scenario (see results). Once we established that each species responded differently over the stress gradient and to competition (see results), we proceeded to a second phase of analysis (Fig. 2), in which we modeled each species separately. To establish if ND had an effect we again constructed two different scenarios: with and without taking neighbor identity into account. As before, the most complex model for each scenario was built and then dredged. In some instances, the difference in AICc between the best model in each of the two scenarios was <2 (Table 2). Thus, based on the available evidence there is no reason to favor any scenario, and the proper procedure is to consider both alternatives in light of the support each has. A way to do this is by means of model averaging using Akaike weights (Johnson and Omland 2004). From the averaged model, we obtained three performance curves for each species over the soil depth gradient: (i) growing alone (N = 0), (ii) growing with conspecifics (N = 10) and (iii) with growing with heterospecifics (N = 10). We used these curves to obtain interpretable measures of environmental constraints, competition and niche differences over the depth gradient as follows: let f 0 (d) be the biomass (or fecundity) of a plant grown alone in depth d, f c (d) the biomass of a plant grown with conspecifics, f h (d) the corresponding figure for heterospecifics. Environmental constraints, E(d), along the soil depth gradient were measured as E(d) = 1− f0(d) max(f0(d)) We used f 0 (d) because we wanted a measure of how each species responds to the environment without taking biotic interactions into account. Large values indicate strong environmental constraint, and if E(d) = 1 the species is completely excluded from a soil depth. In this case the environmental constraint becomes an environmental filter (Kraft et al. 2015). Intraspecific competition, C c (d), was defined as: Cc(d) = 1− fc(d) f0(d) If C c (d) = 0 there is no interaction, and values of 1 indicate competitive exclusion. Negative values indicate facilitation. A similar formula can be derived for interspecific competition, C h (d). As a measure of the magnitude of niche differences we used S(d) = fh(d) fc(d) − 1 S(d) is positive when there are niche differences and negative if interspecific competition is stronger than intraspecific competition. RESULTS In the first phase of the analysis where we constructed two different model scenarios (with or without focal species) to see if all species responded equally or if each one had a unique response to the soil depth gradient, we found that the best model included focal species in interaction with soil depth and N (Appendix S3). Among the models that did not include interactions with focal species, the best one had a ∆AICc = 4.21 and was the 17th best model (Appendix S3). Therefore, the evidence shows that each species has a unique response over the stress gradient, justifying the analysis of the responses of each one separately. In the second phase of the analysis we did not always find enough evidence to conclude that we should or should not include neighbor identity (conspecifics or heterospecifics) in the analyses. In about half of the cases the difference in AICc between the best model in each of the two scenarios was smaller than 2 (Table 2). Only one species had a ∆ > 2 favoring one scenario over the other for both growth and fecundity. This species was Crusea and the best model for growth and fecundity supported the scenario that includes neighbor identity. The rest of the species did not have enough evidence to support one scenario over the other. Thus, we used Akaike weighted model averaging to assess species’ response to environment and competition. D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGYRESEARCH ARTICLE 189JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 B io m a ss Depth sp 1 sp 2 B io m a ss Depth sp 1 sp 2? P h a se 1 : S p e ci e s re sp o n se s d iff e r o ve r th e g ra d ie n t? YES: include sta"s"cal interac"ons with depth NO: without sta"s"cal interac"ons with depth B io m a ss o f sp 1 Depth alone w heterospecifics w conspecifics Joint analysis (Not applicable) Analyze species separately P h a se 2 : H o w d o e n v ir o n m e n t a n d in te ra c" o n s a ff e ct b io m a ss ? D o e s n e ig h b o r id e n " ty m a # e r? B io m a ss o f sp 1 Depth alone w indifferent neighbors Allowing effects of neighbors to differ Neighbor iden"ty is irrelevant Model averaging (Akaike-weighted) B io m a ss o f sp 1 Depth alone w heterospecifics w conspecifics E n v ir o n m e n ta l co n st ra in t E (d ) Depth S tr e n g th Depth Intraspecific compe""on Cc(d) Stabilizing mechanisms S(d) Interspecificc compe""on Ch(d) C a lc u la " o n o $ h e v a ri a b le s o f in te re st Figure 2: Flow chart of our statistical analysis. In Phase 1 we tested different models to establish if there are niche differences over the stress gradient. On the left, the two species respond differently to depth, so there is a statistical interaction between focal species and soil depth. Because this was the case, we modeled each species’ biomass separately, which was Phase 2. Again, we had models corresponding to two different scenarios. On the left, biomass response to soil depth was predicted for number of neighbors and neighbor identity. This allows for niche differences if conspecifics have a more negative effect on biomass than heterospecifics. On the right are models where biomass response to soil depth is the same regardless of neighbor identity, and there are no niche differences. Finally, the Akaike-averaged model was used to calculate changes in environmental constraints, intra- and interspecific competition, and magnitude of niche differences. Table 2: ∆AIC for the models for growth and fecundity of each species Focal species Growth Fecundity Without neighbor identity With neighbor identity Without neighbor identity With neighbor identity Aristida 0 1.94 0 0.18 Conyza 0 1.43 — — Crusea 6.28 0 3.58 0 Plantago 0 0.76 0 2.98 Porophyllum 0.53 0 2.89 0 Tagetes 2.14 0 1.29 0 For each of the 12 species/performance combination, two possible scenarios are compared: one without neighbor identity, in which intra and interspecific effects are identical, and one with neighbor identity, in which the model supports differences in the effects of con- and heterospecific neighbors. ∆AIC = 0 means that this scenario included the best model for a species. Conyza did not reproduce in our experiment. D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGY RESEARCH ARTICLE 190 JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 0 0. 1 0. 2 0. 3 0 4 8 12 0. 00 0. 05 0. 10 0. 15 0 0. 1 0. 2 0. 3 0 20 40 60 0 0. 1 0. 2 0 2 4 6 0 0. 1 0. 2 0 2 4 6 5 10 15 20 25 30 0 0. 4 0. 8 5 10 15 20 25 30 0 20 40 60 Soil depth (cm) A b o v e g ro u n d b io m a s s ( g ) N u m b e r o f in fl o re s c e n s e s aditsirAaditsirA Conyza aesurCaesurC Plantago Plantago mullyhporoPmullyhporoP setegaTsetegaT Figure 3: Performance curves for each species along the soil depth gradient obtained from the average growth model (left column) and fecundity model (right column). Black line growing alone (N = 0), medium gray line growing with conspecifics (N = 10) and dotted light gray line growing with heterospecifics (N = 10). Dots represent the raw data. D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGYRESEARCH ARTICLE 191JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 Our E(d) values indicate that shallow soil reduced biomass up to 70% and fecundity up to 50%. Thus it can be a strong environmental constraint for some species like Aristida and Porophyllum (Fig.  4). However, shallow soils did not always affect species in a negative way, and the opposite was observed for four out of the six species studied. Such is the case for Crusea and Plantago (Fig. 4). Depth did not act as an environmental filter because species were never completely excluded at any soil depth. In general, performance diminished with the presence of other individuals indicating that negative interactions prevail in our system. These were especially strong when the interacting individuals were conspecific (Fig. 3). Intraspecific competition increased with soil depth in two species (Aristida and Porophyllum), diminished in two more (Crusea and Plantago) and remained nearly constant in the remaining species (Conyza and Tagetes) (Fig. 5). However, these changes were less pronounced than those in interspecific competition. As a result, in most species, interspecific competition intensity became much smaller than the intraspecific competition at some extreme of the gradient. Only two species experienced facilitation (Conyza and Plantago), and it was only found at the deep-soil end of the gradient (Fig. 5). However, this result must be taken with caution because we have few data regarding performance with heterospecifics in deep soil (Fig. 4). We found stronger intra-relative to interspecific competition indicating ND at least at some portion of the soil depth gradient for most of the studied species. The sole exception was Porophyllum (Fig. 5) that was outcompeted by the other species along the entire gradient. Some species experienced greater niche differences in shallow soils, as was the case for Aristida, Crusea, Tagetes and (for growth only) Plantago (Fig. 5). Conyza and Plantago showed greater niche differences in deep soils (Fig. 5). These two species were the same ones that could have experienced facilitation there (Conyza and Plantago). DISCUSSION Our findings suggest that the effects of drought stress depend on species identity and, in consequence, that each species performs 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 5 10 15 20 25 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 5 10 15 20 25 Soil depth (cm) E n v ir o n m e n ta l c o n s tr a in t E (d ) Aristida ogatnalPaesurC setegaTmullyhporoP Conyza Figure 4: Environmental constraints (E(d)) along the soil depth gradient for each species. Solid line is (E(d)) on growth. Dotted line is (E(d)) on fecundity. D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGY RESEARCH ARTICLE 192 JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 better at a different portion of the soil depth gradient. Hydric stress represented a strong environmental constraint in a few species that reduced biomass substantially. Unexpectedly, some species had poor performance in deep soils. This could be caused by the reduction of natural enemies in the dry soil (de Vries et al. 2012; Kardol et al. 2010; van der Putten et  al. 2016). Facilitation never occurred in shallow -1 .0 0. 0 1. 0 2. 0 -1 .0 0. 0 1. 0 2. 0 -1 .0 0. 0 1. 0 2. 0 -1 .0 0. 0 1. 0 2. 0 -1 .0 0. 0 1. 0 2. 0 5 10 15 20 25 -1 .0 0. 0 1. 0 2. 0 5 10 15 20 25 Soil depth (cm) M a g n itu d e aditsirAaditsirA Conyza aesurCaesurC Plantago Plantago mullyhporoPmullyhporoP setegaTsetegaT Figure 5: Magnitude of intraspecific competition (gray line), interspecific competition (dotted line) and niche differences (black line) along the soil depth gradient calculated from the growth (left column) and fecundity (right column) models. D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGYRESEARCH ARTICLE 193JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 soils, supporting the idea that competition predominates under extremely harsh conditions when stress is caused by the lack of a resource (Maestre et al. 2009). We found evidence for greater intra- than interspecific competition, i.e. for ND in five out of the six species, but it was usually restricted to some portion of the depth gradient. Because soil depth is heterogeneous, and species have ND only at some portions of the soil depth gradient, we believe that different species may persist in different soil depths. Thus, even if some species cannot coexist through ND at a given depth, we may find them co-occuring in a larger, heterogeneous area that comprises the depths that are appropriate for each of them. As expected, ND was slightly more frequent and stronger in shallow soils, apparently because niche overlap between species became small with stress. Thus, niche differences may be stabilizing coexistence at the stressful end of the gradient. However, this pattern was not observed in all species. Environmental constraints Soil depth frequently affected species performance (i.e. the biomass and fecundity of four out of our six species was affected by soil depth when growing alone), causing substantial reductions of performance in shallow soils in two species (Aristida and Porophyllum). This is in agreement with previous findings at the study site (Martorell et al. 2015; Villarreal-Barajas and Martorell 2009). However, these species were not completely excluded from shallow soils when growing alone, so depth was not an environmental filter sensu Kraft et al. (2015). Unexpectedly, four species had a better performance in shallow soils. This may be due to interactions with soil biota, rather than environmental constraints related to hydric stress. Drought may reduce the abundance of edaphic organisms such as microbial pathogens (de Vries et al. 2012; Kardol et al. 2010; van der Putten et al. 2016). Thus, plants may be favored in shallow soil with reduced natural-enemy abundance. Species interactions There is still a debate on interaction strength under extreme stress. The SGH states that competition should be weak under such conditions (Brooker and Callaghan 1998) while recent work suggests that competition may be strong at both ends of the stress gradient (Maestre and Cortina 2004; Maestre et al. 2009; Tielbörger and Kadmon 2000). Our results partially support the SGH, as competition tended to be stronger at the benign extreme of the gradient in four species (Aristida, Porophyllum, Crusea and Tagetes). However, the opposite was true for the other two species (Plantago and Conyza) and is in agreement with strong competition at both ends of a stress gradient. This may be expected when the stressor is the lack of a resource (Maestre and Cortina 2004; Maestre et al. 2009; Tielbörger and Kadmon 2000), which is water in our study. In these cases facilitation should be restricted to intermediate conditions (Maestre and Cortina 2004; Tielbörger and Kadmon 2000). This could be the case of Plantago and Conyza. The increasing trend in facilitation from shallow to intermediate depth soils is well supported by the data, but there were nearly no observations for facilitation at deep soils, so the conclusion that facilitation is greater when stress is lowest may be an unjustified extrapolation. Niche differentiation There was strong evidence for two forms of differentiation between species. One is local, and was observed wherever intraspecific competition was larger than interspecific competition, i.e. where there was positive S(d) in Fig.  5. The other form occurred at larger spatial scales, and was clearly related to soil depth heterogeneity. Because depth is closely related with water availability (Villarreal-Barajas and Martorell 2009), this suggests some form of hydrological-ND. This was the case, for instance, of Tagetes and Conyza: each species experiences greater intra- than interspecific competition in a different section of the depth gradient. This suggests that Conyza could experience stable coexistence with the other species in the community in deep soils, while Tagetes would be able to persist where the soil is shallow. In fact, previous work has found that species are segregated over the soil depth gradient at the study site due to substantial ND (Martorell et al. 2015). Please see online Appendix S1 for an in depth discussion of the use of ND as evidence for stable coexistence. We found no ND for Porophyllum anywhere along the soil depth gradient. Although this does not necessarily imply that stable coexistence is impossible (Broekman et al. 2019), the conditions for this to occur are more restrictive. Porophyllum is scarce: in 2012, we found only 5 individuals of Porophyllum in 3320 dm2 of the semiarid grassland, whereas the other five species were more abundant (average: 1985, range: 33–4705 individuals in 3320 dm2; Martínez-Blancas et al. 2018). However Porophyllum is abundant in isolated, tiny patches of soil where other species are absent. We believe that sink–source dynamics may be relevant for Porophyllum such that these isolated, competition-free patches are functioning as a source of seeds for the continuous grassland that acts like a sink (Mouquet and Loreau 2002, 2003; Pulliam 1988). As expected, ND tended to increase toward the shallow-soil end of the gradient (Aristida, Crusea and Tagetes, plus Plantago when growth was analyzed; Fig. 5), although exceptions were also observed. This suggests larger ND under stressful conditions. This pattern is consistent with some observations in nature, such as the greater divergence in traits related to water usage and growth-form diversity as aridity increases (Bernard- Verdier et al. 2012; Cody 1991). Nevertheless, in those studies, different species occur over the stress gradient (but see Wainwright et al. 2018). This is not our case, where the same species occurred everywhere. The fact that the same pattern is observed when comparing different species or the same species under different conditions suggests that ND is somehow promoted under stress, be it either through plasticity (Pérez- Ramos et al. 2019) or by speciation processes. The changes in ND over our gradient reflect changes in interaction strengths with stress. When ND became larger in shallow soils, it was because interspecific competition decreased more than the intraspecific one. Stronger reductions in interspecific competition, such as the ones we observed, may be expected if niche breath diminishes in stressful environments, reducing niche overlap. This is in line with the idea that traits that allow for survival in stressful environments come at the cost of a poor ability to deal with a broad range of conditions (Boulangeat et al. 2012; Crisp et al. 2009), and with the evidence for greater niche specialization as conditions become stressful (Bonetti and Wiens 2018; Boulangeat et al. 2012; Entling et al. 2007; Thuiller et al. 2004). In our study, it is likely that only specialized genotypes with narrow niches survived in shallow soils, while individuals with any niche breadth survived under more benign conditions. Abiotic factors such as water availability may structure communities by acting as environmental constraints or by altering species interactions. ND occurred at different scales: (i) a local one where spatial heterogeneity plays no role; and (ii) a larger one involving the soil depth heterogeneity, where different sets of species occur in neighboring patches that differ in soil depth. Our results suggest that ND are somewhat more important under stressful conditions. It seems that this is the result of a trade-off between niche breadth and the capacity to endure extreme stress. However, the number of species that we studied was small and there were exceptions to the general pattern, so it is still unclear if this result is widespread in nature. Supplementary material Supplementary material is available at Journal of Plant Ecology online. Appendix S1: Stabilizing niche differences and their relevance for species coexistence Appendix S2: Model selection Appendix S3: Models used to estimate biomass as a function of depth and interactions D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 JOURNAL OF PLANT ECOLOGY RESEARCH ARTICLE 194 JOURNAL OF PLANT ECOLOGY | VOL 13 | April 2020 | 185–194 Funding This research was supported by the Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México (IN212618). Conflict of interest statement. None declared. Acknowledgements We thank D.  García-Meza, M.  Romero, V.  Zepeda, A.  Cervantes, I.  Belastegui, C.  Vazquez, R.  Díaz-Talamantes and P.  Trujano who helped in the field and lab. J. A. Meave and K. Boege provided valuable feedback during the conception and realization of the study. This work is presented in partial fulfillment toward A.M.-B’s doctoral degree in the Posgrado en Ciencias Biológicas at UNAM. Agradecemos a la comunidad de Concepción Buenavista por su amistad y apoyo. REFERENCES Adler  PB, Fajardo  A, Kleinhesselink  AR, et  al. (2013) Trait-based tests of coexistence mechanisms. Ecol Lett 16:1294–306. Anderson HM, Hutson V, Law R (1992) On the conditions for permanence of species in ecological communities. Am Nat 139:663–8. Bartoń K (2018) MuMIn: Multi-Model Inference. R Package Version 1.42.1. https:// CRAN.R-project.org/package=MuMIn (1 April 2019, date last accessed) Bates D, Mächler M, Bolker BM, et al. (2015) Fitting linear mixed-effects models using lme4. J Stat Softw 67:1–48. Bernard-Verdier M, Navas ML, Vellend M, et  al. (2012) Community assembly along a soil depth gradient: contrasting patterns of plant trait convergence and divergence in a Mediterranean rangeland. J Ecol 100:1422–33. Bimler MD, Stouffer DB, Lai HR, et al. (2018) Accurate predictions of coexistence in natural systems require the inclusion of facilitative interactions and environmental dependency. J Ecol 106:1839–52. Bonetti  MF, Wiens  JJ (2018) Evolution of climatic niche specialization : a phylogenetic analysis in amphibians. Proc R Soc B Biol Sci 281:1–9. Boulangeat  I, Lavergne  S, Van  Es  J, et  al. (2012) Niche breadth, rarity and ecological characteristics within a regional flora spanning large environmental gradients. J Biogeogr 39:204–14. Broekman MJE, Muller-Landau HC, Visser MD, et al. (2019) Signs of stabilisation and stable coexistence. Ecol Lett 22:1957–75. Brooker  RW, Callaghan  TV (1998) Environmental gradients: a model. Oikos 81:196–207. Brooker  RW, Maestre  FT, Callaway  RM, et  al. (2008) Facilitation in plant communities: the past, the present, and the future. J Ecol 96:18–34. Brooks M, Kristensen K, Van Benthem KJ, et al. (2017) glmmTMB balances speed and flexibility among packages for Zero-inflated Generalized Linear Mixed Modeling. R J 9:378–400. Callaway RM (2007) Positive Interactions and Interdependence in Plant Communities, 1st edn. Dordrecht, Netherlands: Springer. Chase  JM (2007) Drought mediates the importance of stochastic community assembly. Proc Natl Acad Sci U S A 104:17430–4. Chesson P (2000) Mechanisms of maintenance of species diversity. Annu Rev Ecol Syst 31:343–66. Chisholm RA, Pacala SW (2011) Theory predicts a rapid transition from niche- structured to neutral biodiversity patterns across a speciation-rate gradient. Theor Ecol 4:195–200. Cody ML (1991) Niche theory and plant growth form. Vegetatio 97:39–55. Crisp MD, Arroyo MT, Cook LG, et al. (2009) Phylogenetic biome conservatism on a global scale. Nature 458:754–6. de Vries FT, Liiri ME, Bjørnlund L, et al. (2012) Legacy effects of drought on plant growth and the soil food web. Oecologia 170:821–33. Entling  W, Schmidt  MH, Bacher  S, et  al. (2007) Niche properties of Central European spiders: shading, moisture and the evolution of the habitat niche. Glob Ecol Biogeogr 16:440–8. Goldberg  DE, Scheiner  SM (2001) ANOVA and ANCOVA: field competition experiments. In Scheiner SM, Gurevich J (eds). Design and Analysis of Ecological Experiments. New York, United States of America: Oxford University Press, 77–98. Gravel D, Canham CD, Beaudet M, et al. (2006) Reconciling niche and neutrality: the continuum hypothesis. Ecol Lett 9:399–409. Grime JP (1977) Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory. Am Nat 111:1169–94. Halloy S (1990) A morphological classification of plants, with special reference to the New Zealand alpine flora. J Veg Sci 1:291–304. Hart SP, Marshall DJ (2013) Environmental stress, facilitation, competition, and coexistence. Ecology 94:2719–31. Hubbell  SP, Foster RB (1986) Biology, chance and history and the structure of tropical rain forest tree communities. In Diamond JM, Case TJ (eds). Community Ecology. New York, United States of America: Harper & Row, 314–29. Johnson JB, Omland KS (2004) Model selection in ecology and evolution. Trends Ecol Evol 19:101–8. Kardol P, Cregger MA, Campany CE, et al. (2010) Soil ecosystem functioning under climate change: plant species and community effects. Ecology 91:767–81. Kraft NJB, Adler PB, Godoy O, et al. (2015) Community assembly, coexistence and the environmental filtering metaphor. Funct Ecol 29:592–9. Levine  JM, HilleRisLambers  J (2009) The importance of niches for the maintenance of species diversity. Nature 461:254–7. Lieberman M, Lieberman D (2007) Nearest-neighbor tree species combinations in tropical forest: the role of chance, and some consequences of high diversity. Oikos 116:377–86. Maestre  FT, Callaway  RM, Valladares  F, et  al. (2009) Refining the stress-gradient hypothesis for competition and facilitation in plant communities. J Ecol 97:199–205. Maestre FT, Cortina J (2004) Do positive interactions increase with abiotic stress? A test from a semi-arid steppe. Proc Biol Sci 271:S331–3. Martinez-Blancas A, Paz H, Salazar GA, et al. (2018) Related plant species respond similarly to chronic anthropogenic disturbance: implications for conservation decision-making. J Appl Ecol 55:1860–70. Martorell C, Almanza-Celis CAI, Pérez-García EA, et al. (2015) Co-existence in a species-rich grassland: competition, facilitation and niche structure over a soil depth gradient. J Veg Sci 26:674–85. Martorell C, Martínez-López M (2014) Informed dispersal in plants: Heterosperma pinnatum (Asteraceae) adjusts its dispersal mode to escape from competition and water stress. Oikos 123:225–31. Martorell C, Zepeda V, Martínez-Blancas A, et al. (2017) A diversity world record in a grassland at Oaxaca, México. Bot Sci 95:1–7. Matías L, Godoy O, Gómez-Aparicio L, et al. (2018) An experimental extreme drought reduces the likelihood of species to coexist despite increasing intransitivity in competitive networks. J Ecol 106:826–37. Mouquet  N, Loreau  M (2002) Coexistence in metacommunities: the regional similarity hypothesis. Am Nat 159:420–6. Mouquet  N, Loreau  M (2003) Community patterns in source-sink metacommunities. Am Nat 162:544–57. Pérez-Ramos IM, Matías L, Gómez-Aparicio L, et al. (2019) Functional traits and phenotypic plasticity modulate species coexistence across contrasting climatic conditions. Nat Commun 10:2555. Pulliam HR (1988) Sources, sinks, and population regulation. Am Nat 132:652–61. R Development Core Team (2018) R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. http:// www.R-project.org (29 November 2019, date last accessed) Silvertown J, Araya Y, Gowing D (2015) Hydrological niches in terrestrial plant communities: a review. J Ecol 103:93–108. Thuiller W, Lavorel S, Midgley G, et al. (2004) Relating plant traits and species distributions along bioclimatic gradients for 88 Leucadendron taxa. Ecology 85:1688–99. Tielbörger K, Kadmon R (2000) Temporal environmental variation tips the balance between facilitation and interference in desert plants. Ecology 81:1544–53. van  der  Putten WH, Bradford MA, Pernilla Brinkman E, et  al. (2016) Where, when and how plant–soil feedback matters in a changing world. Funct Ecol 30:1109–21. Villarreal-Barajas  T, Martorell  C (2009) Species-specific disturbance tolerance, competition and positive interactions along an anthropogenic disturbance gradient. J Veg Sci 20:1027–40. Wainwright CE, HilleRisLambers  J, Lai HR, et  al. (2018) Distinct responses of niche and fitness differences to water availability underlie variable coexistence outcomes in semi-arid annual plant communities. J Ecol 107:293–306. Wang X, Wiegand T, Kraft NJ, et al. (2016) Stochastic dilution effects weaken deterministic effects of niche-based processes in species rich forests. Ecology 97:347–60. Zambrano J, Marchand P, Swenson NG (2017) Local neighbourhood and regional climatic contexts interact to explain tree performance. Proc R Soc B Biol Sci 284:1–9. D o w n lo a d e d fro m h ttp s ://a c a d e m ic .o u p .c o m /jp e /a rtic le -a b s tra c t/1 3 /2 /1 8 5 /5 6 9 9 9 2 2 b y U n iv e rs id a d N a c io n a l A u to n o m a d e M e x ic o u s e r o n 2 2 M a y 2 0 2 0 37 Appendix S1 Stabilizing niche differences and their relevance for species coexistence Modern coexistence theory leads our understanding of diversity maintenance in communities (Chesson 2000, Adler et al. 2007). This theory recognizes that both niche and neutral processes act simultaneously on communities. Niche differences allow coexistence because species’ requirements may be different enough to preclude competitive exclusion. The signature of niche differentiation is a stronger intra- than interspecific competition, making species limit their own populations more than they are limited by others (Levine and HilleRisLambers 2009). Under such conditions, overall intraspecific competition diminishes with the species’ abundance (Anderson et al. 1992, Chesson 2000), allowing the population’s growth rate to increase. Such “increase when rare” process actively precludes extinction and is thus a stabilizing mechanism (Chesson 2000). In contrast, equalizing mechanisms allow coexistence because species have similar competitive abilities, so no one tends to exclude others from the community (Chesson 2000, Adler et al. 2007). Two species with slightly asymmetric competitive abilities need only a small difference in their niches to coexist. If one species is competitively superior, stabilizing mechanisms need to be stronger (niche differences must be larger) for them to coexist (Chesson 2000, Adler et al. 2007). Niche differences relevant for the stabilization of coexistence are those that affect population growth rates (Chesson 2000). However, we measured performance as biomass, an attribute of the individual. It would thus be important to assess whether biomass reflect population dynamics. To do so, we calculated the Spearman correlation between mean biomass in each plot and the number of inflorescences produced per seed sown. The number of inflorescences is highly correlated with the number of seeds 38 produced (Primack, 1978; Zimmerman & Gross, 1984; Devlin, 1989). Using data from Zepeda & Martorell (2019a) we found a high correlation between inflorescence and seed numbers in the two species shared with our study, Aristida (r = 0.87) and Tagetes (r = 0.94). The correlation is high (mean r = 0.92) for all of the other annual species in Zepeda & Martorell (2019). Thus, the number of inflorescences produced per seed sown is a close proxy for the number of seeds produced per seed sown, i.e., the net reproductive rate. In classic population ecology, the net reproductive rate indicates the extent of growth that annual plant populations experience in a year (Leverich & Levin, 1979; Caswell 2001). Because the number of seeds is proportional to the number of inflorescences (Zimmerman & Gross 1984; Devlin 1989), our procedure would have no effect on Spearman’s r, which is insensitive to measurement scale. The biomass of every species correlated strongly with the net reproductive rate (Table S1.1). Thus, it is a good proxy of population-level performance, and as such is useful for measuring niche differences mechanisms and proceeding with the analysis. Stronger intra- than interspecific competition, making species limit their own populations more than others limit them is a necessary (although not sufficient) diagnostic for niche differentiation and a condition for stable coexistence in simple situations. However, how our measures of performance, even if they correlate with population level performance, are useful for interpreting stabilizing niche differences is not straightforward. We therefore present the following considerations. 39 Table S1.1. Spearman correlation between biomass and net reproductive rate p rho Aristida adscensionis < 0.001 0.73 Crusea divaricata < 0.001 0.97 Plantago nivea < 0.001 0.94 Porophyllum linaria < 0.001 0.95 Tagetes micrantha < 0.001 0.84 Consider a system with two species, i and j. Assuming a Lotka-Volterra model for competition, biomass can be modeled as , (1) where in the biomass of species i at time t in the intraspecific competition treatment, is the initial biomass of species i, ri is the growth rate, is the intraspecific competition coefficient, and Ni is the number of conspecifics in the plot. The Lotka- Volterra model requires biomass to be replaced by population numbers, but by changing the units of r, we can transform size into numbers. In our system where biomass is closely correlated with population growth rates this poses no problem. When competition is interspecific, , (2) where is the interspecific effect of j on i, and Nj is the number of species j individuals in the plot. For species j the respective equations are: and (3) w i ,t intra = w i ,0 e r i 1−α ii N i( ) w i ,t intra w i ,0 α ii w i ,t inter = w i ,0 e r i 1−α ij N j( ) α ij w j ,t intra = w j ,0 e r i 1−α jj N j( ) 40 . (4) When estimating competition and niche differentiation we kept Ni = Nj = 10. Our estimate of frequency dependence (see main text) for species i was therefore , (5) which indicates negative frequency dependence if aii > aij, i.e., if the effect of i on itself is greater that the effect of j on i. A similar equation can be derived for species j. For there to be stable coexistence between i and j, the following must hold (Broekman et al., 2019): . (6) If Si(d) and Sj(d) are both negative, then aij / aii and aji / ajj are both < 1 and so is their product, which is the figure in the square root in eq. (6). This means that stable coexistence can occur if there is negative frequency dependence for both species in a two species system. However, we do not have data on pairwise effects, but instead measures of how five different heterospecific individuals affect each focal species. This does not necessarily imply that stabilization occurs between each possible pair of species (Chesson 2000). As a final note of caution, our experiment was conducted over a single year, so the effects of interannual variability were absent. Thus, we missed contributions to stable coexistence that are promoted by temporal fluctuations (e.g., storage effect; Chesson 2000). However, fluctuation-dependent stabilization plays a minor role in coexistence at the study site (Zepeda and Martorell 2019b) w j ,t inter = w j ,0 e r i 1−α ji N i( ) S i (d ) =1− w i ,t inter w i ,t intra =1− e r i 1−α ij 10( ) e r i 1−α ii 10( ) α ij α ji α ii α jj <1 41 References Adler PB, HilleRislambers J, Levine JM (2007) A niche for neutrality. Ecol Lett 10:95– 104. Anderson HM, Hutson V, Law R (1992) On the conditions for permanence of species in ecological communities. Am Nat 139:663–668. Broekman MJE, Muller-Landau HC, Visser MD, et al. (2019) Signs of stabilisation and stable coexistence. Ecol Lett 22:1957–1975. Chesson P (2000) Mechanisms of maintenance of species diversity. Annu Rev Ecol Syst 31:343–366. Caswell H (2001) Matrix Population Models. Sunderland, MA: Sinnauer Associates. Devlin B (1989) Components of seed and pollen yield of Lobelia cardinalis: Variation and Correlations. Am J Bot 76:204–2014 Leverich WJ, Levin DA (1979) Age-specific survivorship and reproduction in Phlox drummondii. Am Nat 113:881–903. Levine JM, HilleRisLambers J (2009) The importance of niches for the maintenance of species diversity. Nature 461:254–257. Primack R (1978) Regulation of Seed Yield in Plantago. J Ecol 66:836–847. Zepeda V, Martorell C (2019a) Seed mass equalises the strength of positive and negative plant-plant interactions in a semi-arid grassland. Oecologia 190: 287–296. Zepeda V, Martorell C (2019b) Fluctuation-independent niche differentiation and relative non-linearity drive coexistence in a species-rich grassland. Ecology 100: e02726. Zimmerman M, Gross RS (1984) The relationship between flowering phenology and seed set in an herbaceous perennial plant, Polemonium folisissimum Gray. Am Midl Nat 42 111:185–191. 43 Appendix S2 Model selection We did not have any a priori reason to suspect that the effects of depth and interactions on species growth should be linear or not, or to prefer any specific form non-linear response. More complex functions (high-degree polynomials, for instance) have better fits, but this may be the result of overfitting, i.e., the inclusion of experimental error into the model. Because experimental error is random and not expected to be replicated, the predictive power (the ability to predict new observations) of an overfitted model should be low. This principle has been used, for instance, to determine the appropriate level of model complexity in generalized additive models (Wood 2008). To determine the appropriate wriggliness that our model should allow to describe the responses of plant biomass, we used cross validation to measure the predictive power of different models that included quadratic and cubic terms. We first fitted all of the nested models within a model that took into account quadratic and cubic forms of our variables of interest. That is, a model with the total number of neighbors (same species and other species neighbors grouped together) and soil depth as quadratic and cubic terms and we also included all possible interactions between linear, quadratic, and cubic neighbor and soil depth terms: Where g(B) is a function (the link function) of biomass, d is soil depth, N is total number of neighbors and ε is the random effect of the block. We applied the function MuMIn to this model to fit all of the nested models. We did this for each one of the six species in our experiment separately. We tried identity, log, and inverse link functions. g(B) = β 0 +β 1 d +β 2 N +β 3 d 2 +β 4 N 2 +β 5 d 3 +β 6 N 3 +β 7 dN +β 8 dN 2 +β 9 dN 3 +β 10 d 2N +β 11 d 2N 2 +β 12 d 2N 3 +β 13 d 3N +β 14 d 3N 2 +β 15 d 3N 3 +ε k 44 Once all of the nested models for each species had been fitted, we performed cross validation on five of them. The models that we selected for cross validation were: 1) the model with the lowest AICc, 2) the model with the lowest AICc and no cubic terms, 3) the model with the lowest AICc, no cubic terms, and no interaction between quadratic terms, 4) the model with the lowest AICc, no cubic terms, and no N2, and 5) the model with the lowest AICc, no cubic terms, and no d2. We performed leave-one-out cross validation on these models by leaving out one data point each time we fitted them until each of the data points had been excluded once from the fitted models. We then estimated the predictive error as , where Bi is the biomass of the plant that was left out and is the biomass estimated by the model fitted to the respective reduced dataset. Finally, we determined which type of model had the lowest predictive error for most of the species. We generally found that cubic terms did not increase the predictive power of our models (Table S2.1). We also found that the identity link increased the models’ predictive power, and that models with inverse link error distribution tended to not converge (Table S2.1). Thus, we decided to use models with identity and to include only quadratic terms in the models used to test our predictions. Table S2.1. Predictive model errors from the leave-one-out runs for each of the models we selected 1) through 5). The numbers in bold are the smallest model errors. More than one model may have the same error because the best model in the different scenarios was the same. Some models did not converge with certain links. B i − B̂ i( ) 2 ∑ B̂ i 45 Species Link functions 1) 2) 3) 4) 5) Aristida adscensionis Identity 0.1516 0.1504 0.1554 0.1559 0.1505 Log 0.53 0.1717 0.3751 0.1583 0.2685 Inverse Models did not converge Conyza filaginoides Identity 0.1168 0.1081 0.1081 0.1081 0.1272 Log 0.1256 0.1276 0.1276 0.1276 0.1253 Inverse Models did not converge Crusea diversifolia Identity 0.4820 0.4629 0.4629 0.4629 0.4624 Log 0.5223 0.4814 0.4814 0.4814 0.4814 Inverse 0.5089 0.5114 0.5114 0.5114 0.5114 Plantago nivea Identity 0.27 0.2582 0.2582 0.2582 0.2582 Log 0.2760 0.2884 0.2870 0.2860 0.2842 Inverse 0.2797 0.2858 0.2858 0.2858 0.2858 Portulaca pilosa Identity 0.1297 0.1014 0.1187 0.1014 0.0981 Log 0.1356 0.1249 Model did not converge 0.1250 0.1250 Inverse 0.1285 0.1274 0.1274 0.1274 0.1274 Tagetes micrantha Identity 4.281 4.51 4.61 4.03 4.49 Log 4.87 4.71 4.71 3.96 4.71 Inverse 4.718 4.718 5.09 Models did not converge 46 Appendix S3 Models used to estimate biomass as a function of depth and interactions The phases of analysis described here are the same shown in figure 1 of the main text. Phase one of the analyses We fit all of the possible nested models by applying the function “dredge” from the R package MuMIn (R Core Team, 2017) on the most complex model of each of our scenarios. This function uses a model (in this the most complex one for each scenario) as its main argument; it then fits all the possible nested models within the model it was given. It returns an object listing each model’s parameter estimations and AICc. The models we used to run the dredge function for each scenario are the following. The model corresponding to the same response scenario (no focal species interaction with depth or neighbors) was: Where B is biomass, d is soil depth, Nc is the density of conspecifics, Nh is the density of heterospecifics, s is a categorical variable regarding focal species identity and εk is the random effect of the kth block. Note that this model fits two surfaces: one relating B to d and N for conspecific neighbors, and one for heterospecifics. Both surfaces coincide when Nh and Nc are equal to zero (biomass is the same because the focal plant is alone) but may differ elsewhere. The identity of the focal species is included here as an effect on the intercept because there is no reason to expect that the biomass of different species is the same. The model corresponding to the species-specific response scenario was: B = β 0 + s + β 1 d + β 2 Nc + β 3 Nh + β 4 d 2 + β 5 Nc 2 + β 6 Nh 2 + β 7 dNc + β 8 dNh + β 9 d 2 Nc + β 10 d 2 Nh + β 11 dNc 2 + β 12 dNh 2 + β 13 d 2 Nc 2 + β 14 d 2 Nh 2 + ε k 47 In total we fitted 33271 models and then determined by AICc which one of these models predicted final biomass better (Δ AIC < 2). If the best model included an effect of the different focal species other than a change in the intercept, then running specific models for each species may be justified. If, on the contrary, the best model by two AIC points only included a change in the intercept by the focal species, then all species behave the same way and using species identity for each model is not justified. Because the best model contained focal species in interaction with N and soil depth, we ran separate models for each species in Phase 2 of the analysis. Phase two of the analyses After determining that each species’ fitness responded differently to neighbors and the soil depth gradient, we obtained models for each of the six species separately. Again we had two groups of models, one group included the effect of conspecific and heterospecific neighbors pooled together (no stabilizing mechanism) and the other group had models that included the different effects of heterospecific and conspecific neighbors (allowing for stabilizing mechanisms). The most complex model for the scenario without stabilizing mechanism was: B = β 0 + s+β 1 ds+β 2 Ncs+β 3 Nhs+β 4 d 2 s+β 5 Nc 2 s+β 6 Nh 2 s+β 7 dNcs+β 8 dNhs+β 9 d 2 Ncs+β 10 d 2 Nhs+β 11 dNc 2 s+β 12 dNh 2 s+β 13 d 2 Nc 2 s+β 14 d 2 Nh 2 s+ε k B = β 0 +β 1 d +β 2 N +β 3 d 2 +β 4 N 2 +β 5 dN +β 6 d 2 N +β 7 dN 2 +β 8 d 2 N 2 +ε k 48 where N is the number of neighbors, regardless of their identity. The most complex model for the scenario with possible stabilizing mechanisms was: As before, we applied the dredge function on the most complex models pertaining to each scenario and obtained all the possible nested models. In total, we fitted 16636 models for each species and calculated the ΔAICc for each one. We then averaged the predicted values of each model using their Akaike weights and obtained performance curves for each species over the soil depth gradient. These curves were used to obtain easily interpretable measures of environmental constraints, competition and stabilizing niche differences for Phase 3 of the analysis. References R Core Team 2017. R: A language and environment for statistical computing. – R Foundation for Statistical Computing. Wood, S. N. 2008. Fast stable direct fitting and smoothness election for generalized additive models. - J. Am. Stat. Assoc. 70: 495–518 g = β 0 +β 1 d +β 2 c+β 3 h+β 4 d 2 +β 5 c2 +β 6 h2 +β 7 dc+β 8 dh+β 9 d 2c+β 10 d 2h+β 11 dc2 +β 12 dh2 +β 8 d 2c2 +β 13 d 2h2 + Zu B = β 0 +β 1 d +β 2 Nc+β 3 Nh+β 4 d 2 +β 5 Nc 2 +β 6 Nh 2 +β 7 dNc+β 8 dNh+β 9 d 2 Nc+β 10 d 2 Nh+β 11 dNc 2 +β 12 dNh 2 +β 13 d 2 Nc 2 +β 14 d 2 Nh 2 +ε k 49 Capítulo III. Species alliances and hidden niche dimensions drive species clustering along a hydric gradient in a semiarid grassland A. Martínez-Blancas, I. X. Belaustegui y C. Martorell. 2022. Ecology Letters. 25:2651-2662. Ecology Letters. 2022;25:2651–2662. wileyonlinelibrary.com/journal/ele | 2651© 2022 John Wiley & Sons Ltd. INTRODUCTION Niche differentiation is usually thought of as a requi- site for species coexistence. Nevertheless, a number of theoretical studies and empirical data find a ten- dency for groups of species in the community to resem- ble each other in their niches (D'Andrea et al.,  2020; Scheffer et al., 2015; Scheffer & van Nes, 2006a; Segura et al., 2011; Vergnon et al., 2009). This raises the ques- tion of how diversity is maintained. When species in the same community occupy different niches they compete more strongly with their conspecifics than with hetero- specifics. Thus, when species become less abundant and conspecific interactions become rare, competition is al- leviated and the population grows, avoiding extinction. Without such mechanism, it remains unclear how groups of species with similar niches may persist in the commu- nity (Broekman et al., 2019; Chesson, 2000). The emergent neutrality model builds on the fact that coexistence does not depend exclusively on the niche, but on fitness differences that arise from the species com- petitive abilities. Species will coexist only if niche dif- ferences overcome fitness differences (Adler et al., 2007; Chesson,  2000). Thus, coexistence is possible despite minimal niche differentiation if fitness differences are small. Under emergent neutrality (Figure 1b), species self- organise into regularly spaced clusters of two or more ecologically similar species along a niche axis (water availability, Figure 1) (Scheffer & Van Nes, 2006b). Such pattern contrasts with the expectation of species occur- ring evenly over the niche axis (Figure  1a) as a result of competitive exclusion between species with similar niches. Nevertheless, in the emergent neutrality model, coexistence within clusters is made possible by small fitness differences, and coexistence between clusters by large niche differences (Vergnon et al., 2009). Emergent L E T T E R Species alliances and hidden niche dimensions drive species clustering along a hydric gradient in a semiarid grassland Alejandra Martínez- Blancas1,2 | Ian Xul Belaustegui1 | Carlos Martorell1 Received: 10 May 2022 | Revised: 12 September 2022 | Accepted: 13 September 2022 DOI: 10.1111/ele.14122 1Facultad de Ciencias, Departamento de Ecología y Recursos Naturales, Universidad Nacional Autónoma de México, Coyoacán, Ciudad de México, Mexico 2Posgrado en Ciencias Biológicas, Universidad Nacional Autónoma de México, Coyoacán, Ciudad de México, Mexico Correspondence Carlos Martorell, Facultad de Ciencias, Universidad Nacional Autónoma de México, Circuito Exterior S/N, Cd. Universitaria 04510, CDMX, Mexico. Email: martorell@ciencias.unam.mx Funding information Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México, Grant/Award Number: PAPIIT IN- 212618 Editor: Eric Seabloom Abstract Clustering of species with similar niches or traits occurs in communities, but the mechanisms behind this pattern are still unclear. In the emergent neutrality model, species with similar niches and competitive ability self- organise into clusters. In the hidden- niche model, unaccounted- for niche differences stabilise coexistence within clusters. Finally, clustering may occur through alliances of species that facilitate each other. We tested these hypotheses using population- growth models that consider interspecific interactions parameterised for 35 species using field data. We simulated the expected community dynamics under different species- interaction scenarios. Interspecific competition was weaker within rather than between clusters, suggesting that differences in unmeasured niche axes stabilise coexistence within clusters. Direct facilitation did not drive clustering. In contrast, indirect facilitation seemingly promoted species alliances in clusters whose members suppressed common competitors in other clusters. Such alliances have been overlooked in the literature on clustering, but may arise easily when within cluster competition is weak. K E Y W O R D S coexistence, competition, emergent neutrality, environmental gradients, hidden niches, hydric stress, indirect facilitation, indirect interactions, species clumping, species interactions 2652 | DRIVERS OF SPECIES CLUSTERING neutrality has been observed in theoretical and empirical studies. Several computational simulations of communi- ties structured by competition result in self- organised clusters of species (Bonsall et al.,  2004; Doebeli & Dieckmann, 2003; Roelke & Eldridge, 2008; Scheffer & van Nes, 2006a; Scranton & Vasseur, 2016). In plankton communities, groups of similar- sized species emerge be- cause of shared resource use (Segura et al., 2011; Vergnon et al., 2009). Clusters of species of similar body size have also been observed in communities of beetles (Scheffer et al., 2015) and birds (Thibault et al., 2011). In a trop- ical forest, groups of plants that resemble each other in traits associated with competition for light have recently been observed (D'Andrea et al., 2020). Using life- history traits, Herault and Thoen (2009) found groups of simi- lar plants in open habitats. However, only one study has explored if clusters of species are maintained by similari- ties in their niche and fitness (Segura et al., 2011). The problem with the emergent neutrality model is that small fitness differences can easily overcome the tiny niche differences between species in the same clus- ter, making long- term coexistence unlikely. In studies, only transient clusters are observed (Scheffer & Van Nes, 2006b), and at least some degree of niche differenti- ation is required for them to persist (Barabás et al., 2013). Thus, as an alternative to the emergent neutrality model, it has been proposed that clusters arise when differen- tiation over unobserved niche axes, known as hidden niches (Figure 1c), stabilise the coexistence between sub- sets of species in the community (Barabás et al., 2013; D'Andrea et al., 2020). Differentiation over hidden niche axes should reduce competition between species in the same cluster (Figure 1c), whereas the huge niche overlap within clusters in the emergent neutrality model should result in an extremely strong within- cluster competition (Figure 1b). Finally, direct and indirect positive interactions may also produce species clustering. Indirect effects occur when a species exerts an effect on a second one by affect- ing the abundance of a third species that interacts with both (Worthen & Moore, 1991). For example, if a species constrains a competitively dominant species, it may indi- rectly facilitate the weaker competitors that are hindered by the competitively dominant species (Levine,  1999). Such indirect facilitation has been observed between a few herbs (Levine, 1999). These positive interactions may produce species alliances, that is, groups of species that facilitate each other (Figure 1d). Such facilitation may be direct, but also indirect by either constraining a com- mon antagonist or facilitating a common benefactor. The term alliance was introduced by Wilson  (1986) for a necrophagous arthropod association involving beetles and mites. Beetles directly benefit mites by transporting them to new habitats, while mites indirectly enhance bee- tle survival by feeding on the eggs of its competitors. In consequence, other mite species who also use beetles for transportation are indirectly facilitated (Wilson,  1983; Wilson & Knollenberg, 1987). We tested whether the species of a semiarid grassland assemble into clusters over a hydrological niche axis. We F I G U R E 1 Expected patterns under different species clustering mechanisms. each line represents a species. Under emergent neutrality (panel b) there is a complete overlap of niches and thus competition is strong within clusters but weak between them. Coexistence is possible due to equal competitive abilities. In the hidden niche scenario (panel c) competition is weak because niche differentiation between species occurs over unobserved axes. Species that differ in their niches over unobserved axes are represented in different colours. Under the species alliance hypotheses (panel d), species may become clustered because they facilitate each other directly (purple) or indirectly by suppressing strong common competitors from other clusters (black). In all cases, removal of interactions weakens clustering (panel a). If there are species alliances due to direct facilitation, the sole removal of facilitative interactions should suffice to cause such weakening. Furthermore, the removal of all but one species in a cluster should cause the abundance of the single remaining species to increase strongly in the emergent neutrality scenario, weakly in the hidden niche scenario, and to decrease if there are species alliances. (a) (b) (c) (d) | 2653MARTÍNEZ- BLANCAS ET AL. then analysed if there was evidence for emergent neutral- ity or if rather hidden niches or species alliances are be- hind species clustering. We chose to use a hydrological axis because its importance component for plant niches (Silvertown,  2004) and because of its particular impor- tance in determining niche differentiation in the studied grassland (Martínez- Blancas & Martorell, 2020; Martorell et al., 2015). Because all mechanisms of clustering depend on interactions, we first parameterised population- growth models for 35 species using field data to estimate pairwise competition and facilitation, which is frequent at the study site (Martorell & Freckleton,  2014; Villarreal- Barajas & Martorell, 2009). These models were then used to simulate the expected community dynamics under different scenar- ios of species interactions and community membership. To test if clusters arise from species interactions, as expected from all three hypotheses, we compared the clustering predicted by the model with and without interspecific in- teractions. To establish if clusters were neutral or if they involve hidden niches, we compared the strength of within- and between cluster competition. We would expect within cluster competition to be greater than that between clus- ters if there is emergent neutrality, but the opposite under a hidden- niche scenario (Figure 1). To determine whether clusters arise from alliances caused by direct facilitation, we assessed whether clusters disappear when facilitation is removed from the model and whether pairwise facilitation was stronger within- than between- clusters (Figure 1). To test if species alliances result from indirect interactions, we assessed if species were favoured by the joint occurrence of the other species in its cluster. If this is the case, species abundance would increase when the whole set of its allies is present. M ATERI A L A N D M ETHODS Using field data, we parameterised 35 population growth models, one for each of the most abundant species (Table 1); the fitted parameters were the intrinsic popu- lation growth rates and interaction (competition and fa- cilitation) coefficients between each species pair along a hydric gradient. With these models, we simulated the abundance of each species while interacting with each other over the gradient. We then determined whether species clustering occurred in the simulated community, either considering all interactions, none of them or only competition. The interaction strength within and among species clusters was assessed, and the effects of the joint presence of the species in a cluster had on its own mem- bers were measured. Study site and data collection Our data come from a semi- arid grassland in Concepción Buenavista, Oaxaca, southern Mexico, at an elevation of 2275 m a.s.l. Annual precipitation is 530.3 mm and mean temperature is 16.3°C. Over 600 species occur in this grassland, about half of which are annual (Martorell et al., 2022). Up to 25 species and 135 individuals may be found in a single 0.1 × 0.1 m square, which is likely the scale over which interactions occur between the tiny plants (<5  cm tall; Martorell & Freckleton,  2014; Martorell et al.,  2017). The site has thin soils rarely deeper than 30 cm, with depth changing rapidly over distances of a few centimetres (Martorell et al.,  2015; Villarreal- Barajas & Martorell,  2009). The soil depth gradient is related to water availability, where thin soils are drier than deeper ones. Soil with a 2 cm depth has a water potential three times more negative than 30 cm soil (Martorell & Martínez- López,  2014; Villarreal- Barajas & Martorell, 2009). Other variables, such as nu- trient availability may concomitantly change with soil depth but they have not been analysed. The soil depth gradient is related to niche differentiation in this com- munity (Martínez- Blancas & Martorell, 2020; Martorell et al., 2015). The data come from 25 0.5 ha sites, with the whole depth gradient occurring in most. Eight 1 × 1 m quad- rats were randomly chosen within each site, and in each quadrat we randomly selected 20 0.1 × 0.1 m squares. We quantified the abundance of all vascular plants within the same squares yearly between 2001 and 2017 (some sites were only recorded since 2008 or 2012) towards the end of the rainy season in October. The average number of observations per species was 3031 ranging from 15,676 to 326. Presence/absence data instead of abundance was recorded for two grass species (Bouteloua chondrosiodes and Hilaria cenchroides) because they are cespitose, making it difficult to distinguish individuals. Soil depth was measured at each 0.1 × 0.1 m square by inserting a metal rod into the ground until it reached the bedrock. Population growth model To determine how the species' populations change as they interact with each other over the gradient, we used a discrete- time model that describes the annual changes in their density. Considering that interactions have posi- tive and negative effects simultaneously (Brooker & Callaghan, 1998; Callaway & Walker, 1997), the popula- tion growth model for each species is a modified version of the Beverton– Holt model: where N j,t is the observed abundance of species j in year t in a 0.1 × 0.1 m square, the subindex i identifies the focal species, D i is the proportion of individuals that die each (1)Ni,t+1= ( 1−Di ) Ni,t+Ni,t Bi,t ∏ j≠i ( 1+Nj,t )! ij 1+ ∑ j"ijNj,t 2654 | DRIVERS OF SPECIES CLUSTERING year and was fixed to one for annual species but allowed to vary for perennials, B i,t is the number of new individ- uals of the focal species produced each year t in the ab- sence of interactions, and α ij and β ij are measures of the competitive and facilitative effects, respectively, of the as- sociated species j on the focal species. Both were assumed to be constant through time. We assumed that no intraspe- cific facilitation occurs, as is common in plants (Adler et al., 2018; Bonanomi et al., 2010). Seed banks at the study site are small and seed survival is low (Zepeda, 2021), so no seed bank was considered. Because individual recruit- ment and interactions at our study site may change along the hydric gradient (Martínez- Blancas & Martorell, 2020), we modelled B i,t as a function of soil depth d as and α ij as (2)Bi,t = eai,t+bi,td+ci,td 2 , TA B L E 1 Species included in the study, their life cycle and their ability to persist in the community in three different scenarios: (a) with species interactions, (b) without species interactions and, (c) without facilitation. An, indicates annual life cycle, pe perennial lifecycle and Fe indicates facultative. X indicates scenarios where the species did not go extinct Species Family Life cycle Scenario (a) (b) (c) Aristida adscensionis L. Poaceae Fe X X X Aristida divaricate Humb. & Bonpl. ex Willd. Poaceae Pe X X X Bouteloua chondrosioides (Kunth) Benth. ex S. Watson Poaceae Pe X Bouteloua hirsute Lag. Poaceae Pe X X X Bouteloua scorpioides Lag. Poaceae Pe X X X Crusea diversifolia W.R. Anderson Rubiaceae An Cyperus seslerioides Kunth Cyperaceae Pe X Dalea humilis G. Don Fabaceae An X Euphorbia mendezii Boiss. Euphorbiaceae An X X X Evolvulus sericeus Sw. Convolvulaceae Pe X X X Florestina purpurea (Brandegee) Rydb. Asteraceae An X X X Heliotropium foliosissimum J.F. Macbr. Boraginaceae Pe Helianthemum glomeratum (Lag.) Lag. Cistaceae Pe Heterosperma pinnatum Cav. Asteraceae An X X Hilaria cenchroides Kunth Poaceae Pe X X Microchloa kunthii Desv. Poaceae Pe X X Milla biflora Cav. Asparagaceae Pe X X X Muhlenbergia peruviana (P. Beauv.) Steud. Poaceae An X X X Muhlenbergia phalaroides (Kunth) P.M. Peterson Poaceae Pe X X X Muhlenbergia rigida (Kunth) Kunth Poaceae Pe X X X Oxalis lunulate Zucc. Oxalidaceae Pe X X X Phemeranthus oligospermus (Brandegee) Ocampo Montiaceae Pe X X X Plantago nivea Kunth Plantaginaceae Fe X X X Portulaca pilosa L. Portulacaceae Fe X X X Richardia tricocca (Torr. & A. Grey) Standl. Rubiaceae Pe X X X Sanvitalia procumbens Lam. Asteraceae An X X X Sedum oteroi Moran Crassulaceae Pe X X X Sporobolus tenuissimus (Mart. ex Schrank) Kuntze Poaceae Pe X X X Stenandrium dulce (Cav.) Nees Acanthaceae Pe X X X Stevia ephemera Grashoff Asteraceae An X Tagetes micrantha Cav. Asteraceae An X X Thymophylla aurantiaca (Brandegee) Rydb. Asteraceae An X X X Tridax coronopifolia (Kunth) Hemsl. Asteraceae Fe Tripogandra purpurascens Commelinaceae An X Tripogon spicatus (S. Schauer) Handlos Poaceae Pe X | 2655MARTÍNEZ- BLANCAS ET AL. Facilitation was not explicitly modelled as a function of soil depth. See Appendix S1 for details on the ratio- nale of model. Modelling the interactions between all possible pairs of species in our study would have been imposible given the complexity of the model. Thus, we performed sparse modelling, in which we only estimated interactions be- tween species pairs for which there was sufficient sta- tistical evidence, which were selected using Bayesian Variable Selection (BVS; Garcia- Donato & Forte, 2018; Appendix  S2). The remaining species were added to- gether and allowed to interact with the focal species as a single ‘multispecies’ for which interaction parameters were estimated. This was observed to improve the mod- el's AIC (Appendix S1). All the parameters in equations (1– 3) were estimated through maximum likelihood using non- linear regres- sion of the observed N i,t+1 on the abundances of all spe- cies selected from BVS using package TMB (Kristensen et al.,  2016) for R (R Core Team,  2021). Because the movement of seeds between experimental units can bias parameter estimation, we explicitly accounted for dis- persal during model fitting (Appendices S1 and S3). Community simulations and cluster identification To determine whether there are species clusters and their underlying causes, we simulated the abundance of all species in our study for different soil depths (from 3 to 28 cm every 0.1 cm) using the fitted models (Equation 1) for all species simultaneously. The system was iterated until species abundances stabilised (Appendix S4). We conducted three sets of simulations. In the first, all in- teractions (including those with our ‘multispecies’) were allowed to occur. To test whether direct facilitation in- duced species to cluster over the depth gradient, we per- formed the same simulations but setting all β to zero. To evaluate whether interactions are the cause of clustering, we set all interspecific α and β to zero (Appendix S4). We used the k- means algorithm implemented by D'Andrea et al.  (2019, 2020) to find species clusters. K- means assigns species to a predetermined number of clusters k so that within- cluster variation is minimised. We tried all possible values of k between 1 and 8, and calculated the gap statistic (G k ), which is the difference in within- cluster variation between the observed and null communities (D'Andrea et al., 2020). Positive G k val- ues indicate tighter than expected clustering in the data. Selecting the number of clusters that maximised G k re- sulted in very large numbers of clusters in our dataset, so we selected the smallest value of k for which G k was significant. The k- means algorithm was originally implemented for species traits, which are usually described with a single value per species. However, we were interested in clustering of hydrological niches, which are reflected in species abundance over the depth gradient. Thus, rather than a single value for each species, we need to consider the abundance at each depth. Because k- means minimises the differences between species in a cluster, we calculated the differences between each species pair comparing their abundances over the depth gradient using Hellinger distances, which measure how different two curves are. We then used polar ordination to reduce the Hellinger's distance between all species pairs to a two- dimensional Euclidean space (which captured 67% of the variation, see Appendix S4). This allowed the con- struction of an appropriate null model for the k- means al- gorithm on two dimensions. The k- means algorithm was applied to the community simulations with and without all interactions, and without facilitation (Appendix S4). We also applied k- means to our original field data and compared the resulting clusters to those obtained using our model using the Sørensen similarity index. We used a null model to assess the significance of these compari- sons (Appendix S4). Strength of interactions within and between clusters To test whether clustering was consistent with emergent neutrality, hidden niches or facilitation, we compared the average strength of interactions between and within clusters using the following: where !within is the average of all interspecific competition coefficients (α) between species that belong to the same cluster, and !between is the respective average consider- ing species pairs that do not belong to the same cluster. Positive !d values indicate that within- cluster competition is stronger than between- cluster competition, as expected from the emergent neutrality hypothesis. Negative values indicate that competition is stronger between rather than within clusters and suggests hidden niches. Because com- petition changes with soil depth in our model, we calcu- lated !d for four soil depths representative of each cluster. To compare within and between cluster facilitation, we cal- culated ! using the same procedure as above. Facilitation does not change with soil depth in our model, so we ob- tained a single value of!. To assess significance, we com- pared observed !d and ! with null values obtained from 1000 communities generated by shuffling species interac- tion coefficients while maintaining cluster size and num- ber constant. To measure the indirect effects that species in the same cluster had on each other, we simulated the (3)!ij = egi+hid . (4)!d = !within − !between, 2656 | DRIVERS OF SPECIES CLUSTERING abundance of each species without the presence of the other species within its cluster over the depth gra- dient. These abundances were compared with our original simulation with all species and all interac- tions using the log response ratio (Choler et al., 2001) LRR i,d  = ln(NC i,d / NWC i,d ); where NC is the equilibrial abundance of the focal species i at soil depth d simu- lated with all the species in the community, and NWC is the equilibrial abundance of species i at soil depth d simulated without the other species in its cluster. Positive LRR values indicate facilitation, as expected when there are species alliances. If, in addition, ! ≤ 0, it can be assumed that positive LRR values are due to indirect facilitation. RESU LTS Species clustering In nature and in modelled scenarios, most species per- sisted along the entire soil depth gradient (Figure 2). Despite this, with field data we found significant clus- tering with k = 7 (Figure 3), indicating that there are seven clusters. However, nearly half of them contained one or two species. A better classification (Figure 2) was obtained using five clusters, which was marginally significant (p  =  0.062, Figure  3). We found four sig- nificant clusters in all simulated scenarios (p = 0.002 with all interactions, p = 0.011 without facilitation, and P  =  0.048 with no interactions, Figure  3). However, when we eliminated interactions, there was one large cluster comprising most species, and three clusters with just one or two species, implying that no proper clus- tering occurs without interactions. When we included interactions, some species were driven into shallow soils, forming a new cluster (Figure 4), and two species moved to deeper soils, forming a new cluster with some species that already occurred there. The inclusion of interactions resulted in a greater spread of species over the gradient (standard deviation of species distribu- tion's centroids = 4.7) than without interactions (3.8). Species density decreased on average 77% with interac- tions compared to when interactions were eliminated (Figure 2). The clusters obtained with all interactions and without facilitation comprised approximately the same species (Appendix  S5). From now on, we will refer to each of the clusters that were observed when all interactions were considered in the model as clusters 1 to 4 in order of increasing soil depth such that clusters 1 and 2 are the shallow soil clusters and 3 and 4 are the deep soil clusters. F I G U R E 2 Simulated species abundances from (a) simulations with all interactions (b) simulations without facilitation and (c) without interactions. Same- colour lines represent abundances of species within the same cluster. In (d), we show the species clusters obtained from real data (identified with letters in red and black lines) and simulations with all interactions (coloured ovals), and their approximate location over the soil depth gradient. The first two letters of each species and its genus are shown. (a) (b) (c) (d) | 2657MARTÍNEZ- BLANCAS ET AL. When all interactions were included in the model, species distribution over the gradient resembled that of species in nature: the soil depth where species attained their maximum density in the simulation was correlated with that observed in the field by Pedraza and Martorell (2019) (r = 0.74, p < 0.001). When we eliminated all inter- actions, this correlation became weaker but remained significant (r = 0.44, p = 0.029). Clusters obtained from our simulation also resembled nature. We found that three of them were more similar than expected by chance to their counterparts obtained from field data (the exception being cluster 3), with an averaged similarity of 63% (Table S5.1). Nearly all differences between clusters obtained from real and simulated data can be attributed to a single cluster of species that in nature tend to occur in deep soil (Cluster C in Figure 2). In our model, very few species displayed this behaviour (Figure 2), and thus the species in the deep soil cluster were assigned to other clusters. In our simulations with all interactions, 27 spe- cies persisted in at least a portion of the depth gradient while eight species went extinct. When we eliminated facilitation from our simulations, only 24 species per- sisted (Table 1), and when we eliminated all interactions a slightly different set of 27 species occurred. Direct and indirect interactions between and within clusters There were 20% (241) meaningful (≠ 0) pairwise interspe- cific interactions between our 35 species. Intraspecific competition was substantially stronger than interspecific competition (Appendix S4). When competition changed with depth, it frequently was stronger in shallow soils (h i was negative in 58% of species pairs; Appendix S4). The interactions with the ‘multispecies’ were negligible. Competition was stronger between clusters rather than within them at all soil depths examined but was only significant in shallow soils (clusters 1 and 2 with p values of 0.040 and 0.048) and not in deep soils (clusters 3 and 4 with p values of 0.098 and 0.519). This was tested for the average soil depth where species in the cluster occur pondered by their abundance over the gradient: F I G U R E 3 Gap statistic G against potential number of clusters under three scenarios and field data: (a) with all interactions present, (b) without facilitation (c) without species interactions and (d) field data. Red line indicates 0.95 quantile of the expected gap statistic obtained from 1000 null sets. The black line is our observed gap statistic G k . If the black line is above the red one, data support k clustering. (a) (b) (c) (d) F I G U R E 4 Species position along the soil depth gradient with and without species interactions. Same- colour circles represent species in the same cluster. Species position was calculated as the average soil depth where the species is present weighted by its abundance. 2658 | DRIVERS OF SPECIES CLUSTERING 9.9 cm (cluster 1), 15.48 cm (cluster 2), 19.9 cm (cluster 3) and 27.18 cm (cluster 4). Average values of the log response ratio (LLR) were positive within every cluster except for the third one (Figure 5). Given the fact that within- cluster direct fa- cilitation was not particularly strong (! ≈ 0), such posi- tive LRR values suggest positive indirect effects from the other species in the cluster. Almost all LRR for cluster 2 were positive values. Clusters 1 and 4 had several neg- ative LRR values, but the majority were positive along the portion of the gradient where their species were most abundant. Negative values of LRR were observed in cluster 3 were near zero along the soil depth interval where its species are abundant. DISCUSSION We found that species cluster along the hydric gradient with or without species interactions. However, clusters were poorly structured when interspecific interactions were excluded from the model. This suggests that spe- cies interactions drive clustering. Interspecific compe- tition was weaker within clusters rather than between them, suggesting that differences in unmeasured niche dimensions stabilise coexistence within clusters. Direct facilitation does not seem to drive clustering, as clus- ters persisted after removing it and the intensities of within- and between- cluster facilitation were similar. In contrast, indirect facilitation seems to promote species alliances in three clusters whose members increased each other's abundances when they are present together. The role of interactions in species clustering Removing species interactions from our model yielded poor clustering, with all but one cluster containing one or two species. This supports previous arguments that clustering in communities is driven by species interac- tions (D'Andrea et al., 2020; Segura et al., 2011; Vergnon et al.,  2009). Direct facilitation did not play a funda- mental role in such clustering. Even though facilitation caused niches to shift, as previously suggested (Bulleri et al., 2016), when it was eliminated clusters only changed slightly. Thus, the emergence of groups of species could not be attributed to facilitated species cooccurring with their benefactors over the depth gradient (Filazzola et al., 2018). In contrast, competition was a major driver of clustering. As expected if interspecific competition causes spatial segregation (Shigesada et al., 1979), species spread out over the soil depth gradient when interactions were included in the model (Figure 4). This is in line with previous findings that competition reduces the over- lap of hydrological niches at the study site (Martorell F I G U R E 5 Log- response ratio (LRR) for each clusters 1 through 4 along the soil depth gradient. The black line represents the cluster median while the blue ribbon represents the interquartile range. The dark blue bar at the bottom is the cluster position calculated as the interval comprising the average soil depth where each species in the cluster is present weighted by its abundance | 2659MARTÍNEZ- BLANCAS ET AL. et al.,  2015). The species that competition forced into shallow soil might persist there because they tolerate dry conditions where their competitors are absent. These tol- erant species store water in their tissues (Phemeranthus and Portulaca) or have annual life cycles, so they need not to endure the driest months (Heterosperma). Species that were driven into deep soil may be escaping from their competitors by tolerating higher quantities of path- ogens (negative plant– soil feedbacks) present in deeper soil (Martorell et al., 2021), but evidence for this is yet scarce. As a result of greater species spread, species at the extremes of the gradient formed clusters that were clearly differentiated from the rest. Mechanisms of species clustering There is currently a debate on the mechanisms that drive species clustering in nature (Barabás et al., 2013; Vergnon et al., 2013). While earlier studies favour emer- gent neutrality (Sakavara et al., 2018; Scheffer et al., 2015; Scheffer & van Nes,  2006a; Segura et al.,  2011, 2013; Vergnon et al., 2009, 2012), some authors have recently advocated for hidden niches (Barabás et al.,  2013; D'Andrea et al., 2020; Segura et al., 2011). This debate is mostly based on the outcomes of computational models or theoretical arguments. Besides the present work there is, to the best of our knowledge, only one more study pro- viding empirical data on clustering mechanisms (Segura et al.,  2011), and both support different mechanisms for clustering. While our results suggest hidden niches, emergent neutrality seems to drive clustering in a natural phytoplankton community because most clusters occur on portions of a gradient where species' fitness is high and equal. Nevertheless, hidden niches may still stabilise coexistence within some clusters in this aquatic system (Segura et al., 2011). Our results point to hidden niches rather than emer- gent neutrality. An important issue here is that stron- ger intra- than interspecific competition is a signature of niche differentiation (Barabás et al., 2013; Broekman et al., 2019). This is what we found, so our results are inconsistent with a strictly neutral scenario where niches are identical. Furthermore, if, as required by the neutral model, species in the same cluster have similar niches, the difference between intra- and interspecific competi- tion would be smaller when considering species in the same cluster rather than in different ones. This would be observed if interspecific competition within clusters were stronger— and closer to the strong intraspecific competition— than between clusters. Nevertheless, the opposite was observed, indicating lower niche overlap within clusters, as expected if there were differentiation over unmeasured niche dimensions. Likely candidates for such unmeasured dimensions are seed predation and plant– soil feedbacks. At the study site, the probability that an ant removes a seed increases with seed density. (García- Meza et al., 2021). This may enhance negative density dependence and thus increase the difference between intra- and interspecific regulation. Strong plant– soil feedbacks can mediate niche differentiation. They are important in this grass- land, especially in shallow soils (Martorell et al., 2021) where the greatest values of !d were observed. The study site is also nutrient poor (particularly in phosphorus; Martorell et al., 2021) such that variation in nutrient lim- itation along the soil depth gradient may also provide ad- ditional axes for niche differentiation. Light availability gradients (which may be related to clustering, D'Andrea et al., 2020) do not seem to play a role in niche differenti- ation. Plants are unlikely to shade each other because of their tiny size, and thus do not seem to be light- limited. Temporal fluctuations do not contribute to niche differ- entiation either (Zepeda & Martorell, 2019). In contrast with the general trend for interactions to reduce strongly population densities (77% on average), the abundances of the species in clusters 1, 2 and 4 were increased by the joint occurrence of the other species in the same cluster. This evidence, altogether with clus- ter persistence after removing direct facilitation and the observed ! ≈ 0, suggests that clustering arises from species alliances through indirect facilitation. Within trophic levels, indirect facilitation occurs when a spe- cies facilitates another one by suppressing its competi- tor (Levine, 1999; Wilson, 1986). This scenario is more likely to occur when species require different resources, so that competition between the benefactor and the fa- cilitated species is weaker than the indirect facilitative effect (Davidson, 1985; Levine, 1999). This was our case: indirect facilitation mediated by the suppression of spe- cies in other clusters (because the !between values are large) could predominate over the weak within- cluster com- petition (small !within). In fact, whenever hidden niches drive clustering, species alliances may arise due to the larger between than within cluster competition. Thus, hidden niches and species alliances may jointly reinforce species clustering. Intransitive competition could also generate indirect facilitation (Figure S6.1), but we found no evidence for intransitivity (Appendix S6). Model considerations There are a few ways in which our data or model could have biased our results. Studies, like ours, based on ob- servational rather than experimental data, tend to find weaker inter- than intra- specific competition (Adler et al., 2018; Tuck et al., 2018), possibly biasing the results towards greater niche differentiation. However, obser- vational, rather than experimental data, involve larger timescales that allow processes that affect species in- teractions, such as plant– soil feedbacks (important in this grassland, Martorell et al., 2021), to play out (Adler et al., 2018). Unmeasured environmental variables may 2660 | DRIVERS OF SPECIES CLUSTERING also bias interaction estimates. For instance, if a spe- cies is abundant under certain conditions where others are absent, we may wrongly infer competition. We con- trolled for unmeasured environmental factors by includ- ing quadrats as random effects in our model. Data spanning many years and even decades, may be essential to capture rare demographic events on which population persistence may depend (Félix- Burruel et al., 2021). We used the birth rates (B) from only 7 years for our simulations. Unobserved variation over such pe- riod may explain why a few common species at the study site oddly became extinct in our simulations. Other spe- cies have been declining, so the model seems to correctly predict their extinction. Despite these considerations, our model reproduced distribution of species maximal abundance along the soil depth gradient. Also, the tendency for competition to increase in shallow soils is in agreement with exper- imental studies in this grassland (Martínez- Blancas & Martorell, 2020; Villarreal- Barajas & Martorell, 2009), and in line with the hypothesis that competition in- creases with stress related to lack of resources (Maestre et al.,  2009). With the exception of cluster 3, all the clusters predicted by our simulations were similar to the ones observed in nature. Nevertheless, some spe- cies that achieved their maximal observed density in deeps soils were incorrectly modelled to have their maxima elsewhere. Thus, they misclassified in clusters 1– 3 instead of cluster 4, which was the one occurring in deep soil. The density of some of these species had a multimodal distribution over the depth gradient. In several of them, the maximum predicted by our sim- ulations corresponded to a secondary maximum ob- served in the real data, suggesting that the simulated clusters reflect the true distribution of the misclassi- fied species to some extent. It is interesting that only cluster 3, which had no counterpart in the real data, lacked evidence for species alliances. Concluding remarks Our study highlights the need for quantitative stud- ies with empirical data to disentangle the effects of species interactions on the structure and composi- tion of communities in nature. As models suggest, interspecific interactions were key for species cluster- ing. Differentiation of the niche over unmeasured di- mensions was important for clustering over the stress gradient because it allowed for the stabilisation of co- existence, although we can only speculate that such stabilisation was strong enough to overcome fitness differences (Broekman et al.,  2019). Formerly over- looked indirect interactions may have an important role in species clustering, and may easily arise from within- cluster stabilising mechanisms such as hidden niches. Under such scenario, clusters are in fact species alliances that indirectly facilitate each other by sup- pressing common competitors in other clusters. AU T HOR CON T R I BU T IONS AM- B and CM designed the study and performed the research, AM- B, CM and IXB collected data, performed the modelling work and analysed the data. AM- B and CM wrote the first draft of the manuscript, and all au- thors contributed substantially to revisions. ACK NOW LEDGEM EN TS LFVV Boullosa, A Cervantes, D García- Meza, F Herce, G Martínez, M Martínez, F Pedraza, C Vázquez, V Zepeda assisted us in the field. MA Romero offered com- putational support. EJ González, K Boege and J Meave discussed with us some ideas. CONACyT provided sti- pend for AMB during her PhD studies in the Posgrado de Ciencias Biológicas, UNAM. PAPIIT- UNAM (IN225511, IN220514 and IN212618) funded the research. Agradecemos a la comunidad de Concepción Buenavista por permitirnos realizar este estudio en su territorio. F U N DI NG I N FOR M AT ION Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México, Grant/ Award Number: PAPIIT IN- 212618 PEER R EV I EW The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ele.14122. DATA AVA I LA BI LI T Y STAT EM EN T The data and code used for analyses and simulations are available at GitHub https://zenodo.org/badge/ lates tdoi/46592 6013 and figshare: https://doi.org/10.6084/ m9.figsh are.20468 535.v3 ORCI D Alejandra Martínez- Blancas  https://orcid. org/0000-0003-4669-0257 Carlos Martorell  https://orcid. org/0000-0002-0758-3953 R E F ER E NC E S Adler, P.B., HilleRislambers, J. & Levine, J.M. (2007) A niche for neu- trality. Ecology Letters, 10, 95– 104. Adler, P.B., Smull, D., Beard, K.H., Choi, R.T., Furniss, T., Kulmatiski, A. et al. (2018) Competition and coexistence in plant communi- ties: intraspecific competition is stronger than interspecific com- petition. Ecology Letters, 21, 1319– 1329. Barabás, G., D'Andrea, R., Rael, R., Meszéna, G. & Ostling, A. (2013) Emergent neutrality or hidden niches? Oikos, 122, 1565– 1572. Bonanomi, G., Incerti, G., Capodilupo, M. & Mazzoleni, S. (2010) Rare self- facilitation in terrestrial plants as compared to aquatic sessile organisms: empirical evidences and causal mechanisms. Community Ecology, 11, 148– 159. Bonsall, M.B., Jansen, V.A.A. & Hassell, M.P. (2004) Life history trade- offs assemble ecological guilds. Science (80- . )., 306, 111– 114. | 2661MARTÍNEZ- BLANCAS ET AL. Broekman, M.J.E., Muller- Landau, H.C., Visser, M.D., Jongejans, E., Wright, S.J. & de Kroon, H. (2019) Signs of stabilisation and sta- ble coexistence. Ecology Letters, 22, 1957– 1975. Brooker, R.W. & Callaghan, T.V. (1998) The balance between positive and negative plant interactions and its relationship to environ- mental gradients: a model. Oikos, 81, 196– 207. Bulleri, F., Bruno, J.F., Silliman, B.R. & Stachowicz, J.J. (2016) Facilitation and the niche: implications for coexistence, range shifts and ecosystem functioning. Functional Ecology, 30, 70– 78. Callaway, R.M. & Walker, L.R. (1997) Competition and facilitation: a synthetic approach to interactions in plant communities. Ecology, 78, 1958– 1965. Chesson, P. (2000) Mechanisms of maintenance of species diversity. Annual Review of Ecology and Systematics, 31, 343– 366. Choler, P., Michalet, R. & Callaway, R.M. (2001) Facilitation and competition on gradients in alpine plant communities. Ecology, 82, 3295– 3308. D'Andrea, R., Guittar, J., O'Dwyer, J.P., Figueroa, H., Wright, S.J., Condit, R. et al. (2020) Counting niches: abundance- by- trait pat- terns reveal niche partitioning in a neotropical forest. Ecology, 101, e03019. D'Andrea, R., Riolo, M. & Ostling, A.M. (2019) Generalizing clusters of similar species as a signature of coexistence under competi- tion. PLoS Computational Biology, 15, 1– 19. Davidson, D.W. (1985) An experimental study of diffuse competition in harvester ants. The American Naturalist, 125, 500– 506. Doebeli, M. & Dieckmann, U. (2003) Speciation along environmental gradients. Nature, 421, 259– 264. Félix- Burruel, R.E., Larios, E., González, E.J. & Búrquez, A. (2021) Episodic recruitment in the saguaro cactus is driven by mul- tidecadal periodicities. Ecology, 102, 1– 11. Filazzola, A., Sotomayor, D.A. & Lortie, C.J. (2018) Modelling the niche space of desert annuals needs to include positive interac- tions. Oikos, 127, 264– 273. Garcia- Donato, G. & Forte, A. (2018) Bayesian testing, variable se- lection and model averaging in linear models using R with BayesVarSel. R J., 10, 155– 174. García- Meza, D., Andresen, E., Ríos- Casanova, L. & Martorell, C. (2021) Seed density in monospecific and mixed patches affects individual and collective foraging in ants. Insectes Sociaux, 68, 81– 92. Herault, B. & Thoen, D. (2009) How habitat area, local and regional factors shape plant assemblages in isolated closed depressions. Acta Oecologica, 35, 385– 392. Kristensen, K., Nielsen, A., Berg, C.W., Skaug, H. & Bell, B. (2016) TMB: automatic differentiation and Laplace approximation. Journal of Statistical Software, 70, 1– 21. Levine, J.M. (1999) Indirect facilitation: evidence and predictions from a riparian community. Ecology, 80, 1762– 1769. Maestre, F.T., Callaway, R.M., Valladares, F. & Lortie, C.J. (2009) Refining the stress- gradient hypothesis for competition and facilitation in plant communities. Journal of Ecology, 97, 199– 205. Martínez- Blancas, A. & Martorell, C. (2020) Changes in niche differ- entiation and environmental filtering over a hydric stress gradi- ent. Journal of Plant Ecology, 13, 185– 194. Martorell, C., Almanza- Celis, C.A.I., Pérez- García, E.A., Sánchez- Ken, J.G., Perez- Garcia, E.A. & Sanchez- Ken, J.G. (2015) Co- existence in a species- rich grassland: competition, facilita- tion and niche structure over a soil depth gradient. Journal of Vegetation Science, 26, 674– 685. Martorell, C. & Freckleton, R.P. (2014) Testing the roles of competi- tion, facilitation and stochasticity on community structure in a species- rich assemblage. Journal of Ecology, 102, 74– 85. Martorell, C., García- Meza, D., Martínez- Blancas, A., Zepeda, V. & Vazquez- Ribera, C. (2022) Pastizales de la región Chocholteca: un despreciado récord mundial de diversidad vegetal. In: Cruz Angón, A., Nájero Cordero, K.C., Cruz Medina, J. & Melgarejo, E.D. (Eds.) La biodiversidad en Oaxaca estudio de esttao. Ciudad de México: Comisión Nacional para el Conocimiento y Uso de la Biodiversidad, pp. 290– 301. Martorell, C., Martínez- Blancas, A. & García- Meza, D. (2021) Plant– soil feedbacks depend on drought stress, functional group, and evolu- tionary relatedness in a semiarid grassland. Ecology, 0, e03499. Martorell, C. & Martínez- López, M. (2014) Informed dispersal in plants: heterosperma pinnatum (Asteraceae) adjusts its dispersal mode to escape from competition and water stress. Oikos, 123, 225– 231. Martorell, C., Zepeda, V., Martínez- Blancas, A., García- Meza, D. & Pedraza, F. (2017) A diversity world record in a grassland at Oaxaca, México. Botanical Sciences, 95, 1– 7. Pedraza, F. & Martorell, C. (2019) Allocating species in Grime’s strat- egy space: an alternative to trait-based approaches. Botanical Sciences, 97, 649– 660. R Core Team. (2021). R: a language and environment for statistical computing. Roelke, D.L. & Eldridge, P.M. (2008) Mixing of supersaturated as- semblages and the precipitous loss of species. The American Naturalist, 171, 162– 175. Sakavara, A., Tsirtsis, G., Roelke, D.L., Mancy, R. & Spatharis, S. (2018) Lumpy species coexistence arises robustly in fluctuating resource environments, 115. Scheffer, M. & van Nes, E.H. (2006a) Self- organized similarity, the evolutionary emergence of groups of similar species. Proceedings of the National Academy of Sciences, 103, 6230– 6235. Scheffer, M. & Van Nes, E.H. (2006b) Species of groups of similar species. Proceedings of the National Academy of Sciences, 103, 6230– 6235. Scheffer, M., Vergnon, R., Van Nes, E.H., Cuppen, J.G.M., Peeters, E.T.H.M., Leijs, R. et al. (2015) The evolution of function- ally redundant species; Evidence from beetles. PLoS One, 10, e0137974. Scranton, K. & Vasseur, D.A. (2016) Coexistence and emergent neu- trality generate synchrony among competitors in fluctuating en- vironments. Theoretical Ecology, 9, 353– 363. Segura, A.M., Calliari, D., Kruk, C., Conde, D., Bonilla, S. & Fort, H. (2011) Emergent neutrality drives phytoplankton species coex- istence. Proceedings of the Royal Society B: Biological Sciences, 278, 2355– 2361. Segura, A.M., Kruk, C., Calliari, D., García- Rodriguez, F., Conde, D., Widdicombe, C.E. et al. (2013) Competition drives clumpy species coexistence in estuarine phytoplankton. Scientific Reports, 3, 1– 6. Shigesada, N., Kawasaki, K. & Teramoto, E. (1979) Spatial segregation of interacting species. Journal of Theoretical Biology, 79, 83– 99. Silvertown, J. (2004) Plant coexistence and the niche. Trends in Ecology & Evolution, 19, 605– 611. Thibault, K.M., White, E.P., Hurlbert, A.H. & Ernest, S.K.M. (2011) Multimodality in the individual size distributions of bird com- munities. Global Ecology and Biogeography, 20, 145– 153. Tuck, S.L., Porter, J., Rees, M. & Turnbull, L.A. (2018) Strong re- sponses from weakly interacting species. Ecology Letters, 22, 1845– 1852. Vergnon, R., Dulvy, N.K. & Freckleton, R.P. (2009) Niches versus neutrality: uncovering the drivers of diversity in a species- rich community. Ecology Letters, 12, 1079– 1090. Vergnon, R., van Nes, E.H. & Scheffer, M. (2012) Emergent neutral- ity leads to multimodal species abundance distributions. Nature Communications, 3, 1– 6. Vergnon, R., van Nes, E.H. & Scheffer, M. (2013). Interpretation and predictions of the emergent neutrality model: a reply to Barabás et al. Oikos, 122, 1573– 1575. Villarreal- Barajas, T. & Martorell, C. (2009) Species- specific distur- bance tolerance, competition and positive interactions along an anthropogenic disturbance gradient. Journal of Vegetation Science, 20, 1027– 1040. 2662 | DRIVERS OF SPECIES CLUSTERING Wilson, D.S. (1983) The effect of population structure on the evolu- tion of mutualism: a field test involving burying beetles and their phoretic mites. The American Naturalist, 121, 851– 870. Wilson, D.S. (1986) Adaptive indirect effects. In: Diamond, J. & Case, T.J. (Eds.) Community ecology. New York: Harper & Row, pp. 437– 444. Wilson, D.S. & Knollenberg, W.G. (1987) Adaptive indirect effects: the fitness of burying beetles with and without their phoretic mites. Evolutionary Ecology, 1, 139– 159. Worthen, W.B. & Moore, J.L. (1991) Higher- order interactions and indirect effects: a resolution using laboratory drosophila com- munities. The American Naturalist, 138, 1092– 1104. Zepeda, V. (2021). Variabilidad climática y mantenimiento de la diver- sidad vegetal en un pastizal rico en especies. PhD dissertation. Mexico City: Universidad Nacional Autónoma de México. Zepeda, V. & Martorell, C. (2019) Fluctuation- independent niche differentiation and relative non- linearity drive coexistence in a species- rich grassland. Ecology, 100, e02726. SU PPORT I NG I N FOR M AT ION Additional supporting information can be found online in the Supporting Information section at the end of this article. How to cite this article: Martínez- Blancas, A., Belaustegui, I.X. & Martorell, C. (2022) Species alliances and hidden niche dimensions drive species clustering along a hydric gradient in a semiarid grassland. Ecology Letters, 25, 2651–2662. Available from: https://doi.org/10.1111/ele.14122 62 Appendix 1. Model and Model fitting Population growth model Our population growth model includes two components. The first one (1-Di , where Di is the probability that an individual dies; the subindex i indicates the focal species) corresponds to the individuals that survive in a given 0.1 ´ 0.1 m square, while the second one represents the individuals that arrive and get established at that square. Arriving individuals may be from the same or nearby squares. Individuals that survive from one time step to the next are generally large and not susceptible to interspecific interactions (Martorell & Freckleton 2014), so Di does not depend on interactions in our model. However, we assume that fecundity is affected by interactions and environmental stress. Soil depth is related to water availability at our study site. Deep soils have greater water availability than shallow ones, which are more stressful, and individuals grow larger in deeper soils. Given the relationship between size and reproduction (Lacey 1986), in our model the intrinsic birth rate Bi,t(d) depends on soil depth d. Note that the intrinsic birth rate is also year-specific (t) as population growth has been observed to change with time in previous studies at the study site (Martorell & Freckleton 2014). Birth rate responds to species interactions, which are both competitive (depending on a function C(Nt), where Nt is a vector containing the number of individuals of all species in the square at time t) and facilitative (represented by the function F(Nt)) at our study site (Villarreal-Barajas & Martorell 2009; Martorell & Freckleton 2014). In our model, facilitation is only interspecific as intraspecific facilitation is not common in plants (Bonanomi et al. 2010; Adler et al. 2018). Note that the model does not take the seed 63 bank into account because it is small and short lived in our grassland (Zepeda 2021). Thus, the general form of the model is: !!,#$% = #1 − &! + (!,#(*) &!((") *!((",+) , !!,# (S1.1) The model was parameterized using data collected at 25 0.5 ha sites at the semi-arid grassland. Eight 1 × 1 m quadrats were randomly chosen within each site, and in each one we randomly selected 20 0.1 × 0.1 m squares. Of the 25 sites, eight of them were followed for 15 years (14 transitions), eight more of them were followed for 7 years (6 transitions) and nine of them were followed for 6 years (5 transitions). Population growth and competition as a function of soil depth. We tested different models that incorporate Bi,t as a function of soil depth (d) and because we had no prior reasons to suspect any particular form of the curve. We tested some models in which birth rates changed monotonically with soil depth, Bi,t = ai,t +bi,t d (S1.2) (!,# = .,!,"$-!,"+ (S1.3) (!,# = .,!,"$-!,"./+ (S1.4) (!,# = ,!,"+ -!,"$+ (S1.5) However, it is likely that birth rate has an optimum at intermediate depths, as it is common in the response of species to environmental factors (Fraaije et al. 2015). Thus, also tested some models with quadratic terms: Bi,t = ai,t +bi,t d + ci,t d2 (S1.6) (!,# = .,!,"$-!,"+$0!,"+$ (S1.7) 64 (!,# = .,!,"$-!,"./+$0!,"./$+ (S1.8) where a, b and c are parameters that relate the soil depth to recruitment in species i and year t. Equations S1.4-5 and S1.7-8 have the property that they only produce positive values, as it should be the case for numbers of recruits. Equation S1.5 reaches an upper asymptote. We also tested models in which some parameters did not change annually, but they provided poorer fits than the ones shown. For the interactions, we used the function reported in Martorell & Freckleton (2014), which have proven to be appropriate for the study system: /! = 1 +01!,1(*)!1,# 2 13% (S1.9) 2! =341 + !1,#54!,% 2 15! (S1.10) where a is competition and b is facilitation. Note that if 0 < a and 0 ≤ c < 1, the competition coefficient overwhelms facilitation at high densities of the associated species, representing the inevitable outcome that when the associate species is very abundant it will displace the focal species as it monopolizes all the available space. This also prevent the occurrence of an “orgy of mutual benefits” (May & McLean 2007), where a couple of species help each other to grow to infinity. To make 1!,1 dependent of soil depth we used the equation 1!,1 = .6!,%$7!,%+ (S1.11) 65 Where, g and h are parameters to be estimated from the data. The exponential grants that alphas are positive, so that they correspond to competition and facilitation. Facilitation is not explicitly modeled as a function of soil depth even though it would also be expected to change with hydric stress (Brooker & Callaghan 1998; Pugnaire et al. 2004). However, because B, a and b are statistically correlated, trying to estimate the three of them as functions of d caused model identifiability problems. We chose to allow competition to change with soil depth and keep facilitation constant, which allows the balance between both interactions to vary over the hydric gradient. We assessed the fit of all the possible models for Bi,t by plugging equations S1.2 through 8 into our population model (S1.1). In these preliminary fits we used data for seven annual and perennial species (Table S1.1) for which we had large amounts of data and and we assumed that competition did not change with depth. The models were fitted assuming a negative binomial distribution, which describes appropriately our data (Martorell & Freckleton 2014; Zepeda & Martorell 2019). A function evaluating the likelihood of each model was coded in C++ and integrated into the R environment (R Development Core Team, 2010) using TMB (Kristensen et al. 2016). We used nlminb (R Development Core Team, 2010) to find the solution of maximum likelihood given the function coded in C++ and a set of starting values for the parameters. Each model was fitted 100 times using different initial values for the parameters to reduce the chances of getting estimates for local minima. The initial values for the parameters were obtained using Latin hypercube sampling with the function geneticLHS in the package lhs (Carnell 2019). From these model fits, the model with the largest likelihood was selected (Bolker 2007). We used the Akaike information criterion (AIC) to compare between models 66 (Table S1.1). These models were fitted for a subset of species in the community. Each of the species were modelled as the focal species i with the rest of the species j in the set for interspecific interactions. Table 1. DAIC values for the population growth models. Each model includes a different way of modelling new individuals as a function of soil depth. The best model or models for each species are in bold. Eq. (1.2) Eq. (1.3) Eq. (1.4) Eq. (1.5) Eq. (1.6) Eq (1.7) Eq. (1.8) Aristida adscensionis 152.46 20.22 94.45 509.69 181.68 0 156.46 Florestina purpurea 78.16 39.36 40.87 383.01 112.69 0 82.16 Plantago nivea 367.2 123.18 101.06 1828.18 910.9 0 371.2 Microchloa kunthii 179.82 97.52 81.38 2154.18 447.26 0 184.46 Stevia ephemera 395.33 47.18 158.98 674.86 20971.4 0 399.32 Tagetes micrantha 127.44 0 27.74 638.7 372.1 1.74 131.44 Thymophyla aurantiaca 478.48 155.94 204.42 741.78 579.42 0 482.48 Because Eq. 1.7 had DAIC = 0 in the majority of the species we chose used this equation to model every species’ population. Thus, the function used to model population growth was: !!,#$% = 71 − &! + .,!$-!+$0!+$ ∏ 9%$:%,"; &!,%' %(! %$∑ =)!,%*+!,%,:%," ' %-. 8!!,# (S1.12) For the two species that had presence-absence data, this equation had to be slightly modified. Because we are dealing with such data, Ni,t is a number between zero and one, 67 that, as in the models for counts correspond to the average density in a 0.1 ´ 0.1 m square. From the equation of the mathematical expectation, the average number of occupied squares in time t + 1 is equal to the probability that a square is occupied at time t (which is equal to Ni,t) multiplied by the probability that it remains occupied plus the probability that a square is not occupied multiplied by the probability that it is colonized. Assuming that the number of new recruits is, as before, (!,#(*) &!((") *!((",+) , and that it follows a Poisson distribution, the probability that none of them gets established is .>?!,"(+) /!(1")3!(1", ). Thus !!,#$% = !!,# 91 − &!.>?!,"(+) /!(1") 3!(1", ): + 41 − !!,#5 91 − .>?!,"(+) /!(1") 3!(1", ): (S1.13) In the first term of the sum, a Di is required because we are considering the scenario where the individual in the occupied square dies and it is not replaced, while in the second term the square is already unoccupied at time t. This equation can be rearranged to !!,#$% = 1 − 41 − (1 − &!)!!,#5.>?!,"(+) /!(1") 3!(1", ) (S1.14) The functions for the birth rate, competition and facilitation are S1.7, S1.9 and S1.10, just as in the case for species with count data. Model parametrization Because of the difficulty to fit the entire set of parameters in the model, we fit it in two phases (Fournier et al. 2012). Because dispersal between squares can bias the competition estimates (Freckleton & Watkinson 2001), we had to consider it when fitting the model 68 (Martorell & Freckleton 2014; Zepeda & Martorell 2019), and the first phase was used mainly to obtain the parameters for the dispersal kernel. Once these parameters were estimated, they were maintained fixed while estimating the rest of the parameters in the model during phase 2. Phase 1: The number of new individuals (Bi,t) and competition (1!,1) were held constant along the soil depth, but was allowed to vary for each year. Only interspecific interactions were included in the model fitting. As before, we assumed a negative binomial distribution, except for species with presence/absence data, for which we used a Bernoulli distribution. As before, we evaluated the likelihood in C++, and optimized it using TMB, nlminb, and 500 sets of starting values for the parameters obtained using Latin hypercube sampling In this phase, the population growth model was the following: !!,@,#$% = (1 − &!)!!,@,# + ∑ ?A!," %$B!,!:!,4," !!,C,#<@,CC (S1.15) where Ni,k,t is the observed abundance of the focal species i in year t in the kth 0.1 × 0.1 m square, (=!,#is the number of new individuals of the focal species at each year t in the absence of interactions assuming independence from the soil depth, and kk,l is the probability that a seed produced in square l reaches square k. This probability was assumed to follow the distribution proposed by Clark et al. (1999) (See appendix 3). Note that the summation goes through all squares, including the focal one, k. Di was fixed to one for annual species (Table 1). For the species with presence – absence data, we used instead: !!,@,#$% = 1 − 41 − (1 − &!)!!,@,#5.>∑ 56!," .*7!,!8!,4," :!,4,"D9,44 (S1.16) 69 Phase 2: In this phase the model included recruitment and competition as functions of soil depth. We also included the 1 × 1 m quadrats as a random effect. It also included the interaction coefficient (a) for the “multispecies”, which encompasses all species that were not selected as part of the model using Bayesian variable selection (See appendix 2). The abundances of those species add up to the abundance of the multispecies. We included the multispecies because we found that it improved model AIC in all of the 16 species where we tested its inclusion. We found problems of convergence during this phase because Bi,t(d), a and b were very correlated and depended on soil depth, resulting in identifiability problems. This was aggravated if we made b also dependent on depth. Thus, we kept b constant over depth, and, to simplify the number of parameters that needed to be calculated, we used the from phase 1. To do so we assumed that the average birth rate over the depth gradient was equal to : (=!,# = ∑ .,!,"$-!,"+/%F$0!,"(+/%F)$GHF +3IF 251 (the division by ten in the equation is required because the average was obtained using data calculated every 0.1 mm). The summation was preformed from d = 30 to d = 280 because soil depths in our study ranged from 30 mm to 280 mm. This was because plants growing in soil depths < 30 mm are extremely rare while soil depths > 280 mm are also rare and do not provide plants with extra moisture. The division by 251 was performed because there are 251 millimeters between 30 and 280 cm. Rearranging equation 2 we obtain: (=!,# = .,!," ∑ .-!,"+$0!,"+$!3% By applying a natural logarithm and rearranging we obtain: Bi ,t Bi ,t 70 @!,#3 ./?A!," ./∑ =:!,",*;!,", $ !-" (S1.17) Plugging equation 11 into our model we obtain the model that we used to obtain the rest of the parameters !!,#$% = (1 − &!)!!,# +∑ = <=56!," <=∑ ? :!,",*;!,", $ !-" *:!,",*;!,", $ ∏ 9%$:%,"; &!% %(! %$∑ =)!*+!,:%,"% C !!,#<@,C (S1.18) The parameters were estimated through maximum likelihood using non-linear regression of the observed Ni,t+1 on Ni,t. The facilitation (bij) and mortality (Di) parameters were bound between 0 and 1 using a logistic (inverse log-odds) transformation. The equivalent expression for species with presence-absence data was obtained from the same assumptions on ai,t (Equation S1.17), and the software for model fitting was the same used above. The models were fitted 500 times using different initial values from a Latin hypercube as before. References Adler, P.B., Smull, D., Beard, K.H., Choi, R.T., Furniss, T., Kulmatiski, A., et al. (2018). Competition and coexistence in plant communities: intraspecific competition is stronger than interspecific competition. Ecol. Lett., 21, 1319–1329. Bolker, B.M. (2007). Ecological Models and Data in R. Princeton University Press, Princeton. Bonanomi, G., Incerti, G., Capodilupo, M. & Mazzoleni, S. (2010). Rare self-facilitation in terrestrial plants as compared to aquatic sessile organisms: empirical evidences and causal mechanisms. Community Ecol., 11, 148–159. Brooker, R.W. & Callaghan, T. V. (1998). The balance between positive and negative 71 plant interactions and its relationship to environmental gradients: a model. Oikos, 81, 196–207. Carnell, R. (2021). lhs: Latin Hypercube Samples. R package version 1.1.3. https://CRAN.R-project.org/package=lhs Clark, J.S., Silman, M., Kern, R., Macklin, E. & HilleRisLambers, J. (1999). Seed dispersal near and far: Patterns across temperate and tropical forests. Ecology, 80, 1475–1494. Fournier, D.A., Skaug, H.J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M.N., et al. (2012). AD Model Builder: Using automatic differentiation for statistical inference of highly parameterized complex nonlinear models. Optim. Methods Softw., 27, 233–249. Fraaije, R.G.A., ter Braak, C.J.F., Verduyn, B., Breeman, L.B.S., Verhoeven, J.T.A. & Soons, M.B. (2015). Early plant recruitment stages set the template for the development of vegetation patterns along a hydrological gradient. Funct. Ecol., 29, 971–980. Freckleton, R.P. & Watkinson, a R. (2001). Asymetric competition between plant species. Funct. Ecol., 15, 615–623. Kristensen, K., Nielsen, A., Berg, C.W., Skaug, H. & Bell, B. (2016). TMB: Automatic differentiation and laplace approximation. J. Stat. Softw., 70, 1–21. Lacey, E.P. (1986). Onset of Reproductioin in Plants: Size-versus age-dependency. Trends Ecol. Evol., 1, 72–75. Martorell, C. & Freckleton, R.P. (2014). Testing the roles of competition, facilitation and stochasticity on community structure in a species-rich assemblage. J. Ecol., 102, 74– 72 85. May, R. & McLean, A. (2007). Theoretical Ecology: Principles and Applications. Third Edit. Oxford University Press. Pugnaire, F.I., Armas, C. & Valladares, F. (2004). Soil as a mediator in plant-plant interactions in a semi-arid community. J. Veg. Sci., 15, 85–92. Villarreal-Barajas, T. & Martorell, C. (2009). Species-specific disturbance tolerance, competition and positive interactions along an anthropogenic disturbance gradient. J. Veg. Sci., 20, 1027–1040. Zepeda, V. (2021). Variabilidad climática y mantenimiento de la diversidad vegetal en un pastizal rico en especies. Zepeda, V. & Martorell, C. (2019). Fluctuation-independent niche differentiation and relative non-linearity drive coexistence in a species-rich grassland. Ecology, 100, e02726. 73 Appendix 2. Bayesian Variable Selection To reduce the computational time required to adjust our models and avoid overfitting we limited the number of interacting species. To do so in an unbiased way, for each of the focal species models, we used Bayesian variable selection (BVS) to determine if companion species had a strong enough competitive or facilitative effect to be included in the model (Ovaskainen et al. 2017). We also used this method to determine if soil depth had an effect on focal species birth or species interactions. BVS only works with linear models, and our model (see equation 1 in the main text) cannot be linearized. Thus, as a proxy, we used the modified Gompertz population growth model used by (Ovaskainen et al. 2017), !!,#3!!,#.JK!9%>∑B!%./:!%,";L (S2.1) where Ni,t is the observed abundance of the focal species i in year t and αij is a measure of the competitive effect. We modified such model to allow for soil depth to affect the population growth rate (r =a + bd + cd2) and species interactions (α = g + hd), where d is soil depth. Then, we expanded the equation and applied the natural logarithm to both sides of the equation to eliminate the exponentials and obtain the following: ln!!,#$% = ln!!,# + @ + C* + D*G − ∑ 4E1 + F1* + G1*G +H1*I5ln!1,#1 (S2.2) where s = ag, u = ah + bg, v = bh + cg and w =ch. We used this equation with the R function GibbsBvs in the package BayesVarSel (Garcia-Donato & Forte 2018) to run model selection. The most basic model visited was ln!!,#$% = ln!!,# + @. Every time BVS found high inclusion probability for any parameter interacting with a companion, it was justified to incorporate this species into our population growth 74 model. Similarly, high inclusion probabilities of the soil depth parameter (d) suggest the inclusion of soil depth in the birth rates of our original model. A high inclusion probability for cubic terms means that including an effect of soil depth on companion species interactions in our population growth model may be justified. Validation of Bayesian Variable Selection Before going forth with BVS, it was necessary to validate whether the Gompertz model can be used to draw conclusions on the parameters of our model. Thus, we applied BVS to simulated data obtained using our population model. In other words, we simulated population size in time t+1 (Nt+1) in different scenarios using the equation: !!,#$% = (1 − &!)!!,# + =@!*:!,*;!, $ :!," %$∑ 96!,%$7!,%+;% :%," (S2.3) Where Di is mortality. Note that we have removed the exponential from the interaction parameter I!,1 + ℎ!,1* to allow the interaction coefficient to be negative, which in this simplified version of the model corresponds to facilitation rather than competition. It was important for us to see if facilitative interactions were also detectable. We simulated sets with 300 data. The values of parameters D, a, b and c were fixed for all simulations to 0.6, 1, 0.25 and -0.005 respectively. In all cases we included five species in the simulation: the focal species (species 1), an associated species that interacted with the focal one (species 2), and three (species 3–5) that did not interact to test whether BVS was able to screen the species that should be included in the model from accompanying species that have no role in it. One of these three species (species 5) had its densities correlated with soil depth in some simulations. This is a likely scenario which may lead to the detection of spurious interactions. For instance, if two species occur at different soil depth, BVS may incorrectly suggest a competitive interaction between them because 75 the presence of the associated species is correlated with a poor performance of the focal species. Because we were mainly interested in knowing if BVS correctly identified which interacting species to include in the model and how soil depth affected these interactions, we conducted simulations using different values for g and h. We then selected 300 values for d and the number Nj,t of individuals of the five species using a uniform distribution. Soil depth spanned an interval of 0 to 30 cm as occurs at our study site. Population sizes were simulated to be between 0 and 5 as frequently occurs within our study unit (10 × 10 cm squares). We then obtained the respective 300 Ni,t+1 using equation (S2.3). We incorporated noise to the data by adding a random number (normally distributed with mean = 0 and a standard deviation = 0.1.) to the values of Ni,t+1 obtained. We ran BVS on the simulated data using the linearized Gompertz model and observed if the relevant variables were identified. If for example gi,j was very small or gi,j = 0 we expected BVS to return small inclusion probabilities for selecting these variables. Detection of interacting species To establish how detectable species interactions were, we simulated Ni,t+1 with different values for the interaction coefficient of species 2 (g2). Focal species and associate species two, three, four and five interactions coefficients were held constant (g1 = 0.3, g3 = 0, g4 = 0, g5 = 0). Parameter h = 0 for all species in this simulation and noise was not included. BVS correctly identified interactions with species 2 with high posterior probabilities regardless of their intensity or sign (Table S2.1), suggesting that both competition and facilitation were detectable by BVS even when interaction intensity was small. None of the parameters set to zero in the simulation were selected by BVS. 76 Table S2.1. Detection of interactions differing in intensity in data without noise. Data correspond to the inclusion probabilities obtained from Bayesian variable selection. Values above 0.5 are shown in bold. Inclusion probability g2 = 0.01 g2= 0.1 g2 = 0.25 g2 = 0.5 g2 = -0.5 g2 = -0.1 g2 = -0.01 d 1.00 1.00 1.00 1.00 1.00 1.00 1.00 d 2 1.00 1.00 1.00 1.00 1.00 1.00 1.00 N1,t 1.00 1.00 1.00 1.00 1.00 1.00 1.00 N2,t 1.00 1.00 1.00 1.00 1.00 1.00 1.00 N3,t 0.01 0.01 0.01 0.01 0.00 0.01 0.01 N4,t 0.02 0.02 0.02 0.02 0.02 0.02 0.02 N5,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d × N1,t 0.05 0.05 0.15 0.05 0.05 0.05 0.06 d2 × N1,t 0.03 0.03 0.15 0.03 0.02 0.03 0.03 d 3 × N1,t 0.01 0.01 0.14 0.02 0.01 0.02 0.02 d × N1,t 0.17 0.17 0.01 0.16 0.17 0.16 0.17 d 2 × N1,t 0.06 0.06 0.02 0.06 0.06 0.06 0.06 d 3 × N1,t 0.04 0.04 0.02 0.05 0.05 0.05 0.04 d × N2,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d 2 × N2,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d 3 × N2,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d × N3,t 0.01 0.01 0.02 0.01 0.01 0.01 0.01 d 2 × N3,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d 3 × N3,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d × N4,t 0.05 0.05 0.01 0.05 0.05 0.05 0.05 d 2 × N4,t 0.04 0.04 0.01 0.04 0.04 0.04 0.04 d 3 × N4,t 0.03 0.03 0.01 0.03 0.03 0.03 0.04 d × N5,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d 2 × N5,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 d 3 × N5,t 0.01 0.01 0.01 0.01 0.01 0.01 0.01 77 Next, added noise to Ni,t+1, as our field-collected data are necessarily noisy using the same parameters as above. When noise was added the interaction with species 2 was always detectable, although inclusion probabilities were not too high when the interaction was weak (g2 £ 0.1). The intraspecific interaction was missed when was very small (- 0.1 < g1 < 0.1 ), but appeared in interaction with soil depth with high probabilities. BVS also included some parameters (interactions with soil depth) that could incorrectly suggest that h = 0 (Table S2.2). Table S2.2. Detection of interactions differing in intensity in data with noise. Data correspond to the inclusion probabilities obtained from Bayesian variable selection. Values above 0.5 are shown in bold. Inclusion probability g2 = 0.01 g2= 0.1 g2 = 0.25 g2 = 0.5 g2 = -0.5 g2 = -0.1 g2 = -0.01 d 1.00 1.00 1.00 1.00 1.00 1.00 1.00 d 2 1.00 1.00 1.00 1.00 1.00 1.00 1.00 N1,t 0.02 0.02 1.00 1.00 1.00 1.00 0.01 N2,t 0.54 0.54 1.00 1.00 1.00 1.00 1.00 N3,t 0.01 0.01 0.15 0.01 0.01 0.02 0.01 N4,t 0.02 0.02 0.01 0.01 0.01 0.01 0.01 N5,t 0.02 0.02 0.04 0.01 0.01 0.03 0.03 d × N1,t 0.92 0.92 0.97 0.95 0.95 1.00 0.54 d2 × N1,t 0.52 0.52 0.50 0.46 0.46 1.00 0.44 d 3 × N1,t 0.41 0.41 0.23 0.15 0.15 1.00 0.06 d × N1,t 0.01 0.01 0.01 0.01 0.01 0.02 0.01 d 2 × N1,t 0.01 0.01 0.01 0.01 0.01 0.02 0.01 d 3 × N1,t 0.02 0.02 0.01 0.01 0.01 0.02 0.01 d × N2,t 0.50 0.50 0.06 0.02 0.02 0.64 0.18 d 2 × N2,t 0.07 0.07 0.04 0.02 0.02 0.25 0.24 d 3 × N2,t 0.05 0.05 0.03 0.03 0.03 0.08 0.20 78 d × N3,t 0.02 0.02 0.10 0.01 0.01 0.01 0.01 d 2 × N3,t 0.02 0.02 0.06 0.01 0.01 0.01 0.01 d 3 × N3,t 0.01 0.01 0.03 0.01 0.01 0.02 0.01 d × N4,t 0.02 0.02 0.01 0.01 0.01 0.01 0.02 d 2 × N4,t 0.02 0.02 0.01 0.01 0.01 0.02 0.01 d 3 × N4,t 0.02 0.02 0.01 0.01 0.01 0.02 0.01 d × N5,t 0.02 0.02 0.03 0.01 0.01 0.03 0.04 d 2 × N5,t 0.02 0.02 0.02 0.01 0.01 0.03 0.04 d 3 × N5,t 0.02 0.02 0.02 0.01 0.01 0.02 0.05 Spurious interactions due to density-soil depth correlations. To simulate the situation where associated species density is correlated with soil depth, we added 1- d/k to the abundance of companion species 5 (N5), such that it became more abundant at shallow soils. We varied k to change the intensity of the species’ correlation to soil depth (k = 2, k = 5, k = 10 and k = 20), making the correlation stronger when values of k are smaller. Parameters g1 and g2 where set to 0.3, while g3 –g5 and all h were set to zero. Noise was not included in this simulation. We found that BVS did not include species 5 as an interacting species except when its correlation with soil depth was very strong (k = 2; Table S2.3). As before, in some cases intraspecific interaction was detected in interactions with soil depth, but not on its own. We then ran a similar simulation to understand if a correlation of species abundance with stress was further affected by noise, with similar results (Table S2.4). Table S2.3. Detection of spurious correlations due to correlations between soil depth and density of associated species. These results do not include noise. Data correspond to the 79 inclusion probabilities obtained from Bayesian variable selection. Values above 0.5 are shown in bold. Inclusion probability k = 2 k = 5 k = 10 k = 20 d 1.00 1.00 1.00 1.00 d 2 0.01 0.99 1.00 1.00 N1,t 0.02 0.27 0.42 0.08 N2,t 1.00 0.72 1.00 1.00 N3,t 0.01 0.01 0.01 0.03 N4,t 0.03 0.01 0.01 0.01 N5,t 0.01 0.02 0.02 0.01 d × N1,t 1.00 0.05 1.00 0.05 d2 × N1,t 1.00 0.09 0.69 0.84 d 3 × N1,t 1.00 0.11 0.32 0.85 d × N1,t 0.01 0.22 0.05 0.02 d 2 × N1,t 0.02 0.11 0.03 0.01 d 3 × N1,t 0.02 0.10 0.03 0.01 d × N2,t 0.08 0.93 0.02 0.01 d 2 × N2,t 0.20 0.76 0.01 0.01 d 3 × N2,t 0.22 0.35 0.01 0.01 d × N3,t 0.01 0.02 0.01 0.01 d 2 × N3,t 0.01 0.02 0.01 0.01 d 3 × N3,t 0.01 0.01 0.01 0.01 d × N4,t 0.01 0.01 0.01 0.01 d 2 × N4,t 0.01 0.01 0.02 0.00 d 3 × N4,t 0.01 0.01 0.02 0.01 d × N5,t 0.01 0.03 0.01 0.01 d 2 × N5,t 0.01 0.04 0.01 0.01 d 3 × N5,t 1.00 0.46 0.01 0.01 Table S2.4. Detection of spurious correlations due to correlations between soil depth and density of associated species. These results include noise. Data correspond to the 80 inclusion probabilities obtained from Bayesian variable selection. Values above 0.5 are shown in bold. Inclusion probability k = 2 k = 5 k = 10 k = 20 d 1.00 1.00 1.00 1.00 d 2 0.33 1.00 1.00 1.00 N1,t 0.80 0.02 0.06 0.01 N2,t 1.00 1.00 1.00 1.00 N3,t 0.01 0.01 0.00 0.00 N4,t 0.00 0.00 0.01 0.00 N5,t 0.04 0.01 0.00 0.01 d × N1,t 0.02 0.01 0.09 0.93 d2 × N1,t 0.03 0.00 0.16 0.08 d 3 × N1,t 0.05 0.00 0.43 0.04 d × N1,t 0.14 0.01 0.04 0.00 d 2 × N1,t 0.07 0.01 0.03 0.00 d 3 × N1,t 0.09 0.01 0.02 0.00 d × N2,t 0.05 0.24 0.01 0.01 d 2 × N2,t 0.02 0.26 0.01 0.01 d 3 × N2,t 0.01 0.19 0.01 0.01 d × N3,t 0.01 0.01 0.00 0.00 d 2 × N3,t 0.01 0.01 0.00 0.00 d 3 × N3,t 0.01 0.01 0.00 0.00 d × N4,t 0.00 0.00 0.00 0.00 d 2 × N4,t 0.00 0.01 0.00 0.00 d 3 × N4,t 0.00 0.00 0.00 0.00 d × N5,t 0.01 0.41 0.00 0.00 d 2 × N5,t 0.69 0.28 0.00 0.01 d 3 × N5,t 0.07 0.29 0.01 0.02 Simulation of a soil depth effect on species interactions 81 Soil depth was expected to have an effect not only on birth but also on the strength of species interactions. We tested if this effect was detectable using simulations where h was different from zero. As before, g1 and g2 where set to 0.3, while g3 –g5 were set to zero. We first modified h1 to vary the strength of the intraspecific interactions (focal species) as a function of soil depth. Noise was included in the data. If BVS finds high inclusion probability for cubic terms, then it is unequivocally detecting the effect of soil depth on species interactions that we simulated because w in equation S2.2 can only be different from zero if h is different from zero. Using this criterion, BVS correctly identified that h1 should be included in the model except for h1 = 0.05, which was nevertheless close to the critical probability of 0.5 (Table S2.5). In these analyses BVS indicated incorrectly that interactions with other species should be included in the model. Table S2.5. Detection changes of intraspecific interactions with soil depth. These results include noise. Data correspond to the inclusion probabilities obtained from Bayesian variable selection. Values above 0.5 are shown in bold. Inclusion probability h1 = 0.5 h1 = 0.25 h1 = 0.1 h1 = 0.05 h1 = 0.01 d 1.00 1.00 1.00 1.00 1.00 d 2 1.00 1.00 1.00 0.79 1.00 N1,t 0.61 0.99 0.82 0.35 0.22 N2,t 1.00 0.13 0.95 0.28 0.61 N3,t 0.12 0.21 0.21 0.22 0.12 N4,t 0.06 0.12 0.10 0.23 0.12 N5,t 0.12 0.10 0.07 0.14 0.44 d × N1,t 1.00 1.00 1.00 0.98 1.00 d2 × N1,t 1.00 1.00 1.00 0.47 1.00 d 3 × N1,t 1.00 1.00 1.00 0.45 1.00 d × N1,t 0.14 0.15 0.59 0.24 1.00 d 2 × N1,t 0.10 0.13 0.34 0.28 1.00 d 3 × N1,t 0.09 0.13 0.29 0.30 1.00 82 d × N2,t 1.00 0.96 0.50 0.33 0.63 d 2 × N2,t 1.00 1.00 0.34 0.36 0.61 d 3 × N2,t 1.00 1.00 0.32 0.35 0.55 d × N3,t 0.08 0.10 0.10 0.36 0.17 d 2 × N3,t 0.08 0.09 0.09 0.80 0.17 d 3 × N3,t 0.09 0.08 0.08 0.90 0.17 d × N4,t 0.08 0.07 0.08 0.28 0.13 d 2 × N4,t 0.08 0.07 0.09 0.26 0.14 d 3 × N4,t 0.08 0.08 0.09 0.25 0.14 Then, we simulated an effect of soil depth on the interaction species 2 by including non- zero values of h2. In this case, the cubic term was only selected by BVS when h was quite large, although the respective quadratic or linear terms were selected except for the lowest values of h (h2 = 0.01, Table S2.6). In contrast with the scenario in which intraspecific interactions changed with soil depth, spurious interactions with species 3 – 5 were not detected in these simulations. Table S2.6. Detection changes of interspecific interactions with soil depth. These results include noise. Data correspond to the inclusion probabilities obtained from Bayesian variable selection. Values above 0.5 are shown in bold. Inclusion probability h2 = 0.5 h2 = 0.25 h2 = 0.1 h2 = 0.05 h2 = 0.01 d 1.00 1.00 1.00 1.00 1.00 d 2 1.00 1.00 1.00 1.00 1.00 N1,t 0.93 1.00 1.00 1.00 1.00 N2,t 0.07 0.23 1.00 0.68 1.00 N3,t 0.04 0.09 0.02 0.02 0.01 N4,t 0.03 0.02 0.04 0.01 0.01 N5,t 0.03 0.02 0.02 0.10 0.01 d × N1,t 0.52 0.06 0.60 0.47 0.01 d2 × N1,t 0.65 0.08 0.27 0.22 0.01 d 3 × N1,t 0.78 0.09 0.15 0.13 0.01 d × N1,t 1.00 1.00 1.00 0.53 0.01 83 d 2 × N1,t 0.97 0.76 0.41 0.27 0.01 d 3 × N1,t 0.73 0.24 0.15 0.34 0.01 d × N2,t 0.04 0.16 0.03 0.32 0.01 d 2 × N2,t 0.04 0.10 0.02 0.08 0.01 d 3 × N2,t 0.05 0.08 0.02 0.05 0.01 d × N3,t 0.05 0.04 0.02 0.02 0.01 d 2 × N3,t 0.05 0.03 0.02 0.01 0.01 d 3 × N3,t 0.05 0.03 0.02 0.01 0.01 d × N4,t 0.03 0.02 0.06 0.01 0.01 d 2 × N4,t 0.03 0.01 0.08 0.01 0.01 d 3 × N4,t 0.04 0.02 0.11 0.01 0.00 Conclusions With the exception of cubic terms of associated species, BVS always identified correctly which variables should be included in the model, although it frequently selected additional variables that did not play a role in the generation of the artificial data. The inclusion of extra parameters is not ideal for computation times, but it is likely that these parameters will later be estimated as ~ 0 when fitting our model to the data. In fact, this was the case for several of the parameters chosen by BVS. Because of this and the difficulty to detect interaction changes based on cubic terms, when fitting our model we preferred to include a parameter h for all species in which any interaction between density and soil depth was found. We also included both competition and facilitation (a and b in equation 1 in the main text) for all associated species that were selected by BVS. Intraspecific competition was always included in the model based on the finding of previous studies that it has strong effects at the study site (Martorell & Freckleton 2014; Zepeda & Martorell 2019). Associated species correlation with soil depth only affected model selection when correlation was very strong. This seems to be relatively unusual in our study site but may still have an effect on our results. 84 References Garcia-Donato, G. & Forte, A. (2018). Bayesian testing, variable selection and model averaging in linear models using R with BayesVarSel. R J., 10, 155–174. Martorell, C. & Freckleton, R.P. (2014). Testing the roles of competition, facilitation and stochasticity on community structure in a species-rich assemblage. J. Ecol., 102, 74– 85. Ovaskainen, O., Tikhonov, G., Dunson, D., Grøtan, V., Engen, S., Sæther, B., et al. (2017). How are species interactions structured in species-rich communities ? A new method for analysing time-series data. Proc. R. Soc. B Biol. Sci., 284, 1–7. Zepeda, V. & Martorell, C. (2019). Fluctuation-independent niche differentiation and relative non-linearity drive coexistence in a species-rich grassland. Ecology, 100, e02726. Appendix 3. Dispersal Kernel Development of the model for dispersal estimation In models where interactions are estimated from time series in natural systems, dispersal may result in overestimated intraspecific competition strengths (Freckleton & Watkinson 2000). To avoid this problem, we introduced dispersal in the model fitting procedure, as suggested by Martorell & Freckleton (2014). Given the computationally demanding process involved in estimating the dispersal parameters, these were estimated using a simplified, single species version of the model shown in equation (1) of the main text that comprised birth (B) and death (D) rates, as well as intraspecific competition. Note that in this appendix we are using a slightly different notation as only one species is involved and spatial terms need to be handled more explicitly than in the main text and other appendices. In the model, birth rate varied annually due to previous evidence (Martorell & Freckleton 2014; Zepeda & Martorell 2019). The dispersal parameters estimated using this model where then used during the fitting of the multispecies model with soil depth and facilitation (See Appendix S1). Such multi-phase approach, where some parameters are fitted in simple models and are later fed into more complicated ones, is used when there is high dimensionality (Fournier et al. 2012). To model seed dispersal we used a 2-D dispersal kernel originally proposed by Clark et al. (1999). This model involves a parameter related to mean dispersal distance η, and another parameter that determines the kernel shape c. The kernel is given by the mathematical expression: !(#) = !" e#$ ! "%#, (S3.1) where r is dispersal distance and M is a normalizing term which ensures that the volume under the whole surface is equal to 1 (making this function into a probability density distribution), and is given by: ( = &'($)$$#%* , (S3.2) where Γ is the gamma function. This kernel is radially symmetric, has adaptable kurtosis, and generalizes a few common probability density functions: * = 1 gives an exponential distribution, while * = 2 gives a Gaussian distribution (Figure S3.1, Clark et al 1999). If the mother plant is located at coordinates x’ and y’, equation (S3.1) can also be expressed in cartesian coordinates as a function ,(-, /): ,(-, /) = !" e #+%&'(')*$+&,(,)*$" , # (S3.3) Figure S3.1 Two contrasting shapes for the dispersal kernel. The function ! with parameter * = 1 is shown in the left panel and with parameter * = 2 in the right. The horizontal axes correspond to spatial coordinates, and the vertical axis is the probability density that a seed reaches a specific position in space. The mother plant is located at the center of the plot, where the probability density is maximal. Our data consisted of measurements of population sizes for all species occurring in 0.1 ⨉ 0.1 m subsquares of 1 ⨉ 1 m squares. Using the above kernel (equation S3.1) we can calculate the proportion of seeds produced in one subsquare that arrive to another. The full 1 ⨉ 1 m square consists of 100 subsquares which can be numbered with two indices running between 1 and 10, and where the top-left subsquare is the (1,1)-subsquare and the bottom-right subsquare is the (10,10)-subsquare. Assuming that seeds disperse from the center of each square, we get the following expression for the proportion of seeds k produced1 in from the (i, j)-subsquare that disperse into the (k, l)-subsquare: 2((3, 4), (5, 6)) = ∫ ∫ ,(-, /) 8- 8/-#./0.2-#.#0.23#4/0.23#4#0.2 . (S3.4) Using this result, we modeled the population size in each subsquare as a function of time in a way similar to the population growth model shown in equation (1) in the main text: 9-,3(: + 1) = (1 − =)9-,3(:) + ∑ ∑ ?.,4(:)!046!!0.6! 2@(3, 4), (5, 6)A, (S3.5) where Ni,j(t) is the population size at the (i, j)-square at time : , = is the death rate, B is the intraspecific competition coefficient, and ?.,4(:) = 7-∙9.,0(;)!/=9.,0(;) (where Bt is the birth rate at time t) represents the number of individuals that are produced in the (i, j)-subsquare. This equation tells us that the number of individuals in subsquare (k, l) at time : + 1 is equal to the proportion of the individuals in that square at time : that survived, plus the sum of all the individuals that dispersed into the subsquare from all subsquares (including the same one). 1 It must be noted that when we refer to seed numbers we are only considering those that will get established, not the total numbers of seeds produced, dispersed or arriving to a subsquare However, a few problems had to be resolved before applying the above approach. First, because we only had data from 20 out of the 100 subsquares in each square, we could not perform the sum over all subsquares in equation S3.5. To solve this problem, we performed a linear interpolation of values for subsquares that were missing. To do so, we built a Delaunay triangulation of the 1 m2 square taking the subsquares for which we had data as vertices (Lee & Schachter 1980). For the triangulation to cover the whole square, its four corners need to be considered as vertices in the triangulation regardless of whether they were in the sample of 20 subsquares or not. We obtained the function that described the unique plane passing through the zi,j(t) values corresponding to the three vertices of the triangle. Then, for each subsquare for which we had no data within a triangle we calculated the respective z value on such plane (Figure S3.2). For corners for which we had no data and thus a zi,j(t) value could not be calculated, we assigned the zi,j(t) value of the closest measured subsquare. If two or more subsquares were at the same distance from the corner, we used the average of their zi,j(t) values. Figure S3.2. Interpolation of seed production of subsquares for which there were no data (gray). In this example the tree subsquares shown in red have data and are the corners or one of the triangles obtained from the Delauney triangulation. Seed production is assumed to change linearly over the triangle. Second, because f always decreases with the radius, and since subsquares near the edge of the square are on average farther from the rest of the subsquares, these subsquares receive, on average, a smaller portion of the seeds dispersed within the square. This does not correspond to the situation in nature, where the 1 ⨉ 1 m square was just a portion of the grassland delimited by imaginary lines. To address this problem, we obtained the matrix A which elements are C-,3 = ∑ ∑ 2@(3, 4), (5, 6)A4. (S3.6) The matrix A has larger values towards its center and smaller values towards the edges due to the artificial lack of seed donors near the border of the square. We then calculated the matrix A’ whose elements are C-,3> = ?@A(B)C1,2 (S3.7) A’ is an edge correcting matrix which has entries equal to 1 at locations that receive the maximal proportion of seeds, and grows towards the edges, where the proportion is smaller. We then calculate sk,j(t), which is the number of seeds that arrive to subsquare (k,j)at time t corrected for the bias experienced by subsquares near the border of the square DD,E(E) = ∑ ∑ ?.,4(:)!046!!0.6! 2@(3, 4), (5, 6)AC-,3> . (S3.8) Finally, in the above model there would be a difference between the sum of all sk,j(t) values and the sum of all zi,j(t) values calculated for the subsquares. This difference, D(t), would represent the number of successfully dispersed individuals which are lost to the exterior of the square. In a more natural setting, assuming that the density of the population in the vicinity of the square is similar to that inside it, this loss to the exterior would be compensated by a comparable arrival of seeds from the vicinity into the square. Since we had no information about the population outside the square, the most parsimonious solution was to distribute D(t) between all the 100 subsquares equally. Thus, equation (S3.5) after correcting for edge effects and dispersal out of the square becomes: 9-,3(: + 1) = (1 − =)9-,3(:) + F-,3(:) + F(;)!00 . (S3.9) Implementation In this work a maximum likelihood approach was taken to estimate the parameters of the model (S3.9) for each individual species. For a given species and a point in parameter space, the likelihood was calculated in the following way: for each subsquare sampled and year, the estimated population size for the next year was calculated using equation (S3.9), and this value was used as the mean of a negative binomial distribution (the variance parameter q of the distribution was estimated during model fitting); then we obtained the probability of the observed number of individuals in the subsquare in the next year from this distribution. Assuming independence we have that the probability of all the events is equal to the product of the probabilities of the individual events. From this we obtain the following expression for the log likelihood of the parameters given the data: lnJℒ(L ∣∣ NGHI;/! )O = 6P(∏89R(9GHI;/!, S, 9JI;;/!)) = ∑ ln(89R(9GHI;/!, S, 9JI;;/!)). (S3.10) Where L represents a particular set of parameters, NGHI;/! is the total data for a given species, 89R is the negative binomial distribution, 9GHI;/! and 9JI;;/! are the observed and estimated population sizes at time : + 1 respectively, and the sum is over all the available valid data from field observations. Finally, the nlminb optimization routine and the TMB (Kristensen et al. 2016) package in R were used to search for parameters that maximized the likelihood (See Appendix S1 for details). Note: in the TMB package the 89R distribution is implemented using size and prob parameters, so in the code we use the F3?T = q and the equality U#VW = KK/934--+5. Estimation of parameters via maximum likelihood is a process that requires the calculation of equation (S3.9) several thousand times. It was thus important to optimize the code so that the calculations ran as fast as possible. This was done by simplifying the interpolation of zi,j(t) and the calculation of the dispersal probabilities. Interpolation of zi,j(t) values: This process requires the application of Delaunay’s triangulation and the interpolation of zi,j(t). Much of these calculations can be done as part of the pre-processing of the data before searching for the maximum likelihood parameters. To do so, first we made the interpolation for each square using the deldir (Turner 2021) package, and from which we then created a 100 × 20 matrix that codifies the linear interpolation of the zi,j(t) values. To understand how this was done, let [! = J3!, 4!, ?.5,45(:)O , [& = J3&, 4&, ?.$,4$(:)O , [L = J3L, 4L, ?.6,46(:)O be three vectors whose i and j values correspond to the vertices of a minimum triangle in the triangulation. If we take a point (i*, j*) inside the triangle, there is only one value z* such that the vector P*= (i*, j*, ?.∗,4∗(:)) falls on the plane that passes through all three vectors [!, [&, [L. This vector P* must satisfy @([& − [!) × ([L − [!)A ⋅ ([∗ − [!) = 0. (S3.11) Where × represents the cross product of vectors in space. The intuition behind the equation above is that vectors [& − [! and [& − [! represent two sides of the triangle defined by the three points, so those vectors lie parallel to the plane that passes through the points. The cross product of two non-zero vectors result in a vector which is normal to the plane defined by the vectors, which means ([& − [!) × ([L − [!) is normal to the plane we are looking for. We then observe that if [∗ lies on this plane then ([∗ − [!) is parallel to the plane, and thereby perpendicular to ([& − [!) × ([L − [!). Finally, two vectors are perpendicular in space if and only if their dot product is equal to zero. Thus (S3.11) correctly characterizes the points on the plane that pass through the three points. Solving the equality (S3.11) for ? gives an expression of the following form: ?.∗,4∗(:) = b!?.5,45(:) + b&?.$,4$(:) + bL?.6,46(:) (S3.12) Where b! = 1 − b& − bL , b& = .64∗#46.∗4$.6#.$46 and bL = 4$.∗#.$4∗4$.6#.$46 . Which means we can express the value of ?.∗,4∗(:) as a linear combination of the values ?.5,45(:), ?.$,4$(:), ?.6,46(:) with coefficients which are fully determined by the geometry of the triangulation and are independent of the parameters that we need to estimate via maximum likelihood. This means that we can avoid calculating these values in the likelihood maximization algorithm by instead supplying a 100 ´ 20 matrix B for each 1 m2 square such that each row corresponding to one subsquare and contains the b coefficients values required to calculate its zi,j(t) value. The columns in this matrix correspond to the 20 subsquares for which we have data. Given a subsquare with coordinates (3, 4) there are four cases. If (3, 4) is a vertex in the triangulation (i.e., one of the measured subsquares) then its corresponding row in the ^ matrix will have only a unique value of 1 in the column that corresponds to itself. If (3, 4) is an unmeasured vertex in one of the corners of the 1 m2 square, it will have entries equal to !N# in the columns corresponding to the P* closest subsquares (if there is only one closest subsquare P* = 1). Another case is when (3, 4) falls within a triangle whose vertices are all measured subsquares. In this case we put the coefficients b!, b&, bL defined above at the corresponding position for the vertices. The final case is when (3, 4) falls inside of a triangle whose vertices are not all measured subsquares (this is the case when at least one is an unmeasured corner). When this happens instead of just putting value b- at the corresponding vertex position, we put O1N#. That is, we distribute the value b- evenly over the measured subsquares from which the value of the corner subsquare is defined. Constructing B in this way it follows that, if _̀ is a 20 element vector containing the zi,j(t) values calculated for each of the subsquares for which we have data, then ^ ⋅ _̀ is a 100 element vector with the interpolated values for all 100 subsquares. There is still however one detail which was obviated in the discussion above. That is how to test whether a particular point (3, 4) falls within a triangle if we are given its vertices. The idea for this was based on the following observations: if we take a line segment between two vertices S! and S&of the triangulation, we can easily define a transformation ab (the composition of a translation and a rotation) that takes S! to (0,0) and places S& on the positive - axis. Now, for every other point q* we say that it is on the outer side of the line segment if ab(S∗) is in the upper semi-plane (its / coordinate is positive), and we say it is on the inner side of the line segment otherwise. To see the connection with the problem at hand notice that, given a triangle with vertices S!, S&, SL, a point S∗ will be in the interior of the triangle if and only if it is on the same side of all three line segments U!U&, U&UL, ULU!. The way this is implemented in the code is that we fill a matrix sideValMat with values 1 or −1 depending on whether a point is on the inner or outer side respectively of each side of a triangle in the triangulation (indices of all 100 subsquares on the vertical and indices of sides on the horizontal). Then we find the triangle each point belongs to by checking, for each triangle the sum of these values for each point. A point will be inside the triangle if and only if the sum of the three values is exactly equal to 3 or −3. In the code, dd is a matrix which contains all the ^ matrices for all distinct squares stacked on top of each other. The matrix sq contains the information to map each set of 20 measured subsquares to the index of the square which corresponds to it. The matrix yearmat has as many rows as datapoints and as many columns as measured years, and each row has entries all equal to 0 except for a 1 in the column number which corresponds to the year when the measurement was taken, and is used to apply the correct birth rate which depends on the year. All of these are supplied as parameters to the TMB optimizer, and are programmed in R. Preprocessing of the data includes removal of measurements from squares that were totally empty in a certain year, as well as removal of data points which were known to be incorrect or not useful. This is performed by the CleanDT function. Calculation of dispersal probabilities: In practice it was not possible to find an analytic expression for the dispersal probabilities k (equation S3.4), so numeric approximation was used. This is a computationally intensive procedure. To reduce the computation time required, rather than solving equation S3.4 for each of the 5,050 possible pairs of subsquares in the 1 ⨉ 1 m square, we used the property that dispersal between subsquares depends only on the relative positions to each other. The maximal distance between two squares happens when they are in opposite corners, and thus are separated by nine squares in both x and y axes. Therefore, all the possible relative positions between pairs of squares can be represented in a 19 ⨉ 19 matrix M in which the mother plant is located at the center of the matrix, i.e, in the cell (10,10) (figure S3.3). If the elements of M are defined as e-,3 = 2((10,10), (5, 6)), the matrix contains all possible values of k given parameters f and *. Furthermore, due to radial symmetry only a subset of the entries of M must be calculated, and the rest can be filled in by some reflection (figure S3.3). This has the computational advantage that we no longer have to calculate all 361 entries of the matrix solving the double integral numerically, but only do this 55 times. In the C++ file the matrix M is called rainMtrxFull. Having obtained M, 10 ⨉ 10 submatrices are taken from it (as shown in figure S3.3) representing dispersal from each of the 100 subsquares and stored in a 10 ⨉ 10 ⨉ 100 array sResTot. Adding up the 100 layers of sResTot entry-by-entry we obtain the matrix A mentioned above; and the dot product of the full (^ ⋅ _̀) value vector with each of the 100-long columns of sResTot(keeping the first and second dimensions fixed) gives the corresponding total dispersal into each subsquare. Figure S3.3 Dispersal probability matrix M. Light colors represent greater probabilities. The region enclosed in red represents the 55 entries of the matrix from which all others can be obtained by simple reflection. The blue square encloses a region which represents the submatrix we would take to obtain the dispersal matrix from subsquare (8,6). For each subsquare position in the square a submatrix can be obtained in a similar way. The last part of the C++ code used for rapid optimization with TMB package is just a straightforward calculation of the likelihood as described above. References Clark, J.S., Silman, M., Kern, R., Macklin, E. & HilleRisLambers, J. (1999). Seed dispersal near and far: Patterns across temperate and tropical forests. Ecology, 80, 1475–1494. Fournier, D.A., Skaug, H.J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M.N., et al. (2012). AD Model Builder: Using automatic differentiation for statistical inference of highly parameterized complex nonlinear models. Optim. Methods Softw., 27, 233–249. Freckleton, R.P. & Watkinson, A.R. (2000). On detecting and measuring competition in spatially structured plant communities. Ecol. Lett., 3, 423–432. Kristensen, K., Nielsen, A., Berg, C.W., Skaug, H. & Bell, B. (2016). TMB: Automatic differentiation and laplace approximation. J. Stat. Softw., 70, 1–21. Lee, D.T. & Schachter, B.J. (1980). Two algorithms for constructing a Delaunay triangulation. Int. J. Comput. Inf. Sci., 9, 219–242. Martorell, C. & Freckleton, R.P. (2014). Testing the roles of competition, facilitation and stochasticity on community structure in a species-rich assemblage. J. Ecol., 102, 74–85. Turner, R. (2021). deldir: Delaunay Triangulation and Dirichlet (Voronoi) Tessellation. R package version 1.0-6. https://CRAN.R-project.org/package=deldir Zepeda, V. & Martorell, C. (2019). Fluctuation-independent niche differentiation and relative non-linearity drive coexistence in a species-rich grassland. Ecology, 100, e02726. 96 Appendix 4. Community simulations Using our fitted models for each species we simulated our community along the soil depth gradient. This was done by iterating equation (1) in the main text. In each iteration we selected randomly one year between 2010 and 2016, and used the birth rates estimated for that year. Because this is a stochastic model, species densities do not converge to constant a value, but instead to a stationary distribution. It appeared that 100 iterations sufficed to reach the stationary state. The model then was iterated more times, and the average abundance of each species over time was estimated. In some preliminary simulations we tried 1000, 5000, 10 000 iterations to calculate the mean density over time using different initial densities from a uniform distribution between 0.001 and 1. The result was nearly the same using 5000 and 10 000 iterations, suggesting that the system has a single stationary state and that 10 000 iterations provide a reasonable estimate of its average. Nevertheless, the system was iterated 100 000 times to reach a more precise estimate of the species abundances. This procedure was conducted individually for every depth between 3 and 28 cm every 0.1 cm. To test our hypotheses, it was necessary to eliminate facilitation and all interactions from our models. To eliminate facilitation we set bi,j = 0. Thus, our model without facilitation was !!,#$% = (1 − &!)!!,# + !!,# ?!," %$∑ B!%:%,"% . (S4.1) To eliminate all interactions simulations were not necessary because abundances do not depend on the abundances of other species. Under these conditions, the stable state N* can be solved analytically as !∗ = ?>N B!!N . (S4.2) 97 Comparing hydrological niches Hellinger’s distance H is an equivalent of the Euclidean distance that is used to compare two probability density functions P and Q as K(L, N) = % √G O∑4PQ! − PR!5G. It has a maximum value of one if both distributions are completely different, i.e., if P has nonzero values for every q = 0 and vice versa. The main advantage of H over Euclidean distance is that in the latter the niches of two species that never co-occur over the depth gradient can be smaller than the distance between species that co-occur at some places (Conde & Domínguez 2018). To calculate Hellinger’s distances we first standardized the simulated abundance of each species over the depth gradient to sum one. Hellinger’s distances satisfy the triangle inequality H(P,Q) £ H(P,R) + H(Q,R). This makes it appropriate for use in polar ordination, which is based on triangular relationships between each species and (typically) the two that are most distant to each other (Bray & Curtis 1957). We used polar ordination to reduce the Hellinger’s distances between all species pairs to a two-dimensional Euclidean space where all species can be placed. This greatly simplified the construction of an appropriate null model for the k-means algorithm. Polar ordination was conducted using package natto (Oksanen 2021) for R. The first two polar ordination axes combined captured > 60 % of the variation between species. Implementing the k-means algorithm and its null model The k-means algorithm searches for the number of clusters that minimize variation within groups for an a priori defined number of groups k. The within-group variation is reduced as k is increased. Thus, the optimum value of k is arguably equal to the number of species because then each group contains one species and there is no within-group variance. To 98 solve this problem, D’Andrea et al. (2020) propose the use of the gap statistic. The value of k is selected as the number of clusters that maximizes the differences in within-cluster dispersion between observed and null communities. We applied this procedure to the results of the polar ordination (Fig. S1) using k values between 1 and 8. To do so we used function kmeans in R. We found that the repeated application of this function using the default parameters to the same dataset produced different clusters, so we set the attributes iter.max = 10 000 000 and nstart = 100 which resulted consistent results and the succesfull minimization of within-group variation. This algorithm was then applied to 1000 null communities and define the gap index for k clusters as S@ = % %FFF ∑ log(&@∗) − log (&@)PQCCR (S4.3) where &@∗ is the within cluster dispersion in a given null community and Dk is the dispersion in the observed data. The gap index measures the degree to which the observed niche dispersion differs from what is expected by chance. Positive values of Gk indicate tighter than expected clustering in the data for k clusters. This procedure was performed for cluster number ranging from one to eight. We did not consider more than eight clusters because, having only 27 species in the analysis, more than eight clusters would imply that most clusters would include only one or two species and precludes the examination of emergent neutrality, hidden niches or species alliances. In its original formulation, the algorithm selects the number of clusters that maximizes Gk. However, this resulted in very large numbers of clusters in our dataset, so we selected the smallest value of k for which Dk was significantly larger than &@∗. The p value was calculated is the proportion of null communities with a larger gap statistic than the observed community. 99 To see if our results differed from random, we applied this algorithm to 1000 null communities. To recreate the distribution that the species show in our polar ordination (Fig. S1), were most species occur near the center of the plot and become sparser towards the periphery, we used a bivariate normal distribution with parameters resembling those of our polar ordination to create our null communities. To do so, we estimated the means and variance-covariance matrix of the bivariate normal distribution that best fit the data in the polar ordination. We then used these parameters and function rmvnorm from package mvtnorm (Genz et al. 2021) to produce a random set of species to which the k-means algorithm was applied. We then compared these results using the null communities to our observed results. 100 Figure S4.1. Polar ordination of species using the Hellinger distance matrix. The species abundance data in this example come from the simulations with all interactions. See full species names in Table 1. -0.2 0.0 0.2 0.4 0.6 0.8 1.0 0 .0 0 .2 0 .4 0 .6 0 .8 Axis 1 A x is 2 A. adscensionis A. divaricata. BouhirB. scorpioides E. mendezii E. sericeus F. purpurea H. pinnattum M. kunthii M. biflora M. peruviana M. phalarioides M. rigida O. lunulata P. oligospermus P. nivea P. pilosa R. tricocca S.procumbens S. oteroi S. tenu folia S. dulce T. micrantha T. aura thiac T. purpurascens B. chondrosioides H. cenchroides 101 Table S4.1. Parameter a for each species obtained from the fitted models and used for our simulations. Year 2001 2002 2003 2004 2005 2006 2007 2010 2011 2012 2013 2014 2015 2016 Aristida adscensionis 4.92 0.41 0.75 2.59 3.39 0.40 2.05 1.36 -1.31 0.51 2.37 2.28 0.72 2.75 Aristida divaricata 0.20 2.03 0.89 1.55 -0.35 1.37 0.88 Bouteloua hirsuta 1.23 2.30 -0.54 0.71 -1.61 1.63 -1.68 0.79 0.36 -0.01 0.17 0.00 -0.14 0.04 Bouteloua scorpioides 1.55 1.70 0.88 2.86 1.04 0.53 1.98 Crusea diversifolia 0.33 1.42 -0.93 0.73 0.41 -2.39 2.15 0.10 0.24 -0.30 0.60 -1.60 0.13 0.94 Cyperus seslerioides -0.14 0.34 -0.42 -0.55 -1.77 0.60 1.26 Dalea humilis -0.01 0.17 -0.95 -1.08 -0.82 -0.61 0.87 Euphorbia mendezii 0.08 -1.66 0.81 1.65 -0.39 -0.30 1.81 Evolvolus sericeus -10.28 -0.73 0.17 -3.59 -0.13 -3.06 -1.19 -1.38 -1.86 -1.09 -1.84 Florestina purpurea -0.08 -0.27 0.12 -1.16 1.02 0.73 2.10 0.60 0.31 0.42 0.97 0.27 0.24 1.18 Heliotropium foliossisimum -3.75 -1.25 -0.98 -2.70 -1.59 -3.93 -1.47 Helianthemum glomeratum -1.47 -0.21 -3.36 0.68 -1.09 -1.70 0.47 Heterosperma pinnatum -0.58 -1.33 -0.16 -2.45 -1.44 -0.37 1.83 -0.29 -0.13 0.59 0.53 -0.44 0.29 1.15 Microchloa kunthii 0.55 0.30 -0.71 0.51 0.11 -0.47 -0.22 0.50 0.12 -0.02 0.23 -0.16 -0.50 0.41 Milla biflora 0.72 0.27 -0.17 0.31 -0.11 0.11 0.46 Muhlenbergia peruviana 1.85 0.74 0.16 0.85 -1.47 -1.38 3.04 1.53 1.20 -0.19 1.79 -0.69 -2.35 0.44 Muhlenbergia phalaroides -1.51 -0.24 -1.86 -0.10 0.85 0.67 -0.05 Muhlenbergia rigida -0.88 -1.89 -1.43 -2.83 -1.65 -2.96 -4.86 Oxalis lunulata -0.04 1.11 1.17 -0.21 -0.26 1.71 0.73 Phemeranthus oligospermus 0.10 2.43 -2.16 2.41 -0.14 -0.37 -0.13 -0.04 -2.15 1.09 1.16 Plantago nivea -0.51 -0.47 -0.40 -1.59 -1.56 -1.50 1.84 -0.03 -0.58 -1.48 -0.90 -1.44 -2.23 0.22 Porophyllum pilosa 2.18 1.12 -1.65 -2.31 -0.55 -2.84 -0.35 -1.75 -1.52 0.17 0.66 -0.70 -0.88 0.01 Richardia tricocca -0.46 -0.01 -2.65 -2.91 -0.62 -0.56 0.26 0.14 -0.38 -3.82 0.17 -1.55 -1.24 -0.01 102 Sanvitalia procumbens 1.14 2.00 -1.40 -0.85 0.29 1.70 3.27 -0.20 0.18 -0.54 1.58 0.29 1.18 1.50 Sedum oteroi 0.39 0.79 -0.01 0.63 1.18 -0.15 0.64 Sporobolus tenuissimus 1.63 1.23 2.22 1.11 2.17 -1.00 1.01 0.40 1.45 0.13 1.25 0.44 0.49 1.60 Stenandrium dulce -0.41 -0.35 -0.81 -1.33 -0.40 -2.85 -0.45 Stevia ephemera -2.95 3.54 0.93 -3.33 0.21 -2.88 1.28 0.60 0.23 0.25 0.53 -0.87 -0.45 0.91 Tagetes micrantha -0.38 0.02 -0.89 -1.41 0.33 -3.73 -0.58 0.59 0.39 -0.91 0.05 -2.18 -0.58 -0.65 Thymophylla aurantiaca -0.06 0.35 0.78 0.00 0.42 0.35 2.96 -0.58 0.80 -0.74 1.37 -0.23 0.93 0.56 Tridax coronopifolia -1.17 -0.57 -1.34 -0.49 1.60 0.62 1.22 -1.69 0.12 -0.91 -0.42 -0.33 0.57 1.05 Tripogandra purpurascens 0.10 -0.29 -2.39 -6.62 3.36 -14.15 -0.55 2.01 -1.56 -2.63 -0.44 -3.15 -0.80 -3.38 Tripogon spicatus 0.10 0.64 -0.51 -0.04 -1.32 -0.11 -0.09 Bouteloua chondrosioides -1.86 -2.62 -0.95 -2.39 -1.13 -2.68 0.15 -0.71 -0.76 -0.62 -2.83 -0.70 -1.54 -1.84 Hilaria cenchroides -5.07 -23.96 -18.94 -20.43 -15.47 0.24 -23.51 0.41 2.19 -1.11 -3.17 -1.49 -3.19 -2.72 103 Table S4.2. Parameter b for each species obtained from the fitted models and used for our simulations. Year 2000 2001 2002 2003 2004 2005 2006 2010 2011 2012 2013 2014 2015 2016 Aristida adscensionis -0.08 0.37 0.52 -0.24 -0.24 0.32 -0.05 0.32 0.63 0.86 0.00 0.10 0.44 -0.30 Aristida divaricata 0.46 0.27 0.27 0.06 0.60 0.31 0.26 Bouteloua hirsuta -0.36 -0.35 0.83 -0.18 0.73 -0.38 0.98 0.31 -0.08 0.39 0.15 0.08 0.39 0.10 Bouteloua scorpioides 0.59 -0.01 0.68 -0.20 0.60 0.72 -0.02 Crusea diversifolia 0.11 -0.23 0.78 -0.13 0.01 0.33 -0.13 0.02 0.22 0.14 0.05 0.52 0.15 0.05 Cyperus seslerioides 0.35 0.21 0.54 0.27 0.37 0.58 -0.05 Dalea humilis 0.16 -0.10 0.48 0.64 0.29 0.50 0.14 Euphorbia mendezii 0.13 0.59 0.34 0.06 0.38 0.32 0.08 Evolvolus sericeus 3.32 0.40 0.54 1.09 0.40 1.05 0.32 0.48 0.78 0.34 0.94 Florestina purpurea 0.81 0.30 -0.31 0.68 0.03 0.45 -0.12 0.34 0.21 0.24 0.56 -0.15 0.78 0.24 Heliotropium foliossisimum 0.93 0.20 0.32 0.38 0.23 0.89 0.44 Helianthemum glomeratum -0.06 -0.02 0.75 -0.09 0.48 0.76 0.03 Heterosperma pinnatum 1.04 0.63 0.77 0.69 1.56 0.22 0.51 0.79 0.57 0.39 0.43 0.28 0.59 0.06 Microchloa kunthii 0.03 0.57 0.55 0.01 0.37 0.55 0.39 0.17 0.25 0.27 0.22 0.27 0.61 0.01 Milla biflora -0.12 -0.20 -0.09 0.18 0.67 0.24 0.26 Muhlenbergia peruviana -0.25 0.76 0.65 -0.10 0.14 30.47 -0.86 0.15 0.16 -0.41 0.26 0.63 0.65 0.34 Muhlenbergia phalaroides 0.54 0.58 0.48 0.37 0.16 0.49 0.13 Muhlenbergia rigida 0.04 0.66 0.07 0.93 0.28 0.61 1.26 Oxalis lunulata 0.45 0.22 0.21 0.42 0.06 0.36 0.38 Phemeranthus oligospermus -0.69 -0.72 0.52 -1.69 -1.48 0.74 0.16 -0.92 0.33 0.33 -0.15 Plantago nivea 0.36 0.59 0.55 0.75 0.63 0.65 -0.10 0.28 0.44 0.81 0.51 0.22 1.10 0.04 Porophyllum pilosa -81.75 0.29 1.14 1.96 0.64 222.23 0.78 0.22 0.42 0.38 -0.14 0.46 0.72 0.85 104 Richardia tricocca 0.21 0.26 1.00 1.10 0.42 0.19 0.48 0.06 0.59 1.46 0.29 0.91 0.73 0.72 Sanvitalia procumbens 0.15 -0.50 1.18 0.44 0.53 -0.07 -0.11 0.69 0.68 0.93 0.38 -0.08 0.47 0.00 Sedum oteroi -1.58 0.13 -0.30 -0.26 -0.46 0.63 -1.77 Sporobolus tenuissimus -0.19 0.29 -1.25 0.43 -0.64 0.29 0.08 0.27 -0.07 -0.03 -0.03 0.24 0.32 0.01 Stenandrium dulce 0.15 0.12 0.27 0.83 0.35 1.19 0.11 Stevia ephemera 0.70 0.48 0.01 1.01 0.19 0.41 0.15 0.51 0.51 0.12 0.49 0.09 0.39 -0.07 Tagetes micrantha 0.16 0.39 0.81 0.44 0.33 0.50 1.24 0.13 0.23 0.21 0.21 0.08 0.27 0.48 Thymophylla aurantiaca 0.93 0.68 0.18 0.27 0.18 0.10 -0.15 0.66 0.05 0.86 0.42 -0.06 0.53 0.22 Tridax coronopifolia 1.09 0.82 0.67 0.24 0.02 -0.05 0.08 1.03 0.08 0.94 0.86 -0.15 0.05 0.37 Tripogandra purpurascens -0.06 0.57 0.33 2.74 -0.43 78.40 0.75 0.00 0.79 0.14 1.10 0.03 -0.09 1.77 Tripogon spicatus -0.12 0.07 0.32 -0.02 0.54 0.22 0.03 Bouteloua chondrosioides 0.42 0.87 0.28 0.71 0.38 0.96 0.01 0.36 0.41 0.30 0.68 0.30 0.64 0.51 Hilaria cenchroides 0.85 0.20 -0.04 0.03 -0.95 -0.01 -0.95 -0.02 -1.19 0.42 0.79 0.53 0.98 0.89 105 Table S4.3. Parameter c for each species obtained from the fitted models and used for our simulations. Year 2001 2002 2003 2004 2005 2006 2007 2010 2011 2012 2013 2014 2015 2016 Aristida adscensionis 28.60 -24.90 29.40 12.90 -1.20 -25.50 -1.50 25.20 -27.80 39.30 -20.00 0.50 0.83 -1.66 Aristida divaricata 0.50 0.48 -1.24 -0.99 -0.06 -0.42 -0.62 Bouteloua hirsuta 24.43 -35.53 32.31 -27.88 31.18 -34.25 25.01 -0.28 -2.35 0.98 -36.92 -1.59 0.30 -1.70 Bouteloua scorpioides 3.01 -5.26 5.02 -6.60 3.16 1.42 -1.81 Crusea diversifolia -47.63 -72.52 29.41 -66.69 -43.79 -40.90 -46.81 -45.63 -0.95 -27.18 -33.03 0.52 -1.80 -0.92 Cyperus seslerioides 1.04 -0.52 1.75 -0.59 -0.16 0.38 -2.33 Dalea humilis -0.50 -0.08 -1.04 1.26 0.12 0.50 -0.02 Euphorbia mendezii 0.99 0.03 1.25 -0.82 -1.07 0.15 0.54 Evolvolus sericeus -0.24 0.77 2.14 0.43 0.13 0.72 -2.00 -1.33 1.21 -4.71 3.14 Florestina purpurea 37.85 -43.86 43.70 -1.60 -42.92 43.77 -40.22 0.27 36.43 0.55 2.61 -38.44 44.53 0.34 Heliotropium foliossisimum 0.52 -0.33 0.79 -0.98 -0.14 -0.06 0.51 Helianthemum glomeratum -290.73 -233.12 278.06 -253.17 1.27 192.34 -205.07 Heterosperma pinnatum 13.62 -20.15 19.66 -1.01 9.16 -26.10 21.05 17.35 17.54 -0.34 -0.45 -0.66 1.36 -17.30 Microchloa kunthii -0.87 2.60 1.74 -8.72 0.18 12.73 1.29 -1.24 0.57 0.03 -0.91 -0.29 11.52 -1.32 Milla biflora 2.10 -4.57 -5.31 0.95 4.57 0.98 -2.95 Muhlenbergia peruviana -26.03 19.40 24.91 22.26 22.42 -24.78 22.57 -20.77 -30.68 -31.93 22.31 1.78 -24.34 -2.42 Muhlenbergia phalaroides -1.19 0.49 -1.89 -1.59 -0.95 2.90 0.23 Muhlenbergia rigida -9.58 22.40 4.96 16.51 -1.21 2.30 14.33 Oxalis lunulata 0.08 -0.60 -0.54 0.36 -1.23 0.65 0.24 Phemeranthus oligospermus -0.07 1.05 0.31 -0.39 0.68 0.50 0.39 -0.24 -0.58 -0.30 0.26 Plantago nivea 0.22 0.57 34.58 28.57 1.19 1.75 -22.88 1.00 -0.17 47.30 -0.23 -2.52 30.53 -31.27 Porophyllum pilosa -20.85 35.80 97.24 97.94 -19.93 -97.04 18.93 -70.24 1.13 58.23 42.32 111.30 120.39 64.93 Richardia tricocca -16.77 -19.69 20.02 1.55 -1.06 -15.77 -0.33 -15.98 0.59 11.08 -0.95 16.07 0.89 18.23 Sanvitalia procumbens -1.72 -21.98 34.99 -24.92 29.93 -31.71 -26.18 19.85 19.26 35.44 -0.96 -31.37 37.82 -38.41 Sedum oteroi -0.88 0.39 -0.77 0.11 0.91 0.11 -0.84 106 Sporobolus tenuissimus 90.67 -89.92 -104.79 74.77 -54.65 -0.62 90.16 -113.61 -1.04 -101.73 -2.28 -0.56 77.96 -0.89 Stenandrium dulce -30.61 -56.66 -1.66 52.72 0.17 39.96 -50.64 Stevia ephemera 31.82 26.63 -0.87 0.84 -24.91 -25.02 -27.68 2.19 2.65 -27.37 0.18 -1.79 -1.08 -17.77 Tagetes micrantha -1.32 -0.95 13.69 -8.25 -0.80 -15.33 11.55 -9.54 -11.45 -5.47 -10.04 -14.22 -1.10 -0.05 Thymophylla aurantiaca 10.69 22.79 -18.50 0.39 12.88 0.24 -24.69 0.66 -0.91 22.64 0.35 15.01 25.26 -1.85 Tridax coronopifolia 43.01 3.18 44.64 -46.13 -43.13 -37.76 -40.08 2.53 -43.55 2.49 43.24 -40.27 -42.06 1.98 Tripogandra purpurascens -48.30 -47.30 84.20 37.00 -56.10 -45.90 -0.06 -110.00 -0.14 -104.00 137.00 -86.70 -109.00 1.11 Tripogon spicatus 0.66 -0.78 0.32 -0.75 0.55 1.24 -0.05 Bouteloua chondrosioides 0.16 0.55 0.04 -0.81 0.82 0.92 0.68 -0.47 -0.10 0.78 -0.65 -0.32 1.00 0.24 Hilaria cenchroides -1.09 -0.83 -0.46 0.63 -0.82 -0.12 -0.25 -0.42 0.75 0.39 -1.04 -0.46 0.23 0.65 107 Table S4.4. Parameter D (mortality) for each species obtained from the fitted models and used for our simulations. Species D Aristida adscensionis 0.80 Aristida divaricata 0.43 Bouteloua hirsuta 0.50 Bouteloua scorpioides 0.42 Crusea diversifolia 1.00 Cyperus seslerioides 0.73 Dalea humilis 1.00 Euphorbia mendezii 1.00 Evolvolus sericeus 0.65 Florestina purpurea 1.00 Heliotropium foliossisimum 0.43 Helianthemum glomeratum 0.75 Heterosperma pinnatum 1.00 Microchloa kunthii 0.68 Milla biflora 0.63 Muhlenbergia peruviana 1.00 Muhlenbergia phalaroides 0.32 Muhlenbergia rigida 0.28 Oxalis lunulata 0.81 Phemeranthus oligospermus 0.80 Plantago nivea 0.81 Porophyllum pilosa 0.93 Richardia tricocca 0.71 Sanvitalia procumbens 1.00 Sedum oteroi 0.34 Sporobolus tenuissimus 0.66 Stenandrium dulce 0.52 Stevia ephemera 1.00 Tagetes micrantha 1.00 Thymophylla aurantiaca 1.00 Tridax coronopifolia 0.91 Tripogandra purpurascens 1.00 Tripogon spicatus 0.65 Bouteloua chondrosioides 0.26 Hilaria cenchroides 0.54 108 Table S4.5. Parameter g for each species obtained from the fitted models and used for our simulations. A. adscensionis A. divaricata B. hirsuta B. scorpioides C. diversifolia C. seslerioides D. humilis E. mendezii E. sericeus F. purpurea H. foliossisimum H. glomeratum A. adscensionis 2.25 -1.31 0.49 A. divaricata 3.74 B. hirsuta 0.75 -1.54 B. scorpioides -1.17 3.58 C. diversifolia -0.68 1.86 -0.31 -0.22 C. seslerioides 1.10 1.15 -0.54 D. humilis 0.79 0.09 E. mendezii 0.08 -1.36 2.19 E. sericeus 1.22 F. purpurea -0.06 0.95 H. foliossisimum 0.84 H. glomeratum 4.25 -1.55 H. pinnatum -0.20 -33.68 1.35 -34.61 1.10 M. kunthii -18.44 -19.62 -9.71 M. biflora M. peruviana M. phalaroides M. rigida O. lunulata 0.48 P. oligospermus -0.45 P. nivea -61.23 -1.97 P. pilosa R. tricocca S. procumbens -1.25 S. oteroi S. tenuissimus -86.95 S. dulce S. ephemera -0.30 12.84 1.66 1.25 T. micrantha -18.81 -21.90 3.18 T. aurantiaca -0.45 -0.95 T. coronopifolia 1.61 T. purpurascens 1.21 T. spicatus -1.39 0.40 B. chondrosioides 0.26 0.81 0.39 -0.72 0.53 -0.78 H. cenchroides 0.91 0.96 109 H. pinnatum M. kunthii M. biflora M. peruviana M. phalaroides M. rigida O. lunulata P. oligospermus P. nivea P. pilosa R. tricocca S. procumbens A. adscensionis -30.30 -22.90 -31.30 -0.31 A. divaricata 0.42 B. hirsuta -27.80 B. scorpioides C. diversifolia -1.78 -0.70 C. seslerioides 0.65 0.30 -0.71 D. humilis -1.48 -0.44 0.54 E. mendezii -1.15 E. sericeus F. purpurea -50.84 -2.32 H. foliossisimum H. glomeratum -0.41 H. pinnatum 0.38 0.24 -26.18 M. kunthii -0.29 -16.52 -20.24 M. biflora 2.07 M. peruviana 1.12 -0.07 M. phalaroides 2.01 M. rigida -0.55 O. lunulata 0.67 1.05 -0.81 1.78 0.71 P. oligospermus 0.88 -0.67 -0.47 P. nivea -86.31 -0.90 -75.24 P. pilosa 0.35 R. tricocca 0.54 S. procumbens 1.24 S. oteroi S. tenuissimus -132.77 -1.32 -2.34 S. dulce -30.35 0.20 S. ephemera -0.36 1.75 1.43 -24.58 -29.69 -24.84 -2.47 T. micrantha -1.28 -17.14 0.08 -19.01 T. aurantiaca -22.26 -1.67 -2.41 -23.37 T. coronopifolia T. purpurascens 0.76 0.80 T. spicatus -0.32 0.45 B. chondrosioides 0.10 -0.03 -0.20 H. cenchroides 1.07 110 S. oteroi S. tenuissimus S. dulce S. ephemera T. micrantha T. aurantiaca T. coronopifolia T. purpurascens T. spicatus B. chondrosioides H. cenchroides Multispecies A. adscensionis -40.60 -33.90 -0.39 -20.60 -32.70 A. divaricata -0.30 B. hirsuta 0.09 -37.53 B. scorpioides -0.22 -5.02 C. diversifolia -0.86 -49.15 3.46 0.82 -49.67 C. seslerioides -0.15 -1.35 -0.32 -0.57 0.03 -3.21 D. humilis -0.41 0.30 -1.87 E. mendezii 1.49 1.37 -0.61 -1.95 E. sericeus -3.18 F. purpurea -51.14 -2.85 -45.87 -43.10 H. foliossisimum -0.19 -1.44 H. glomeratum 0.20 H. pinnatum -0.93 -28.53 -4.01 -8.40 -0.14 -22.32 M. kunthii -20.21 -20.69 -21.12 -19.37 -19.94 M. biflora -0.75 -5.72 M. peruviana -24.84 -3.86 M. phalaroides -4.35 M. rigida 0.66 O. lunulata 0.84 -0.14 -0.93 P. oligospermus 0.16 P. nivea -1.69 -2.00 -1.26 -2.70 -77.65 P. pilosa -2.58 R. tricocca 0.33 -3.38 -21.51 -15.77 -9.31 S. procumbens -1.40 -41.80 S. oteroi 2.71 -0.11 S. tenuissimus 1.09 -1.33 -1.17 -127.69 S. dulce 1.60 -0.52 S. ephemera 1.04 -28.42 -43.41 -0.99 -2.78 -28.60 1.77 -28.12 T. micrantha 0.33 -0.07 -1.11 -13.79 -16.42 T. aurantiaca -24.57 -0.94 0.75 0.43 1.59 10.25 -24.01 T. coronopifolia 1.52 -0.11 -55.79 -55.00 T. purpurascens -114.00 -0.23 0.23 T. spicatus -0.05 0.91 1.83 0.92 -0.10 B. chondrosioides 0.39 0.24 -0.79 -1.19 H. cenchroides 1.52 -0.82 111 Table S4.6. Parameter h for each species obtained from the fitted models and used for our simulations. A. adscensionis A. divaricata B. hirsuta B. scorpioides C. diversifolia C. seslerioides D. humilis E. mendezii E. sericeus F. purpurea H. foliossisimum H. glomeratum A. adscensionis 0.00 0.00 A. divaricata B. hirsuta 0.00 B. scorpioides 0.00 C. diversifolia -0.12 0.00 0.00 C. seslerioides 0.00 0.14 D. humilis 0.00 E. mendezii 0.00 0.00 E. sericeus F. purpurea 0.00 H. foliossisimum H. glomeratum -0.12 H. pinnatum 0.00 0.00 -0.05 0.00 0.00 M. kunthii 0.00 0.00 0.00 M. biflora M. peruviana M. phalaroides M. rigida O. lunulata 0.00 P. oligospermus 0.00 P. nivea 0.00 0.00 P. pilosa R. tricocca S. procumbens 0.00 S. oteroi S. tenuissimus 0.00 S. dulce S. ephemera 0.00 -2.31 0.00 0.00 T. micrantha 0.00 0.00 -0.10 T. aurantiaca 0.00 0.00 T. coronopifolia 0.00 T. purpurascens 0.00 T. spicatus 0.00 0.00 B. chondrosioides 0.00 0.00 0.00 0.00 0.00 0.00 H. cenchroides 0.00 0.00 112 H. pinnatum M. kunthii M. biflora M. peruviana M. phalaroides M. rigida O. lunulata P. oligospermus P. nivea P. pilosa R. tricocca S. procumbens A. adscensionis 0.00 0.00 0.00 0.00 A. divaricata 0.00 B. hirsuta -7.95 B. scorpioides C. diversifolia 0.00 0.12 C. seslerioides 0.00 -0.09 0.00 D. humilis 0.00 0.00 0.00 E. mendezii 0.00 E. sericeus F. purpurea 0.00 0.00 H. foliossisimum H. glomeratum -0.04 H. pinnatum 0.05 0.00 M. kunthii 0.00 0.00 M. biflora M. peruviana 0.00 M. phalaroides M. rigida O. lunulata -0.15 0.00 0.00 0.00 P. oligospermus 0.00 0.00 P. nivea -4.68 0.00 P. pilosa R. tricocca S. procumbens S. oteroi S. tenuissimus 0.00 0.00 0.00 S. dulce -20.47 0.00 S. ephemera -0.18 -0.11 0.00 0.00 0.00 0.00 0.00 T. micrantha 0.15 0.00 0.00 0.00 T. aurantiaca 0.00 0.00 0.00 0.00 T. coronopifolia T. purpurascens 0.00 0.00 T. spicatus 0.00 0.00 B. chondrosioides 0.00 0.00 0.00 H. cenchroides 0.00 113 S. oteroi S. tenuissimus S. dulce S. ephemera T. micrantha T. aurantiaca T. coronopifolia T. purpurascens T. spicatus B. chondrosioides H. cenchroides A. adscensionis 0.00 0.00 0.00 0.00 A. divaricata B. hirsuta 0.00 B. scorpioides 0.00 C. diversifolia 0.03 0.00 -0.40 0.00 C. seslerioides 0.00 0.00 0.00 0.00 0.00 D. humilis 0.00 0.00 E. mendezii -0.24 0.00 0.00 E. sericeus F. purpurea 0.00 0.07 0.00 H. foliossisimum 0.00 H. glomeratum H. pinnatum 0.00 0.00 0.13 0.23 0.00 M. kunthii 0.00 0.00 0.00 0.00 M. biflora 0.00 M. peruviana 0.00 M. phalaroides M. rigida O. lunulata 0.00 0.00 P. oligospermus P. nivea -0.01 0.00 0.00 0.15 P. pilosa R. tricocca 0.00 0.00 0.29 0.68 S. procumbens 0.00 S. oteroi S. tenuissimus 0.00 0.00 S. dulce S. ephemera 0.00 0.00 0.00 0.00 0.00 0.00 T. micrantha 0.03 0.00 0.00 T. aurantiaca 0.00 -0.03 0.00 0.00 -4.40 T. coronopifolia 0.00 0.00 T. purpurascens 0.00 T. spicatus 0.00 -0.15 0.00 B. chondrosioides 0.00 0.00 H. cenchroides 114 Table S4.6. Parameter ß for each species obtained from the fitted models and used for our simulations. A. adscensionis A. divaricata B. hirsuta B. scorpioides C. diversifolia C. seslerioides D. humilis E. mendezii E. sericeus F. purpurea H. foliossisimum H. glomeratum A. adscensionis -1.18 -19.60 A. divaricata B. hirsuta -2.28 B. scorpioides 4.62 C. diversifolia -1.94 0.44 -69.85 C. seslerioides -2.09 -1.08 D. humilis -1.68 E. mendezii -2.89 0.13 E. sericeus F. purpurea -1.81 H. foliossisimum H. glomeratum 1.55 H. pinnatum -1.18 -1.87 -23.72 -1.75 28.35 M. kunthii -2.17 -14.34 -17.10 M. biflora M. peruviana M. phalaroides M. rigida O. lunulata -0.92 P. oligospermus -1.41 P. nivea -1.37 -37.11 P. pilosa R. tricocca S. procumbens -1.76 S. oteroi S. tenuissimus -2.23 S. dulce S. ephemera -24.91 -24.14 -20.32 -23.85 T. micrantha -1.94 -14.42 1.73 T. aurantiaca -19.23 -23.17 T. coronopifolia 45.17 T. purpurascens -81.70 T. spicatus -0.79 -0.45 B. chondrosioides 0.11 -0.50 -1.02 -0.56 -1.07 -1.57 H. cenchroides -0.87 -0.45 115 H. pinnatum M. kunthii M. biflora M. peruviana M. phalaroides M. rigida O. lunulata P. oligospermus P. nivea P. pilosa R. tricocca S. procumbens A. adscensionis -2.33 -23.90 -1.88 0.57 A. divaricata 0.27 B. hirsuta -1.58 B. scorpioides C. diversifolia -0.04 55.68 C. seslerioides -0.90 -1.85 -1.29 D. humilis -1.68 -0.45 -0.39 E. mendezii -1.07 E. sericeus F. purpurea -0.95 -3.60 H. foliossisimum H. glomeratum 0.01 H. pinnatum -29.76 -3.11 M. kunthii -16.81 -2.90 M. biflora M. peruviana -0.51 M. phalaroides M. rigida O. lunulata -1.99 -1.57 -0.48 -1.26 P. oligospermus 0.23 -0.66 P. nivea -1.99 -1.51 P. pilosa R. tricocca S. procumbens S. oteroi S. tenuissimus -1.82 0.60 -1.28 S. dulce -1.41 0.38 S. ephemera -27.70 1.26 -24.33 -3.15 -0.89 -0.34 -24.76 T. micrantha 0.57 -2.99 0.56 -17.19 T. aurantiaca -0.87 -4.38 -18.48 -3.21 T. coronopifolia T. purpurascens -1.72 126.00 T. spicatus -1.84 -2.13 B. chondrosioides 0.63 0.04 -1.06 H. cenchroides -0.96 116 S. oteroi S. tenuissimus S. dulce S. ephemera T. micrantha T. aurantiaca T. coronopifolia T. purpurascens T. spicatus B. chondrosioides H. cenchroides Multispecies A. adscensionis -2.47 -2.27 -24.60 -21.40 -3.19 A. divaricata -0.45 B. hirsuta -46.75 -2.98 B. scorpioides 2.95 -5.26 C. diversifolia -1.12 -0.81 1.08 -0.92 -85.29 C. seslerioides -0.42 -1.74 -1.39 1.23 -1.83 -2.20 D. humilis -1.12 -1.75 -1.38 E. mendezii -2.24 -2.70 -1.20 -0.31 E. sericeus -3.60 F. purpurea -2.20 -2.96 -34.72 -3.62 H. foliossisimum 0.46 -0.68 H. glomeratum 1.68 H. pinnatum -32.40 -1.44 -24.21 -2.44 -3.49 -29.06 M. kunthii -1.83 -2.36 -1.12 -15.16 -14.19 M. biflora -2.46 -1.69 M. peruviana -1.62 -3.41 M. phalaroides -2.56 M. rigida 22.78 O. lunulata 0.55 -0.63 -1.57 P. oligospermus -0.01 P. nivea -44.45 -39.97 -0.73 -0.37 -59.67 P. pilosa -1.74 R. tricocca -25.43 -19.59 0.28 0.62 -2.81 S. procumbens -27.53 -30.59 S. oteroi -0.44 S. tenuissimus -1.72 -79.79 -101.33 S. dulce -0.08 S. ephemera -2.07 -24.37 -1.00 -5.17 -2.69 23.25 -5.08 T. micrantha 0.56 -16.83 -13.73 -1.20 T. aurantiaca -1.29 -21.85 0.48 -0.16 -5.55 -2.75 T. coronopifolia 1.45 -0.76 -26.57 T. purpurascens -0.50 -0.07 T. spicatus 0.64 0.04 0.81 -0.75 B. chondrosioides -0.56 -0.15 -0.72 H. cenchroides -0.44 129 References Bray, J.R. & Curtis, J.T. (1957). An Ordination of the Upland Forest Communities of Soutther Wisconsin. Ecol. Monogr., 27, 325–349. Conde, A. & Domínguez, J. (2018). Scaling the chord and Hellinger distances in the range [0,1]: An option to consider. J. Asia-Pacific Biodivers., 11, 161–166. D’Andrea, R., Guittar, J., O’Dwyer, J.P., Figueroa, H., Wright, S.J., Condit, R., et al. (2020). Counting niches: Abundance-by-trait patterns reveal niche partitioning in a Neotropical forest. Ecology, 101, e03019. Genz, A., Bretz, F., Miwa, T., Mi, X., Leisch, F., Scheipl, F., et al (2021). mvtnorm: Multivariate Normal and t Distributions. R package version 1.1-2. URL http://CRAN.R-project.org/package=mvtnorm Oksanen, J. (2021). natto: An Extreme 'vegan' Package of Experimental Code. R package version 0.2. https://github.com/jarioksa/natto 130 Appendix 5. Species clusters in the three scenarios Figure S5.1. Species clusters along the soil depth gradient. Blue, green, orange, and pink ovals represent clusters present in the scenario with all interactions. Thick border ovals represent clusters present in the scenario without facilitation. Species with same- colored fonts represent species clusters when interactions are absent from our model. Note that species with black font are those that went extinct when we eliminated interactions and species outside of circles went extinct when all interactions were present or when facilitation was eliminated. Soil Depth (cm) Heterosperma pinnatum Tagetes micrantha Phemeranthus oligospermum Portulaca pilosa Euphorbia mendezii Florestina purpurea Aristida adscensionis Aristida divaricata Bouteloua hirsuta Bouteloua scorpioides Micrchloa kunthii Sedum oteroii Sporobolus tenuissiums Thymophylla aurantiaca Muhlenbergia peruviana Muhlenbergia phalaroides Richardia tricocca Milla biflora Sanvitalia procumbes Evolvulus sericeus Muhlenbergia rigida Stenandrium dulce Plantago nivea Hillaria cenchroides Bouteloua chondrosiodes Tripogandra purpurascens Oxalis lunulata 3 288 13 2118 Stevia ephemera Cyperus seslerioides Dalea humilis Tripogon spicatus 131 Appendix 5. Species clusters in the three scenarios and in the original field data Species clusters in the field data To analyze if clusters are also observed in the original field data, we applied the same k- means algorithm that we applied to our simulations. To obtain the average abundance of each species over the soil depth gradient we fitted general additive models (GAM) using a negative binomial distribution of the error in package mgcv for R (Wood, Pya, and Säfken 2016). The k-means algorithm indicated that seven clusters significantly differed from random in the field data (Fig 3), although many of these clusters contained very few species. A better classification was obtained using five clusters, which was marginally significant (P = 0.0618). To simplify the comparison between these clusters (hereafter referred to as “observed clusters”, OC, and identified with letters) and those obtained from our model with all of the interactions (modelled clusters, MC, which are identified by numbers as in the main text), which had four clusters, we combined two of the five OC. We specifically combined clusters C and D because they were most similar to cluster 3 (Table S5.1). We used Sørensen’s similarity index to compare the OC and MC. To determine whether OC and MC were more similar than expected by chance, we used a null model. In it, we assigned randomly the species in each MC to each of the four OC. This maintained the number of species in MC constant. We then calculated Sørensen similarity indices for the shuffled data. This was done 1000 times. We identified which clusters occurred roughly over the same depths and had similar species, and considered them to be more or less the same clusters. Thus, we considered clusters 1 and A, 2 and B, 3 and C-D, and 4 and E the be the same clusters. Similarity between cluster 3 and the pooled clusters C-D was only 36% and was not significant despite having combined cluster C and D to try to maximize similarity between OC and 132 MC. The three remaining pairs of clusters averaged a similarity of 61% and were all significant (table S5.1). These comparisons were conducted using one-tailed tests because these cluster pairs were expected to be more similar than expected by chance. We also compared the similarity between the remaining pairs of clusters. It was 10 %, and (with only one exception) it did not differ from random. These tests were conducted using one- tailed tests because we expected cluster pairs to be more different than expected by chance (Table S5.1). Table S5.1 Sørensen similarity between clusters. Clusters obtained from our models are represented by columns while clusters obtained from field data are represented by rows. Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster A 50 % (0.0362) 20 % (0.1685) 0 % (0.2340) 0 % (0.5587) Cluster B 0 % (0.3197) 67 % (0.0045) 25 % (0.6477) 0 % (0.5610) Cluster C and D 27 % (0.7312) 37 % (0.5980) 36 % (0.1340) 15 % (0.5995) Cluster E 0 % (0.3189) 0 % (0.0098) 0 % (0.2333) 67 % (0.0024) References Wood, Simon N., Natalya Pya, and Benjamin Säfken. 2016. “Smoothing Parameter and Model Selection for General Smooth Models.” Journal of the American Statistical Association 111 (516): 1548–1563. doi:10.1080/01621459.2016.1180986. 133 Appendix 6. Testing for intransitivity In the model of species alliances presented in the main text, species that belong to different clusters compete more strongly that species that occur in the same cluster (Figure S6.1 A). In this scenario, when all species in a focal cluster are removed except one (the triangle in Figure S6.1), species in other clusters are released from competition and grow more than the sole remaining species in the focal cluster, which was only weakly regulated by its cluster mates. Under such a scenario, the remaining species would diminish in size due to augmented inter-cluster competition, as observed in our results. This favors coexistence between species in a cluster, which indirectly facilitate each other. However, there is another scenario in which species may facilitate each other indirectly: intransitive competition (Figure S6.1 B and C). Under such scenario, strong competitors suppress all the species that are below them in a competitive hierarchy, but there are loops of species that resemble a rock-paper-scissors game. If the intransitive loop comprises species of different clusters we may observe a pattern similar to the one that we observed, but only in some species (Figure S6.1 B). In others, the opposite would occur, depending on their position in the loop (Figure S6.1 C). Intransitive loops within clusters (Figure S6.1 D) may also promote coexistence between species in a cluster, but the removal of all species but one would result in an increase in the abundance of the remaining species. This scenario is not consistent with our results and therefore can be discarded. 134 Figure S6.1. Indirect facilitation and transitive competition scenarios. The pink and the blue ovals represent two different clusters. Each species is represented by a different shape, the direction of the arrows indicates a negative effect on the species towards which the arrow points, and arrow thickness corresponds to interaction strength. The bottom row represents what we would observe if we were to eliminate all the species in the blue cluster except the triangle. Shapes that become larger represent an increase in the species’ population density while shapes that become smaller represent species’ whose populations diminish. The scenario A corresponds to no intransitivity, B and C show two outcomes of transitivity depending on the position of species in an intransitivity loop that involves two clusters, and D corresponds to within-cluster intransitivity. Tests for intransitivity require a single interaction matrix. We had two instead: one for facilitation containing b values and another for competition containing a values (see equation 1 in the main text). We obtained our interaction matrix A by calculating the A B C D A ll s p e ci e s A ft e r re m o v a l 135 derivative of the growth rate of the focal species i with respect to each of the density of the associate species Nj, and evaluating it at the equilibrial density !1∗: @!,1 = 1 !! VI4!! , !1 , … 5V!1 X :% ∗ where g(Ni,Nj,…) is equation 1 in the main text. This gave us a simple measure of how species growth rate changes in response to changes in the density of the associate species. Thus, a positive derivative means that there is facilitation while a negative derivative implies competition in the community at equilibrium. To assess intransitivity, we tested for nestedness in our interaction matrix A. Nestedness reflects orderly sequences of competition as would be expected if there are transitive competitive hierarchies (Ulrich et al. 2014; Laird and Schamp 2006; Ulrich et al. 2018). We used the weighted interactions nestedness estimator (WINE, Galeano, Pastor, and Iriondo 2009) because it considers the strength (or weight) of each interaction, and not just its net effect defined as which species wins or loses in the competitive interaction. Positive WINE values indicate nestedness and thus transitive, whereas negative WINE indicates non-nestedness and therefore intransitivity. We compared this estimator to a null model where we shuffled the species interaction coefficients of each focal species without shuffling the intraspecific interaction coefficient. Because our interactions change with soil depth, the process was repeated for four depths that correspond roughly to the depth were each cluster was observed: 9.9, 15.5, 19.9 and 28 cm. Nestedness did not differ from what was expected by chance in our interaction matrix for any soil depth (Table S6.1), indicating that there is no evidence for either strong transitivity or intransitivity. We also tested for nestedness within and between clusters. It was important to test for nestedness between different clusters because it is the scenario that 136 could explain our results. No matrices were significantly nested, and we only found significant non-nestedness, and thus intransitivity, in shallow soils for species in cluster 3, and in deep soil when considering interactions between species in clusters 2 and 3 (Table S6.1). It is important to note that none of these patterns were found in the range of depths that characterized these clusters, as species in cluster 2 reach their maximum density near 15.5 cm, and species in cluster 3 near 19.9 cm. In consequence, intransitivity between clusters is not likely to be responsible for the patterns that we observed in our data when removing all but one species in each cluster. Thus, intransitivity was not an important factor driving species clustering. Table S6.1 Weighted interactions nestedness estimator (WINE) values at different soil depths for all species considered simultaneously (column All) and between species occurring in different pairs of clusters. Cluster 4 had only two species, and thus there were insufficient data for most analyses that included that cluster. Values in parentheses are two-tailed P values. Depth All Cluster1 Clusters 1 and 2 Clusters 1 and 3 Cluster 2 Clusters 2 and 3 Clusters 2 and 4 Cluster 3 9.9 cm 0.020 (0.7817) -0.429 (0.3001) 0.339 (0.3548) 0.379 (0.4218) 0.044 (0.3056) 0.102 (0.2420) 0.276 (0.5402) -0.109 (0.0379) 15.5 cm 0.077 (0.7575) 0.742 (0.7010) 0.362 (0.7511) 0.354 (0.3934) 0.058 (0.4900) 0.164 (0.8521) 0.283 (0.5854) -0.213 (0.0299) 19.9 cm 0.021 (0.6905) 0.330 (0.4420) 0.338 (0.6460) 0.400 (0.5908) 0.070 (0.6631) 0.156 (0.8577) 0.286 (0.6217) 0.173 (0.2694) 28 cm 0.025 (0.8402) 0.748 (0.9551) 0.295 (0.4522) 0.341 (0.3464) 0.041 (0.1835) -0.099 (0.0171) 0.213 (0.0708) 0.211 (0.6010) References Galeano, Javier, Juan M. Pastor, and Jose M. Iriondo. 2009. “Weighted-Interaction Nestedness Estimator (WINE): A New Estimator to Calculate over Frequency Matrices.” Environmental Modelling and Software 24 (11). Elsevier Ltd: 1342–1346. doi:10.1016/j.envsoft.2009.05.014. Laird, Robert A., and Brandon S. Schamp. 2006. “Competitive Intransitivity Promotes 137 Species Coexistence.” American Naturalist 168 (2): 182–193. doi:10.1086/506259. Ulrich, Werner, Yasuhiro Kubota, Agnieszka Piernik, and Nicholas J. Gotelli. 2018. “Functional Traits and Environmental Characteristics Drive the Degree of Competitive Intransitivity in European Saltmarsh Plant Communities.” Journal of Ecology 106 (3): 865–876. doi:10.1111/1365-2745.12958. Ulrich, Werner, Santiago Soliveres, Wojciech Kryszewski, Fernando T. Maestre, and Nicholas J. Gotelli. 2014. “Matrix Models for Quantifying Competitive Intransitivity from Species Abundance Data.” Oikos 123 (9): 1057–1070. doi:10.1111/oik.01217. 138 Capítulo IV. Discusión general Utilizamos dos enfoques contrastantes, uno experimenta y uno observacional, para estudiar la importancia de un gradiente hídrico en la coexistencia de especies vegetales en un pastizal rico en especies. En el Capítulo II se estudió experimentalmente la diferenciación de nicho a lo largo del gradiente de estrés hídrico para establecer si en algún extremo del gradiente la diferenciación de nicho era mayor. En el Capítulo III, analizamos con datos observacionales de catorce años (Martorell y Freckleton 2014; Martínez-Blancas y Martorell 2018; Zepeda y Martorell 2019) si hay grupos a lo largo del gradiente hídrico, como resultado de las interacciones y exploramos qué mecanismo particular puede generar estos grupos. Ambos acercamientos brindaron información valiosa acerca del papel del ambiente, y en particular del gradiente de estés hídrico, sobre la coexistencia de las especies. A manera de resumen, nuestros hallazgos específicos en el Capítulo II sugieren que los efectos del estrés hídrico dependen de la identidad de la especie. También encontramos reducción de la biomasa de algunas especies en el extremo seco del gradiente, lo que representó una restricción ambiental, señal de estrés hídrico. Sin embargo, el extremo seco del gradiente no representó un filtro ambiental pues no restringió la supervivencia de ninguna especie. No observamos facilitación en suelos delgados. Esto apoya a la hipótesis modificada del gradiente del estrés en donde la competencia (y no la facilitación) predomina en condiciones de estrés intenso generado por la escasez de un recurso (Maestre et al. 2009). También encontramos evidencia de diferenciación de nicho, es decir, mayor competencia intraespecífica que interespecífica, en por lo menos cierta porción del gradiente para cinco de las seis especies utilizadas en nuestro experimento. Esto podría 139 indicar que las especies pueden coexistir en escalas espaciales grandes. Es decir, en áreas heterogéneas que por lo tanto presenten las condiciones adecuadas para cada especie. Sin embargo, también encontramos diferenciación de nicho en porciones similares del gradiente de profundidad de suelo para algunas especies. Esto es evidencia de que algunas combinaciones de especies pueden coexistir a escala local en profundidades de suelo similares. Finalmente, la diferenciación de nicho fue ligeramente más frecuente en suelos delgados con mayor estrés que en suelos profundos, pero esto no se observó en todas las especies. Esto sugiere que la amplitud de los nichos de las especies podría reducirse, reduciendo así el traslape de los mismos y en consecuencia permitiendo la coexistencia. Recapitulando el Capítulo III, encontramos que las especies sí se agrupan a lo largo el gradiente hídrico. Nuestros hallazgos además sugieren que las interacciones son importantes para este agrupamiento. Incluso, la presencia de la competencia parece generar cierta segregación de las especies a lo largo del gradiente hídrico. Esto sugiere que la competencia genera una reducción en el traslape de nichos en el sitio. Adicionalmente, la competencia interespecífica es más débil al interior de los grupos que entre ellos. Esto sugiere que los nichos de estas especies están diferenciados en otros ejes, además del gradiente hídrico, que estabiliza a la coexistencia de especies dentro de los grupos. De otra manera, si estos grupos fueran neutrales, hubiéramos observado una fuerte competencia interespecífica debido a las grandes similitudes entre las especies. La facilitación directa no parece jugar un papel en la formación de grupos, ya que los grupos persisten aun cuando eliminamos a la facilitación de nuestras simulaciones. Además, la intensidad de la facilitación fue similar dentro de los grupos que entre grupos. Sin embargo, la facilitación indirecta parece promover alianzas entre las especies dentro de un mismo grupo. Parece ser que los miembros de estos grupos incrementan las abundancias de otras especies dentro de 140 su mismo grupo cuando estas están presentes en la comunidad. Este fue el caso de tres de cuatro de los grupos que encontramos. Los factores abióticos como la disponibilidad de agua pueden estructurar a las comunidades, ya sea porque representan una restricción ambiental (Silvertown 2004) o por modificar las interacciones entre las especies (Callaway 2007; Brooker et al. 2008). En consecuencia, el gradiente hídrico puede ser un importante eje de diferenciación de nichos para la coexistencia de especies (Silvertown 2004). Nuestros resultados confirman la importancia del gradiente hídrico para la coexistencia, pero también sugieren que otros ejes del nicho son los que permiten la coexistencia a escalas locales. En el Capítulo III encontramos que las especies forman grupos a lo largo del gradiente hídrico. Las especies que pertenecen a diferentes grupos probablemente pueden coexistir gracias a las diferencias en su tolerancia al estrés hídrico. Esto coincide con lo que observamos en el Capítulo II en cuanto a que la coexistencia de algunas especies se puede dar a escalas espaciales grandes donde encontramos diferentes profundidades de suelo. Lo que observamos en el capítulo, II en cuanto a que las especies también pueden coexistir a escalas espaciales pequeñas, coincide con las especies que encontramos que coexisten dentro de un mismo grupo en el Capítulo III. Sin embargo, a esta escala las especies requieren de diferencias en otros ejes de su nicho, además el hídrico, para poder coexistir. Algunos de los ejes del nicho que podrían estabilizar la coexistencia en este pastizal a escalas locales son los enemigos naturales como patógenos del suelo (Martorell et al. 2021) o insectos granívoros (García-Meza et al. 2021). Diferenciación de nichos a lo largo del gradiente de estrés hídrico En cuanto a las preguntas específicas al Capítulo II, nuestros resultados sugieren que la diferenciación de los nichos puede ser más frecuente en ambientes estresantes. Sin embargo, 141 el número de especies en nuestro experimento fue pequeño y dos de las cinco especies que difieren en su nicho mostraron el patrón contrario, pues no tuvieron diferenciación de nicho a mayor estrés. No obstante, es posible que la naturaleza ambigua de nuestros resultados no sea solamente consecuencia de nuestro número pequeño de especies. A la fecha ya hay varios estudios que no encuentran una relación clara entre el estrés y la diferenciación de nichos o la igualación de las capacidades competitivas (Bimler et al. 2018; Germain et al. 2018; Wainwright et al. 2018). Por lo tanto, posiblemente no exista un patrón en la naturaleza respecto a la intensidad del estrés y la diferenciación de nichos y esto puede deberse a que las interacciones intra- e interespecíficas de las especies responden de maneras variadas al estrés. Quizás la teoría del gradiente del estrés y la teoría modificada del gradiente de estrés también sean evidencia de esto. La teoría del gradiente de estrés (Bertness y Callaway 1994; Callaway, 2007; Brooker et al., 2008) establece que predominan las interacciones positivas con estrés. Por otro lado, unos hallazgos contrastantes inspiraron la teoría modificada del gradiente de estrés (Tielbörger y Kadmon 2000; Maestre y Cortina 2004), la cual establece que la competencia se vuelve preponderante a mayor estrés, particularmente cuando la escasez de un recurso es lo que está causando el estrés. Nosotros, tanto en el Capítulo I como en el Capítulo II, encontramos que estos dos escenarios se presentan en nuestro sistema de estudio. Es decir, algunas especies experimentan menor competencia e incluso facilitación a mayor estrés mientras que otras especies experimentan el patrón contrario (mayor competencia a mayor estrés). Esto no permite establecer una direccionalidad general en cuanto a la prevalencia de la diferenciación de nichos (la cual depende del balance entre la competencia intra e interespecífica) en cualquiera de los dos extremos del gradiente de estrés. Asimismo, no encontramos evidencia de que los filtros ambientales estuvieran operando en el sitio de estudio con las especies que utilizamos. A pesar del estrés hídrico 142 intenso que representan los suelos delgados de este pastizal, todas las especies lograron establecerse aun en los suelos más delgados y secos. Sin embargo, los suelos delgados sí produjeron una reducción fuerte en la fecundidad y en la biomasa de la mayoría de las especies. Por lo tanto, los suelos delgados sí representan una restricción en el desempeño de las especies de parte del ambiente. Algunas especies no mostraron una reducción de fecundidad ni biomasa. Esto se puede deber a que algunas especies pueden variar en su susceptibilidad a las diferentes de formas de estrés (Grime 1977). Grupos de especies a lo largo del gradiente de estrés hídrico En el Capítulo III nos preguntamos si hay grupos de especies a lo largo del gradiente hídrico de la comunidad del pastizal y cuáles qué procesos generan estos grupos. Encontramos que sí hay grupos de especies y que las interacciones entre especies propician el agrupamiento de especies. Anteriormente, un gran número de estudios señalaban a la neutralidad emergente como un mecanismo de agrupamiento (Scheffer y van Nes 2006; Vergnon et al. 2009, 2012; Segura et al. 2011; Scheffer et al. 2015; Sakavara et al. 2018). Sin embargo, algunos estudios también señalaban que ciertos nichos ignorados en los estudios de neutralidad emergente, llamados nichos ocultos, podrían estar estabilizando a la coexistencia dentro de los grupos (Barabás et al. 2013; D’Andrea et al. 2020). Nuestros resultados son consistentes con la hipótesis de los nichos ocultos. Además, encontramos que la facilitación indirecta también es un mecanismo que propicia al agrupamiento de las especies en esta comunidad. La facilitación indirecta ocurre cuando las especies dentro de un mismo grupo suprimen competitivamente a las especies de los otros grupos, facilitando así a los miembros de su mismo grupo. Se ha señalado, además, que esto se propicia cuando las especies no compiten por los mismos recursos (Levine 1999), es decir, cuando 143 sus nichos se diferencian en otros ejes, tal como encontramos nosotros. A pesar de que la facilitación indirecta intratrófica se había observado previamente en hierbas (Levine 1999; Callaway y Pennings 2000), nunca se había relacionado con la formación de grupos dentro de las comunidades. Ventajas y desventajas de los datos experimentales y los datos observacionales En esta tesis se evaluó a la coexistencia de especies a lo largo de un gradiente hídrico con datos experimentales (Capítulo II) y con datos observacionales (Capítulo III). Ambos acercamientos nos ayudaron a entender cómo pueden coexistir un número enorme de especies en el pastizal semiárido de Concepción Buenavista (Martorell et al. 2017); sin embargo, cada uno tiene sus ventajas y sus desventajas. En los estudios experimentales podemos manipular algunas de las condiciones ambientales que queremos poner a prueba (Freckleton y Watkinson 2000; Tuck et al. 2018), aislando una o más variables como tratamientos experimentales. Sin embargo, los estudios experimentales suelen ser más breves ya que es difícil mantenerlos por periodos largo. Por esta razón, pueden no reflejar los efectos de largo plazo que experimentan las plantas en condiciones naturales que sí se ven reflejados en estudio observacionales de largo plazo (Diamond 1983; Adler et al. 2018). Un ejemplo de estos procesos de largo plazo son la retroalimentación suelo-planta que se ha visto que tiene una importante repercusión sobre la supervivencia de los individuos de este pastizal (Martorell et al. 2021). Otra desventaja que pueden tener los estudios experimentales es que frecuentemente hay que eliminar parcial o totalmente a la vegetación original del sitio para preparar el experimento. Esto puede alterar la estructura física el suelo (Freckleton y Watkinson 2000), lo que puede modificar los resultados y por lo tanto las conclusiones de 144 los estudios. En contraste, los estudios con datos observacionales no tienen estas desventajas y se ha visto que arrojan buenas estimaciones de las interacciones, siempre y cuando se tome en cuenta el movimiento de semillas y se realicen a la escala espacial correcta (Freckleton y Watkinson 2000). A pesar de esto, también son controversiales, ya que se ha encontrado que pueden llevar a una subestimación de la competencia interespecífica de las especies (Adler et al. 2018; Tuck et al. 2018). En estudios de coexistencia de especies, esto puede representar una sobreestimación de la diferenciación de los nichos. Se ha argumentado que la subestimación de la competencia interespecífica ocurre porque en las comunidades naturales las especies podrían ya estar espacialmente confinadas por la competencia a sitios donde no compiten tan fuertemente con sus vecinos (Tuck et al. 2018). A este fenómeno se le ha llamado el fantasma de la competencia pasada. En particular, nuestro estudio con datos experimentales tuvo la ventaja de que pudimos manipular la profundidad del suelo, y por lo tanto la disponibilidad de agua, al medir la competencia intra- e interespecífica. Esto nos permitió calcular la diferenciación de los nichos entre cada una de las especies y las demás especies del experimento. Casi todas las especies tuvieron diferenciación de su nicho respecto al resto en al menos algún intervalo de profundidad de suelo. Esto demostró que, a pesar de ser un estudio experimental donde la competencia interespecífica suele ser más fuerte, la regulación intraespecífica fue suficientemente fuerte para estabilizar la coexistencia. Queda pendiente establecer si la diferenciación de los nichos de estas especies es suficiente para sobreponerse a las diferencias en las capacidades competitivas entre ellas para poder coexistir establemente. Como fue un experimento de tan solo algunos meses, durante la época de lluvias, también quedó pendiente establecer algunos efectos que se presentan a largo plazo. Un ejemplo de este tipo de efectos es la retroalimentación suelo planta y puede afectar tanto a la diferenciación de los nichos de 145 las plantas como a las diferencias en sus capacidades competitivas (Yan et al. 2022). Sin embargo, el que diferentes especies tuvieran nichos diferentes de las otras en diferentes porciones del gradiente de profundidad de suelo fue un importante antecedente para el capítulo subsecuente, donde nos abocamos a determinar si las especies se agrupaban en diferentes porciones del gradiente como resultado de la competencia. Nuestro estudio con datos observacionales nos permitió trabajar con muchas más especies, lo cual era necesario para establecer si se forman grupos de especies a lo largo del gradiente de profundidad del suelo. Un experimento en el campo para calcular interacciones entre el mismo número de especies hubiera resultado prohibitivo debido al gran número de posibles combinaciones de especies (Diamond 1983). Por esta razón, los datos observacionales son esenciales en estudios donde están involucradas un gran número de especies. Respecto a la desconfianza que existe hacia la estimación de la competencia interespecífica con datos observacionales, nuestros hallazgos son consistentes con estudios previos en el sitio (Villarreal-Barajas y Martorell 2009) y con los resultados experimentales de esta tesis donde frecuentemente en el extremo más estresante del gradiente de profundidad de suelo incrementa la competencia (hipótesis del gradiente de estes modificado, Maestre y Cortina 2004; Tielbörger y Kadmon 2000). Además, nuestras simulaciones se asemejan a la distribución de las especies a lo largo del gradiente de estrés hídrico previamente observado por Pedraza y Martorell (2019). Por lo tanto, nuestros datos observacionales no solo fueron de gran utilidad para entender a la formación de grupos en las comunidades, sino que también parecen ser confiables para la estimación de interacciones. A pesar de que ambos métodos que utilizamos en esta tesis presentan desventajas para medir correctamente la intensidad de las interacciones, las conclusiones a las que llegamos no se contraponen. Con ambos enfoques, la competencia tiende a ser mayor en suelos 146 delgados, un resultado que ya se había observado previamente en otros experimentos en este sitio (Villarreal-Barajas y Martorell, 2009; Martorell et al. 2015). Además, ambos estudios coinciden en que la diferenciación de nichos ocurre a dos escalas: 1) escala local donde la diferenciación de nicho entre las especies parece ocurrir en otros ejes del nicho y 2) a escalas espaciales mayores donde las especies se diferencían en sus requerimientos hídricos. Por lo tanto, haber utilizado dos métodos distintos y haber llegado a conclusiones similares de robustez a los hallazgos de la tesis. Un exceso de diferenciación de los nichos Un exceso de diferenciación de los nichos entre las especies que coexisten en una comunidad es un resultado frecuente en estudios de coexistencia de especies (Turnbull et al. 2005; Levine y HilleRisLambers 2009; Adler et al. 2010; Chu y Adler 2015; Buche et al. 2022). Al decir esto, nos referimos a que la diferenciación de nicho es mucho más de la necesaria para contrarrestar las diferencias en las capacidades competitivas de las especies que corresponden a los mecanismos igualadores (Chu y Adler 2015). Esto significa que la competencia intraespecífica es tanto más fuerte que la interespecífica, que a pesar de que una especie tenga habilidades competitivas muy fuertes, ésta no desplazaría competitivamente al resto de las especies. Nosotros encontramos algo similar a pesar de que no medimos los mecanismos igualadores. En el Capítulo III observamos que es muy probable que otros ejes del nicho, además de la profundidad de suelo, estén dando estabilidad a las especies dentro de un mismo grupo. También se evidenció en el Capítulo II, donde se encontraron diferencias de nicho para cinco de las seis especies en el estudio. Si bien las diferencias de los nichos que encontramos no son suficientes para poder decir que hay coexistencia estable entre dos especies, sí son un requisito para pueda haber coexistencia entre especies (Broekman et al. 147 2019). Un exceso de diferenciación de nichos ya se había observado previamente en la comunidad estudiada (Zepeda y Martorell 2019). Por lo tanto, nuestro trabajo reitera la importancia de la diferenciación de los nichos para la coexistencia en el sitio. Sin embargo, la mayoría de los estudios dentro el marco de teoría moderna de la coexistencia se ha realizado en plantas y particular con hierbas anuales (Buche et al. 2022). Por lo tanto, aunque nosotros incorporamos también hierbas perenes, aún queda mucho que aprender sobre la coexistencia entre otros tipos de organismos. Ya que, al menos en las comunidades vegetales, la existencia de la diferenciación de nicho entre las especies ser generalizada, actualmente existe un esfuerzo para discernir en qué ejes del nicho está ocurriendo esta diferenciación. Nuestro trabajo aunado a otros trabajos (Wainwright et al. 2018; Bimler et al. 2018; Germain, Mayfield y Gilbert 2018) apuntan a que la diferenciación de nicho no está relacionada a la disponibilidad de agua. Tampoco parece depender de las fluctuaciones ambientales (Zepeda y Martorell 2019) ni de los enemigos naturales que se encuentran en el suelo (Yan et al. 2022). Por lo tanto, queda por explorar de donde proviene esta diferenciación de los nichos. La presencia de enemigos naturales, como patógenos del suelo o herbívoros, o la dispersión y movimiento de semillas y propágulos son candidatos interesantes para futuros trabajos que tengan como objetivo averiguar cuáles son los ejes del nicho donde ocurre esta diferenciación entre las especies que permita la coexistencia estable. Conclusiones Aunque no encontramos una tendencia fuerte de mayor diferenciación de nichos a mayor estrés hídrico, el Capítulo II reveló que las especies estudiadas tienen diferencias en sus respectivos nichos que pueden contribuir a la coexistencia estable. Las especies presentaron 148 diferenciación en diferentes porciones del gradiente, pero algunas de ellas contaron con esta diferenciación en porciones similares del gradiente. Esto nos dio pistas de lo que confirmamos en el Capítulo III, en el que encontramos que sí se forman grupos de especies a lo largo del gradiente hídrico. Sin embargo, a pesar de la importancia del gradiente hídrico para la comunidad del pastizal semiárido, nuestras simulaciones apuntan a que este no es el único eje que permite la coexistencia al interior de los grupos. Los grupos que observamos no se mantienen por sus similitudes, como sucedería en una comunidad neutral, sino que más bien otras dimensiones de sus nichos que no fueron abordadas en este estudio son las responsables de estabilizar la coexistencia entre las especies dentro del mismo grupo. Además de esta diferenciación de sus nichos, las especies dentro de un mismo grupo se facilitan indirectamente al suprimir a las especies de otros grupos, reforzando la cohesión de las especies al interior de su grupo. Este tipo de facilitación indirecta puede ser propiciada cuando las especies tiene diferentes requerimientos de recursos, situación que no sucedería si el nicho hidrológico fuera el único factor que estabiliza la coexistencia de esta comunidad. Por lo tanto, podemos concluir que a pesar de la importancia del nicho hidrológico de las plantas y en particular para esta comunidad, la coexistencia estable entre ellas probablemente no sería posible sin las diferencias de sus nichos en ejes adicionales que estabilizan la coexistencia. Nuestro estudio se suma al gran número de trabajos que ya habían encontrado un exceso de diferenciación de los nichos estabilizando a las comunidades de plantas (Levine y HilleRisLambers 2009; Adler et al. 2010; Chu y Adler 2015; Zepeda y Martorell 2019; Buche et al. 2022). Queda pendiente averiguar en qué ejes del nicho está ocurriendo esta diferenciación de nicho entre las especies. El efecto de enemigos naturales podría ser crucial en la diferenciación de nicho (Martorell et al. 2021) por lo que sería util implementar 149 experimentos donde se limite o se excluya la presencia de enemigos naturales incluyendo patógenos en el suelo. En la ecología de comunidades es particularmente dificil establecer patrones ya que las comunidades están compuestas por una multitud de diferentes especies (Lawton 1999). Sin embargo parece ser menos complicado establecer la presencia de ciertos patrones en la ecología de poblaciones a pesar de que existen variaciones entre especies e incluso entre poblaciones de la misma especie (Lawton 1999). La teoría moderna de la coexistencia nos permite estudiar a las comunidades a partir de las dinamicas poblacionales de las diferentes especies que las componen (Chesson 2000). Esto puede representar una ventaja para tratar de buscar grandes patrones en la ecología de la comunidades. La utilización de este marco teórico es cada vez común y han empezado a emergir ciertos patrones en la ecología de comunidades. Por ejemplo, parece ser bastante común un exceso de diferenciación de nichos en las comunidades (Buche et al. 2022). Este trabajo se suma a este patrón habiendo utilizaddo a la teoría moderna de la coexistencia como marco teórico de referencia con datos de campo así como datos experimentales como observacionales. Adicionalmente logramos un buen entendimiento del papel del gradiente hídrico en la coexistencia de especies de un pastizal semiárido. Referencias bibliográficas Adler PB, Ellner SP, Levine JM (2010) Coexistence of perennial plants: an embarrassment of niches. Ecol Lett 13:1019–1029. Adler PB, Smull D, Beard KH, et al. (2018) Competition and coexistence in plant communities: intraspecific competition is stronger than interspecific competition. Ecol Lett 21:1319–1329. 150 Barabás G, D’Andrea R, Rael R, Meszéna G, Ostling A (2013) Emergent neutrality or hidden niches? Oikos 122:1565–1572. Bertness MD, Callaway R (1994) Positive interactions in communities. Trends Ecol Evol 9:191–193. Bimler MD, Stouffer DB, Lai HR, Mayfield MM (2018) Accurate predictions of coexistence in natural systems require the inclusion of facilitative interactions and environmental dependency. J Ecol 106:1839–1852. Buche L, Spaak JW, Jarillo J, De Laender F (2022) Niche differences, not fitness differences, explain predicted coexistence across ecological groups. J Ecol 110:2785– 2796. Broekman MJE, Muller-Landau HC, Visser MD, Jongejans E, Wright SJ, Kroon H de (2019) Signs of stabilisation and stable coexistence. Ecol Lett 22:1957–1975. Brooker RW, Maestre FT, Callaway RM, et al. (2008) Facilitation in plant communities: the past, the present, and the future. J Ecol 96:18–34. Callaway RM (2007) Positive Interactions and Interdependence in Plant Communities. 1st ed. Dordrecht, Netherlands: Springer. Callaway RM, Pennings SC (2000) Facilitation may buffer competitive effects: Indirect and diffuse interactions among salt marsh plants. Am Nat 156:416–424. Chu C, Adler PB (2015) Large niche differences emerge at the recruitment stage to stabilize grassland coexistence. Ecol Monogr 85:373–392. D’Andrea R, Guittar J, O’Dwyer JP, et al. (2020) Counting niches: Abundance-by-trait patterns reveal niche partitioning in a Neotropical forest. Ecology 101:e03019. Diamond JM (1983) Laboratory, field and natural experiments. Nature 304:586–587. Freckleton RP, Watkinson AR (2000) On detecting and measuring competition in spatially 151 structured plant communities. Ecol Lett 3:423–432. García-Meza D, Andresen E, Ríos-Casanova L, Martorell C (2021) Seed density in monospecific and mixed patches affects individual and collective foraging in ants. Insectes Soc 68:81–92. Germain RM, Mayfield MM, Gilbert B (2018) The ‘filtering’ metaphor revisited: competition and environment jointly structure invasibility and coexistence. Biol Lett 14:20180460. Grime JP (1977) Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory. Am Nat 111:1169–1194. Lawton JH (1999) Are there general laws in ecology? Oikos 84:177–192. Levine JM (1999) Indirect facilitation: Evidence and predictions from a riparian community. Ecology 80:1762–1769. Levine JM, HilleRisLambers J (2009) The importance of niches for the maintenance of species diversity. Nature 461:254–257. Maestre FT, Callaway RM, Valladares F, Lortie CJ (2009) Refining the stress-gradient hypothesis for competition and facilitation in plant communities. J Ecol 97:199–205. Maestre FT, Cortina J (2004) Do positive interactions increase with abiotic stress? A test from a semi-arid steppe. Proc R Soc B 271:S331–S333. Martínez-Blancas A, Paz H, Salazar GA, Martorell C (2018). Related plant species respond similarly to chronic anthropogenic disturbance: Implications for conservation decision-making. J Appl Ecol 55:1860–1870. Martorell C, Almanza-Celis CAI, Pérez-García EA, Sánchez-Ken JG, Perez-Garcia EA, Sanchez-Ken JG (2015) Co-existence in a species-rich grassland: Competition, facilitation and niche structure over a soil depth gradient. J Veg Sci 26:674–685. 152 Martorell, C. y Freckleton, R.P. (2014). Testing the roles of competition, facilitation and stochasticity on community structure in a species-rich assemblage. J Ecol 102, 74–85 Martorell C, Martínez‐Blancas A, García‐Meza D (2021) Plant–soil feedbacks depend on drought stress, functional group, and evolutionary relatedness in a semiarid grassland. Ecology 0:e03499. Martorell C, Zepeda V, Martínez-Blancas A, García-Meza D, Pedraza F (2017) A diversity world record in a grassland at Oaxaca, México. Bot Sci 95:1–7. Pedraza F y Martorell C (2019) Allocating species in Grime’s strat- egy space: an alternative to trait-based approaches. Bot Sci 97:649–660. Sakavara A, Tsirtsis G, Roelke DL, Mancy R, Spatharis S (2018) Lumpy species coexistence arises robustly in fluctuating resource environments. 115. Scheffer M, Nes EH van (2006) Self-organized similarity, the evolutionary emergence of groups of similar species. Proc Natl Acad Sci 103:6230–6235. Scheffer M, Vergnon R, Nes EH Van, et al. (2015) The Evolution of Functionally Redundant Species; Evidence from Beetles. PLoS One 10:e0137974. Segura AM, Calliari D, Kruk C, Conde D, Bonilla S, Fort H (2011) Emergent neutrality drives phytoplankton species coexistence. Proc R Soc B Biol Sci 278:2355–2361 10.1098/rspb.2010.2464. Silvertown J (2004) Plant coexistence and the niche. Trends Ecol Evol 19:605–611. Tielbörger K, Kadmon R (2000) Temporal environmental variation tips the balance between facilitation and interference in desert plants. Ecology 81:1544–1553. Tuck SL, Porter J, Rees M, Turnbull LA (2018) Strong responses from weakly interacting species. Ecol Lett 10.1111/ele.13163. Turnbull LA, Manley L, Rees M (2005) Niches , rather than neutrality , structure a 153 grassland pioneer guild. 1357–1364. Vergnon R, Dulvy NK, Freckleton RP (2009) Niches versus neutrality: uncovering the drivers of diversity in a species-rich community. Ecol Lett 12:1079–1090. Vergnon R, Nes EH van, Scheffer M (2012) Emergent neutrality leads to multimodal species abundance distributions. Nat Commun 3:1–6. Villarreal-Barajas T, Martorell C (2009) Species-specific disturbance tolerance, competition and positive interactions along an anthropogenic disturbance gradient. J Veg Sci 20:1027–1040. Wainwright CE, HilleRisLambers J, Lai HR, Loy X, Mayfield MM (2018) Distinct responses of niche and fitness differences to water availability underlie variable coexistence outcomes in semi-arid annual plant communities. J Ecol 107:293–306. Yan X, Levine JM, Kandlikar GS (2022) A quantative synthesis of soil microbial effects on plant species coexistence. Proc Natl Acad Sci 119: e2122088119. Zepeda V, Martorell C (2019) Fluctuation-independent niche differentiation and relative non-linearity drive coexistence in a species-rich grassland. Ecology 100:e02726.