Co-expression and Transcriptome Analysis of Marchantia polymorpha Transcription Factors Supports Class C ARFs as Independent Actors of an Ancient Auxin Regulatory Module

We performed differential gene expression (DGE) and co-expression analyses with genes encoding components of hormonal signaling pathways and the ∼400 annotated transcription factors (TFs) of M. polymorpha across multiple developmental stages of the life cycle. We identify a putative auxin-related co-expression module that has significant overlap with transcripts induced in auxin-treated tissues. Consistent with phylogenetic and functional studies, the class C ARF, MpARF3, is not part of the auxin-related co-expression module and instead is associated with transcripts enriched in gamete-producing gametangiophores. We analyze the Mparf3 and MpmiR160 mutant transcriptomes in the context of coexpression to suggest that MpARF3 may antagonize the reproductive transition via activating the MpMIR11671 and MpMIR529c precursors whose encoded microRNAs target SQUAMOSA-PROMOTER BINDING PROTEIN-LIKE (SPL) transcripts of MpSPL1 and MpSPL2. Both MpSPL genes are part of the MpARF3 co-expression group corroborating their functional significance. We provide evidence of the independence of MpARF3 from the auxin-signaling module and provide new testable hypotheses on the role of auxin-related genes in patterning meristems and differentiation events in liverworts.


