Are Cell Junctions Implicated in the Regulation of Vitellogenin Uptake? Insights from an RNAseq-Based Study in Eel, Anguilla australis

At the onset of puberty, ovarian follicles become competent to incorporate large amounts of vitellogenin (Vtg). Using an RNAseq-based approach, transcriptomes from pre-vitellogenic (PV) and early vitellogenic (EV) ovaries from wild-caught eel, Anguilla australis, were compared to investigate the expression of specific genes encoding cell junction proteins that could be involved in regulating Vtg uptake. Partial support was found for the mechanical barrier hypothesis proposing that the access of Vtg to the oolemma is restricted by a tight junction (TJ) network within the granulosa cell layer, which changes between the PV and EV stage. Among 25 genes encoding TJ-constituting proteins, five were down-regulated and two were up-regulated. A chemical barrier hypothesis stating that gap junctions (GJs) are involved in modulating Vtg uptake was not supported, as only five GJs were found to be expressed in the ovary with no significant changes in expression between stages. Furthermore, the endocytic pathway was found to be up-regulated during the PV-EV transition. Finally, the study showed that gene expression patterns may help identify suitable candidates involved in the regulation of Vtg uptake, and provided novel sequence data for A. australis, including putative Vtg receptors corresponding to Lr8 and Lrp13 members of the low-density lipoprotein receptor family.


Introduction
In order to produce high quality eggs that have the ability to self-sustain the development of viable offspring, female oviparous animals need to actively deposit nutrients into the developing oocytes to form the egg yolk: a mass rich in lipids and proteins that will serve as a food source for the future embryo. The developmental period associated with yolk accrual into the oocyte is known as vitellogenesis (vitellogenic growth)-it is characterised by the incorporation of the bulk of the yolk proteins that are obtained through the uptake of very high-density lipoproteins, so-called vitellogenins (Vtg), from the circulation. During this growth phase, different molecular events occurring within the ovarian follicle can affect Vtg uptake. The Vtg receptor (Vtgr)-a member of the low-density lipoprotein receptor (LDLR) supergene family-is probably the most distinct factor. However, extensive evidence in many teleost species indicates that the vtgr gene is already expressed and translated at early stages of development [1][2][3], and its transcription either decreases [3,4] or remains constant during subsequent oocyte growth [5]. Indeed, the fact that pre-vitellogenic (PV) follicles are not competent to internalise Vtg in vivo nor in vitro [6], even though the receptors are available [1], prompts the notion that other elements must be involved in activating and supporting the oocyte's ability to accumulate great quantities of Vtg [6].
Following estrogen-induced synthesis in the liver, Vtg is transported to the ovary via the bloodstream where it exits the capillaries surrounding the theca cells. Once there, it must make its way through the intercellular spaces among granulosa cells to reach the oocyte membrane [7]. Indeed, open channels between the granulosa cells, associated with oocyte developmental stage and size, were shown to provide access for exogenous proteins to the oolemma in various oviparous animals [8], including Coturnix japonica [9], Fundulus heteroclitus [10], and Xenopus laevis [11,12]. Thus, the structural changes that the granulosa cell layer undergoes during oogenesis may well represent the first barrier to affect the accessibility of Vtg to the oocyte. In this context, tight junctions (TJs) create a paracellular barrier to selectively modulate the transport of molecules and water between neighbouring cells [13]. TJs are, for example, known to effectively regulate the epithelial permeability and paracellular transport in many tissues from teleost fish [14], namely gill [15], blood-brain barrier [16], and blood-testis barrier [17]. Furthermore, data from the zebrafish (Danio rerio) support a role for TJs in follicle development in view of changes in gene expression of cell junction proteins during oogenesis and their response to endocrine regulation [18]. Hence, a stage-specific TJ network could generate a more "compact" or "leaky" tissue, as such regulating the intercellular space between the granulosa cell layer and restricting or allowing Vtg to reach the oocyte membrane, respectively.
Other structural changes pertaining to cell-cell interactions during follicular development may also be linked to the regulation of Vtg accumulation. Notably distinctive is the appearance of interdigitating projections of microvilli from both granulosa cells and the oocyte, and the assembly of gap junctions (GJs) that connect them, occurring before and during vitellogenic growth (see [19] and references therein). Additionally, the subsequent retraction of the granulosa microvilli during oocyte maturation, (after termination of Vtg uptake and completion of oocyte growth) is well-documented ( [19] and references therein). Due to this strong relationship between the granulosa cells and the oocyte during vitellogenic growth, cell junctions involved in cell-cell transport and communication, i.e., GJs, conceivably also play a role in stimulating and maintaining Vtg accrual. Indeed, gap junctional intercellular communication is required for uptake of Vtg in X. laevis oocytes [20], associated with the sharing of signalling molecules such as cyclic adenosine monophosphate (cAMP) [21].
In order to investigate the involvement of cell junction proteins in the regulation of Vtg uptake, the following hypotheses were tested: (i) a mechanical barrier hypothesis, in which Vtg uptake is regulated by a network of TJs acting as a physical barrier restricting Vtg accessibility to the oolemma, and (ii) a chemical barrier hypothesis, in which cell junctions involved in cell-cell transport act as a chemical barrier of an external signalling molecule required for Vtg uptake. Consequently, a transcriptomic approach was used to compare the expression of genes encoding cell junction-constituting proteins between ovaries in the PV stage (previous to Vtg incorporation) and the early vitellogenic stage (EV stage, right after Vtg uptake has started). Additionally, the expression of genes likely involved in the molecular pathway that dictates the fate of Vtg after reaching the oolemma was examined. This included the need to identify the sequences corresponding to Vtgr candidates in the species studied, the New Zealand shortfinned eel, Anguilla australis [22]. Like other temperate anguillid eels, this species spends most of its life in freshwater systems until it starts a long oceanic migration towards its breeding grounds. Accordingly, both immature eels and early maturing, migratory eels can be captured from the wild [23], representing the PV and the EV ovarian stages, respectively. Although artificial reproduction of eels has been achieved by applying intense hormonal treatments [24,25], egg hatchability and larval survival rates remain a challenge for a sustainable production at commercial level, which still relies heavily on natural stocks. Elucidation of important mechanisms that assure the production of high-quality eggs, notably Vtg uptake, may help improve freshwater eels' artificial reproduction and conservation.

