Genecology and ecophysiology of the maintenance of foliar phenotypic polymorphisms of Leptospermum recurvum (Myrtaceae) under oscillating atmospheric desiccation in the tropical-subalpine zone of Mount Kinabalu, Borneo

We investigated genecology and ecophysiological mechanisms of the polymorphism of leaf trichome density of Leptospermum recurvum Hook. f. (Myrtaceae) in the deglaciated summit zone above 3,000 m a.s.l. of Mt. Kinabalu, Borneo. Various phenotypes with variable foliar trichome densi-ties occurred sympatrically in the same population, and the composition of coexisting phenotypes varied substantially among populations. We conducted a common garden experiment by sowing seeds from multiple maternal trees of different leaf trichome densities. We found a significant relation between pubescence of maternal trees and offspring, which indicated that leaf trichome density had a genetic basis. Microsatellite analysis revealed that there was no barrier to gene flows among phenotypes or among populations, and very low neutral genetic differentiation among populations with high gene flows for both directions of phenotypes. The soils in the sites dominated by pubescent trees were significantly more desiccated than in the sites dominated by glabrous trees during a short drought. Glabrous trees had a significantly greater mortality rate than pubescent trees after an intensive El Niño drought (13.7 vs. 3.9%)


| INTRODUCTION
Remarkable adaptive radiation occurs often in isolated environments such as oceanic islands or summit zones of tropical mountains. The most remarkable example is the silver-sword alliance radiated from a single ancestor to several distinct genera in Hawaii (Carr, 1987). Radiated species are genetically isolated from each other. By contrast, some species demonstrate extreme intraspecific phenotypic polymorphism within/among populations, which may be still at an incipient stage of adaptive radiation, as seen in Metrosideros polymorpha Gaudich. (Myrtaceae) of Hawaii (Izuno et al., 2017;Kitayama, Mueller-Dombois, Pattison, & Cordell-Hart, 1997;Tsujii, Onoda, Izuno, Isagi, & Kitayama, 2015). Genetically based intraspecific phenotypic polymorphism may be a precursor to local adaptation. To what extent genetically based phenotypes are differentiated between populations depends on the balance between selection and gene flow (Lenormand, 2002;Mckay & Latta, 2002). In the case of a restricted gene flow, the alleles related to adaptive traits can be effectively selected, and genetic drift can also promote genetic differentiation. Conversely, because a gene flow has homogenizing effects on genetic structure, maintaining differentiation of phenotypes under a strong gene flow requires ongoing strong divergent selection; alternatively, the phenotypes may be neutral to the selection. Mechanisms of phenotypic differentiation driven by ongoing natural selection have not been well elucidated except for plant species occurring in very distinct edaphic and biotic conditions (Kittelson & Maron, 2001;Linhart & Grant, 1996;Sambatti & Rice, 2006).
Leptospermum recurvum Hook. f. (Myrtaceae) presents another example of phenotypic polymorphism presumably at an incipient radiation stage on an isolated tropical mountain. Here, we discuss genecology and ecophysiological mechanisms of foliar polymorphisms using this species as a model plant. L. recurvum is a woody species, endemic to Mt. Kinabalu (4,095 m), Borneo, and occurs primarily on deglaciated rock faces from 2,700 to 4,095 m a.s.l. It has been derived from its widespread lower montane counterpart of Leptospermum flavescens Sm., and radiated into Pleistocene postglacial vacant niches of the tropical alpine zone (Lee & Lowry, 1980). L. recurvum demonstrates a continuum of phenotypic polymorphisms in leaf trichome density, that is, ranging from glabrous to pubescent phenotypes (see Figure S1). Peculiarly, in this species, various phenotypes often occur sympatrically in the same population (Lee & Lowry, 1980), and the composition of coexisting phenotypes varies among populations. Some populations consist of entirely glabrous phenotypes, some with entirely pubescent ones and the rest with variously mixed phenotypes (Lowry, 1969). It has long been speculated that the magnitude of phenotypic variations of leaf trichomes of this species within/among populations reflects the degree of soil desiccation, because leaf trichomes might be adaptive to the avoidance of desiccation stress (Lee & Lowry, 1980;Lowry, 1969). It is known that intermittent atmospheric desiccation and irregular El Niño droughts characterize the alpine zone of Kinabalu (see study site), which causes a variation of desiccation stress in association with variable soil depth and soil waterholding capacity. L. recurvum occurs primarily on denuded granitic rock faces above 3,000 m as a dwarf shrub or a cushion plant (Kitayama, 1992;Smith, 1980), but extends down to 2,700 m on deeper soils of oligotrophic ultrabasic substrate as a small tree (Aiba & Kitayama, 1999), again suggesting that this species varies in life form according to soil depth. L. recurvum is primarily monoclinous often with andromonoecy flowers, and appears to be pollinated by flies and bees (Kudo & Suzuki, 2004), suggesting the presence of rampant gene flows among populations.
The differentiation of phenotypes in L. recurvum may not be complete, so that it is predicted that gene flow and divergent selection are involved in this species. Leaf trichome density has been considered one of xeromorphic characters because leaf trichomes increase boundary layer resistance of transpiration (Amada, Onoda, Ichiei, & Kitayama, 2017;Ripley, Pammenter, & Smith, 1999;Wuenscher, 1970). If pubescent phenotypes represent adaptive variations to desiccation stress, El Niño droughts may cause difference in mortality among the phenotypes (i.e., variation in leaf trichomes) and then act as one end of divergent selection in populations dominated by pubescent phenotypes (see Section 2). Lee and Lowry (1980) suggested lower mortality of pubescent phenotypes than glabrous ones after an El Niño drought by a field survey, but they investigated 52 individuals only. In addition, although an allozyme analysis was conducted in their study, few polymorphic loci were detected, and the pattern of gene flow remained unveiled. Leaf trichome density is probably genetically based for this species, as seen in model plant Arabidopsis thaliana (L.) Heynh (Pattanaik, Patra, Singh, & Yuan, 2014) or in Hawaiian M. polymorpha (Izuno et al., 2017;Kitayama et al., 1997;Tsujii et al., 2015), but it is unknown to date.
In our study, we ask what evolutional processes (selection, gene flow or drift) shape phenotypic polymorphisms of L. recurvum, and how the droughts associated with El Niño influence the ecological and genetic process in the subalpine zone of Mt. Kinabalu. For the first question, we examined the two following hypotheses: (a) leaf trichome density is a genetically based trait and (b) there is no barrier (or a small if any) to gene flow among populations where different phenotypes are dominant. To examine these hypotheses, we performed a common garden experiment and a population genetic analysis with 10 polymorphic microsatellite markers. For the second question, we examined the following hypotheses: (c) soil moisture of the site dominated by pubescent phenotypes is more susceptible to droughts than those dominated by glabrous phenotypes, (d) the mortality of pubescent phenotypes in a given site after an El Niño drought is lower than that of glabrous phenotypes in the same site and (e) pubescent phenotypes have a physiological adaptation to desiccation. To examine these hypotheses, we conducted field survey and physiological experiments. Through these examinations, we will discuss the influence of El Niño droughts on the maintenance of phenotypic polymorphism of this species.
The term "drought" means two phenomena. One is that rainfall is very little during a certain period and the other is that water content of air, soil or tissues is low. We use "drought" as the former context (e.g., El Niño droughts), whereas the latter context is expressed by "desiccation" (e.g., soil or atmospheric desiccation).