INTRODUCTION
Transcriptome studies in model systems provide a foundation to characterize genetic pathways in the absence of mutational studies. Transcriptome co-expression studies can highlight how specification of complex phenotypes depends on the activity of coordinated batteries of regulatory genes. In plants, many co-expression analyses focused on secondary metabolite production (Ruprecht et al., 2011;Ruprecht and Persson, 2012;Cassan-Wang et al., 2013;Ferreira et al., 2016;Turco et al., 2017), although co-expressed modules for hormonal genes (Ruprecht et al., 2017b;Singh et al., 2017), or developmental factors have also been described (Ichihashi et al., 2014). The functional significance of co-expression modules has been tested by further differential gene expression (DGE), protein-protein interaction (Kemmeren et al., 2002;Ichihashi et al., 2014), epigenetic (Turco et al., 2017) or functional analyses (Yokoyama et al., 2007;Ruprecht et al., 2011;Cassan-Wang et al., 2013;Ichihashi et al., 2014). Furthermore, studies have corroborated conservation of sets of co-expressed genes in closely related species (Ficklin and Feltus, 2011;Ichihashi et al., 2014;Netotea et al., 2014;Ferreira et al., 2016) as well as in deep evolutionary time (Stuart et al., 2003;Gerstein et al., 2014;Ruprecht et al., 2017a). From a simple ancestral genetic network, gene/genome duplication events could trigger the formation of novel co-expression clusters that can be re-wired and co-opted to pattern novel biological processes (Ruprecht et al., 2017b).
With the establishment of model bryophytes, an increasing number of transcriptome datasets describing multiple developmental and environmental conditions in mosses (Devos et al., 2016;Ortiz-Ramírez et al., 2016;Perroud et al., 2018) and liverworts (Alaba et al., 2015;Frank and Scanlon, 2015;Higo et al., 2016), are providing novel avenues to perform reverse genetic studies in non-vascular plant lineages. Marchantia polymorpha is a dioecious liverwort model system  amenable to gene editing (Sugano et al., 2014) and gene silencing (Flores-Sandoval et al., 2016) techniques. Regulatory genes encoded in the M. polymorpha genome exhibit little genetic redundancy, with many transcription factor (TF) families represented by single paralogs (Bowman et al., 2017), making it an attractive system to characterize gene families that display genetic redundancy in other systems. The life cycle of M. polymorpha is haploid dominant (Figure 1A), with the diploid generation dependent upon the maternal haploid plant. Haploid spores germinate producing a sporeling, a developmental stage with active cell proliferation that forms a simple polarity between anchoring rhizoids and photosynthetic cells (Shimamura, 2016). In accordance with previous work (Bowman et al., 2017;Boehm et al., 2017), day 1 sporelings are single celled but actively differentiate chloroplasts, day 2 sporelings initiate a rhizoid cell, in day 3 sporelings the rhizoid cell elongates and the first divisions of the shoot occur, clearly establishing anatomical apical-basal polarity. By day 4, sporelings undergo proliferation of photosynthetic shoot tissue and further elongation of the rhizoid ( Figure 1A). Using light cues, sporelings transition into prothalli, a two-dimensional heart-shaped body in which a meristematic zone with an apical cell producing derivatives in two planes is established . When the apical cell ( Figure 1B) transitions to producing derivatives in four planes (top, bottom and lateral planes) a three-dimensional vegetative thallus ( Figure 1C) is formed. Growth in thalli occurs at apical regions with differentiated tissues produced from apical cell derivatives and the apical meristem dichotomously branching ( Figure 1C) at regular intervals (Shimamura, 2016;Solly et al., 2017). Upon inductive far-red to red light ratios (Kubota et al., 2014;Inoue et al., 2016;Linde et al., 2017) apices differentiate umbrella-like sexual gametangiophores (Figures 1D,E) wherein gametes (Figures 1F,G) are produced (Higo et al., 2016;Koi et al., 2016;Rövekamp et al., 2016;Yamaoka et al., 2018). In an moist environment, a sperm cell swims and fertilizes an egg cell, thus forming a zygote that undergoes several rounds of cell proliferation ( Figure 1H) before differentiating sporogenous ( Figure 1I) tissue (Shimamura, 2016).
Auxin is a key phytohormone that triggers multiple aspects of land plant physiology and development in a context dependent fashion. Comparative analysis between charophycean algae and land plant genomes and transcriptomes indicates that the canonical F-box mediated auxin signaling pathway evolved in the land plant ancestor (Bowman et al., 2017). Surprisingly, charophycean algae are able to transcriptionally respond to exogenous auxin (Hori et al., 2014;Ohtaka et al., 2017;Romani, 2017;Mutte et al., 2018) even in the absence of specific TIR1 F-box receptors, AUX/IAAs with auxin-sensitive degron (II) domains or distinct class A and B ARFs (Bowman et al., 2017;Flores-Sandoval et al., 2018;Mutte et al., 2018). M. polymorpha possesses single paralogs of all the canonical auxin-signaling pathway genes, retaining the predicted ancestral embryophyte state Kato et al., 2015;Bowman et al., 2017;Mutte et al., 2018). In M. polymorpha, the class A ARF (MpARF1) is necessary to elicit physiological and transcriptional responses to exogenous auxin, acting as a transcriptional activator (Kato et al., 2017). The class B ARF (MpARF2) has not yet been characterized but is regarded as a transcriptional repressor (Kato et al., 2015). The class C ARF (MpARF3) is a target of the embryophyte-specific miRNA160 and antagonizes differentiation events in multiple developmental contexts with loss-of-function Mparf3 alleles able to respond to exogenous auxin (Flores-Sandoval et al., 2018;Mutte et al., 2018), and strong gain-of-function MpARF3 alleles resembling auxin-depleted mutants Flores-Sandoval et al., 2018). Phylogenetic analysis indicates Class C ARFs evolved in a charophycean algal ancestor prior to the emergence of the auxin-signaling module (Flores-Sandoval et al., 2018;Mutte et al., 2018). Two additional genes, the NON-CANONICAL ARF (MpNCARF) and NON-CANONICAL IAA (MpNCIAA), are predicted to form part of the auxin-mediated transcriptional response (Flores-Sandoval et al., 2018;Mutte et al., 2018). NCARFs evolved from a class A-like ARF in the embryophyte ancestor, and in M. polymorpha act synergistically with MpARF1 to promote auxin signaling despite lacking a B3-DNA binding domain (Flores-Sandoval et al., 2018;Mutte et al., 2018). NCIAA genes form a sister clade to the AUX/IAA auxin signaling repressor genes, and as with the case of the class C ARFs, evolved prior to the origin of land plants (Flores-Sandoval et al., 2018;Mutte et al., 2018). Mpnciaa mutants are also sensitive to exogenous auxin (Mutte et al., 2018) suggesting they are not essential for auxin signaling. Despite conservation of these components across land plants, it is not known whether these genes act in common co-expression groups and whether their co-expression reflects functional dependency. Furthermore, the low number of annotated TFs, aMIR precursors and hormonerelated genes in M. polymorpha allow accessible computation of a regulatory co-expression matrix, which perhaps is difficult in angiosperm model systems due to their high genetic redundancy. In this study we create a preliminary expression landscape of TFs defining developmental transitions in M. polymorpha based on available transcriptomic data. With co-expression analysis, we define gene clusters related to different tissues and Young wild-type sporophyte at similar to the sequenced 13-day sporophytes from Frank and Scanlon (2015). (I) Mature wild-type sporophytes protruding from archegoniophore. Asterisks = apical notches and gc = gemmae cups. developmental stages. Moreover, we reveal a cluster of auxinrelated genes that are validated by significant overlap with auxininducible transcripts. We further show that the class C auxin response factor is expressed independently of this cluster and may negatively influence reproductive transitions in M. polymorpha.

Co-expression Analysis
Pearson's correlations were calculated for all Marchantia TFs and for genes with putative connection with auxin biology: correlation coefficients were calculated amongst TFs and whole genome co-expression partners were calculated only for auxinrelated genes. Gene sets were generated considering those with correlation > 0.8 and P-value < 0.001 as significantly coexpressed. Venn and Euler plots were generated using R package VennDiagram as well as VENNY (Oliveros, 2007). UpSet plots were generated using its respective R package (Lex et al., 2014). Heatmaps for TF co-expression matrixes were generated using Heatmapper (Babicki et al., 2016) by applying average linkage clustering to rows and columns and Euclidian distance measurements.