Animals and Sample Collection
The study was approved by the University of Otago Animal Ethics Committee and was conducted in accordance with the guidelines of the Australian & New Zealand Council for the Care of Animals in Research and Teaching. Wild shortfinned eels were caught from Lake Ellesmere (South Island, New Zealand; year 2019) using fyke nets. Based on morphological features [26], PV and EV eels (n = 6 per stage) were identified, euthanised with an overdose of benzocaine (0.3 g/L, Sigma-Aldrich, Merck KGaA, Darmstadt, Germany), and dissected at a local makeshift field station. Total body weight was measured and after dissection, ovary and liver were weighed to calculate somatic indices: gonadosomatic index (GSI) and hepatosomatic index (HSI), respectively. Ovarian fragments were fixed in 4% paraformaldehyde for histological analysis (oocyte diameter measurement and confirmation of gonadal stage, Section 2.2), and flash frozen for RNA extraction and next-generation sequencing (Section 2.3).

Histological Analysis
After fixation, ovarian fragments were dehydrated and embedded in Technovit 7100 resin (Heraeus Kulzer GmbH and Co., Hanau, Germany) following the manufacturers' instructions. Sections were cut at 2 µm using a Reichert-Jung microtome (Model 2050, Cambridge Instruments GmbH, Nussloch, Germany) and stained with Polychrome I (11.3% methylene blue and 0.2% azure II) and Polychrome II (0.2% basic fuchsin) stains. Microscope images were captured with an Olympus camera SC100 and an Olympus adaptor U-TVO.5XC-3 attached to an Olympus microscope BX51 (Olympus Corporation, Tokyo, Japan). Oocyte diameter (OD) was calculated from captured images using ImageJ v1.52a (National Institutes of Health, Bethesda, MD, USA). Oocytes with a visible nucleus were measured according to Lokman et al. [27], and data from the 50% highest-ranked oocytes were retained to calculate the average OD, and to eliminate those oocytes that may have been sectioned off-center.