| Study site and meteorological conditions
Study site is the tropical subalpine zone above 3,000-3,800 m a.s.l. on Mount Kinabalu. Kinabalu is the highest mountain in Southeast Asia and represents tropical high mountains of this region. Rainfall around Mt. Kinabalu is usually ample throughout the year with more than 2,500 mm per year in the lowland to the subalpine zone (Kitayama, Lakim, & Wahab, 1999). However, the subalpine zone of Mt. Kinabalu is subjected to intermittent atmospheric desiccation due to the interplays of ascending radiation-driven air during a day and descending air during a night, which causes rapid oscillation of daily saturation deficits at high altitudes above 3,000 m (Kitayama, 1996). Therefore, subalpine (and alpine) plants are exposed to rapidly alternating atmospheric desiccations, although annual rainfall is ample (Kitayama, Ando, Repin, & Nais, 2014). Moreover, an extensive drought occurs in the subalpine zone during El Niño events. While severe droughts can occur also in the lowland during El Niño to the extent that causes higher tree mortality (Aiba & Kitayama, 2002;Potts, 2003), the magnitude of droughts is greater above about 3,000 m than at lower altitudes with less precipitation and a higher evaporation demand Kitayama et al., 2014). This is because the Walker circulation, a synoptic air circulation in the equatorial Pacific, becomes weakened during an El Niño event and brings about a stronger air subsidence with dry conditions at higher altitudes in the western Pacific (Kitayama, 1996;Kitayama et al., 2014;Nullet & Juvik, 1994). For example, in the El Niño drought in 2010 prior to our study, the total precipitation during the peak driest 50 days at 3270 m was only 6.1 mm (unpublished data). This intermittent and/or irregular occurrence of desiccation would result in a specific adaptation which is different from the adaptation under chronic desiccation (e.g., Holmgren et al., 2006;Stenseth et al., 2002).
The large parts of the region above 3,000 m of Mt. Kinabalu are covered with granite bed rocks, and therefore plant communities are assembled only in places where glacial deposits and debris are deposited, for example, gentle slope ridges, knick line (concave), depression of bed rock and rock crevices. Various types of vegetation are shaped according to soil depth, from a forest with closed canopy to scrub or sparse communities on rock crevices, which are distributed in a mosaic manner. The magnitude of the tree mortality varied among plant communities (Kudo & Kitayama, 1999;Smith, 1980). Aiba and Kitayama (2002) and Kitayama et al. (2014) suggested that the alpine and subalpine vegetation is controlled by droughts based on the study of tree mortality after a drought and vegetation physiognomy.