Enrichment Analysis
Comparison of co-expression groups with differentially expressed genes were performed by enrichment analysis using Fisher's exact test. For the auxin enrichment analysis we manually selected auxin related co-expression groups (and unrelated genes as negative controls) as a bait against: Auxin inducible genes from Physcomitrella patens (p.adjust < 0.05; Lavy et al., 2016), Arabidopsis thaliana (p.adjust < 0.01; Omelyanchuk et al., 2017) and M. polymorpha (p.adjust < 0.05; Mutte et al., 2018). For the first two species listed we determined orthologs using Phytozome 1 . −log 10 P-values were reported in each Figure.

GO and Protein Family Enrichment
Gene sets were functionally characterized via GO term enrichment analysis using an in-house Fisher's exact test algorithm implemented in R (v3.3.2). We compiled Biological Process related GO term annotation for M. polymorpha genes using annotations available in Phytozome for M. polymorpha and A. thaliana and P. patens orthologs. We included in our GO term database both one-to-many and many-to-one relationships in order to obtain a more accurate annotation. For protein families, we used the same algorithm based on a different database. In this case, we annotated protein families using the HMMscan algorithm implemented in HMMER3 (Wheeler and Eddy, 2013) over the M. polymorpha proteome with default parameters and using Pfam profiles 2 .

Transcription Factors Controlling Developmental Transitions in M. polymorpha
We performed DGE analysis on publicly available RNA-Seq libraries to (1) identify transcription factors defining or influencing Marchantia developmental transitions and (2) use the obtained datasets as references for enrichment analysis (Frank and Scanlon, 2015;Higo et al., 2016;Bowman et al., 2017). Datasets span multiple tissues throughout the life cycle, including intensively sampled sporelings during the first 5 days of development (Bowman et al., 2017), apical cells (Frank and Scanlon, 2015), gametangiophores, antheridia (Higo et al., 2016), and 13-day old sporophytic micro-dissections (Frank and Scanlon, 2015). In our analysis, 82% of differentially expressed TFs with P < 0.01 (edgeR) have fold changes (FC) below 4× in most tissues except sporophytes (Figure 2A and Supplementary Figure S1A), suggesting feedback regulation preventing drastic changes in regulatory gene expression. To account for both read abundance and fold-changes, we used the product of logFC and logCPM (counts per million) to categorize the transcripts enriched in a particular tissue. Of the 405 TFs in the M. polymorpha genome (Supplementary Table S1), 45 were not differentially expressed in any of the performed pair-wise comparisons (Supplementary Figure S1B).
In order to find TFs specific to a developmental transition, we compared all upregulated TFs in each tissue using a cut off of logFC and logCPM > 0 ( Figure 2B and Supplementary Table S2). Mature sporophyte-enriched TFs sequenced by JGI (Bowman et al., 2017) were included in this comparison to validate the 13-day microdissected sporophytes with both datasets showing a high degree of overlap ( Figure 2B). Mature sporophytes have the largest number of uniquely upregulated TFs (Figure 2B), followed by 24-h germinated spores, gametangiophores and thalli. The apical cell and antheridia had a single uniquely upregulated TF each ( Figure 2B). As expected, antheridia and antheridiophores, male and female gametangiophores and thalli/24-h sporelings show the highest degree of similarity. Surprisingly, 13-day sporophytes, but not mature sporophytes, had common and exclusive upregulation with archegoniophores and antheridia. Finally, antheridia, antheridiophores and archegoniophores formed a large group of uniquely upregulated TFs ( Figure 2B). The TFs enriched in each tissue library will be described in detail in the following paragraphs.

Reproductive Transition
A total of 126 TFs (87 upregulated, Supplementary Table S9) are shared with similar dynamics in thalli vs. archegoniophore and thalli vs. antheridiophore comparisons ( Figure 4C). Nine of these genes were enriched in no other tissue but gametangiophores albeit being mildly upregulated ( Figure 2B and Supplementary Table S2). Thirteen genes were exclusive to gametangiophores and antheridia, corroborating their specificity in reproductive roles ( Figure 2B and Supplementary Table S2). MpBNB is in this former group, consistent with its role in both gametangiophore and gamete formation (Yamaoka et al., 2018). Additional genes were enriched but not specific to gametangiophores (Supplementary Table S2). A majority of these TFs ( Figure 4C) have 1:1 ratios of logCPM between male and female tissues with only MpR2R3-MYB1 being higher in females (logCPM female/male = 1.3) and MpBZR3 and MpR2R3-MYB6 preferentially expressed in males (logCPM male/female = 2.5 in both cases). Despite their high expression in males, complete MpR2R3-MYB1 transcripts were exclusively found in females, suggesting male expression may be represented by incomplete transcripts.
There are additionally seven genes uniquely upregulated in 13-day sporophytes (Figure 2B)