Total RNA Extraction, cDNA Library Preparation, and High-Throughput Sequencing
Ovarian total RNA was extracted using NucleoSpin RNA Kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) following manufacturer's instructions. The total RNA concentration was measured with the Qubit RNA Assay Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA). The relative purity level was determined using a DeNovix spectrophotometer DS-11 (DeNovix Inc., Wilmington, DE, USA), which gave a 260/280 nm absorbance ratio of about 2. The quality and integrity of total RNA was further analysed using an Agilent 5300 Fragment Analyzer (Otago Genomics Services, University of Otago, Dunedin, New Zealand) obtaining the following RNA integrity numbers (RIN): 7.1 ± 1.5 (PV stage, n = 6), and 10 (EV stage, n = 6).
The extracted total RNA was sent frozen using an Ultra-freeze MAX Bio-bottle transport system (Bio-bottle Ltd., Auckland, New Zealand) to the sequencing contractor (Beijing Genomics Institute, BGI, Hong Kong). Once there, a second quality control was performed using an Agilent 2100 Bioanalyzer yielding RIN values of 3.5 ± 0.9, and 9.6 ± 0.2 for the PV and EV stages, respectively ( Figure S1, see Section 4).
The strand-specific cDNA library construction (MGIEasy RNA Directional Library Prep Set V2.1, MGI Tech Co. Ltd., a subsidiary of the BGI Group) and subsequent highthroughput sequencing were performed by BGI. Briefly, samples were first enriched for mRNA using oligo-dT selection (Poly(A) mRNA enriched by oligo(dT) beads), and then fragmented and reverse-transcribed. After adaptor ligation and PCR amplification, cDNAs were compacted into DNA nanoballs in order to sequence 100 bp paired-end reads on a DNBseq platform.
Initial clustering was performed by visualising Euclidean distance between vsttransformed counts (variance stabilising transformation) using the pheatmap v1.0.12 [38] and the DeSeq2 v1.22.2 packages [39] in R v3.5 [40]. Differential expression (DE) analysis was performed using the edgeR 3.24.3 [41,42] and the limma package v3.38.3 [43]. Genes with small counts were filtered out using the filterByExpr function from the edgeR package. The voom function [44] was then used to transform count data to log 2 counts per million (logCPM), to normalise counts using quantile normalisation, and to estimate the meanvariance relationship before linear modelling using lmFit and eBayes functions [45]. A gene ontology (GO) enrichment analysis carried out with Goseq v1.22.0 [46] was used as an exploratory analysis. All genes with an adjusted p value after FDR correction (q value) of 0.05 [47] were considered. Later, an exhaustive search of target genes was conducted (Section 2.5) where DEGs were identified with a q value < 0.05 and a threshold of 2-fold change (2-FC, i.e., −1 ≥ log 2 FC ≥ 1).

Target Genes
Target genes encoding proteins constituting cell junctions, i.e., GJs and TJs (Table 1), as well as genes involved in clathrin-mediated endocytosis, intracellular vesicle trafficking, and Vtg processing (Table S2), were searched on the annotated transcriptome, and their deduced protein sequences were retrieved. Gene nomenclature was based on the corresponding A. anguilla top hit of a BLASTp search against the nr protein database from the National Centre for Biotechnology Information (NCBI) [48], except for the Vtgr candidates found. The A. anguilla sequences on NCBI correspond to RefSeq GCF_013347855.1, Annotation Release 100 (accessed on 26 August 2021). Then, the target genes were checked for differential expression using the previously established cut-off. Only the genes passing the filtering were used (a higher number of genes encoding cell junctions, i.e., claudins and GJs, were found in the transcriptome with very low expression to be used in DE analysis). When multiple gene IDs were found to be associated with the same functional annotation, only one gene ID was used based on their highest read counts, unless it was confirmed that they correspond to different gene sequences after nucleotide and protein alignments. Additionally, a heatmap of genes involved in the pathway associated with Vtg uptake was visualised using the R package pheatmap v1.0.12 [38] after transformation of the counts with the DeSeq2 v1.26.0 implemented regularised log transformation [39].
The protein sequences corresponding to the Vtgr candidates, deduced from complete open reading frames, were further analysed using the CD-Search tool [49] from the Conserved Domain Database on NCBI [50], and they were subjected to phylogenetic analysis. To construct a phylogenetic tree, amino acid sequences corresponding to members of the LDLR family (low-density lipoprotein receptor 8-LR8-, low-density lipoprotein receptor-related protein 13-Lrp13-, and low-density lipoprotein receptor-related protein 4-LRP4-) from different species were retrieved from NCBI (see species and accession numbers in Figure S2). After alignment of sequences using the ClustalW algorithm [51], the tree was constructed in MEGA v7 [52] using the Maximum Likelihood method based on the JTT matrix-based model [53] applying 1000 bootstrap replicates. A discrete Gamma distribution was used to model evolutionary rate differences among sites (2 categories, +G parameter = 1.7003). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 8.87% sites). All positions with less than 95% site coverage were eliminated. The LRP4 members were used as out-group. Table 1. Expression of genes encoding cell junctions during progression from the PV stage to the EV stage in A. australis ovarian tissue. Genes encoding tight junction proteins as well as gap junction (gj) proteins are shown. Both differentially and not differentially expressed genes are included, along with their respective q values and log 2 FC. Fold change in expression is represented by log 2 FC and it is relative to the PV stage. A positive log 2 FC indicates up-regulation in the PV stage and down-regulation in the EV stage, while a negative log 2 FC means up-regulation in the EV stage and down-regulation in the PV stage. Significant DEGs were considered when showing q value < 0.05 and −1 ≥ log 2 FC ≥ 1, representing a 2-FC. Gene names are based on the A. anguilla top hit (A. anguilla genome, NCBI RefSeq: GCF_013347855.1, Annotation Release 100) from the BLASTp search of deduced protein sequences retrieved from the transcriptome. The read counts (mean ± SEM) are shown for each stage (n = 6 per stage).