| Sampling site and phenotypic composition of populations
Population study was conducted between 3,100 and 3,800 m on the southern slope of Mt. Kinabalu (Figure 1), where granite bedrocks occurred and a mosaic of vegetation with various structures was formed. Since those L. recurvum populations which exceeded 2 m in canopy height tended to occupy a large area (each population more than several hundreds of meter square), we ensured an enough sampling size in investigating phenotypic composition and edaphic environment, and in applying a population genetic analysis at these sites. These populations demonstrated either a predominantly pubescent population or a predominantly glabrous population. Therefore, the sites for these investigations were divided into two population (site) types; sites dominated by glabrous Leptospermum (glabrous sites, GS) and sites dominated by pubescent Leptospermum (pubescent sites, PS). In our preliminary observation, the species richness and understory density of the plant communities where these populations occur were each higher in GS than in PS.
On the other hand, for investigating mortality rate and leaf traits, we selected the sites where populations had a lower canopy height than 2 m and mixture of various phenotypes, in order to evaluate the differences between phenotypes at the same site and to eliminate environmental effects among sites.
Composition of phenotypes by leaf trichome density was investigated in 12 sites (populations) including eight GS and four PS (Figure 1). Site names in Figure 1 are encoded by the following codes; first two digits indicate the altitude of a site followed by its predominant leaf pubescence (e.g., 33P indicating a PS at 3300 m). Because multiple sites for GS were selected from the same altitude, site names were annotated with e (east) or w (west) according to their relative positions. We evaluated the leaf trichome density of 20-60 individuals in each site, using six classes as described below. The leaf trichome density was semi-quantitatively evaluated by visual observation and divided into the following six classes (0-5): 0, most leaves have no trichome on both sides; 1, some leaves have no trichome but others have trichome only on a main vein or a leaf base on abaxial side; 2, most leaves are covered thinly by trichome on abaxial side; 3, most leaves are covered thickly on abaxial side; 4, most leaves are covered thickly on abaxial side and covered thinly on adaxial side and 5, some leaves are like Class 4 but others are covered thickly on both sides ( Figure S1).

| Common garden experiments
To examine the genetic basis of pubescence, we carried out a common garden experiment and compared pubescence between parents and offspring. In June 2011, approximately 10 capsules of L. recurvum were collected from four or five trees per site: 33P, 35Ge, 37P and another site on rock-crevice at 3,800 m (38R). The pubescence of sampled maternal trees is shown in Table S1. Twenty to thirty ripe seeds per maternal tree were sown immediately, and germinated and grown in a greenhouse exposed to open air but protected from rain at the nursery of Kinabalu Park (1,610 m, mean annual temperature 18.1 C, relative light intensity 27%) for 3 months. In September 2011, randomly selected 8-20 healthy seedlings per maternal tree were replanted into 6-cm cubed pots filled with a 1:2:2 mixture of river sand, peat and forest topsoils. Forest topsoils were collected in a forest stand nearby the greenhouse. After germinated, pots were watered regularly, and given fertilizer and insecticide when required. Finally, in July 2012, we obtained a total of 133 living seedlings. Mortality (for period between sowing and measurements) was 51% with no differences among maternal pubescence or among sampling sites. Total branch length and pubescence were measured for all seedlings. For evaluating offspring pubescence, we used two newly defined dichotomous categories (0 0 , all leaves had no trichomes on any sides; 1 0 , some leaves had trichomes at least on any part of leaves). This is because many seedlings had some juvenile leaves that were long, thin and glabrous at the time when the measurements were conducted. To test the effects of maternal pubescence on offspring pubescence, we used a generalized F I G U R E 1 Map of the study sites. Glabrous-type-dominant sites (GS) are marked with open circles, pubescent-type-dominant sites (PS) with closed circles for investigation of phenotypic composition, edaphic environment and genetic analysis. Site names were defined by the following codes; first two digits indicate the altitude of a site followed by its predominant leaf pubescence (e.g., 33P indicating a PS at 3,300 m). Because some sites for GS were selected from the same altitude, site names were annotated with e (east) or w (west) according to their relative positions. Mortality census was conducted around the two areas with closed squares. Leaf traits measurements were conducted using population at the site with open squares. The closed triangle indicates Low's Peak, the summit of Mt. Kinabalu linear mixed model (GLMM) with a binomial error distribution and a logistic link function ("lmer" function in package "lme4" (Bates, Maechler, & Bolker, 2012)) in R, with offspring pubescence as binomial dependent variable (1 0 , 0 0 ). The maternal pubescence (0-5 as continuous variable), the offspring total branch length and the seedsampling sites were fitted as fixed terms, and the identity of maternal tree was fitted as a random term. The reason why offspring total branch length was used as a fixed term is that the ratio of mature leaves to juvenile leaves appeared to increase with seedling size and that we considered this parameter as covariate correcting the effects of immaturity.