Marchantia Regulatory Genes Are Organized in Co-expression Groups
In order to predict functional roles of clusters of TFs sharing similar expression dynamics, we created a co-expression matrix comparing averages of normalized reads per kilobase per million (RPKMs) across 11 different tissues. Pearson coefficients were calculated for each of the annotated Marchantia TFs, microRNA (miR) precursors as well as hormone biosynthetic and signaling genes. Two stages of sporeling development (48 and 96 h), multiple thallus RNAs, archegoniophores, antheridiophores, antheridia and older sporophytes (JGI) were used in this analysis. Using a minimal RPKM value of 5 as a threshold, a Pearson coefficient heatmap ( Figure 6A) suggests at least 9 highly correlated groups, some of which show a greater degree of correlation than others. A distance dendrogram (RPKM < 1) supported the occurrence of multiple groups (Supplementary Figure S8) and we focused on describing nine of the most conspicuous in basis of their tissue-specificity or putative functional relationships (Supplementary Figure S9). The antheridia (group 1) and sporophyte-enriched genes (group 3) combined comprise ∼1/2 of the total sample. This is in agreement with these tissues being the most divergent in terms of form and physiology.

Testing Co-expression Group Robustness in Thalloid Haploid Tissues
To examine co-expression group robustness and mitigate the biasing effects elicited by the non-thalloid sporelings, antheridia, and diploid sporophytic tissues, we generated a coexpression matrix without these three groups (Supplementary Figure S11). In this matrix, MpARF1 forms a co-expression group with MpGRF, MpMADS2, MpCOI, and MIR11707, an amiR targeting MpAGO1 (Tsuzuki et al., 2015;Lin et al., 2016), that remains correlated with the auxin group. As expected, the auxin group now includes MpTAA and MpNCIAA but excludes MpNCARF and MpAUX/LAX. The biggest coexpression groups in this matrix are formed by reproductive genes or antheridiophore-specific genes (Supplementary Figure S11). Within the reproductive group, former groups 2 and 7 have merged to include MpARF3, MpSPL1-3, MpC3HDZ, MpCAMTA, MpNCARF, and MpTPL, which was previously associated with antheridia. MpAUX/LAX is grouped with the light signaling group, which has gained MpLUX, formerly associated with sporophytes. MpWIP, MpKAN, and MpRR-A form their own group, which has overlap with the auxin group.

Co-expression Groups Are Dynamic at the Earliest Stages of Development
To examine how co-expression groups form at the earliest developmental transitions, we generated a matrix exclusive to sporeling RNA libraries. Multiple sporeling co-expression groups are formed that have different affinities with general coexpression groups (Supplementary Figures S2, S3, S5, S12). Auxin-related genes in particular show independent expression patterns prior to their clustering in later developmental stages. For example, the auxin signaling repressor MpIAA forms a wellsupported group (Pearson > 0.9) with genes that putatively promote cell proliferation such as the class I TCP gene MpTCP1