Statistics
Data from somatic indices and OD were tested for normality and equal variances using the Shapiro-Wilk [54] and Levene's [55] tests, respectively, and arcsine transformed when necessary (GSI and HSI). Thereafter, the data were subjected to a t-test using InfoStat v.2018 [56] to compare PV and EV eels. All gene sets studied were evaluated under the Fisher's exact test in R v3.5 [40] to check their enrichment in the DEGs group.

Characterisation of Ovarian Developmental Stage
When progressing from the PV to the EV stage, significant increases were seen in the GSI (from 0.89 ± 0.17% to 3.28 ± 0.69%, t = 7, df = 10, p value < 0.0001), the HSI (from 0.57 ± 0.04% to 1.04 ± 0.13%, t = 7.11, df = 10, p value < 0.0001), and the OD (from 109 ± 6 µm to 271 ± 6 µm, t = 14.56, df = 10, p value < 0.0001), respectively. Further histological analysis of ovaries confirmed the developmental stages assigned using morphological traits while sampling. An evident increase in oocyte size due to great accumulation of lipid droplets and yolk proteins (Vtg incorporation) was seen in EV oocytes when compared to PV oocytes, which only showed a few lipid droplets and no signs of Vtg accumulation in their cytoplasm (Figure 1).
the JTT matrix-based model [53] applying 1000 bootstrap replicates. A discrete Gamma distribution was used to model evolutionary rate differences among sites (2 categories, +G parameter = 1.7003). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 8.87% sites). All positions with less than 95% site coverage were eliminated. The LRP4 members were used as out-group.

Statistics
Data from somatic indices and OD were tested for normality and equal variances using the Shapiro-Wilk [54] and Levene's [55] tests, respectively, and arcsine transformed when necessary (GSI and HSI). Thereafter, the data were subjected to a t-test using In-foStat v.2018 [56] to compare PV and EV eels. All gene sets studied were evaluated under the Fisher's exact test in R v3.5 [40] to check their enrichment in the DEGs group.

Characterisation of Ovarian Developmental Stage
When progressing from the PV to the EV stage, significant increases were seen in the GSI (from 0.89 ± 0.17% to 3.28 ± 0.69%, t = 7, df = 10, p value < 0.0001), the HSI (from 0.57 ± 0.04% to 1.04 ± 0.13%, t = 7.11, df = 10, p value < 0.0001), and the OD (from 109 ± 6 μm to 271 ± 6 μm, t = 14.56, df = 10, p value < 0.0001), respectively. Further histological analysis of ovaries confirmed the developmental stages assigned using morphological traits while sampling. An evident increase in oocyte size due to great accumulation of lipid droplets and yolk proteins (Vtg incorporation) was seen in EV oocytes when compared to PV oocytes, which only showed a few lipid droplets and no signs of Vtg accumulation in their cytoplasm ( Figure 1).

De Novo Transcriptome and Downstream Analysis
An average of 14.25 million 100 bp paired-end reads per sample were sequenced (min: 12.73, max: 14.65). The de novo assembly generated 282,469 transcripts for 171,620 putative genes, of which 43,975 (25.6%) were functionally annotated. After filtering out

De Novo Transcriptome and Downstream Analysis
An average of 14.25 million 100 bp paired-end reads per sample were sequenced (min: 12.73, max: 14.65). The de novo assembly generated 282,469 transcripts for 171,620 putative genes, of which 43,975 (25.6%) were functionally annotated. After filtering out lowly expressed genes, 32,447 genes were kept, of which 20,810 (64.1%) were functionally annotated. Initial clustering performed using the Euclidean distance between samples was visualised on a heatmap to assess the variability within the dataset (Figure 2). The samples clustered into two distinct phenotypic groups, PV and EV, indicating that the between-group variance was greater than the within-group variance.  regulated BP terms included endocytosis (GO:0006897, q value 7.6 × 10 −6 ) and membrane organisation (GO:0061024, q value 0.02) (Figure 2, Table S1), representative genes of which are also denoted in Section 3.4. The third BP term (regulation of cation transmembrane transport, GO:1904062, q value 0.03), along with the MF term (potassium channel inhibitor activity, GO:0019870, q value 0.04), included a group of genes involved in modulating the frequency, rate or extent of cation transmembrane transport corresponding to the serine/threonine protein kinases WNK1, WNK2, and WNK3 ( Figure 2, Table S1). Among the most down-regulated GO terms during the PV-EV transition, integral component of membrane (GO:0016021, q value 3.3 × 10 −8 , cellular component -CC-category) and identical protein binding (GO:0042802, q value 3.4 × 10 −5 , MF category) (Figure 2, Table S1) included genes encoding TJ-constituting proteins denoted in Section 3.3. Finally, 4878 genes were differentially expressed with a q value < 0.05 and a threshold of 2-FC; 2851 genes were up-regulated and 2027 genes were down-regulated when progressing from the PV to the EV stage (Figure 3).