| Population genetic analysis
We analyzed 12 L. recurvum populations at the sites where phenotypic composition was investigated, using 10 microsatellite markers developed by Ando, Kaneko, Isagi, Repin, and Kitayama (2013); Lrec005, Lrec134, Lrec139, Lrec206, Lrec273, Lrec264, Lrec296, Lrec349, Lrec368 and Lrec475. For each population, leaf samples were taken from 19 to 27 adult trees, which were widely spaced in between without canopy overlaps, to avoid clonal plants. We placed these samples immediately in silica gel for later DNA extraction. DNA extraction, PCR amplification of 10 microsatellite loci and genotyping of each individual were carried out following the procedure of Ando et al. (2013).
The mean number of alleles per locus, allelic richness, expected heterozygosity (H e ) (Nei, 1987), inbreeding coefficient (F IS ) and deviation from Hardy-Weinberg equilibrium were calculated over all loci in each population based on 2,400 randomizations using FSTAT (Goudet, 1995(Goudet, , 2001. To determine genetic variance partitioned into three levels, that is, among site types, among populations within site types, and within populations, we conducted an analysis of molecular variance (AMOVA) (Excoffier, Smouse, & Quattro, 1992). Fixation indices were computed and tested using 20,000 permutations for each level of genetic structure: F CT for variation among site types, F SC for variation among populations within site types and F ST for variation within populations. Genetic variance and fixation indices were both computed by ARLEQUIN (Excoffier, Laval, & Schneider, 2005). In addition, the number of migrants per generation for diploid organisms, 4Nm (four times the product of effective population size and immigration rate) among site types and the mutation-scaled effective population size Θ, which is the product of effective population size and mutation rate, for each site type were estimated with a coalescent approach using the Bayesian implementations available in MIGRATE-N (Beerli, 2006;Beerli & Felsenstein, 2001). For this analysis, allelic data were pooled for each site type. The Bayesian inference was run with one long chain, an increment of 20, a sampling of 100,000 and a burn-in of 50,000. The prior distribution was uniform, and the priors were 0.0005-200 and 0-7,000 for Θ and 4Nm, respectively.

| Characterizing the edaphic environment
To reveal the site condition of each population type, we evaluated soil chemical and physical properties. We selected three sites each for GS and PS (33Ge, 35Ge and 37Ge for GS; 33P, 35P and 37P for PS; Figure 1), and randomly collected three 10-cm-deep soil core samples beneath the litter layer in each site, and combined them as one sample for each site. Soil bulk density and soil water content were measured. In addition, soil total carbon and total nitrogen were determined using a CN analyzer (JM1000CN, J SCIENCE, Kyoto, Japan). To examine the role of El Niño as the selective pressure, we needed to evaluate soil moisture during El Niño droughts, but El Niño did not occur during the period of our soil-moisture survey. However, instead of El Niño droughts, we encountered a short drought event in which the cumulative precipitation was only 5.2 mm from June 14 to July 1, 2012. On the last day of the short drought, predawn (0300-0600 hr) xylem water potential (Ψ X ) was measured as an index of "realized desiccation" (Richter, 1997; but see Donovan, Linton, & Richards, 2001) at two sites for each site type (33Ge and 35Ge for GS, and 33P and 35P for PS; Figure 1). Ψ X was measured for randomly selected six or seven 1 m-height saplings per site, using a pressure chamber (model 1003, PMS Instruments, New York, NY). To eliminate nighttime water loss of leaves, the measured branches were enveloped with plastic bags overnight. Because we could not find enough L. recurvum saplings in understory, we also used saplings of several other species in all sites.
For the data derived from soil core, we compared the differences of bulk density and water content between site types using a t test in R 2.15.1 (R Core Team, 2012). The effect of site type on predawn Ψ X was analyzed in R using a two-way analysis of variance (ANOVA) with site types (GS or PS) and site altitudes (3,300 m or 3,500 m) as explanatory terms.

| Mortality rate in the field
From January to March 2010, Mt. Kinabalu was affected with a drought associated with El Niño (see introduction).
Many L. recurvum died back after the drought, especially in populations with lower canopy height than 1 m such as crevice vegetation and small sandy patchy scrubs. In September 2010, to clarify the relation between phenotypes and mortality, at around 3,300 and 3,500 m, we inventoried randomly selected 296 individuals with the height ≧10 cm at randomly selected five sites whose populations demonstrated a mixture of phenotypes (smaller than 5 × 15 m 2 ) with canopy height less than 1 m. We checked the condition (dead/alive) of the selected individuals and measured the stem basal diameter and pubescence of them. We judged those individuals to be dead when they had no green leaves on any branches. Even on dead individuals, many reddish brown leaves remained and we could evaluate pubescence of dead individuals based on the remaining leaves. We also found dead individuals with no leaves, but such individuals tended to have white and decayed bark with some lichen, compared to dead individuals with leaves remained. Therefore, such individuals were considered to be killed before the drought, and we excluded such individuals from the census. To test the effects of pubescence on mortality, we used a GLMM with a binomial error distribution and a logistic link function ("lmer" function) in R, with the condition of individuals as a binomial dependent variable (alive 1, dead 0). The pubescence of individuals (0-5 as continuous variables) and the stem basal diameter were fitted as fixed terms, and the identity of the site was fitted as a random term.

| Ecophysiological leaf-trait measurements
To study difference in leaf-level adaptation between phenotypes, we evaluated various leaf traits using a leaf gas exchange measurement, a leaf chemical analysis and a leaf-drying curve analysis (L-D analysis, see Cape & Percy, 1996). All plant materials used here were collected from individuals, which were widely spaced in between without canopy overlaps, in a single population at 3260 m, where various phenotypes occurred sympatrically. The population was a pure stand of L. recurvum near the mountain hut (Waras Hut) and composed of individuals of similar heights (1-2 m). They demonstrated a range of pubescence classes and appeared to be randomly distributed within a small patch of less than 20-m diameter on a homogenous soil substrate. We, therefore, assumed that the site was nearly homogeneous in environmental conditions and that the individuals grown here experienced similar environmental effects like common gardens. Therefore, the differences in leaf traits obtained here were considered to reflect the differences of genetic effects among variations.
The photosynthetic gas exchange of sun leaves was measured in situ with a portable infrared gas analyzer system with a leaf chamber (LI-6400; LI-COR, Lincoln, NE). An intact shoot with multiple leaves was inserted to the leaf chamber. A special care was taken to avoid leaf overlaps inside the chamber while measuring photosynthesis. After the measurement, leaves were excised from each shoot and scanned for determining leaf areas, which were later used for correcting photosynthetic gas exchange rates. The following parameters were calculated: assimilation rate per unit leaf area (A area ) and per unit leaf mass (A mass ), nitrogen content per leaf area (N area ), nitrogen-use efficiency (quotient A area divided by N area ), stomatal conductance ( g s ), transpiration rate (Tr), ratio of internal and external CO 2 concentration (C i /C a ), internal conductance ( g i ) (quotient A area divided by C i ) and respiration rate per unit leaf area (R area ) and per unit leaf mass (R mass ). Foliar nitrogen content was determined on ovendried leaves after field tests (see below). We compared these parameters between the extremes of phenotypes (i.e., pubescence Classes 0 and 5) using six individuals for each phenotype and two replications for each individual. The measurements were taken between 0800 and 1400 in situ, and the leaf temperature, vapor pressure deficit and ambient CO 2 concentration were adjusted to 20 C, 0.7-1.2 kPa and 370 ppm, respectively. Light levels were maintained at 2500 μmol m −2 s −1 using LED lights (model 6400-02B; LI-COR), which was higher than a light saturation point but lower than photoinhibition levels for both glabrous and pubescent phenotypes confirmed by preliminary experiments.
After the gas exchange measurements, the measured leaves were dried in an oven at 60 C, and leaf mass per area (LMA) was calculated. Total leaf nitrogen and carbon concentrations (C mass and N mass ) and the carbon isotope compositions of the leaves that were used in the gas exchange measurements were determined using an elemental analyzer (EA1180; Fisons Instruments, Milan, Italy) connected to a mass spectrometer (Delta S, Finnigan MAT, Bremen, Germany) with an interface (Conflo II, Finnigan MAT). We corrected data with international standard substances (Tayasu, Hirasawa, Ogawa, Ohkouchi, & Yamada, 2011). The δ 13 C signature was calculated versus the Vienna Pee Dee Belemnite standard. The standard deviation of δ 13 C in standard samples was 0.09‰.
The ability of preventing water loss at minimum stomatal openness was evaluated using L-D curves following Cape and Percy (1996). Comparing the extremes of phenotypes (pubescence Classes 0 and 5), we selected six individuals for each phenotype. From the top of canopy for each individual, we collected a 40-cm-long branch, which was subsequently re-cut under water immediately.
In the laboratory, the cut ends were soaked in water, and the 40-cm-long branches were covered with black plastic bags for a day to be saturated with water. Subsequently, three pieces of branches (0.2-1.1 g saturated weight) were taken from each 40-cm-long branch (i.e., three replications for each individual) and air-dried on a laboratory bench with the cut ends sealed with PTFE tape in an airconditioned room (23 C, 75% relative humidity), and weighed 17 times at regular intervals for 25 hr. After the experiments, the branches were dried in an oven at 60 C for 48 hr and weighed as oven dry weight. All measured weights were converted into relative water content (RWC, %) calculated as: RWC = fresh weight −oven dry weight ð Þ = saturated weight − oven dry weight ð Þ × 100: An example of the relation between drying time and RWC is shown in Figure S2. The obtained data were fitted for each branch using the following equation by Cape and Percy (1996), where RWC t is the RWC at time t, a is the RWC at equilibrium with the atmosphere, b is the intercept of the curve at time 0 and k is the constant with units of (time) −1 . A larger k means a greater reduction rate of RWC, that is, the lower ability of preventing water loss at a minimum stomatal openness. For fitting our data to the equation, we omitted the first 3-hours data because stomata might have not reached a minimum openness ( Figure S2). R 2 for fits ranged from .9947 to .9999. To examine the difference of the leaf traits between the phenotypes (Classes 0 and 5), we made a linear mixed model (LMM) for each trait using "lme" function in package "nlme" (Pinheiro et al., 2012) in R with pubescence (as categorical data) as a fixed term and identity of individuals as a random term. In addition, we examined the effects of phenotypes on photosynthetic water-use efficiency (WUE, quotient A area divided by Tr) using LMM in R with A area as a dependent variable, Tr and pubescence (as category) as fixed terms, and identity of individuals as a random term.

| Common garden experiments
Maternal pubescence significantly affected offspring pubescence, but sites of the origin of seeds did not (Table 1). For the offspring from pubescent maternal trees (Classes 2-5) 31.4% of 70 offspring were pubescent, whereas only 3.2% of 66 offspring from glabrous maternal trees (Classes 0 and 1) were pubescent. Offspring pubescence was also affected by their own total branch length showing their degree of growth. However, total branch length was not correlated to maternal pubescence (r = .04, p > .2, Pearson's correlation coefficients).

| Population genetic analysis
In contrast to the results from the isozyme analysis by Lee and Lowry (1980), which showed genetic uniformity among plants, all populations showed at least nine polymorphic loci out of 10 microsatellite loci analyzed here, and the number of alleles per locus per population ranged from one to 18, which are informative enough to estimate the genetic structure and gene flow of L. recurvum. The value of allelic richness and expected heterozygosity (H e ) of each population are shown in Table 2, indicating similar values among populations regardless of altitude or the dominant type of phenotypes. Inbreeding coefficients (F IS ) of all populations were at a low level and no population significantly deviated from Hardy-Weinberg equilibrium after a Bonferroni correction (Table 2).
AMOVA revealed that most of the total molecular variance was distributed within population (Table 3), whereas the differentiation among site types and the differentiation among populations within site type accounted for less than 2% of the total variance. In addition, F SC and F ST were significantly deviated from zero, but all the fixation indices were very low. These results suggested that the level of neutral genetic differentiation was very low across the whole populations of L. recurvum investigated in this study.
The mutation-scaled effective population size (Θ) and gene flow (4Nm) estimated with MIGRATE-N are shown in Table 4. The 95% confidential intervals (CI) of both Θ and 4Nm were highly overlapped between the two site types, indicating similar population sizes and a nearly symmetric gene flow. Medians of 4Nm for both directions were high (4Nm > 50), but 95% CI included zero. Median of 4Nm from GS to PS was nearly twice as that from PS to GS.

| Population phenotypic composition
Clear differentiation of phenotypic composition was found between GS and PS ( Figure 2). Dominant phenotypes in GS were Classes 0 and 1, whereas Classes 2-5 T A B L E 1 A binomial GLMM analysis explaining offspring pubescence of the common garden experiment (pubescent [1 0 ] 1, glabrous [0 0 ] 0). Here, offspring pubescence was classified to two new classes, 1 0 (some leaves had trichomes at least on any part of leaves) and 0 0 (all leaves had no trichomes on any sides). SE means standard error of each estimate. The identity of maternal tree was fitted as a random term. The number of total offspring was 133. See the text and Table S1 for seed sampling sites and maternal trees

T A B L E 2
Values of genetic characters and the SE estimated at 10 microsatellite loci for each population. Allelic richness was based on minimum sample size of 18 diploid individuals. H S , expected heterozygosities; F IS , inbreeding coefficient; HWE, proportion of randomizations that gave a larger F IS than the observed based on 2,400 randomizations. Although the HWE value of 37Ge was lower than 0.05, this was not significant level after Bonferroni correction. See Figure 1 for site name of populations were dominant in PS. For GS, the dominance of extreme phenotype, that is, Class 0, was very high, but for PS, Class 5 showed moderate dominance and various phenotypes rather cooccurred. GS and PS were distributed in a mosaic-like manner and within a small area (see Figure 1). For example, the distance between 35P and 35 Ge was less than 100 m.

| Characterizing the edaphic environment and water potential
The soil chemical composition, bulk density and soil water content of saturated soil are shown in Table S2. Both total carbon and total nitrogen concentration were higher in GS than PS, but C/N ratio was not significantly different between site types. Soil bulk density was higher and soil water content was lower in PS than in GS. The pattern of these soil characters of PS corresponds to the lower accumulation of soil organic matter. Community level-predawn xylem water potential (Ψ X ) during the short drought was significantly different between site types (PS vs. GS) and between altitudes of the sites in ANOVA, and the proportion of variance component of site types was much higher than that of altitudes (Table 5). Mean ± SE of Ψ X was −0.21 ± 0.03 MPa and −0.44 ± 0.04 MPa in GS and PS, respectively. Mean of Ψ X in PS was approximately twice as low as in GS suggesting a much greater desiccation stress in PS.

| Mortality rate
A GLMM analysis indicated that pubescent phenotypes had a significantly lower mortality than glabrous phenotypes at the same sites in the census conducted 6 months after an El Niño drought ( Table 6). The mean mortality of all sites was 3.9 and 13.7% for pubescent individuals (Classes 2-5) (n = 128) and glabrous individuals (Classes 0 and 1) (n = 168), respectively. Basal diameter did not significantly influenced mortality.

| Ecophysiological leaf traits
Mean ± SE of all leaf traits for each phenotype obtained in the field are shown in Table 7 with the p value of the effect of pubescence in an analysis using LMM. For gas exchange measurements, almost all traits were not significantly different between phenotypes, but g i of Class 5 (most pubescent) was significantly larger than that of Class 0 (glabrous). The relation between A area and Tr was significantly different between phenotypes (Figure 3). A area of Class 5 was higher at a given Tr than that of Class 0 (significant difference in intercept, p < .01), suggesting that Class 5 (i.e., most pubescent) had a higher WUE than Class 0 did.
Other traits of leaves used in gas exchange measurements were also different between phenotypes. LMA and leaf area of individual leaves of Class 5 were found to be larger than that of Class 0. For leaf chemical composition, C mass of Class 0 was significantly higher than that of Class 5. N mass was not significantly different between phenotypes, although Class 5 showed slightly higher N mass than Class 0 (Table 7). Reflecting the difference of LMA, however, higher N area was shown in Class 5 than in Class 0. C/N ratio was higher in Class 0 than in Class 5. δ 13 C was not significantly different between phenotypes.
T A B L E 4 95% confidential intervals and median of mutation-scaled effective population size (Θ) for each site type and the number of migrants per generation (4Nm) among site types for each direction estimated by MIGRATE. Θ = 4N e μ, where N e and μ are effective population size and mutation rate, respectively. For site types, GS and PS mean glabrous site and pubescent site, respectively 2.50% Median 97.50% The ability of preventing water loss at minimum stomatal openness evaluated using L-D analysis was not significantly different between phenotypes (Table 7), although pubescent leaves were expected to have a higher boundary resistance than glabrous leaves.

| Genetic process of shaping polymorphisms
In our common garden experiment, a significant relation between pubescence of maternal trees and offspring was found, and this indicated that leaf trichome density had a genetic basis as we hypothesized (Hypothesis a). Although most measured seedlings were likely immature, size effects were removed by including seedling size in the GLMM analysis. Furthermore, as our field survey demonstrated, the sympatric occurrence of various phenotypes within a site (Figure 2) suggested that leaf trichome density was not determined by environmental effects only. However, we conducted the common garden experiment under one environment only, and the pattern and extent of plasticity in leaf trichome density remain unknown. Since the effect of site where seeds were collected was not significant in the GLMM analysis, it was unlikely that the environment that affected maternal trees influenced offspring pubescence.
Microsatellite analysis revealed that there was no barrier to gene flows among phenotypes and among site types in agreement with our hypothesis (Hypothesis b). AMOVA results showed very low neutral genetic differentiation between site types, and medians of gene flow (4Nm) calculated by MIGRATE were high for both directions. Although 95% CI of 4Nm included zero, MIGRATE results indicated high levels of gene flow between site types (AMOVA). In addition, we found significant but very low F ST and no population deviated from Hardy-Weinberg equilibrium significantly. These results mean that strong gene flows occur across whole L. recurvum populations investigated in this study. This may be due to flowering synchrony after El Niño droughts (Kudo & Suzuki, 2004) and gregarious visits of flies and bees (personal observation). On the other hand, we found clear differentiation of phenotypic composition between GS and PS ( Figure 2). These results suggest that the genetic processes shaping the phenotypic polymorphisms involve strong gene flows combined with ongoing divergent selection. Effects of genetic drift must be negligible.

| Selective agent and adaptation
The greater desiccation of soil moisture in PS during the short drought than in GS suggested the susceptibility of PS to droughts in agreement with our hypothesis (Hypothesis c). We measured the predawn Ψ X as "realized desiccation" between site types during a 17-day short drought. As a result, the values of Ψ X were not so low to cause tree dieback. However, our results showed that mean of Ψ X in PS were approximately twice as low as in GS. Then, it is inferred that the difference of soil moisture between site types even becomes larger during longer El Niño droughts that usually last for several months, and the decline of soil moisture in PS would result in the lethal desiccation stress for glabrous phenotypes. It is T A B L E 5 ANOVA results of community level predawn xylem water potential as an index of soil moisture

T A B L E 6
A binomial GLMM analysis explaining mortality after the 2010 El Niño drought in five sites with admixture of phenotypes. The condition of individuals was binomial dependent variable (alive 1, dead 0). SE means standard error of each estimate. The identity of investigated sites was fitted as a random term. The number of total investigated individuals was 296 suggested that the decline of soil moisture during El Niño drought acts as the selective agent against glabrous phenotypes, given that higher mortality of glabrous phenotypes than pubescent phenotypes after an El Niño drought. We evaluated "realized desiccation" from predawn Ψ X , assuming an equilibrium between Ψ X and root surface soil water potential. Although some species may not satisfy this assumption (Donovan et al., 2001), measuring multiple species will overcome this problem as we did. Therefore, our results conservatively indicate a greater decline of soil moisture in PS than in GS during the 17-day short drought period.
The primary cause of the different soil-moisture status during a drought between site types must be a difference in water-holding capacity of soils, which covary with topography. Most PS were found on depressions or near knick lines (concave) of granite bed rocks, which were more fragmented and had shallower soils, whereas GS were frequently found in large and continuous forests on gentle slope ridges with deep soils. On the other hand, narrow but deep rock crevices may also provide a prolonged water supply even during drought periods, and many glabrous phenotypes were also found on rock crevices (Lowry, 1969). In addition to topography, the lower accumulation of soil organic matter in PS may influence water-holding capacity of the soils.
Mortality was lower in pubescent phenotypes than in glabrous ones as we hypothesized (Hypothesis d), which suggested higher desiccation tolerance of pubescent phenotypes. In our measurements, pubescent phenotypes of Class 5 demonstrated a higher WUE (A area / Tr) than glabrous phenotypes did (Figure 3), which indicated that at least leaf-level physiological adaptation contributed to the higher desiccation tolerance of pubescent phenotypes. Gutschick (1999) showed two physiological mechanisms of enhancing WUE, increasing g i and decreasing g s . Both phenotypes of L. recurvum T A B L E 7 Mean and SE of all leaf traits for each phenotype. The identity of individuals was fitted as a random term. See Section 2 for explanations of the traits and the pubescence classes showed a similar g s range, but pubescent phenotypes showed higher g i than glabrous phenotypes did ( Table 7). The g i , that is, A area /C i , includes mesophyll conductance and photosynthetic capacity. We did not measure mesophyll conductance, but our results imply a higher photosynthetic capacity of Class 5. Class 5 showed a similar level of N mass to Class 0 but higher LMA and thereby higher N area than those of Class 0. Therefore, these patterns imply that higher g i results from higher nitrogen investment to photosynthesis (Wright, Reich, & Westoby, 2003) and that leaves of Class 5 pack more mesophyll and chloroplasts in unit area and increase photosynthetic capacity per unit area. The fact that leaf-trichome density reflects distribution of intraspecific variations of L. recurvum means that leaf trichome density is related to adaptation. However, we could not find direct evidence of adaptive function of leaf trichomes per se. Although saving water loss was primarily expected as a function of leaf trichomes for desiccation tolerance (Ripley et al., 1999;Wuenscher, 1970), neither transpiration rate in the gas exchange measurements nor the ability of preventing water loss in the L-D analysis was significantly different. Amada et al. (2017) demonstrated a similar result with us in that H 2 O flux resistance by a leaf trichome layer was negligible in Hawaiian M. polymorpha leaves; rather a suite of traits including greater Rubisco was related to the adaptation of pubescent phenotypes to arid environments.
Another function of leaf trichome is related to the absorption and/or reflectance of sun light (Galmés, Medrano, & Flexas, 2007;Skelton, Midgley, Nyaga, Johnson, & Cramer, 2012). In fact, Lee and Lowry (1980) compared the reflectance spectra of the adaxial face of pubescent and glabrous leaves, and suggested that the pubescent leaves reflected 10% higher light above 300 nm than glabrous ones did. The higher reflectance may effectively prevent from photo-damage especially under intermittent severe desiccation stresses (or alternating sunny and misty conditions), because stomatal closure increases the risk of photoinhibition (Galmés et al., 2007).
We could not reveal selective forces that remove pubescent phenotypes in GS. Species richness was greater and understory was denser (personal observation) in GS than in PS, which implied a relatively stable and stressless environment in GS. In such a situation, intraspecific competition may become important as selective pressure than physical stress such as droughts, given that L. recurvum regenerates simultaneously under a large gap (personal observation). Although we found no evidence that glabrous phenotypes have higher fitness under competition than pubescent ones, leaf characters of glabrous phenotypes imply higher competitive ability in terms of light acquisition. First, glabrous phenotypes showed lower LMA than those of pubescent ones (Table 7), and lower LMA can increase leaf area index and light interception in a given carbon resource (Schieving & Poorter, 1999). Second, glabrous leaves can more effectively absorb photosynthetic active radiation than pubescent leaves do (Lee & Lowry, 1980), and also allow direct and scattered light into the canopy at a standing community level. In addition, no investment to leaf trichomes may allow a higher growth rate. If competition and desiccation stress form two ends of divergent selection, populations consisting of a continuum of phenotypes might be formed under an environment where both pressures are not crucial (Kitayama et al., 1997). This is particularly so when populations are young to receive selective forces.

| Environmental factors of shaping phenotypic variation
Divergent selection that leads to phenotypic variation in L. recurvum appears to be related to the interaction between spatial heterogeneity of topographic conditions and temporal fluctuation of rainfall. First, the alpine zone of Mt. Kinabalu consists of highly heterogeneous topographic characters resulting from geological youthfulness F I G U R E 3 A relationship between A area (assimilation rate per unit leaf area) and Tr (transpiration rate) among sampled individuals. Each point indicates one replication (two replications per individual). Closed and open circles indicate phenotypes of Classes 5 (the most pubescent) and 0 (glabrous), respectively. Solid and dashed lines are significant regression lines (p < .01) by a linear mixed model analysis with A area as a dependent variable, and Tr and pubescence as independent variables with identity of individuals as a random term. Pubescence was fitted as a binomial fixed term, that is, Classes 0 (glabrous) and 5 (the most pubescent). Significant difference in intercept was detected at p < .01 and a slow soil accumulation. These conditions are likely to shape large difference of soil water-holding capacity between sites. Furthermore, the fluctuation of rainfall associated with El Niño Southern Oscillation adds to the heterogeneous topography as a factor of shaping an environmental gradient. If rainfall of subalpine zone of Mt. Kinabalu was low throughout a year and below a potential evapotranspiration, most part of this zone would chronically suffer from a desiccation stress. In contrast, if there was no occasional drought with stable and ample rainfall, desiccation stress would rarely occur at any sites. However, in reality, droughts occur rather regularly at 5-6 year intervals on a long-time scale (Jaksic, 2001;Kitayama et al., 2014), which translates to the difference in the level of desiccation among sites depending on water storage determined by topographic conditions. In sites with poor water-holding capacity, strong desiccation stress could occur even during droughts, whereas in sites with sufficient water-holding capacity soil moisture would not decline to a crucial level during droughts. Thus, the spatial and temporal complexity of the subalpine zone of Mt. Kinabalu may contribute to shaping the divergent selection as a proximate factor.
The magnitude of gene flow may also add to the desiccation-related spatial and temporal patterns of sympatric phenotypic variations. Much greater median gene flow 4Nm from GS to PS may suggest that migration of genes from glabrous populations is ample to supplement glabrous phenotypes that are prone to desiccation.
The interaction between El Niño droughts and topographic heterogeneity in tree mortality was studied in Bornean lowland rainforest (e.g., Itoh et al., 2012;Potts, 2003). However, the interaction may be more important in the subalpine region than in the lowland. That is because magnitude of drought was higher in the subalpine zone than in the lower vegetation zone on Mt. Kinabalu, and because the topographic heterogeneity was high in the subalpine zone reflecting a slow soil accumulation and previous glacial erosion in the Last Glacial Maximum (Hope, 2004).

| CONCLUSIONS
Leaf trichome density of L. recurvum was a genetically based trait, and that there was a substantial level of gene flows between the populations of different phenotypic compositions. These results suggested that the process shaping the phenotypic polymorphisms involved strong ongoing divergent selection. Furthermore, physiological measurements suggested that pubescent phenotypes had higher photosynthetic WUE than glabrous phenotypes, and our field census revealed a difference between phenotypes in mortality after an El Niño drought. These results together suggested that pubescent phenotypes were adaptive variation to the environment subject to droughts at least at a leaf level. In addition, soil moisture in the sites dominated by pubescent type (PS) was lower than in the sites dominated by glabrous type (GS) during a short drought. This implied that El Niño droughts could cause large difference in soil moisture among sites and that a greater desiccation stress removed glabrous phenotypes as natural selection in PS leading to the dominance of pubescent phenotypes. Our study, however, could not elucidate physiological mechanisms of the other end of divergent selection which could select for glabrous phenotypes.