UpSet diagrams (Supplementary
Consistent with these results, the MpARF1, MpARF2, MpNCARF, and MpIAA groups share six GO-terms by biological process (Supplementary Figure S14A) and four GO-terms (biological process) are shared between MpARF1, MpARF2 and MpIAA. In contrast, the MpARF3 group showed 12 specific GOterms (Supplementary Figure S14B) while the MpNCIAA group had 11 specific GO-terms by biological process (Supplementary  Figure S14C).

Analysis of MpARF3 Gain and Loss-of-Function Mutant Transcriptomes Confirms Validity of Co-expression Groups
We performed DGE analysis (edgeR and DESeq2) in previously isolated Mparf3/Mpmir160 CRISPR mutants (Flores-Sandoval et al., 2018) to test whether: (1) loss of MpARF3 disrupts the auxin co-expression group, (2) genes co-expressed with MpARF3 are altered in mutant alleles, and (3) Mparf3 or Mpmir160 transcriptomes could facilitate identification of genes responsible for developmental aspects of their respective mutant phenotypes. The edgeR analysis identified 2874 differentially expressed genes in Mparf3 compared to the wild type (P < 0.01, 1524 upregulated). Similarly, 2864 differentially expressed genes were identified in comparisons of MpmiR160 mutants with the wild type (P < 0.01, 1572 upregulated). To filter out steady-state false positives and highlight genes dependent on MpARF3 activity, we searched for genes with converse behaviors in the Mparf3 and Mpmir160 transcriptomes. Genes downregulated in Mparf3 and upregulated in MpmiR160 mutants were considered to be activated, either directly or indirectly, by MpARF3 (Category 1, N = 112; Supplementary Table S18). Conversely, genes upregulated in Mparf3 and down-regulated in Mpmir160 mutants were considered repressed, either directly or indirectly, by MpARF3 (Category 2, N = 115; Supplementary Table S19). Strikingly, genes with similar dynamics, i.e., upregulated (Category 3, N = 501; Supplementary Table S20) or down-regulated (Category 4, N = 470; Supplementary Table S21) in both mutants were ∼4× more abundant, suggesting a large cohort of genes affected by secondary feedback effects independent of the MpmiR160/MpARF3 regulatory module ( Figure 8A). The most prevalent protein family found in Category 1 was peroxidases (PF00141, Supplementary  tandem arrays in scaffold 118 and with high resemblance to fungal lectin genes (Bowman et al., 2017).
Although no TFs were identified as activated by MpARF3 (Category 1), two miR precursors fell into this category, with MpMIR160 (Mapoly0002s0211) and MpMIR11671 (Mapoly0239s0004) supported by both edgeR and DEseq2 analyses ( Figure 8A). Meanwhile, a third miR precursor, MpMIR529c (Mapoly0006s0020) was also identified as activated by MpARF3 using DESeq2 (Supplementary Figure S14D). We had previously suggested that MpARF3 forms a negative feedback loop with MpmiR160 to control both developmental transitions and differentiation vs. totipotency in a context-dependent fashion (Flores-Sandoval et al., 2018).
Three of the downregulated TFs (MpSPL1, 2 and MpR2R3-MYB1) as well as MpARF3 itself were classified as reproductive genes in our wild-type DGE analysis (Figure 5). Furthermore, we had previously characterized MpARF3 as an antagonist of reproductive transitions, with Mpmir160 alleles insensitive to gamentangiophore-inducing far-red light treatment and Mparf3 alleles as hypersensitive to far-red light treatment (forming more gametangiophores per area than wild type; Flores-Sandoval et al., 2018). We therefore measured the overlap between genes enriched in gametangiophores (either male or female) and the Mparf3 or MpmiR160 mutants. A significant number of genes (N = 304, P = 1.5 × 10 −82 , Hypergeometric test) were upregulated in both Mparf3 mutants and wild-type archegoniophores ( Figure 8B). Consistently, a significant number of genes were upregulated in MpmiR160 mutants and downregulated in archegoniophores (N = 307, P = 3.7 × 10 −52 , Hypergeometric test). Enrichment analysis (Fisher's exact test) using the datasets obtained in our DGE analysis (edgeR, P < 0.01) show that overall transcripts activated by MpARF3 (category 1) are similar to those downregulated in wild type gametangiophores (P < 0.01, Figure 8C). Conversely, transcripts repressed by MpARF3 (category 2) are significantly upregulated in wild type gametangiophores, in particular the subset of reproductive transition genes and not the female or male-enriched genes ( Figure 8C). Additionally, 13-day sporophytes also have transcripts that seem to be controlled by MpARF3 although with an ambiguous trend unlike that observed for reproductive tissues. Although some auxin-related or auxin-induced genes were differentially expressed in Mparf3 transcriptomes (Flores-Sandoval et al., 2018), not a single auxininduced gene had converse expression patterns in the Mparf3 and Mpmir160 transcriptomes ( Figure 8A). Finally, the expression dynamics of the candidate TF/MIR pairs shows that MpARF3, MpSPL1 and 2 expression patterns are positively correlated with peaks of expression in gametangiophores (Figure 8D), while the MpMIR11671 and MpMIR529c precursors are negatively correlated with MpARF3 and show down-regulation in gametangiophores. Thus, an inverse relationship exists between MpARF3 activity and expression of the MpSPL1-2/MpMIR529c/MpMIR11671 modules in M. polymorpha (Figures 8D, 9B).

DISCUSSION
This study provides a first attempt at characterizing DGE patterns and co-expression groups throughout the life cycle of M. polymorpha (Figure 1). Given their hierarchal regulatory character, TFs were chosen as a proxy for the thousands of genes differentially expressed between cell and tissue types. DGE analysis reveals that only a minority of TFs show dramatic foldchanges between developmental transitions in M. polymorpha. However, this figure is qualified by the limited number of tissues and environmental conditions sampled in our analysis (Figure 2). These differentially expressed TFs are candidates for functional analysis to assess whether they determine cell identity or broader processes required for tissue patterning, i.e., meiosis, cell division, expansion, stress-responses.
The relatively small number of annotated TFs (∼400) in M. polymorpha (Bowman et al., 2017) facilitated construction of a co-expression matrix using 11 tissues, including sporelings, thallus, gametangiophores, antheridia and sporophytes. This matrix resolved multiple discrete groups with variable levels of overlap and putative interactions ( Figure 6). As expected, some groups correspond to tissue-specific TF enrichment (e.g., sporophytes, antheridia, archegoniophores or regenerating tissue). Surprisingly, a second kind of group that did not show enrichment in any one type of library (Supplementary Figures S9E,G,I,L) was also identified. These groups include the partners of MpARF1 or MpTOC1, which are involved in auxin and light/circadian responses, respectively. Members of these co-expression groups were differentially expressed in multiple stages but did not change dramatically across our pair-wise comparisons, likely due to lack of inductive environmental conditions or because they are under tight regulatory feedbacks. Thus, our data suggests developmental transitions are also defined by the additive influence of multiple members of a co-expression group.

TF Functions Predicted by Our Analysis
Totipotency A set of TFs that could represent regulating factors controlling aspects of totipotency include the LAS ortholog MpGRAS4, MpANT, and MpARF2. MpGRAS4 is enriched in the apical cell and both 24 and 48 h sporelings. It is further nested in a co-expression group enriched in 24-h regenerating tissue following wounding (Supplementary Figure S9F). Thus, MpGRAS4 is expressed in tissues united by active cell proliferation. Second, the ANT/PLT/BBM ortholog, MpANT, is enriched exclusively in apical cells and 96 h sporelings, raising the possibility that its action is downstream of MpGRAS4. MpANT is also co-expressed in the auxin cluster (Supplementary Figure S10), suggesting that it may be auxin-inducible, as occurs in angiosperm model systems (Galinha et al., 2007;Ding and Friml, 2010). Consistent with this scenario, MpARF2 is also an apical-cell enriched gene whose expression steadily increases from 72-h sporelings (in parallel with MpARF1, Supplementary Figure S4A). It is plausible, given that MpARF1 and MpARF2 act as respective activators or repressors of transcription (Kato et al., 2015), that antagonism between MpARF1 and MpARF2 at 72 h stabilizes MpANT expression in 96 h sporelings, subsequently forming a co-expression group throughout development. Given that LAS has key roles specifying axillary meristems in angiosperms, and that ANT/PLT/BBM orthologs act in meristems in both angiosperms (Aida et al., 2004) and mosses (Aoyama et al., 2012), there might be more overlap between haploid and diploid meristems than previously acknowledged (Frank and Scanlon, 2015).

Cell Proliferation and Auxin
Additional genes may play a role in cellular processes that contribute to, but do not specify, meristematic tissues. For example, a cluster represented by the 24-h sporelings involving MpTCP1, MpGRF, MpC4HDZ and MpIAA genes could represent a cell proliferation group. MpGRF in particular is highly expressed in the first day of sporeling germination, but it is also co-expressed with auxin-related genes throughout the life cycle ( Figure 6A and Supplementary Figure S11). Auxin signaling in 24-h sporelings may be constrained by the presence of the known MpIAA auxin signaling repressor Kato et al., 2015). Given that key aspects of auxin physiology in M. polymorpha involve organ differentiation Flores-Sandoval et al., 2015;Kato et al., 2015), it is tempting to speculate that there may be an antagonism between auxin-dependent differentiation and cell proliferation elicited by MpGRF and other TFs in M. polymorpha.

Differentiation Factors
An auxin-inducible gene (Mutte et al., 2018) characterized in M. polymorpha is the IDDL-like Zinc Finger MpWIP TF, which is essential for air pore differentiation (Jones and Dolan, 2017). Consistently, MpWIP is upregulated in 96-h sporelings at the time of photosynthetic tissue proliferation (Supplementary Figure S5). MpWIP is co-expressed (Supplementary Figure S9H) with the GARP TF MpKAN, whose orthologs influence organ polarity in angiosperms (Eshed et al., 2001). MpKAN and its co-expression partner MpWRKY8 are in turn clearly enriched in 0-h cut thalli that lack apical notches and haven't initiated regeneration. One hypothesis is that members of this group may regulate key aspects of differentiation possibly downstream of auxin and light signaling, given their overlap in expression patterns (MpKAN and MpARF1 have a Pearson coefficient of 0.76, while MpPIF and MpKAN have a coefficient of 0.8).

Reproductive Genes in Marchantia
Comparisons between gametangiophore datasets facilitated distinctions between reproductive transition TFs (enriched in both sexes) vs. female or male-enriched TFs. MpBNB served as a marker to validate this group given its role in promoting formation of gametangiophores irrespective of sex (Yamaoka et al., 2018). A co-expression group including MpSPL1, MpSPL2, MpC3HDZIP, MpABI3B, MpARF3, among others correlates with the vegetative to reproductive transition and is supported by mutant transcriptomes (see below section). Interestingly, auxin-repressed genes are significantly represented in archegoniophore and antheridiophore-enriched transcripts (Supplementary Figure S4B), suggesting a repressor(s) of auxin signaling is a member of the reproductive transition group. Maleenriched TFs were refined by the presence of antheridia libraries, which suggest a robust set of male gamete-enriched TFs and these have been extensively described in previous studies (Higo et al., 2016). MpRWP2/MpMID could provide an interesting candidate for an antheridia specific factor given that its orthologs specify minus gametes in green algae (Ferris and Goodenough, 1997).

Extended Parental to Zygotic Transition in M. polymorpha
The identification of female and antheridia-enriched TFs allowed comparisons with upregulated TFs in 13-day sporophytes, consisting of ∼10 cells and at which point sporogenous tissue differentiation has not commenced (Frank and Scanlon, 2015). Female-enriched TFs shared with young sporophytes include members of 10 TF families. Male-enriched TFs shared with young sporophytes involve members of 7 TF families, although MpBELL4 is the only identified TF exclusively upregulated in antheridia and young sporophytes ( Figure 2B). Figure S6) suggests that only archegoniophore-enriched genes are significantly represented in 13-day sporophyte transcriptomes. Thus, the possibility of RNA contamination from maternal tissues still remains as a plausible explanation to this phenomenon.