Expression of Genes Related to Receptor-Mediated Endocytosis Machinery and Vtg Processing during the PV-EV Transition
Two genes encoding putative Vtgrs were found in the A. australis de novo transcriptome based on their annotation and sequence similarity to other vertebrate Vtgrs (analysis of conserved domains from deduced protein sequences and alignment with corresponding LDLR family members [57], and phylogenetic tree constructed using LR8 and Lrp13 members from several species along with LRP4 members as the out-group, Figure S2). The gene TRINITY_DN701_c1_g1 was designated as lr8, while the gene TRINITY_DN2157_c0_g1 functionally annotated as an LRP4 member was actually designated as lrp13.

Gene Expression of Cell Junctions during the PV-EV Transition.
The group of candidate TJ-constituting proteins was not significantly enrich among the DEGs (Fisher's exact test, p value = 0.07). Among the 25 identified genes e coding different types of TJ-constituting proteins, seven were differentially expresse only two genes were up-regulated during the PV-EV transition (genes encoding claud 12 and MARVEL domain containing protein 2-like-MarvelD2-like/tricellulin-like-

Discussion
A de novo transcriptome from A. australis ovarian tissue was successfully constructed, and an RNAseq approach was used to compare the transcriptomes from PV and EV ovaries. After histological analysis, the sampled groups of wild-caught eels (n = 6 per stage) were a clear representation of these developmental stages, characterising the progression into the early period of vitellogenic growth, when Vtg uptake and accumulation starts to be evident. Great ultrastructural changes occur to the ovarian tissue during development, which are also reflected in its RNA profiling. Normally in teleost fish, ovarian RNA composition changes from a low amount of 18S and 28S rRNA and high amount of 5S rRNA at early stages of development, to higher abundance of 18S and 28S rRNA in more advanced stages [58][59][60][61]. A similar pattern was found in A. australis total RNA profiles from both stages, reflecting the differences in the RIN values ( Figure S1). The GO enrichment analysis also reflected the differences in transcription between PV and EV stages. While only 3 BP terms (endocytosis GO:0006897, membrane organization GO:0061024, and regulation of cation transmembrane transport GO:1904062) and 1 MF term (potassium channel inhibitory activity GO:0019870) were up-regulated during the PV-EV transition, more than 15 terms were instead down-regulated in each GO category, mainly related to RNA metabolism and egg coat formation.
The RNAseq data was further used to investigate which genes may be involved in promoting competence to actively accumulate Vtg in this species. Indeed, the possible involvement of cell junctions in regulating Vtg uptake was evaluated by testing two hypotheses (Section 4.1 and 4.2). Notably, the down-regulated GO terms integral component of membrane (GO:0016021) and identical protein binding (GO:0042802) were found to contain genes encoding TJ-constituting proteins used in hypothesis 1. Similarly, both the endocytosis and membrane organisation GO terms contained genes evaluated further in Section 4.3, where the expression of genes related to receptor-mediated endocytosis machinery and Vtg processing was examined. Ovarian gene expression, represented by colour intensity, is compared between the PV stage and the EV stage. Colour scale gradient from red to blue represents levels of gene expression ranging from large to small, respectively. Up-and down-regulated genes during the PV-EV transition are shown. Each row represents the expression of a gene across all biological replicates from the PV stage or the EV stage (columns).

Discussion
A de novo transcriptome from A. australis ovarian tissue was successfully constructed, and an RNAseq approach was used to compare the transcriptomes from PV and EV ovaries. After histological analysis, the sampled groups of wild-caught eels (n = 6 per stage) were a clear representation of these developmental stages, characterising the progression into the early period of vitellogenic growth, when Vtg uptake and accumulation starts to be evident. Great ultrastructural changes occur to the ovarian tissue during development, which are also reflected in its RNA profiling. Normally in teleost fish, ovarian RNA composition changes from a low amount of 18S and 28S rRNA and high amount of 5S rRNA at early stages of development, to higher abundance of 18S and 28S rRNA in more advanced stages [58][59][60][61]. A similar pattern was found in A. australis total RNA profiles from both stages, reflecting the differences in the RIN values ( Figure S1). The GO enrichment analysis also reflected the differences in transcription between PV and EV stages. While only 3 BP terms (endocytosis GO:0006897, membrane organization GO:0061024, and regulation of cation transmembrane transport GO:1904062) and 1 MF term (potassium channel inhibitory activity GO:0019870) were up-regulated during the PV-EV transition, more than 15 terms were instead down-regulated in each GO category, mainly related to RNA metabolism and egg coat formation.
The RNAseq data was further used to investigate which genes may be involved in promoting competence to actively accumulate Vtg in this species. Indeed, the possible involvement of cell junctions in regulating Vtg uptake was evaluated by testing two hypotheses (Sections 4.1 and 4.2). Notably, the down-regulated GO terms integral component of membrane (GO:0016021) and identical protein binding (GO:0042802) were found to contain genes encoding TJ-constituting proteins used in hypothesis 1. Similarly, both the endocytosis and membrane organisation GO terms contained genes evaluated further in Section 4.3, where the expression of genes related to receptor-mediated endocytosis machinery and Vtg processing was examined.

Mechanical Barrier Hypothesis and Expression of TJ-Constituting Proteins during the PV-EV Transition
The TJ network consists of a group of transmembrane proteins, including claudins, junctional adhesion proteins (JAMs), angulin proteins (LSR, ILDR1, ILDR2), and MARVEL domain proteins (occludin, MARVELD2 or tricellulin, MARVELD3), that are externally associated with their counterparts on an adjacent cell [62]. Intracellularly, they bind to peripheral membrane proteins (tight junction ZO-1, ZO-2, ZO-3) which in turn make contact with the actin cytoskeleton [63]. The main functions of TJs are to maintain tissue integrity and cell polarity, as well as to modulate the permeability by restricting paracellular transport [13]. Evidently, macromolecules like Vtg are able to diffuse across the granulosa cell layer and reach the oocyte surface during active vitellogenic growth. The results obtained in this study showed that among 25 genes encoding TJ-constituting proteins, five genes were down-regulated compared to the two genes up-regulated during the PV-EV transition.
Between neighbouring cells, specific TJ-constituting proteins tend to localise along bicellular contacts, creating so-called bicellular TJs (bTJs), of which claudins are the major components determining the paracellular permeability [13]. The evident down-regulation in gene expression of TJ-encoding elements during the PV-EV transition in A. australis (claudin-like protein ZF-A9, claudin-3-like, MarvelD3, claudin-like protein ZF-A89, and Jam2a) matches the current evidence available from other teleost fish reinforcing the idea that stage-specific changes in TJ-constituting protein abundance may contribute to modulating Vtg paracellular transport. In D. rerio, from a total of 20 genes studied (18 claudins and 2 occludins), six claudins were down-regulated from the PV stage to mid-/late-vitellogenic stage [18], whereas in Gadus morhua, more than 10 genes involved in the TJ pathway were down-regulated from the PV stage to the EV stage [64]. However, the expression of different types of TJ-constituting proteins varies between the species studied, suggesting that species-specific interactions may be ruling this phenomenon. For instance, through observing a stage-specific down-regulation of occludin, Schuster et al. [65] suggested that in avian follicles TJs prevent the access of Vtg to the oocyte surface during early development. Yet, in the present study, the occludin-like gene was not differentially expressed during the PV-EV transition, indicating it may not play a key role during these stages, or it may carry out other functions in A. australis ovarian tissue.
In mammals, a mechanism through which an epithelium can allow the passage of macromolecules without compromising tissue integrity has been proposed. At tricellular contacts, macromolecules are able to diffuse through a central tube formed by bTJs and tricellular TJs (tTJs) [13,66,67]. Tricellulin and angulins are mainly localised at these contact points and may be responsible for sealing this paracellular pathway [66,68,69]. Furthermore, a similar mechanism appears to control paracellular transport through the follicle epithelium during oogenesis in insects (Drosophila sp.) [70], indicating it could be well-conserved across distant taxa. Based on this, we were expecting PV ovaries to highly express tTJ-constituting proteins in order to establish a barrier for macromolecules. However, the genes encoding angulins were not differentially expressed between stages, and unexpectedly, a gene encoding a tricellulin-like protein (MarvelD2-like) was up-regulated in EV ovaries. Then, is it possible that the down-regulation of TJ-constituting proteins is not the only mechanism driving a change in the TJ network promoting permeability? Indeed, there are contradicting findings on how tricellulin expression may affect the macromolecular pathway. Although there is evidence indicating that at low tricellulin expression a pathway for macromolecules is formed in mammalian cells [71], when studying an EpH4 mouse mammary epithelial cell line and comparing wild-type versus tricellulin knockdown cells (by RNA interference), Ikenouchi et al. [66] discovered that the flux of 250 Kda FITC-dextran did not differ between cell types. Accordingly, the increase in expression of the tricellulin-like coding gene in A. australis EV ovary may not be indicative of a decrease in macromolecule flux, as active Vtg accumulation has already started at this stage.
Although genes encoding TJ-constituting proteins were not enriched in the DE group, and the expression patterns of some TJ-encoding proteins were not as expected (between up-regulation and non-differential expression, i.e., tTJs, occludin), the hypothesis should not be rejected. Instead, alternative ways to test it should be explored aiming to elucidate important questions. What is the main force driving tissue permeability? Is it the abundance of TJ-constituting proteins or the interactions between them which creates a tighter epithelium? Could changes in specific TJ-constituting proteins modify the TJ network? Although each TJ-constituting protein may exhibit a unique function, it is likely that the combination of different proteins is what determines the final permeability of an epithelium. Therefore, the differential expression patterns of specific TJ-constituting proteins between the PV stage and the EV stage (including both up-and down-regulations) may create distinct TJ interactions, which in turn could be affecting the macromolecular flux.
Finally, there were putative paralogs (Marvel domain containing proteins and claudin ZF-A89) showing different patterns of expression in terms of regulation (differential expression between stages) and level of expression (number of reads). Intuitively, those genes contributing with a higher number of reads will be the ones contributing more to the pool of translated proteins. However, the gene expression of follicle cells could be diluted or masked by oocyte-specific transcripts, creating difficulties when interpreting the data. Therefore, the localisation of the expressed genes and their function at the protein level is warranted to elucidate their contribution to the TJ-network during follicle development. Indeed, further information will be needed to fully support the idea that a change in the TJ network may facilitate the accessibility of Vtg to the oocyte surface, in turn modulating Vtg uptake. Targeted studies using qPCR, in situ hybridization, knock down or knock out techniques, and immunohistochemistry would provide valuable data to continue testing this hypothesis.

Chemical Barrier Hypothesis and Expression of GJ Proteins during the PV-EV Transition
During oogenesis, the intercommunication between the oocyte and the granulosa cells is crucial for follicle development [72]. The exchange of information can be done through GJ coupling to establish the formation of channels between neighbouring cells, maintaining a rapid and direct communication to synchronise metabolic activities [73]. In teleost fish, the temporal correlation between the appearance of GJs and the onset of vitellogenic growth suggests the association between the two phenomena [74]. Using X. laevis as a model, Monaco et al. [20] concluded that direct GJ communication is required for Vtg uptake as it may be involved in transmitting signalling molecules that may modulate such uptake. Luque et al. [21] subsequently demonstrated that cAMP was the messenger that was transmitted through GJs, regulating the onset and ongoing endocytosis of Vtg in X. laevis. Afterwards, Serrano et al. [75] showed that vitellogenic growth in X. laevis is actively regulated by an increase in the intracellular levels of Ca 2+ , presumably from the cAMP-induced activation of the IP 3 / Ca 2+ pathway. Along with this information, extracellular Ca 2+ is also known to be required for Vtg receptor binding and subsequent uptake in Xenopus oocytes [76,77], suggesting that both intracellular and extracellular levels of Ca 2+ may be implicated in modulating vitellogenic growth.
Even though the predicted scenario for A. australis, in which genes encoding GJ proteins are up-regulated during the PV-EV transition, did not occur, a role for Ca 2+ signalling during vitellogenic growth in A. australis is possible. This notion is supported by the up-regulation of the gene coding for claudin-12 and of the cation transmembrane transport pathway (GO:1904062) in the EV stage. Indeed, both claudin-12 and WNK3, a member of the serine-threonine protein kinases group of genes which were enriched in the biological process GO:1904062, are implicated in Ca 2+ transport. While the former has been implicated in forming channels for Ca 2+ transport in mammalian intestinal epithelia [78,79], the latter plays a role as a positive regulator of the transcellular Ca 2+ transport pathway (human gene expressed in X. laevis oocytes [80]). Lastly, both extracellular and intracellular Ca 2+ are also required for TJ integrity and functioning [81], which further reinforces the possibility that Ca 2+ may regulate Vtg uptake. Whether these proteins perform the same function in teleost fish needs to be confirmed; however, our current data seem to support the implication of cell junctions and calcium ions in regulation of Vtg uptake. Indeed, a deeper understanding on GJ functioning in PV and EV follicles is still needed to extract comprehensive conclusions.

Expression of Genes Related to Receptor-Mediated Endocytosis Machinery and Vtg Processing during the PV-EV Transition
Following receptor recognition by specific members of the LDLR family, Vtg is incorporated into the oocytes through clathrin-mediated endocytosis [7]. Once inside the oocyte, Vtg is processed into yolk proteins by proteolysis, and the receptor is recycled back to the membrane surface [7]. Consequently, our RNAseq analysis was also used to examine the expression of genes related to Vtg recognition, the molecular machinery of clathrin-mediated endocytosis and vesicle trafficking, and Vtg proteolysis.
Two genes encoding putative Vtgrs were found in the A. australis ovarian transcriptome. The deduced proteins from both genes clustered into distinct groups corresponding to Lr8 and Lrp13 members from the LDLR family. Both members have been confirmed to bind Vtg in perciform and salmonid species, i.e., white perch, Morone americana [3], rainbow trout, O. mykiss [82] and cutthroat trout, O. clarki [83]. Therefore, they were both considered in this study as putative Vtgrs. In terms of expression, no difference was found between the PV and EV stages for both genes in A. australis. Similar findings for A. australis have been previously found for lr8 expression by Damsteegt et al. [84] and Jéhannet et al. [85], and in A. anguilla by Morini et al. [5], reinforcing the idea that expression of putative Vtgr does not seem to be a limiting factor for Vtg uptake.
Of note, Opresko & Wiley [86] discovered that under hormone stimulation, enhanced uptake of Vtg in Xenopus oocytes was not due to an increase in Vtgr abundance at the oocyte surface, but to an increase of its internalisation rate. Accordingly, other processes downstream of the endocytic pathway of Vtg (clathrin-mediated endocytosis, vesicle trafficking, Vtg proteolysis) may be involved in the regulation of Vtg accumulation. Sure enough, endocytosis (GO:0006897) and membrane organisation (GO:0061024), both containing genes that actively participate in clathrin-mediated endocytosis, were biological processes found to be up-regulated during the PV-EV transition. Even though clathrin-mediated endocytosis and its associated machinery involved in vesicle trafficking are not exclusive to Vtg uptake, and processes like lipid droplet and cortical alveoli formation might overlap, the active accumulation of Vtg can be considered as the main process happening in the ovary during the PV-EV transition. Knowledge of the function of these specific proteins will be necessary to better understand their contribution to Vtg uptake as they may participate in multiple pathways.
Interestingly, the adaptor protein Arhb, which is known to be required for efficient endocytosis of LDLR family members [87], including that of Vtg in X. laevis [88,89], was not differentially expressed between stages, indicating it may not play an essential role during vitellogenic growth in A. australis. However, the up-regulation of endophilin A3-like and the near up-regulation of endophilin-A3 (log 2 FC −0.94) in EV ovaries of A. australis confirms the pattern found in chicken oocytes. Indeed, Hirayama et al. [90] showed compelling evidence suggesting that in chicken, endophilin A3 is an important component of clathrin-mediated endocytosis during oocyte growth. Finally, other elements from the endocytic machinery (Calm, Numb, sorting nexin 17, etc.) showed differential expression patterns during the PV-EV transition in A. australis, differing to what was found in the swordfish Xiphias gladius when comparing immature and matured ovaries [91]. While differential expression patterns were also found in X. gladius, the DEGs related to the endocytic pathway differed between species.
After endocytosis, Vtg is proteolytically cleaved by specific proteases, i.e., cathepsins. Most of the cathepsins found in the transcriptome were not differentially expressed between PV and EV stages; except for the down-regulation in EV ovaries of cathepsin L1 and nothepsin, the latter likely to be a paralog of cathepsin D [92]. Finally, the early expression of genes encoding these proteases agrees with previous data suggesting that cathepsins may be regulated post-transcriptionally, leading to difficulties when trying to use expression data to explain activity [93][94][95][96]. This issue could be problematic especially in the light of our finding that translation-related GO terms have been down-regulated during the PV-EV transition (MF translation regulator activity GO:0045182, BP cytoplasmic translation GO:0002181, and BP translation GO:0006412). Indeed, additional data on gene and protein expression is needed to further support our findings.

Conclusions
The assembly of a de novo ovarian transcriptome of A. australis, followed by an RNAseq-based study comparing PV and EV ovaries, enabled us to obtain novel gene information from candidates that could be involved in Vtg uptake regulation. Two hypotheses that both implicate cell junctions in Vtg uptake were tested. The first hypothesis-the mechanical barrier hypothesis-was supported on the basis of the differential expression of genes encoding TJ-constituting proteins during the PV-EV transition, suggesting a role for modulating Vtg access to the oocyte surface. No compelling support for the second hypothesis-the chemical barrier hypothesis-was obtained. However, genes encoding clathrin-mediated endocytosis and vesicle trafficking elements were also clearly up-regulated during the PV-EV transition, indicating they, too, may be important factors modulating Vtg accumulation once the protein reaches the oolemma. Lastly, two genes encoding putative Vtgrs were found in A. australis, representing Lr8 (including Lr8+/variants) and Lrp13 members of the LDLR supergene family. Comparing transcriptomes from ovarian tissue at different developmental stages remains a challenge when trying to understand physiological pathways due to different cell types contributing to changing RNA profiles. In addition, the absence of a reliable reference genome plus the known retention of paralogues in teleost fish increases the complexity of next generation sequencing data analysis. However, studies based on comparative analysis of transcript abundances may provide valuable insight into the molecular events and physiological pathways underlying oocyte growth, which can be further investigated in future studies. In conclusion, combination of these data, along with protein expression, functioning, and localisation, is warranted to further unravel these pathways.