Class C ARF and MpNCIAA Genes Are Independent of Other Auxin Related Genes
The discovery of a putative auxin co-expression group among transcription factors (Figure 6) coincides with a similar PIN1related co-expression group observed in moss (Ruprecht et al., 2017b), suggesting an ancestral origin for the auxin coexpression network in land plants. The robustness of the auxin TF co-expression group was tested using whole-genome coexpression analysis for candidate genes and comparisons with the only available auxin-inducible transcriptome in M. polymorpha (Mutte et al., 2018). Enrichment analysis supports significant overlap between MpARF1, MpARF2, MpIAA, MpNCARF, MpYUC2, MpSAUR1, MpPIN1, MpGH3A, and MpILR1 coexpression groups and the M. polymorpha auxin-upregulated transcriptome (Figures 7A,C). Other TFs within the auxin coexpression group but strictly not involved in auxin biosynthesis, transport, conjugation or responses include MpANT, MpSHI, MpbHLH10, 43, 45 and 47, MpMIR11707, MpR2R3-MYB8, and MpC2HDZ ( Figure 7C and Supplementary Figure S12). They could represent upstream or downstream elements regulating the module (Mutte et al., 2018) and are suitable candidates for future characterization. Importantly, many of these genes are also part of the auxin-inducible transcriptome of Physcomitrella but not Arabidopsis (Figure 7C), suggesting they were key components of the auxin co-expression network in the ancient bryophyte ancestor. As MpbHLH43/LRL possibly mediates ectopic rhizoid formation (Breuninger et al., 2016) in response to exogenous auxin in M. polymorpha its association was anticipated. Furthermore, 5 of the 14 TFs induced by auxin (Supplementary Table S24) in M. polymorpha are found within the MpARF1 and MpARF2 co-expression groups ( Figure 7B). Notably, the class C ARF, MpARF3, a gene that is not necessary for physiological and transcriptional auxin responses (Flores-Sandoval et al., 2018;Mutte et al., 2018), forms part of an independent co-expression group highly active in reproductive structures (Figure 8 and Supplementary Figure S9). However, MpARF3 still shares members with MpNCARF (72 genes), MpARF1/MpARF2/MpNCARF/MpIAA (21 genes) and MpNCIAA (9 genes, Supplementary Figure S13) co-expression groups. The significance of this connectivity was tested (Supplementary Table S25), providing ways in which MpARF3 competes for MpARF1 and MpARF2 targets in an auxin-independent manner. Importantly, the MpARF3 and MpNCARF co-expression groups are enriched in reproductive transition genes (Figure 7A), opening the possibility of interaction at this developmental stage. MpNCIAA is another PBI-containing factor that forms an independent group with sporophyte-enriched transcripts (Figure 7 and Supplementary Figure S9) and is also not essential for transcriptional auxin responses (Mutte et al., 2018). Consistent with their auxin independence, both class C ARFs and NCIAAs evolved before the origin the canonical auxin-signaling pathway (and land plants), and thus they may act in independent regulatory networks. The putative roles of NCIAA genes in patterning sporophytic development are of interest although these genes are also expressed in haploid tissues in M. polymorpha.

RNA-Seq Data Supports the Role of MpARF3 as an Inhibitor of Reproductive Transitions
We have previously described MpARF3 as a repressor of differentiation in multiple developmental contexts. For example, Mparf3 mutants display developmental transition, defects, such as ectopic differentiation of air pores, scales, pegged rhizoids, and gametophores, and are deficient in forming undifferentiated gemmalings (Flores-Sandoval et al., 2018;Mutte et al., 2018). Meanwhile strong MpARF3 gain-of-function alleles form undifferentiated calli whose cell identity resembles that of young sporelings or prothalli (Flores-Sandoval et al., 2018). Thus, we anticipated the transcriptome data should reflect aspects of these developmental defects. Using gain-and loss-offunction MpARF3 transcriptomes, were able to discern secondary feedback, treatment-dependent and steady state effects and thus assign with more confidence genes genuinely dependent on the MpMIR160/MpARF3 regulatory module (category 1 and 2 genes, Figure 8). Indeed, enrichment analysis demonstrated that MpARF3 represses, directly or indirectly, genes that promote reproductive transitions (Figure 8). Consistently, MpARF3 activates genes that inhibit production of, or at least are downregulated in, sexual organs. We therefore propose the hypothesis (Figure 9A) that MpARF3 inhibits the reproductive transition in M. polymorpha via activation of the MpMIR11671 and MpMIR529c precursors, which in turn target MpSPL1 and MpSPL2. This appears plausible, given the observed MpSPL1 and MpSPL2 RPKM values in Mparf3/mir160 mutant backgrounds ( Figure 9B) and the roles of SPL genes in regulating heteroblastic changes and specifically in promoting reproductive transitions in angiosperms (Huijser and Schmid, 2011) and possibly mosses (Cho et al., 2012).

AUTHOR CONTRIBUTIONS
EF-S generated data for Figures 1-5. EF-S and FR generated data for Figures 6-8. FR created scripts for enrichment analysis and GO-terms pipelines. EF-S designed all figures, except Figures 6, 7, which were jointly designed by EF-S and FR. EF-S and JB co-wrote the manuscript and interpreted the data. All authors analyzed and revised the data.

FUNDING
This work was made possible by funding from the Australian Research Council (DP160100892 and DP170100049 to JB) and by the Agencia Nacional de Promoción Científica y Tecnológica (PICT2013-3285 and PICT2017-1484).

ACKNOWLEDGMENTS
We thank all past and current members of the Bowman Lab for feedback and helpful discussions. The Monash University Bioinformatics platform provided helpful guidance and feedback. The Monash Next Generation sequencing facility (Micromon) sequenced the Mparf3 and Mpmir160 transcriptomes. We thank Sandra K. Floyd, Tom Dierschke, and John Alvarez for providing images for   Table S10) are also statistically represented in 13-day sporophyte transcriptomes. (B) Venn Diagram of shared TFs between antheridia, archegoniophore, 13-day and JGI sporophytes. Putative maternal, paternal and common TFs continuing expression in 13-day sporophytes are annotated. The + symbol indicates upregulation (logFC > 0). FIGURE S7 | UpSet diagram of tissue-specific TF enrichment at logFC > 1. Red asterisk shows eight uniquely upregulated TFs in apical cells at logFC > 1 compared to other tissues.  Table S1) in 11 tissue libraries using RPKM > 1 as a threshold.  Table S16) RPKM values above 0.75 (to include MpANT) in sporeling tissues. Putative groups were delimited by sets of genes sharing coefficients above 0.8. Genes putatively involved in auxin biosynthesis, perception, signaling and transport are mapped, together with accompanying partners in smaller fonts. (B) Distance dendrogram of tissue libraries used in the analysis. Average replicate RPKM values were used per library. Scale for heatmap indicates Pearson coefficients.  Table S17) in thalli and gametangiophores. Putative groups were delimited by sets of genes sharing coefficients above 0.8. Genes putatively involved in auxin biosynthesis, perception, signaling and transport are mapped, together with accompanying partners in smaller fonts. (B) Distance dendrogram of tissue libraries used in the analysis. Average replicate RPKM values were used per library. Scale for heatmap indicates Pearson coefficients.
FIGURE S13 | (A) UpSet diagram showing shared co-expressed genes (cut-offs indicated) between most PBI-containing genes in M. polymorpha.
Numbers of genes exclusive or shared between groups are shown. Co-expression analyses were performed for all expressed genes in our RPKM matrix. FIGURE S14 | (A) Significant GO-terms (Biological Process) found in the MpARF1, MpARF2, MpNCARF, and MpIAA co-expression groups. X-axis indicates P-Values. (B) Significant GO-terms (Biological Process) found in the MpARF3 co-expression group. X-axis indicates P-values. (C) Significant GO-terms (Biological Process) found in the MpNCIAA co-expresion group. X-axis indicates P-values. (D) TFs differentially expressed in Mparf3 and MpmiR160 transcriptomes using DESeq2 (P < 0.01 cut-off). Genes in blue are consistent with MpmiR160-dependent repression of MpARF3.
TABLE S1 | TFs, MIRs and hormonal genes and accession numbers used in this paper.